Loading...
Searching...
No Matches
bsplinebasis.hh
Go to the documentation of this file.
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
#include <array>
#include <numeric>
#include <dune/common/dynmatrix.hh>
#include <dune/localfunctions/common/localbasis.hh>
#include <dune/common/diagonalmatrix.hh>
#include <dune/localfunctions/common/localkey.hh>
#include <dune/localfunctions/common/localfiniteelementtraits.hh>
#include <dune/geometry/type.hh>
#include <dune/functions/functionspacebases/nodes.hh>
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/leafprebasismixin.hh>
namespace Dune
{
namespace Functions {
// A maze of dependencies between the different parts of this. We need a few forward declarations
template<typename GV, typename R>
class BSplineLocalFiniteElement;
template<typename GV>
class BSplinePreBasis;
template<class GV, class R>
class BSplineLocalBasis
{
friend class BSplineLocalFiniteElement<GV,R>;
typedef typename GV::ctype D;
enum {dim = GV::dimension};
public:
typedef LocalBasisTraits<D,dim,FieldVector<D,dim>,R,1,FieldVector<R,1>,
FieldMatrix<R,1,dim> > Traits;
BSplineLocalBasis(const BSplinePreBasis<GV>& preBasis,
const BSplineLocalFiniteElement<GV,R>& lFE)
: preBasis_(preBasis),
lFE_(lFE)
{}
void evaluateFunction (const FieldVector<D,dim>& in,
std::vector<FieldVector<R,1> >& out) const
{
FieldVector<D,dim> globalIn = offset_;
scaling_.umv(in,globalIn);
preBasis_.evaluateFunction(globalIn, out, lFE_.currentKnotSpan_);
}
void evaluateJacobian (const FieldVector<D,dim>& in,
std::vector<FieldMatrix<D,1,dim> >& out) const
{
FieldVector<D,dim> globalIn = offset_;
scaling_.umv(in,globalIn);
preBasis_.evaluateJacobian(globalIn, out, lFE_.currentKnotSpan_);
for (size_t i=0; i<out.size(); i++)
for (int j=0; j<dim; j++)
out[i][0][j] *= scaling_[j][j];
}
template<size_t k>
inline void evaluate (const typename std::array<int,k>& directions,
const typename Traits::DomainType& in,
std::vector<typename Traits::RangeType>& out) const
{
switch(k)
{
case 0:
evaluateFunction(in, out);
break;
case 1:
{
FieldVector<D,dim> globalIn = offset_;
scaling_.umv(in,globalIn);
preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
for (size_t i=0; i<out.size(); i++)
out[i][0] *= scaling_[directions[0]][directions[0]];
break;
}
case 2:
{
FieldVector<D,dim> globalIn = offset_;
scaling_.umv(in,globalIn);
preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
for (size_t i=0; i<out.size(); i++)
out[i][0] *= scaling_[directions[0]][directions[0]]*scaling_[directions[1]][directions[1]];
break;
}
default:
DUNE_THROW(NotImplemented, "B-Spline derivatives of order " << k << " not implemented yet!");
}
}
unsigned int order () const
{
return *std::max_element(preBasis_.order_.begin(), preBasis_.order_.end());
}
std::size_t size() const
{
return lFE_.size();
}
private:
const BSplinePreBasis<GV>& preBasis_;
const BSplineLocalFiniteElement<GV,R>& lFE_;
// 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_;
};
template<int dim>
class BSplineLocalCoefficients
{
// Return i as a d-digit number in the (k+1)-nary system
std::array<unsigned int,dim> multiindex (unsigned int i) const
{
std::array<unsigned int,dim> alpha;
for (int j=0; j<dim; j++)
{
alpha[j] = i % sizes_[j];
i = i/sizes_[j];
}
return alpha;
}
void setup1d(std::vector<unsigned int>& subEntity)
{
if (sizes_[0]==1)
{
subEntity[0] = 0;
return;
}
/* edge and vertex numbering
0----0----1
*/
unsigned lastIndex=0;
subEntity[lastIndex++] = 0; // corner 0
for (unsigned i = 0; i < sizes_[0] - 2; ++i)
subEntity[lastIndex++] = 0; // inner dofs of element (0)
subEntity[lastIndex++] = 1; // corner 1
assert(size()==lastIndex);
}
void setup2d(std::vector<unsigned int>& subEntity)
{
unsigned lastIndex=0;
// LocalKey: entity number , entity codim, dof indices within each entity
/* edge and vertex numbering
2----3----3
| |
| |
0 1
| |
| |
0----2----1
*/
// lower edge (2)
subEntity[lastIndex++] = 0; // corner 0
for (unsigned i = 0; i < sizes_[0]-2; ++i)
subEntity[lastIndex++] = 2; // inner dofs of lower edge (2)
subEntity[lastIndex++] = 1; // corner 1
// iterate from bottom to top over inner edge dofs
for (unsigned e = 0; e < sizes_[1]-2; ++e)
{
subEntity[lastIndex++] = 0; // left edge (0)
for (unsigned i = 0; i < sizes_[0]-2; ++i)
subEntity[lastIndex++] = 0; // face dofs
subEntity[lastIndex++] = 1; // right edge (1)
}
// upper edge (3)
subEntity[lastIndex++] = 2; // corner 2
for (unsigned i = 0; i < sizes_[0]-2; ++i)
subEntity[lastIndex++] = 3; // inner dofs of upper edge (3)
subEntity[lastIndex++] = 3; // corner 3
assert(size()==lastIndex);
}
public:
void init(const std::array<unsigned,dim>& sizes)
{
sizes_ = sizes;
li_.resize(size());
// Set up array of codimension-per-dof-number
std::vector<unsigned int> codim(li_.size());
for (std::size_t i=0; i<codim.size(); i++)
{
codim[i] = 0;
// Codimension gets increased by 1 for each coordinate direction
// where dof is on boundary
std::array<unsigned int,dim> mIdx = multiindex(i);
for (int j=0; j<dim; j++)
if (mIdx[j]==0 or mIdx[j]==sizes[j]-1)
codim[i]++;
}
// Set up index vector (the index of the dof in the set of dofs of a given subentity)
// Algorithm: the 'index' has the same ordering as the dof number 'i'.
// To make it consecutive we interpret 'i' in the (k+1)-adic system, omit all digits
// that correspond to axes where the dof is on the element boundary, and transform the
// rest to the (k-1)-adic system.
std::vector<unsigned int> index(size());
for (std::size_t i=0; i<index.size(); i++)
{
index[i] = 0;
std::array<unsigned int,dim> mIdx = multiindex(i);
for (int j=dim-1; j>=0; j--)
if (mIdx[j]>0 and mIdx[j]<sizes[j]-1)
index[i] = (sizes[j]-1)*index[i] + (mIdx[j]-1);
}
// Set up entity and dof numbers for each (supported) dimension separately
std::vector<unsigned int> subEntity(li_.size());
if (subEntity.size() > 0)
{
if (dim==1) {
setup1d(subEntity);
} else if (dim==2 and sizes_[0]>1 and sizes_[1]>1) {
setup2d(subEntity);
}
}
for (size_t i=0; i<li_.size(); i++)
li_[i] = LocalKey(subEntity[i], codim[i], index[i]);
}
std::size_t size () const
{
return std::accumulate(sizes_.begin(), sizes_.end(), 1, std::multiplies<unsigned int>());
}
const LocalKey& localKey (std::size_t i) const
{
return li_[i];
}
private:
// Number of shape functions on this element per coordinate direction
std::array<unsigned, dim> sizes_;
std::vector<LocalKey> li_;
};
template<int dim, class LB>
class BSplineLocalInterpolation
{
public:
template<typename F, typename C>
void interpolate (const F& f, std::vector<C>& out) const
{
DUNE_THROW(NotImplemented, "BSplineLocalInterpolation::interpolate");
}
};
template<class GV, class R>
class BSplineLocalFiniteElement
{
typedef typename GV::ctype D;
enum {dim = GV::dimension};
friend class BSplineLocalBasis<GV,R>;
public:
typedef LocalFiniteElementTraits<BSplineLocalBasis<GV,R>,
BSplineLocalCoefficients<dim>,
BSplineLocalInterpolation<dim,BSplineLocalBasis<GV,R> > > Traits;
BSplineLocalFiniteElement(const BSplinePreBasis<GV>& preBasis)
: preBasis_(preBasis),
localBasis_(preBasis,*this)
{}
BSplineLocalFiniteElement(const BSplineLocalFiniteElement& other)
: preBasis_(other.preBasis_),
localBasis_(preBasis_,*this)
{}
void bind(const std::array<unsigned,dim>& elementIdx)
{
/* \todo In the long run we need to precompute a table for this */
for (size_t i=0; i<elementIdx.size(); i++)
{
currentKnotSpan_[i] = 0;
// Skip over degenerate knot spans
while (preBasis_.knotVectors_[i][currentKnotSpan_[i]+1] < preBasis_.knotVectors_[i][currentKnotSpan_[i]]+1e-8)
currentKnotSpan_[i]++;
for (size_t j=0; j<elementIdx[i]; j++)
{
currentKnotSpan_[i]++;
// Skip over degenerate knot spans
while (preBasis_.knotVectors_[i][currentKnotSpan_[i]+1] < preBasis_.knotVectors_[i][currentKnotSpan_[i]]+1e-8)
currentKnotSpan_[i]++;
}
// Compute the geometric transformation from knotspan-local to global coordinates
localBasis_.offset_[i] = preBasis_.knotVectors_[i][currentKnotSpan_[i]];
localBasis_.scaling_[i][i] = preBasis_.knotVectors_[i][currentKnotSpan_[i]+1] - preBasis_.knotVectors_[i][currentKnotSpan_[i]];
}
// Set up the LocalCoefficients object
std::array<unsigned int, dim> sizes;
for (size_t i=0; i<dim; i++)
sizes[i] = size(i);
localCoefficients_.init(sizes);
}
const BSplineLocalBasis<GV,R>& localBasis() const
{
return localBasis_;
}
const BSplineLocalCoefficients<dim>& localCoefficients() const
{
return localCoefficients_;
}
const BSplineLocalInterpolation<dim,BSplineLocalBasis<GV,R> >& localInterpolation() const
{
return localInterpolation_;
}
unsigned size () const
{
std::size_t r = 1;
for (int i=0; i<dim; i++)
r *= size(i);
return r;
}
GeometryType type () const
{
return GeometryTypes::cube(dim);
}
//private:
unsigned int size(int i) const
{
const auto& order = preBasis_.order_;
unsigned int r = order[i]+1; // The 'normal' value
if (currentKnotSpan_[i]<order[i]) // Less near the left end of the knot vector
r -= (order[i] - currentKnotSpan_[i]);
if ( order[i] > (preBasis_.knotVectors_[i].size() - currentKnotSpan_[i] - 2) )
r -= order[i] - (preBasis_.knotVectors_[i].size() - currentKnotSpan_[i] - 2);
return r;
}
const BSplinePreBasis<GV>& preBasis_;
BSplineLocalBasis<GV,R> localBasis_;
BSplineLocalCoefficients<dim> localCoefficients_;
BSplineLocalInterpolation<dim,BSplineLocalBasis<GV,R> > localInterpolation_;
// The knot span we are bound to
std::array<unsigned,dim> currentKnotSpan_;
};
template<typename GV>
class BSplineNode;
template<typename GV>
class BSplinePreBasis :
public LeafPreBasisMixin< BSplinePreBasis<GV> >
{
using Base = LeafPreBasisMixin< BSplinePreBasis<GV> >;
static const int dim = GV::dimension;
class MultiDigitCounter
{
public:
MultiDigitCounter(const std::array<unsigned int,dim>& limits)
: limits_(limits)
{
std::fill(counter_.begin(), counter_.end(), 0);
}
MultiDigitCounter& 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];
}
unsigned int cycle() const
{
unsigned int r = 1;
for (int i=0; i<dim; i++)
r *= limits_[i];
return r;
}
private:
const std::array<unsigned int,dim> limits_;
std::array<unsigned int,dim> counter_;
};
public:
using GridView = GV;
using size_type = std::size_t;
using Node = BSplineNode<GV>;
// Type used for function values
using R = double;
BSplinePreBasis(const GridView& gridView,
const std::vector<double>& knotVector,
unsigned int order,
bool makeOpen = true)
: gridView_(gridView)
{
// \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 know the total number of elements.
assert( std::accumulate(elements_.begin(), elements_.end(), 1, std::multiplies<unsigned>()) == gridView_.size(0) );
for (int i=0; i<dim; i++)
{
// Prepend the correct number of additional knots to open the knot vector
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);
}
BSplinePreBasis(const GridView& gridView,
const FieldVector<double,dim>& lowerLeft,
const FieldVector<double,dim>& upperRight,
const std::array<unsigned int,dim>& elements,
unsigned int order,
bool makeOpen = true)
: elements_(elements),
gridView_(gridView)
{
// Mediocre sanity check: we don't know the number of grid elements in each direction.
// but at least we know the total number of elements.
assert( std::accumulate(elements_.begin(), elements_.end(), 1, std::multiplies<unsigned>()) == gridView_.size(0) );
for (int i=0; i<dim; i++)
{
// Prepend the correct number of additional knots to open the knot vector
if (makeOpen)
for (unsigned int j=0; j<order; j++)
knotVectors_[i].push_back(lowerLeft[i]);
// Construct the actual knot vector
for (size_t j=0; j<elements[i]+1; j++)
knotVectors_[i].push_back(lowerLeft[i] + j*(upperRight[i]-lowerLeft[i]) / elements[i]);
if (makeOpen)
for (unsigned int j=0; j<order; j++)
knotVectors_[i].push_back(upperRight[i]);
}
std::fill(order_.begin(), order_.end(), order);
}
void initializeIndices()
{}
const GridView& gridView() const
{
return gridView_;
}
void update(const GridView& gv)
{
gridView_ = gv;
}
Node makeNode() const
{
return Node{this};
}
size_type maxNodeSize() const
{
size_type result = 1;
for (int i=0; i<dim; i++)
result *= order_[i]+1;
return result;
}
template<typename It>
It indices(const Node& node, It it) const
{
// Local degrees of freedom are arranged in a lattice.
// We need the lattice dimensions to be able to compute lattice coordinates from a local index
std::array<unsigned int, dim> localSizes;
for (int i=0; i<dim; i++)
localSizes[i] = node.finiteElement().size(i);
for (size_type i = 0, end = node.size() ; i < end ; ++i, ++it)
{
std::array<unsigned int,dim> localIJK = getIJK(i, localSizes);
const auto currentKnotSpan = node.finiteElement().currentKnotSpan_;
const auto order = order_;
std::array<unsigned int,dim> globalIJK;
for (int i=0; i<dim; i++)
globalIJK[i] = std::max((int)currentKnotSpan[i] - (int)order[i], 0) + localIJK[i]; // needs to be a signed type!
// Make one global flat index from the globalIJK tuple
size_type globalIdx = globalIJK[dim-1];
for (int i=dim-2; i>=0; i--)
globalIdx = globalIdx * size(i) + globalIJK[i];
*it = {{globalIdx}};
}
return it;
}
unsigned int dimension () const
{
unsigned int result = 1;
for (size_t i=0; i<dim; i++)
result *= size(i);
return result;
}
unsigned int size (size_t d) const
{
return knotVectors_[d].size() - order_[d] - 1;
}
using Base::size;
void evaluateFunction (const FieldVector<typename GV::ctype,dim>& in,
std::vector<FieldVector<R,1> >& out,
const std::array<unsigned,dim>& currentKnotSpan) const
{
// Evaluate
std::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] = oneDValues[i].size();
MultiDigitCounter ijkCounter(limits);
out.resize(ijkCounter.cycle());
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]];
}
}
void evaluateJacobian (const FieldVector<typename GV::ctype,dim>& in,
std::vector<FieldMatrix<R,1,dim> >& out,
const std::array<unsigned,dim>& currentKnotSpan) const
{
// How many shape functions to we have in each coordinate direction?
std::array<unsigned int, dim> limits;
for (int i=0; i<dim; i++)
{
limits[i] = order_[i]+1; // The 'standard' value away from the boundaries of the knot vector
if (currentKnotSpan[i]<order_[i])
limits[i] -= (order_[i] - currentKnotSpan[i]);
if ( order_[i] > (knotVectors_[i].size() - currentKnotSpan[i] - 2) )
limits[i] -= order_[i] - (knotVectors_[i].size() - currentKnotSpan[i] - 2);
}
// The lowest knot spans that we need values from
std::array<unsigned int, dim> offset;
for (int i=0; i<dim; i++)
offset[i] = std::max((int)(currentKnotSpan[i] - order_[i]),0);
// Evaluate 1d function values (needed for the product rule)
std::array<std::vector<R>, dim> oneDValues;
// Evaluate 1d function values of one order lower (needed for the derivative formula)
std::array<std::vector<R>, dim> lowOrderOneDValues;
std::array<DynamicMatrix<R>, dim> values;
for (size_t i=0; i<dim; i++)
{
evaluateFunctionFull(in[i], values[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
oneDValues[i].resize(knotVectors_[i].size()-order_[i]-1);
for (size_t j=0; j<oneDValues[i].size(); j++)
oneDValues[i][j] = values[i][order_[i]][j];
if (order_[i]!=0)
{
lowOrderOneDValues[i].resize(knotVectors_[i].size()-(order_[i]-1)-1);
for (size_t j=0; j<lowOrderOneDValues[i].size(); j++)
lowOrderOneDValues[i][j] = values[i][order_[i]-1][j];
}
}
// Evaluate 1d function derivatives
std::array<std::vector<R>, dim> oneDDerivatives;
for (size_t i=0; i<dim; i++)
{
oneDDerivatives[i].resize(limits[i]);
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
{
for (size_t j=offset[i]; j<offset[i]+limits[i]; j++)
{
R derivativeAddend1 = lowOrderOneDValues[i][j] / (knotVectors_[i][j+order_[i]]-knotVectors_[i][j]);
R derivativeAddend2 = lowOrderOneDValues[i][j+1] / (knotVectors_[i][j+order_[i]+1]-knotVectors_[i][j+1]);
// The two previous terms may evaluate as 0/0. This is to be interpreted as 0.
if (std::isnan(derivativeAddend1))
derivativeAddend1 = 0;
if (std::isnan(derivativeAddend2))
derivativeAddend2 = 0;
oneDDerivatives[i][j-offset[i]] = order_[i] * ( derivativeAddend1 - derivativeAddend2 );
}
}
}
// Working towards computing only the parts that we really need:
// Let's copy them out into a separate array
std::array<std::vector<R>, dim> oneDValuesShort;
for (int i=0; i<dim; i++)
{
oneDValuesShort[i].resize(limits[i]);
for (size_t j=0; j<limits[i]; j++)
oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
}
// Set up a multi-index to go from consecutive indices to integer coordinates
MultiDigitCounter ijkCounter(limits);
out.resize(ijkCounter.cycle());
// 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]]
: oneDValuesShort[k][ijkCounter[k]];
}
}
template <size_type k>
void evaluate(const typename std::array<int,k>& directions,
const FieldVector<typename GV::ctype,dim>& in,
std::vector<FieldVector<R,1> >& out,
const std::array<unsigned,dim>& currentKnotSpan) const
{
if (k != 1 && k != 2)
DUNE_THROW(RangeError, "Differentiation order greater than 2 is not supported!");
// Evaluate 1d function values (needed for the product rule)
std::array<std::vector<R>, dim> oneDValues;
std::array<std::vector<R>, dim> oneDDerivatives;
std::array<std::vector<R>, dim> oneDSecondDerivatives;
// Evaluate 1d function derivatives
if (k==1)
for (size_t i=0; i<dim; i++)
evaluateAll(in[i], oneDValues[i], true, oneDDerivatives[i], false, oneDSecondDerivatives[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
else
for (size_t i=0; i<dim; i++)
evaluateAll(in[i], oneDValues[i], true, oneDDerivatives[i], true, oneDSecondDerivatives[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
// The lowest knot spans that we need values from
std::array<unsigned int, dim> offset;
for (int i=0; i<dim; i++)
offset[i] = std::max((int)(currentKnotSpan[i] - order_[i]),0);
// 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++)
{
// In a proper implementation, the following line would do
//limits[i] = oneDValues[i].size();
limits[i] = order_[i]+1; // The 'standard' value away from the boundaries of the knot vector
if (currentKnotSpan[i]<order_[i])
limits[i] -= (order_[i] - currentKnotSpan[i]);
if ( order_[i] > (knotVectors_[i].size() - currentKnotSpan[i] - 2) )
limits[i] -= order_[i] - (knotVectors_[i].size() - currentKnotSpan[i] - 2);
}
// Working towards computing only the parts that we really need:
// Let's copy them out into a separate array
std::array<std::vector<R>, dim> oneDValuesShort;
for (int i=0; i<dim; i++)
{
oneDValuesShort[i].resize(limits[i]);
for (size_t j=0; j<limits[i]; j++)
oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
}
MultiDigitCounter ijkCounter(limits);
out.resize(ijkCounter.cycle());
if (k == 1)
{
// Complete Jacobian is given by the product rule
for (size_t i=0; i<out.size(); i++, ++ijkCounter)
{
out[i][0] = 1.0;
for (int l=0; l<dim; l++)
out[i][0] *= (directions[0]==l) ? oneDDerivatives[l][ijkCounter[l]]
: oneDValuesShort[l][ijkCounter[l]];
}
}
if (k == 2)
{
// Complete derivation by deriving the tensor product
for (size_t i=0; i<out.size(); i++, ++ijkCounter)
{
out[i][0] = 1.0;
for (int j=0; j<dim; j++)
{
if (directions[0] != directions[1]) //derivation in two different variables
if (directions[0] == j || directions[1] == j) //the spline has to be derived (once) in this direction
out[i][0] *= oneDDerivatives[j][ijkCounter[j]];
else //no derivation in this direction
out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
else //spline is derived two times in the same direction
if (directions[0] == j) //the spline is derived two times in this direction
out[i][0] *= oneDSecondDerivatives[j][ijkCounter[j]];
else //no derivation in this direction
out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
}
}
}
}
static std::array<unsigned int,dim> getIJK(typename GridView::IndexSet::IndexType idx, std::array<unsigned int,dim> elements)
{
std::array<unsigned,dim> result;
for (int i=0; i<dim; i++)
{
result[i] = idx%elements[i];
idx /= elements[i];
}
return result;
}
static void evaluateFunction (const typename GV::ctype& in, std::vector<R>& out,
const std::vector<R>& knotVector,
unsigned int order,
unsigned int currentKnotSpan)
{
std::size_t outSize = order+1; // The 'standard' value away from the boundaries of the knot vector
if (currentKnotSpan<order) // Less near the left end of the knot vector
outSize -= (order - currentKnotSpan);
if ( order > (knotVector.size() - currentKnotSpan - 2) )
outSize -= order - (knotVector.size() - currentKnotSpan - 2);
out.resize(outSize);
// 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][std::max((int)(currentKnotSpan - order),0) + i];
}
}
static void evaluateFunctionFull(const typename GV::ctype& in,
DynamicMatrix<R>& out,
const std::vector<R>& knotVector,
unsigned int order,
unsigned int currentKnotSpan)
{
// It's not really a matrix that is needed here, a plain 2d array would do
DynamicMatrix<R>& N = out;
N.resize(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];
}
}
static void evaluateAll(const typename GV::ctype& in,
std::vector<R>& out,
bool evaluateJacobian, std::vector<R>& outJac,
bool evaluateHessian, std::vector<R>& outHess,
const std::vector<R>& knotVector,
unsigned int order,
unsigned int currentKnotSpan)
{
// How many shape functions to we have in each coordinate direction?
unsigned int limit;
limit = order+1; // The 'standard' value away from the boundaries of the knot vector
if (currentKnotSpan<order)
limit -= (order - currentKnotSpan);
if ( order > (knotVector.size() - currentKnotSpan - 2) )
limit -= order - (knotVector.size() - currentKnotSpan - 2);
// The lowest knot spans that we need values from
unsigned int offset;
offset = std::max((int)(currentKnotSpan - order),0);
// Evaluate 1d function values (needed for the product rule)
DynamicMatrix<R> values;
evaluateFunctionFull(in, values, knotVector, order, currentKnotSpan);
out.resize(knotVector.size()-order-1);
for (size_t j=0; j<out.size(); j++)
out[j] = values[order][j];
// Evaluate 1d function values of one order lower (needed for the derivative formula)
std::vector<R> lowOrderOneDValues;
if (order!=0)
{
lowOrderOneDValues.resize(knotVector.size()-(order-1)-1);
for (size_t j=0; j<lowOrderOneDValues.size(); j++)
lowOrderOneDValues[j] = values[order-1][j];
}
// Evaluate 1d function values of two order lower (needed for the (second) derivative formula)
std::vector<R> lowOrderTwoDValues;
if (order>1 && evaluateHessian)
{
lowOrderTwoDValues.resize(knotVector.size()-(order-2)-1);
for (size_t j=0; j<lowOrderTwoDValues.size(); j++)
lowOrderTwoDValues[j] = values[order-2][j];
}
// Evaluate 1d function derivatives
if (evaluateJacobian)
{
outJac.resize(limit);
if (order==0) // order-zero functions are piecewise constant, hence all derivatives are zero
std::fill(outJac.begin(), outJac.end(), R(0.0));
else
{
for (size_t j=offset; j<offset+limit; j++)
{
R derivativeAddend1 = lowOrderOneDValues[j] / (knotVector[j+order]-knotVector[j]);
R derivativeAddend2 = lowOrderOneDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]);
// The two previous terms may evaluate as 0/0. This is to be interpreted as 0.
if (std::isnan(derivativeAddend1))
derivativeAddend1 = 0;
if (std::isnan(derivativeAddend2))
derivativeAddend2 = 0;
outJac[j-offset] = order * ( derivativeAddend1 - derivativeAddend2 );
}
}
}
// Evaluate 1d function second derivatives
if (evaluateHessian)
{
outHess.resize(limit);
if (order<2) // order-zero functions are piecewise constant, hence all derivatives are zero
std::fill(outHess.begin(), outHess.end(), R(0.0));
else
{
for (size_t j=offset; j<offset+limit; j++)
{
assert(j+2 < lowOrderTwoDValues.size());
R derivativeAddend1 = lowOrderTwoDValues[j] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order-1]-knotVector[j]);
R derivativeAddend2 = lowOrderTwoDValues[j+1] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order]-knotVector[j+1]);
R derivativeAddend3 = lowOrderTwoDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+order]-knotVector[j+1]);
R derivativeAddend4 = lowOrderTwoDValues[j+2] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+1+order]-knotVector[j+2]);
// The two previous terms may evaluate as 0/0. This is to be interpreted as 0.
if (std::isnan(derivativeAddend1))
derivativeAddend1 = 0;
if (std::isnan(derivativeAddend2))
derivativeAddend2 = 0;
if (std::isnan(derivativeAddend3))
derivativeAddend3 = 0;
if (std::isnan(derivativeAddend4))
derivativeAddend4 = 0;
outHess[j-offset] = order * (order-1) * ( derivativeAddend1 - derivativeAddend2 -derivativeAddend3 + derivativeAddend4 );
}
}
}
}
std::array<unsigned int, dim> order_;
std::array<std::vector<double>, dim> knotVectors_;
std::array<unsigned,dim> elements_;
GridView gridView_;
};
template<typename GV>
class BSplineNode :
public LeafBasisNode
{
static const int dim = GV::dimension;
public:
using size_type = std::size_t;
using Element = typename GV::template Codim<0>::Entity;
using FiniteElement = BSplineLocalFiniteElement<GV,double>;
BSplineNode(const BSplinePreBasis<GV>* preBasis) :
preBasis_(preBasis),
finiteElement_(*preBasis)
{}
const Element& element() const
{
return element_;
}
const FiniteElement& finiteElement() const
{
return finiteElement_;
}
void bind(const Element& e)
{
element_ = e;
auto elementIndex = preBasis_->gridView().indexSet().index(e);
finiteElement_.bind(preBasis_->getIJK(elementIndex,preBasis_->elements_));
this->setSize(finiteElement_.size());
}
protected:
const BSplinePreBasis<GV>* preBasis_;
FiniteElement finiteElement_;
Element element_;
};
namespace BasisFactory {
inline auto bSpline(const std::vector<double>& knotVector,
unsigned int order,
bool makeOpen = true)
{
return [&knotVector, order, makeOpen](const auto& gridView) {
return BSplinePreBasis<std::decay_t<decltype(gridView)>>(gridView, knotVector, order, makeOpen);
};
}
} // end namespace BasisFactory
// *****************************************************************************
// This is the actual global basis implementation based on the reusable parts.
// *****************************************************************************
template<typename GV>
using BSplineBasis = DefaultGlobalBasis<BSplinePreBasis<GV> >;
} // namespace Functions
} // namespace Dune
#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
nodes.hh
defaultglobalbasis.hh
leafprebasismixin.hh
Dune::Functions::BasisFactory::bSpline
auto bSpline(const std::vector< double > &knotVector, unsigned int order, bool makeOpen=true)
Create a pre-basis factory that can create a B-spline pre-basis.
Definition bsplinebasis.hh:1289
Dune::Functions::BSplineBasis
DefaultGlobalBasis< BSplinePreBasis< GV > > BSplineBasis
A global B-spline basis.
Definition bsplinebasis.hh:1311
Dune
Definition monomialset.hh:19
Dune::Functions
Definition monomialset.hh:19
Dune::Functions::BasisFactory
Definition argyrisbasis.hh:926
Dune::Functions::BSplineLocalFiniteElement
LocalFiniteElement in the sense of dune-localfunctions, for the B-spline basis on tensor-product grid...
Definition bsplinebasis.hh:369
Dune::Functions::BSplineLocalFiniteElement::BSplineLocalFiniteElement
BSplineLocalFiniteElement(const BSplineLocalFiniteElement &other)
Copy constructor.
Definition bsplinebasis.hh:390
Dune::Functions::BSplineLocalFiniteElement< GV, double >::preBasis_
const BSplinePreBasis< GV > & preBasis_
Definition bsplinebasis.hh:481
Dune::Functions::BSplineLocalFiniteElement::localInterpolation
const BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > & localInterpolation() const
Hand out a LocalInterpolation object.
Definition bsplinebasis.hh:446
Dune::Functions::BSplineLocalFiniteElement::Traits
LocalFiniteElementTraits< BSplineLocalBasis< GV, R >, BSplineLocalCoefficients< dim >, BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > > Traits
Export various types related to this LocalFiniteElement.
Definition bsplinebasis.hh:379
Dune::Functions::BSplineLocalFiniteElement< GV, double >::currentKnotSpan_
std::array< unsigned, dim > currentKnotSpan_
Definition bsplinebasis.hh:488
Dune::Functions::BSplineLocalFiniteElement::BSplineLocalFiniteElement
BSplineLocalFiniteElement(const BSplinePreBasis< GV > &preBasis)
Constructor with a given B-spline basis.
Definition bsplinebasis.hh:383
Dune::Functions::BSplineLocalFiniteElement::localCoefficients
const BSplineLocalCoefficients< dim > & localCoefficients() const
Hand out a LocalCoefficients object.
Definition bsplinebasis.hh:440
Dune::Functions::BSplineLocalFiniteElement::bind
void bind(const std::array< unsigned, dim > &elementIdx)
Bind LocalFiniteElement to a specific knot span of the spline patch.
Definition bsplinebasis.hh:401
Dune::Functions::BSplineLocalFiniteElement< GV, double >::localInterpolation_
BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, double > > localInterpolation_
Definition bsplinebasis.hh:485
Dune::Functions::BSplineLocalFiniteElement::type
GeometryType type() const
Return the reference element that the local finite element is defined on (here, a hypercube).
Definition bsplinebasis.hh:462
Dune::Functions::BSplineLocalFiniteElement::size
unsigned size() const
Number of shape functions in this finite element.
Definition bsplinebasis.hh:452
Dune::Functions::BSplineLocalFiniteElement< GV, double >::localCoefficients_
BSplineLocalCoefficients< dim > localCoefficients_
Definition bsplinebasis.hh:484
Dune::Functions::BSplineLocalFiniteElement::size
unsigned int size(int i) const
Number of degrees of freedom for one coordinate direction.
Definition bsplinebasis.hh:470
Dune::Functions::BSplineLocalFiniteElement< GV, double >::localBasis_
BSplineLocalBasis< GV, double > localBasis_
Definition bsplinebasis.hh:483
Dune::Functions::BSplineLocalFiniteElement::localBasis
const BSplineLocalBasis< GV, R > & localBasis() const
Hand out a LocalBasis object.
Definition bsplinebasis.hh:434
Dune::Functions::BSplinePreBasis
Pre-basis for B-spline basis.
Definition bsplinebasis.hh:507
Dune::Functions::BSplinePreBasis::elements_
std::array< unsigned, dim > elements_
Number of grid elements in the different coordinate directions.
Definition bsplinebasis.hh:1223
Dune::Functions::BSplinePreBasis::gridView_
GridView gridView_
Definition bsplinebasis.hh:1225
Dune::Functions::BSplinePreBasis::R
double R
Definition bsplinebasis.hh:576
Dune::Functions::BSplinePreBasis::evaluateFunctionFull
static void evaluateFunctionFull(const typename GV::ctype &in, DynamicMatrix< R > &out, const std::vector< R > &knotVector, unsigned int order, unsigned int currentKnotSpan)
Evaluate all one-dimensional B-spline functions for a given coordinate direction.
Definition bsplinebasis.hh:1069
Dune::Functions::BSplinePreBasis::evaluateFunction
void evaluateFunction(const FieldVector< typename GV::ctype, dim > &in, std::vector< FieldVector< R, 1 > > &out, const std::array< unsigned, dim > &currentKnotSpan) const
Evaluate all B-spline basis functions at a given point.
Definition bsplinebasis.hh:765
Dune::Functions::BSplinePreBasis::order_
std::array< unsigned int, dim > order_
Order of the B-spline for each space dimension.
Definition bsplinebasis.hh:1217
Dune::Functions::BSplinePreBasis::evaluateAll
static void evaluateAll(const typename GV::ctype &in, std::vector< R > &out, bool evaluateJacobian, std::vector< R > &outJac, bool evaluateHessian, std::vector< R > &outHess, const std::vector< R > &knotVector, unsigned int order, unsigned int currentKnotSpan)
Evaluate the second derivatives of all one-dimensional B-spline functions for a given coordinate dire...
Definition bsplinebasis.hh:1111
Dune::Functions::BSplinePreBasis::evaluateFunction
static void evaluateFunction(const typename GV::ctype &in, std::vector< R > &out, const std::vector< R > &knotVector, unsigned int order, unsigned int currentKnotSpan)
Evaluate all one-dimensional B-spline functions for a given coordinate direction.
Definition bsplinebasis.hh:1012
Dune::Functions::BSplinePreBasis::size
unsigned int size(size_t d) const
Number of shape functions in one direction.
Definition bsplinebasis.hh:756
Dune::Functions::BSplinePreBasis::GridView
GV GridView
The grid view that the FE space is defined on.
Definition bsplinebasis.hh:570
Dune::Functions::BSplinePreBasis::indices
It indices(const Node &node, It it) const
Maps from subtree index set [0..size-1] to a globally unique multi index in global basis.
Definition bsplinebasis.hh:717
Dune::Functions::BSplinePreBasis::evaluate
void evaluate(const typename std::array< int, k > &directions, const FieldVector< typename GV::ctype, dim > &in, std::vector< FieldVector< R, 1 > > &out, const std::array< unsigned, dim > &currentKnotSpan) const
Evaluate Derivatives of all B-spline basis functions.
Definition bsplinebasis.hh:897
Dune::Functions::BSplinePreBasis::size_type
std::size_t size_type
Definition bsplinebasis.hh:571
Dune::Functions::BSplinePreBasis::update
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition bsplinebasis.hh:693
Dune::Functions::BSplinePreBasis::dimension
unsigned int dimension() const
Total number of B-spline basis functions.
Definition bsplinebasis.hh:747
Dune::Functions::BSplinePreBasis::initializeIndices
void initializeIndices()
Initialize the global indices.
Definition bsplinebasis.hh:683
Dune::Functions::BSplinePreBasis::Node
BSplineNode< GV > Node
Definition bsplinebasis.hh:573
Dune::Functions::BSplinePreBasis::gridView
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition bsplinebasis.hh:687
Dune::Functions::BSplinePreBasis::getIJK
static std::array< unsigned int, dim > getIJK(typename GridView::IndexSet::IndexType idx, std::array< unsigned int, dim > elements)
Compute integer element coordinates from the element index.
Definition bsplinebasis.hh:993
Dune::Functions::BSplinePreBasis::makeNode
Node makeNode() const
Create tree node.
Definition bsplinebasis.hh:701
Dune::Functions::BSplinePreBasis::BSplinePreBasis
BSplinePreBasis(const GridView &gridView, const std::vector< double > &knotVector, unsigned int order, bool makeOpen=true)
Construct a B-spline basis for a given grid view and set of knot vectors.
Definition bsplinebasis.hh:597
Dune::Functions::BSplinePreBasis::evaluateJacobian
void evaluateJacobian(const FieldVector< typename GV::ctype, dim > &in, std::vector< FieldMatrix< R, 1, dim > > &out, const std::array< unsigned, dim > &currentKnotSpan) const
Evaluate Jacobian of all B-spline basis functions.
Definition bsplinebasis.hh:796
Dune::Functions::BSplinePreBasis::knotVectors_
std::array< std::vector< double >, dim > knotVectors_
The knot vectors, one for each space dimension.
Definition bsplinebasis.hh:1220
Dune::Functions::BSplinePreBasis::size
size_type size() const
Get the total dimension of the space spanned by this basis.
Definition leafprebasismixin.hh:60
Dune::Functions::BSplinePreBasis::maxNodeSize
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition bsplinebasis.hh:707
Dune::Functions::BSplinePreBasis::BSplinePreBasis
BSplinePreBasis(const GridView &gridView, const FieldVector< double, dim > &lowerLeft, const FieldVector< double, dim > &upperRight, const std::array< unsigned int, dim > &elements, unsigned int order, bool makeOpen=true)
Construct a B-spline basis for a given grid view with uniform knot vectors.
Definition bsplinebasis.hh:649
Dune::Functions::BSplineLocalBasis
LocalBasis class in the sense of dune-localfunctions, presenting the restriction of a B-spline patch ...
Definition bsplinebasis.hh:51
Dune::Functions::BSplineLocalBasis::Traits
LocalBasisTraits< D, dim, FieldVector< D, dim >, R, 1, FieldVector< R, 1 >, FieldMatrix< R, 1, dim > > Traits
export type traits for function signature
Definition bsplinebasis.hh:60
Dune::Functions::BSplineLocalBasis::order
unsigned int order() const
Polynomial order of the shape functions.
Definition bsplinebasis.hh:147
Dune::Functions::BSplineLocalBasis::size
std::size_t size() const
Return the number of basis functions on the current knot span.
Definition bsplinebasis.hh:154
Dune::Functions::BSplineLocalBasis::evaluate
void evaluate(const typename std::array< int, k > &directions, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions and derivatives of any order.
Definition bsplinebasis.hh:104
Dune::Functions::BSplineLocalBasis::evaluateFunction
void evaluateFunction(const FieldVector< D, dim > &in, std::vector< FieldVector< R, 1 > > &out) const
Evaluate all shape functions.
Definition bsplinebasis.hh:76
Dune::Functions::BSplineLocalBasis::evaluateJacobian
void evaluateJacobian(const FieldVector< D, dim > &in, std::vector< FieldMatrix< D, 1, dim > > &out) const
Evaluate Jacobian of all shape functions.
Definition bsplinebasis.hh:89
Dune::Functions::BSplineLocalBasis::BSplineLocalBasis
BSplineLocalBasis(const BSplinePreBasis< GV > &preBasis, const BSplineLocalFiniteElement< GV, R > &lFE)
Constructor with a given B-spline patch.
Definition bsplinebasis.hh:66
Dune::Functions::BSplineLocalCoefficients
Attaches a shape function to an entity.
Definition bsplinebasis.hh:185
Dune::Functions::BSplineLocalCoefficients::localKey
const LocalKey & localKey(std::size_t i) const
get i'th index
Definition bsplinebasis.hh:328
Dune::Functions::BSplineLocalCoefficients::init
void init(const std::array< unsigned, dim > &sizes)
Definition bsplinebasis.hh:263
Dune::Functions::BSplineLocalCoefficients::size
std::size_t size() const
number of coefficients
Definition bsplinebasis.hh:322
Dune::Functions::BSplineLocalInterpolation
Local interpolation in the sense of dune-localfunctions, for the B-spline basis on tensor-product gri...
Definition bsplinebasis.hh:347
Dune::Functions::BSplineLocalInterpolation::interpolate
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition bsplinebasis.hh:351
Dune::Functions::BSplineNode
Definition bsplinebasis.hh:1233
Dune::Functions::BSplineNode::finiteElement
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition bsplinebasis.hh:1257
Dune::Functions::BSplineNode::Element
typename GV::template Codim< 0 >::Entity Element
Definition bsplinebasis.hh:1239
Dune::Functions::BSplineNode::preBasis_
const BSplinePreBasis< GV > * preBasis_
Definition bsplinebasis.hh:1273
Dune::Functions::BSplineNode::FiniteElement
BSplineLocalFiniteElement< GV, double > FiniteElement
Definition bsplinebasis.hh:1240
Dune::Functions::BSplineNode::element_
Element element_
Definition bsplinebasis.hh:1276
Dune::Functions::BSplineNode::bind
void bind(const Element &e)
Bind to element.
Definition bsplinebasis.hh:1263
Dune::Functions::BSplineNode::BSplineNode
BSplineNode(const BSplinePreBasis< GV > *preBasis)
Definition bsplinebasis.hh:1242
Dune::Functions::BSplineNode::element
const Element & element() const
Return current element, throw if unbound.
Definition bsplinebasis.hh:1248
Dune::Functions::BSplineNode::finiteElement_
FiniteElement finiteElement_
Definition bsplinebasis.hh:1275
Dune::Functions::BSplineNode::size_type
std::size_t size_type
Definition bsplinebasis.hh:1238
Dune::Functions::DefaultGlobalBasis
Global basis for given pre-basis.
Definition defaultglobalbasis.hh:53
Dune::Functions::LeafPreBasisMixin
A generic MixIn class for PreBasis.
Definition leafprebasismixin.hh:36
Dune::Functions::LeafPreBasisMixin< BSplinePreBasis< GV > >::size
size_type size() const
Definition leafprebasismixin.hh:60
Dune::Functions::BasisNodeMixin::size
size_type size() const
Obtain the number of basis function in the local node.
Definition nodes.hh:162
Dune::Functions::BasisNodeMixin::setSize
void setSize(const size_type size)
Definition nodes.hh:197
Dune::Functions::LeafBasisNode
Definition nodes.hh:218