Skip to content
Snippets Groups Projects
Commit 11b7ed71 authored by Oliver Sander's avatar Oliver Sander
Browse files

Remove the BSplineBasis implementation 

Development continues in dune-functions
parent 9fe2dd26
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/fufem/functionspacebases/CMakeLists.txt
+ 0
− 1
View file @ 11b7ed71
install(FILES install(FILES
bsplinebasis.hh
conformingbasis.hh conformingbasis.hh
dgpqkbasis.hh dgpqkbasis.hh
dofconstraints.hh dofconstraints.hh
......
dune/fufem/functionspacebases/Makefile.am
+ 1
− 2
View file @ 11b7ed71
SUBDIRS = SUBDIRS =
functionspacebasesdir = $(includedir)/dune/fufem/functionspacebases functionspacebasesdir = $(includedir)/dune/fufem/functionspacebases
functionspacebases_HEADERS = bsplinebasis.hh \ functionspacebases_HEADERS = conformingbasis.hh \
conformingbasis.hh \
dofconstraints.hh \ dofconstraints.hh \
dunefunctionsbasis.hh \ dunefunctionsbasis.hh \
extensionbasis.hh \ extensionbasis.hh \
......
dune/fufem/functionspacebases/bsplinebasis.hh deleted 100644 → 0
+ 0
− 662
View file @ 9fe2dd26
#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
dune/fufem/test/CMakeLists.txt
+ 0
− 1
View file @ 11b7ed71
...@@ -5,7 +5,6 @@ set(GRID_BASED_TESTS ...@@ -5,7 +5,6 @@ set(GRID_BASED_TESTS
basisinterpolatortest basisinterpolatortest
basisinterpolationmatrixtest basisinterpolationmatrixtest
boundarypatchtest boundarypatchtest
bsplinebasistest
coarsegridfunctionwrappertest coarsegridfunctionwrappertest
composedfunctiontest composedfunctiontest
functionintegratortest functionintegratortest
......
dune/fufem/test/Makefile.am
+ 0
− 6
View file @ 11b7ed71
...@@ -5,7 +5,6 @@ TESTS = arithmetictest\ ...@@ -5,7 +5,6 @@ TESTS = arithmetictest\
basisinterpolatortest \ basisinterpolatortest \
basisinterpolationmatrixtest \ basisinterpolationmatrixtest \
boundarypatchtest \ boundarypatchtest \
bsplinebasistest \
coarsegridfunctionwrappertest \ coarsegridfunctionwrappertest \
composedfunctiontest \ composedfunctiontest \
constantfunctiontest \ constantfunctiontest \
...@@ -78,11 +77,6 @@ boundarypatchtest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS) ...@@ -78,11 +77,6 @@ boundarypatchtest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS)
boundarypatchtest_LDADD = $(COMMON_LDADD) $(GRID_LDADD) boundarypatchtest_LDADD = $(COMMON_LDADD) $(GRID_LDADD)
boundarypatchtest_LDFLAGS = $(COMMON_LDFLAGS) $(GRID_LDFLAGS) 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_SOURCES = coarsegridfunctionwrappertest.cc
coarsegridfunctionwrappertest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS) coarsegridfunctionwrappertest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS)
coarsegridfunctionwrappertest_LDADD = $(COMMON_LDADD) $(GRID_LDADD) coarsegridfunctionwrappertest_LDADD = $(COMMON_LDADD) $(GRID_LDADD)
......
dune/fufem/test/bsplinebasistest.cc deleted 100644 → 0
+ 0
− 162
View file @ 9fe2dd26
#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;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment