diff --git a/dune/fufem/functionspacebases/CMakeLists.txt b/dune/fufem/functionspacebases/CMakeLists.txt
index b8c414a0d9f953a692993e6637de7c52170f5619..ea603052cf7d5e00faf72d48cef8dc8224a44e79 100644
--- a/dune/fufem/functionspacebases/CMakeLists.txt
+++ b/dune/fufem/functionspacebases/CMakeLists.txt
@@ -1,5 +1,4 @@
install(FILES
- bsplinebasis.hh
conformingbasis.hh
dgpqkbasis.hh
dofconstraints.hh
diff --git a/dune/fufem/functionspacebases/Makefile.am b/dune/fufem/functionspacebases/Makefile.am
index 30ff6057cf7f5e1869d6e1f8f07e7ffebb308d65..7200901722cd1314168ba8fdf7ae335fd3254bb6 100644
--- a/dune/fufem/functionspacebases/Makefile.am
+++ b/dune/fufem/functionspacebases/Makefile.am
@@ -1,8 +1,7 @@
SUBDIRS =
functionspacebasesdir = $(includedir)/dune/fufem/functionspacebases
-functionspacebases_HEADERS = bsplinebasis.hh \
- conformingbasis.hh \
+functionspacebases_HEADERS = conformingbasis.hh \
dofconstraints.hh \
dunefunctionsbasis.hh \
extensionbasis.hh \
diff --git a/dune/fufem/functionspacebases/bsplinebasis.hh b/dune/fufem/functionspacebases/bsplinebasis.hh
deleted file mode 100644
index 2c65677b28df8e2b4ececc6b76858d19059c7c86..0000000000000000000000000000000000000000
--- a/dune/fufem/functionspacebases/bsplinebasis.hh
+++ /dev/null
@@ -1,662 +0,0 @@
-#ifndef DUNE_FUFEM_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
-#define DUNE_FUFEM_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
-
-/** \file
- * \brief The B-spline global function space basis
- */
-
-#include <dune/common/dynmatrix.hh>
-
-#include <dune/fufem/functionspacebases/functionspacebasis.hh>
-
-namespace Dune
-{
- namespace Fufem {
-
- // Forward declaration, to enable friend declarations
- template <class D, class R, int dim>
- class BSplineLocalFiniteElement;
-
- // Forward declaration, to enable friend declarations
- template <class D, class R, int dim>
- class BSplineLocalBasis;
-
- /** \brief Tensor product of 1d B-Spline functions with arbitrary knot vectors
- *
- * \internal This is an internal implementation class, used to implement the BSplineBasis class.
- *
- * \tparam D Number type used for domain coordinates
- * \tparam R Number type used for spline function values
- * \tparam dim Dimension of the patch
- */
- template<class D, class R, int dim>
- class BSplinePatch
- {
- friend class BSplineLocalFiniteElement<D,R,dim>;
- friend class BSplineLocalBasis<D,R,dim>;
-
- /** \brief Simple dim-dimensional multi-index class */
- class MultiIndex
- {
- public:
-
- /** \brief Constructs a new multi-index, and sets all digits to zero
- * \param limits Number of different digit values for each digit, i.e., digit i counts from 0 to limits[i]-1
- */
- MultiIndex(const std::array<unsigned int,dim>& limits)
- : limits_(limits)
- {
- std::fill(counter_.begin(), counter_.end(), 0);
- }
-
- /** \brief Increment the multi-index */
- MultiIndex& operator++()
- {
- for (int i=0; i<dim; i++)
- {
- ++counter_[i];
-
- // no overflow?
- if (counter_[i] < limits_[i])
- break;
-
- counter_[i] = 0;
- }
- return *this;
- }
-
- const unsigned int& operator[](int i) const
- {
- return counter_[i];
- }
-
- private:
-
- /** \brief The number of different digit values for each place */
- const std::array<unsigned int,dim> limits_;
-
- /** \brief The values of the multi-index. Each array entry has one digit */
- std::array<unsigned int,dim> counter_;
-
- };
-
- public:
-
- /** \brief Constructor: same knot vector and order in all space directions
- * \param makeOpen If this is true, then knots are prepended and appended to the knot vector to make the knot vector 'open'.
- * i.e., start and end with 'order+1' identical knots. Basis functions from such knot vectors are interpolatory at
- * the end of the parameter interval.
- */
- BSplinePatch(const std::vector<R>& knotVector,
- unsigned int order,
- bool makeOpen)
- {
- for (int i=0; i<dim; i++)
- {
- // Prepend the correct number of additional knots to open the knot vector
- //! \todo maybe test whether the knot vector is already open?
- if (makeOpen)
- for (unsigned int j=0; j<order; j++)
- knotVectors_[i].push_back(knotVector[0]);
-
- knotVectors_[i].insert(knotVectors_[i].end(), knotVector.begin(), knotVector.end());
-
- if (makeOpen)
- for (unsigned int j=0; j<order; j++)
- knotVectors_[i].push_back(knotVector.back());
- }
-
- std::fill(order_.begin(), order_.end(), order);
- }
-
- //! \brief Total number of B-spline basis functions
- unsigned int size () const
- {
- unsigned int result = 1;
- for (size_t i=0; i<dim; i++)
- result *= size(i);
- return result;
- }
-
- /** \brief Evaluate all B-spline basis functions at a given point
- *
- * \warning This really evaluates _all_ shape functions, even though most of them
- * are zero at any given point.
- */
- void evaluateFunction (const FieldVector<D,dim>& in,
- std::vector<FieldVector<R,1> >& out,
- const std::array<uint,dim>& currentKnotSpan) const
- {
- out.resize(size());
-
- // Evaluate
- Dune::array<std::vector<R>, dim> oneDValues;
-
- for (size_t i=0; i<dim; i++)
- evaluateFunction(in[i], oneDValues[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
-
- std::array<unsigned int, dim> limits;
- for (int i=0; i<dim; i++)
- limits[i] = size(i);
-
- MultiIndex ijkCounter(limits);
-
- for (size_t i=0; i<out.size(); i++, ++ijkCounter)
- {
- out[i] = R(1.0);
- for (size_t j=0; j<dim; j++)
- out[i] *= oneDValues[j][ijkCounter[j]];
- }
- }
-
- //! \brief Evaluate Jacobian of all B-spline basis functions
- void evaluateJacobian (const FieldVector<D,dim>& in,
- std::vector<FieldMatrix<D,1,dim> >& out,
- const std::array<uint,dim>& currentKnotSpan) const
- {
- out.resize(size());
-
- // Evaluate 1d function values (needed for the product rule)
- Dune::array<std::vector<R>, dim> oneDValues;
-
- for (size_t i=0; i<dim; i++)
- evaluateFunction(in[i], oneDValues[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
-
- // Evaluate 1d function values of one order lower (needed for the derivative formula)
- Dune::array<std::vector<R>, dim> lowOrderOneDValues;
-
- /** \todo Calling evaluateFunction again here is a waste: the lower-order values have
- * already been computed during the first call to evaluateFunction. Still, for the time
- * being I leave it as is to have more readable code. */
- for (size_t i=0; i<dim; i++)
- if (order_[i]!=0)
- evaluateFunction(in[i], lowOrderOneDValues[i], knotVectors_[i], order_[i]-1, currentKnotSpan[i]);
-
- // Evaluate 1d function derivatives
- Dune::array<std::vector<R>, dim> oneDDerivatives;
- for (size_t i=0; i<dim; i++)
- {
- oneDDerivatives[i].resize(oneDValues[i].size());
- for (size_t j=0; j<oneDDerivatives[i].size(); j++)
- {
- if (order_[i]==0) // order-zero functions are piecewise constant, hence all derivatives are zero
- std::fill(oneDDerivatives[i].begin(), oneDDerivatives[i].end(), R(0.0));
- else
- oneDDerivatives[i][j] = order_[i] * (lowOrderOneDValues[i][j] / (knotVectors_[i][j+order_[i]]-knotVectors_[i][j])
- - lowOrderOneDValues[i][j+1] / (knotVectors_[i][j+order_[i]+1]-knotVectors_[i][j+1]) );
-
- }
- }
-
- // Set up a multi-index to go from consecutive indices to integer coordinates
- std::array<unsigned int, dim> limits;
- for (int i=0; i<dim; i++)
- limits[i] = size(i);
-
- MultiIndex ijkCounter(limits);
-
- // Complete Jacobian is given by the product rule
- for (size_t i=0; i<out.size(); i++, ++ijkCounter)
- for (int j=0; j<dim; j++)
- {
- out[i][0][j] = 1.0;
- for (int k=0; k<dim; k++)
- out[i][0][j] *= (j==k) ? oneDDerivatives[k][ijkCounter[k]] : oneDValues[k][ijkCounter[k]];
- }
-
- }
-
- /** \brief Polynomial order of the shape functions
- * \todo Only implemented for 1d patches
- */
- unsigned int order () const
- {
- return order_;
- }
-
- private:
-
- //! \brief Number of shape functions in one direction
- unsigned int size (size_t d) const
- {
- return knotVectors_[d].size() - order_[d] - 1;
- }
-
-
- /** \brief Evaluate all one-dimensional B-spline functions for a given coordinate direction
- *
- * This implementations was based on the explanations in the book of
- * Cottrell, Hughes, Bazilevs, "Isogeometric Analysis"
- *
- * \param in Scalar(!) coordinate where to evaluate the functions
- * \param [out] out Vector containing the values of all B-spline functions at 'in'
- */
- static void evaluateFunction (const D& in, std::vector<R>& out,
- const std::vector<R>& knotVector,
- unsigned int order,
- unsigned int currentKnotSpan)
- {
- out.resize(knotVector.size()-order-1);
-
- // It's not really a matrix that is needed here, a plain 2d array would do
- DynamicMatrix<R> N(order+1, knotVector.size()-1);
-
- // The text books on splines use the following geometric condition here to fill the array N
- // (see for example CottreLL, Hughes, Bazilevs, Formula (2.1). However, this condition
- // only works if splines are never evaluated exactly on the knots.
- //
- // for (size_t i=0; i<knotVector.size()-1; i++)
- // N[0][i] = (knotVector[i] <= in) and (in < knotVector[i+1]);
- for (size_t i=0; i<knotVector.size()-1; i++)
- N[0][i] = (i == currentKnotSpan);
-
- for (size_t r=1; r<=order; r++)
- for (size_t i=0; i<knotVector.size()-r-1; i++)
- {
- R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
- ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
- : 0;
- R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
- ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
- : 0;
- N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
- }
-
- for (size_t i=0; i<out.size(); i++) {
- out[i] = N[order][i];
- }
- }
-
- /** \brief Order of the B-spline for each space dimension */
- array<unsigned int, dim> order_;
-
- /** \brief The knot vectors, one for each space dimension */
- array<std::vector<R>, dim> knotVectors_;
- };
-
-
- /** \brief LocalBasis class in the sense of dune-localfunctions, presenting the restriction
- * of a B-spline patch to a knot span
- *
- * \tparam D Number type used for domain coordinates
- * \tparam R Number type used for spline function values
- * \tparam dim Dimension of the patch
- */
- template<class D, class R, int dim>
- class BSplineLocalBasis
- {
- public:
-
- //! \brief export type traits for function signature
- typedef LocalBasisTraits<D,dim,Dune::FieldVector<D,dim>,R,1,Dune::FieldVector<R,1>,
- Dune::FieldMatrix<R,1,dim>, 2> Traits;
-
- /** \brief Constructor with a given B-spline patch
- *
- * The patch object does all the work.
- */
- BSplineLocalBasis(const BSplinePatch<D,R,dim>& patch)
- : patch_(patch)
- {}
-
- /** \brief Evaluate all shape functions
- * \param in Coordinates where to evaluate the functions, in local coordinates of the current knot span
- */
- void evaluateFunction (const FieldVector<D,dim>& in,
- std::vector<FieldVector<R,1> >& out) const
- {
- out.resize(size());
-
- FieldVector<D,dim> globalIn = offset_;
- scaling_.umv(in,globalIn);
-
- patch_.evaluateFunction(globalIn, out, currentKnotSpan_);
- }
-
- /** \brief Evaluate Jacobian of all shape functions
- * \param in Coordinates where to evaluate the Jacobian, in local coordinates of the current knot span
- */
- void evaluateJacobian (const FieldVector<D,dim>& in,
- std::vector<FieldMatrix<D,1,dim> >& out) const
- {
- out.resize(size());
-
- FieldVector<D,dim> globalIn = offset_;
- scaling_.umv(in,globalIn);
-
- patch_.evaluateJacobian(globalIn, out, currentKnotSpan_);
-
- for (size_t i=0; i<out.size(); i++)
- for (int j=0; j<dim; j++)
- out[i][0][j] *= scaling_[j][j];
- }
-
- //! \brief Evaluate all shape functions and derivatives of any order
- template<unsigned int k>
- inline void evaluate (const typename Dune::array<int,k>& directions,
- const typename Traits::DomainType& in,
- std::vector<typename Traits::RangeType>& out) const
- {
- if (k==0)
- evaluateFunction(in, out);
- else
- DUNE_THROW(NotImplemented, "B-Spline derivatives not implemented yet!");
- }
-
-
- /** \brief Polynomial order of the shape functions
- */
- unsigned int order () const
- {
- DUNE_THROW(NotImplemented, "!");
- }
-
- /** \brief Return the number of basis functions on the current knot span
- * \todo Currently, the number of all basis functions on the entire patch is returned.
- * This will include many that are simply constant zero on the current knot span.
- */
- uint size() const
- {
- return patch_.size();
- }
-
- /** \brief Elements are the non-empty knot spans, here we do the renumbering
- */
- void setCurrentKnotSpan(const std::array<uint,dim>& elementIdx)
- {
- /* \todo In the long run we need to precompute a table for this */
- for (int i=0; i<elementIdx.size(); i++)
- {
- currentKnotSpan_[i] = 0;
-
- // Skip over degenerate knot spans
- while (patch_.knotVectors_[i][currentKnotSpan_[i]+1] < patch_.knotVectors_[i][currentKnotSpan_[i]]+1e-8)
- currentKnotSpan_[i]++;
-
- for (int j=0; j<elementIdx[i]; j++)
- {
- currentKnotSpan_[i]++;
-
- // Skip over degenerate knot spans
- while (patch_.knotVectors_[i][currentKnotSpan_[i]+1] < patch_.knotVectors_[i][currentKnotSpan_[i]]+1e-8)
- currentKnotSpan_[i]++;
- }
-
- // Compute the geometric transformation from knotspan-local to global coordinates
- for (int i=0; i<dim; i++)
- {
- offset_[i] = patch_.knotVectors_[i][currentKnotSpan_[i]];
- scaling_[i][i] = patch_.knotVectors_[i][currentKnotSpan_[i]+1] - patch_.knotVectors_[i][currentKnotSpan_[i]];
- }
- }
- }
-
-
- const BSplinePatch<D,R,dim>& patch_;
-
- // Coordinates in a single knot span differ from coordinates on the B-spline patch
- // by an affine transformation. This transformation is stored in offset_ and scaling_.
- FieldVector<D,dim> offset_;
- DiagonalMatrix<D,dim> scaling_;
-
- // The knot span we are bound to
- std::array<uint,dim> currentKnotSpan_;
- };
-
- /** \brief Local coefficients in the sense of dune-localfunctions, for the B-spline basis on tensor-product grids
-
- \implements Dune::LocalCoefficientsVirtualImp
- */
- template <class D, class R, int dim>
- class BSplineLocalCoefficients
- {
- public:
- /** \brief Standard constructor
- * \todo Not implemented yet
- */
- BSplineLocalCoefficients (const BSplinePatch<D,R,dim>& patch)
- : patch_(patch),
- li_(size())
- {
- std::cout << "WARNING: LocalCoefficients array should be initialized here!" << std::endl;
- }
-
- /** \brief Number of coefficients
- * \todo Currently, the number of all basis functions on the entire patch is returned.
- * This will include many that are simply constant zero on the current knot span.
- */
- std::size_t size () const
- {
- return patch_.size();
- }
-
- //! get i'th index
- const LocalKey& localKey (std::size_t i) const
- {
- return li_[i];
- }
-
- private:
- const BSplinePatch<D,R,dim>& patch_;
- std::vector<LocalKey> li_;
- };
-
- /** \brief Local interpolation in the sense of dune-localfunctions, for the B-spline basis on tensor-product grids
- */
- template<int dim, class LB>
- class BSplineLocalInterpolation
- {
- public:
- //! \brief Local interpolation of a function
- template<typename F, typename C>
- void interpolate (const F& f, std::vector<C>& out) const
- {
- DUNE_THROW(NotImplemented, "BSplineLocalInterpolation::interpolate");
- }
-
- };
-
- /** \brief LocalFiniteElement in the sense of dune-localfunctions, for the B-spline basis on tensor-product grids
- *
- * This class ties together the implementation classes BSplineLocalBasis, BSplineLocalCoefficients, and BSplineLocalInterpolation
- *
- * \tparam D Number type used for domain coordinates
- * \tparam R Number type used for spline function values
- * \tparam dim Dimension of the patch
- */
- template<class D, class R, int dim>
- class BSplineLocalFiniteElement
- {
- public:
-
- /** \brief Export various types related to this LocalFiniteElement
- */
- typedef LocalFiniteElementTraits<BSplineLocalBasis<D,R,dim>,
- BSplineLocalCoefficients<D,R,dim>,
- BSplineLocalInterpolation<dim,BSplineLocalBasis<D,R,dim> > > Traits;
-
- /** \brief Constructor with a given B-spline patch
- */
- BSplineLocalFiniteElement(const BSplinePatch<D,R,dim>& patch)
- : localBasis_(patch),
- localCoefficients_(patch)
- {}
-
- /** \brief Bind LocalFiniteElement to a specific knot span of the spline patch
- * \param ijk Integer coordinates in the tensor product patch
- */
- void bind(const std::array<uint,dim>& elementIdx)
- {
- localBasis_.setCurrentKnotSpan(elementIdx);
- }
-
- /** \brief Hand out a LocalBasis object */
- const BSplineLocalBasis<D,R,dim>& localBasis() const
- {
- return localBasis_;
- }
-
- /** \brief Hand out a LocalCoefficients object */
- const BSplineLocalCoefficients<D,R,dim>& localCoefficients() const
- {
- return localCoefficients_;
- }
-
- /** \brief Hand out a LocalInterpolation object */
- const BSplineLocalInterpolation<dim,BSplineLocalBasis<D,R,dim> >& localInterpolation() const
- {
- return localInterpolation_;
- }
-
- /** \brief Number of shape functions in this finite element */
- uint size () const
- {
- return localBasis_.size();
- }
-
- /** \brief Return the reference element that the local finite element is defined on (here, a hypercube)
- */
- GeometryType type () const
- {
- return GeometryType(GeometryType::cube,dim);
- }
-
- private:
-
- BSplineLocalBasis<D,R,dim> localBasis_;
- BSplineLocalCoefficients<D,R,dim> localCoefficients_;
- BSplineLocalInterpolation<dim,BSplineLocalBasis<D,R,dim> > localInterpolation_;
-
- };
-
- /** \brief A B-spline function space basis on a tensor-product grid
- *
- * \tparam GV GridView, this must match the knot vectors describing the B-spline basis
- * \tparam RT Number type used for function values
- *
- * \todo Various features are not implemented yet:
- * - No multiple knots in a knot vector
- * - No sparsity; currently the implementation pretends that the support of any B-spline basis
- * function covers the entire patch.
- */
- template <class GV, class RT=double>
- class BSplineBasis :
- public FunctionSpaceBasis<
- GV,
- RT,
- BSplineLocalFiniteElement<double,double,GV::dimension> >
- {
- protected:
- typedef typename GV::Grid::ctype ctype;
-
- typedef BSplineLocalFiniteElement<double,double,GV::dimension> LFE;
-
- typedef FunctionSpaceBasis<GV, RT, LFE> Base;
- typedef typename Base::Element Element;
-
- using Base::dim;
- using Base::gridview_;
-
- public:
- typedef typename Base::GridView GridView;
- typedef typename Base::ReturnType ReturnType;
- typedef typename Base::LocalFiniteElement LocalFiniteElement;
- typedef typename Base::LinearCombination LinearCombination;
-
- /** \brief Construct a B-spline basis for a given grid view and set of knot vectors
- *
- * The grid *must* match the knot vectors, i.e.:
- * - The grid must be structured and Cartesian, and have cube elements only
- * - The number of elements in each direction must match the number of knot spans in that direction
- * - In fact, the element spacing in any direction must match the knot spacing in that direction
- * (disregarding knot multiplicities)
- * - When ordering the grid elements according to their indices, the resulting order must
- * be lexicographical, with the x-index increasing fastest.
- *
- * Unfortunately, not all of these conditions can be checked for automatically.
- *
- * \param knotVector A single knot vector, which will be used for all coordinate directions
- * \param order B-spline order, will be used for all coordinate directions
- * \param makeOpen If this is true, then knots are prepended and appended to the knot vector to make the knot vector 'open'.
- * i.e., start and end with 'order+1' identical knots. Basis functions from such knot vectors are interpolatory at
- * the end of the parameter interval.
- */
- BSplineBasis(const GridView& gridview,
- const std::vector<double>& knotVector,
- unsigned int order,
- bool makeOpen = true
- )
- : Base(gridview),
- patch_(knotVector, order, makeOpen),
- localFiniteElement_(patch_)
- {
- // \todo Detection of duplicate knots
- std::fill(elements_.begin(), elements_.end(), knotVector.size()-1);
-
- // Mediocre sanity check: we don't know the number of grid elements in each direction.
- // but at least we now the total number of elements.
- assert( std::accumulate(elements_.begin(), elements_.end(), 1, std::multiplies<uint>()) == gridview_.size(0) );
- }
-
- /** \brief Total number of basis vectors */
- size_t size() const
- {
- return patch_.size();
- }
-
- /** \brief Return LocalFiniteElement implementing the function space on a knot span
- */
- const LocalFiniteElement& getLocalFiniteElement(const Element& e) const
- {
- typename GridView::IndexSet::IndexType elementIndex = gridview_.indexSet().index(e);
- localFiniteElement_.bind(getIJK(elementIndex));
- return localFiniteElement_;
- }
-
- typename GV::IndexSet::IndexType index(const Element& e, const int i) const
- {
- /** \todo Currently, local per-knot-span numbering and global per-patch numbering
- * of dofs coincides. That is not desired -- most degrees of freedom are zero
- * on any given knot span.
- */
- return i;
- }
-
- protected:
-
- /** \brief Compute integer element coordinates from the element index
- * \warning This method makes strong assumptions about the grid, namely that it is
- * structured, and that indices are given in a x-fastest fashion.
- */
- std::array<uint,dim> getIJK(typename GridView::IndexSet::IndexType idx) const
- {
- std::array<uint,dim> result;
- for (int i=0; i<dim; i++)
- {
- result[i] = idx%elements_[i];
- idx /= elements_[i];
- }
- return result;
- }
-
- /** \brief The underlying B-spline basis patch */
- BSplinePatch<double,double,dim> patch_;
-
- /** \brief The LocalFiniteElement presenting views on different knot spans
- *
- * Needs to be mutable, because it must be able to be bound to different elements,
- * while leaving the BSplineBasis object const.
- */
- mutable LocalFiniteElement localFiniteElement_;
-
- /** \brief Number of grid elements in the different coordinate directions */
- std::array<uint,dim> elements_;
-
- };
-
- } // namespace Fufem
-
-} // namespace Dune
-
-#endif
-
diff --git a/dune/fufem/test/CMakeLists.txt b/dune/fufem/test/CMakeLists.txt
index b30e028d7f574e6b8c23f78a016067a3118f6b57..a277183c1cf9abf75172350da0a593ab14847038 100644
--- a/dune/fufem/test/CMakeLists.txt
+++ b/dune/fufem/test/CMakeLists.txt
@@ -5,7 +5,6 @@ set(GRID_BASED_TESTS
basisinterpolatortest
basisinterpolationmatrixtest
boundarypatchtest
- bsplinebasistest
coarsegridfunctionwrappertest
composedfunctiontest
functionintegratortest
diff --git a/dune/fufem/test/Makefile.am b/dune/fufem/test/Makefile.am
index a7473ae0a04292c3d86466052c0049b145debc91..7d66cf8a2c32ff229bd349923eec323766883f41 100644
--- a/dune/fufem/test/Makefile.am
+++ b/dune/fufem/test/Makefile.am
@@ -5,7 +5,6 @@ TESTS = arithmetictest\
basisinterpolatortest \
basisinterpolationmatrixtest \
boundarypatchtest \
- bsplinebasistest \
coarsegridfunctionwrappertest \
composedfunctiontest \
constantfunctiontest \
@@ -78,11 +77,6 @@ boundarypatchtest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS)
boundarypatchtest_LDADD = $(COMMON_LDADD) $(GRID_LDADD)
boundarypatchtest_LDFLAGS = $(COMMON_LDFLAGS) $(GRID_LDFLAGS)
-bsplinebasistest_SOURCES = bsplinebasistest.cc
-bsplinebasistest_CPPFLAGS = $(COMMON_CPPFLAGS)
-bsplinebasistest_LDADD = $(COMMON_LDADD)
-bsplinebasistest_LDFLAGS = $(COMMON_LDFLAGS)
-
coarsegridfunctionwrappertest_SOURCES = coarsegridfunctionwrappertest.cc
coarsegridfunctionwrappertest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS)
coarsegridfunctionwrappertest_LDADD = $(COMMON_LDADD) $(GRID_LDADD)
diff --git a/dune/fufem/test/bsplinebasistest.cc b/dune/fufem/test/bsplinebasistest.cc
deleted file mode 100644
index 3663fdf1f612653732e0c0d4d0e3d34c1a167639..0000000000000000000000000000000000000000
--- a/dune/fufem/test/bsplinebasistest.cc
+++ /dev/null
@@ -1,162 +0,0 @@
-#include "config.h"
-
-/** \file
- * \brief Unit tests for the BSplineBasis class
- */
-
-#include <dune/grid/yaspgrid.hh>
-#include <dune/grid/io/file/vtk.hh>
-
-#include <dune/localfunctions/test/test-localfe.hh>
-
-#include <dune/fufem/functions/vtkbasisgridfunction.hh>
-#include <dune/fufem/functionspacebases/bsplinebasis.hh>
-
-using namespace Dune;
-
-template <class Basis>
-void test(const Basis& basis, bool isPartitionOfUnity)
-{
- typedef typename Basis::GridView GridView;
- GridView gridView = basis.getGridView();
- static const int dim = GridView::dimension;
-
- // Test the LocalFiniteElement
- for (auto it = gridView.template begin<0>(); it!=gridView.template end<0>(); ++it)
- {
- typedef typename Basis::LocalFiniteElement LocalFiniteElementType;
- const LocalFiniteElementType& lFE = basis.getLocalFiniteElement(*it);
-
- // The general LocalFiniteElement unit test from dune/localfunctions/test/test-localfe.hh
- // LocalInterpolation is not implemented yet
- testFE(lFE,DisableLocalInterpolation);
- }
-
- // Make sure the basis is a partition of unity
- if (isPartitionOfUnity)
- {
- for (auto it = gridView.template begin<0>(); it!=gridView.template end<0>(); ++it)
- {
- typedef typename Basis::LocalFiniteElement LocalFiniteElementType;
- const LocalFiniteElementType& lFE = basis.getLocalFiniteElement(*it);
-
- const QuadratureRule<double,dim>& quad = QuadratureRules<double,dim>::rule(it->type(), 3);
- std::vector<FieldVector<double,1> > values;
- for (size_t i=0; i<quad.size(); i++)
- {
- lFE.localBasis().evaluateFunction(quad[i].position(), values);
- double sum = std::accumulate(values.begin(), values.end(), 0.0);
-
- if (std::abs(sum-1.0) > 1e-5)
- DUNE_THROW(Exception, "B-Spline basis is no partition of unity!");
- }
- }
- }
-
- // Write all global basis functions as VTK files, so we can look at them
- for (size_t i=0; i<basis.size(); i++)
- {
- std::vector<FieldVector<double,1> > c(basis.size());
- std::fill(c.begin(), c.end(), 0.0);
- c[i] = 1.0;
-
- SubsamplingVTKWriter<GridView> vtkWriter(gridView,2);
-
- typedef VTKBasisGridFunction<Basis,std::vector<FieldVector<double,1> > > VTKBSplineFunction;
- std::shared_ptr<VTKBSplineFunction> vtkFunction = std::make_shared<VTKBSplineFunction> (basis,c, "value");
- vtkWriter.addVertexData(vtkFunction);
- vtkWriter.write("bsplineglobalbasis_" + std::to_string(i));
-
- }
-}
-
-int main(int argc, char* argv[]) try
-{
- //////////////////////////////////////////////////
- // Test on a 1d grid
- //////////////////////////////////////////////////
- {
- FieldVector<double,1> upper = {1};
- std::array<int,1> elements = {5};
- YaspGrid<1> grid1d(upper,elements);
-
- typedef YaspGrid<1>::LeafGridView GridView;
- GridView gridView = grid1d.leafGridView();
-
- // Create corresponding knot vector
- std::vector<double> knotVector(elements[0]+1);
- for (size_t i=0; i<knotVector.size(); i++)
- knotVector[i] = i * upper[0] / elements[0];
-
- // Create B-spline spaces of different order, and test it
- for (int order=0; order<3; order++)
- {
- typedef Fufem::BSplineBasis<GridView> BasisType;
- BasisType bSplineBasis(gridView, knotVector, order);
-
- // Test it
- test(bSplineBasis,
- true); // This basis is no partition of unity
- }
- }
-
- //////////////////////////////////////////////////
- // Test on a 2d grid
- //////////////////////////////////////////////////
- {
- FieldVector<double,2> upper = {1,1};
- std::array<int,2> elements = {5,5};
- YaspGrid<2> grid2d(upper,elements);
-
- typedef YaspGrid<2>::LeafGridView GridView;
- GridView gridView = grid2d.leafGridView();
-
- // Create corresponding knot vector
- std::vector<double> knotVector(elements[0]+1);
- for (size_t i=0; i<knotVector.size(); i++)
- knotVector[i] = i * upper[0] / elements[0];
-
- // Create B-spline spaces of different order, and test it
- for (int order=0; order<3; order++)
- {
- typedef Fufem::BSplineBasis<GridView> BasisType;
- BasisType bSplineBasis(gridView, knotVector, order);
-
- // Test it
- test(bSplineBasis,
- true); // This basis is a partition of unity
- }
- }
-
- //////////////////////////////////////////////////
- // Test on a 1d grid
- //////////////////////////////////////////////////
- {
- FieldVector<double,3> upper = {1,1,1};
- std::array<int,3> elements = {5,5,5};
- YaspGrid<3> grid3d(upper,elements);
-
- typedef YaspGrid<3>::LeafGridView GridView;
- GridView gridView = grid3d.leafGridView();
-
- // Create corresponding knot vector
- std::vector<double> knotVector(elements[0]+1);
- for (size_t i=0; i<knotVector.size(); i++)
- knotVector[i] = i * upper[0] / elements[0];
-
- // Create B-spline spaces of different order, and test it
- for (int order=0; order<3; order++)
- {
- typedef Fufem::BSplineBasis<GridView> BasisType;
- BasisType bSplineBasis(gridView, knotVector, order);
-
- // Test it
- test(bSplineBasis,
- true); // This basis is a partition of unity
- }
- }
-} catch (Exception e) {
- std::cout << e << std::endl;
-}
-
-