7#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
20#include <dune/localfunctions/common/localbasis.hh>
22#include <dune/localfunctions/common/localkey.hh>
23#include <dune/localfunctions/common/localfiniteelementtraits.hh>
25#include <dune/functions/functionspacebases/nodes.hh>
26#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
27#include <dune/functions/functionspacebases/leafprebasismixin.hh>
34template<
typename GV,
typename R>
35class BSplineLocalFiniteElement;
49template<
class GV,
class R>
54 typedef typename GV::ctype D;
55 enum {dim = GV::dimension};
68 : preBasis_(preBasis),
79 scaling_.
umv(in,globalIn);
81 preBasis_.evaluateFunction(globalIn, out, lFE_.currentKnotSpan_);
91 scaling_.
umv(in,globalIn);
93 preBasis_.evaluateJacobian(globalIn, out, lFE_.currentKnotSpan_);
95 for (
size_t i=0; i<out.size(); i++)
96 for (
int j=0; j<dim; j++)
97 out[i][0][j] *= scaling_[j][j];
102 inline void evaluate (
const typename std::array<int,k>& directions,
104 std::vector<typename Traits::RangeType>& out)
const
114 scaling_.
umv(in,globalIn);
116 preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
118 for (
size_t i=0; i<out.size(); i++)
119 out[i][0] *= scaling_[directions[0]][directions[0]];
125 scaling_.
umv(in,globalIn);
127 preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
129 for (
size_t i=0; i<out.size(); i++)
130 out[i][0] *= scaling_[directions[0]][directions[0]]*scaling_[directions[1]][directions[1]];
147 return *std::max_element(preBasis_.order_.begin(), preBasis_.order_.end());
185 std::array<unsigned int,dim> multiindex (
unsigned int i)
const
187 std::array<unsigned int,dim> alpha;
188 for (
int j=0; j<dim; j++)
190 alpha[j] = i % sizes_[j];
197 void setup1d(std::vector<unsigned int>& subEntity)
208 unsigned lastIndex=0;
209 subEntity[lastIndex++] = 0;
210 for (
unsigned i = 0; i < sizes_[0] - 2; ++i)
211 subEntity[lastIndex++] = 0;
213 subEntity[lastIndex++] = 1;
215 assert(
size()==lastIndex);
218 void setup2d(std::vector<unsigned int>& subEntity)
220 unsigned lastIndex=0;
234 subEntity[lastIndex++] = 0;
235 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
236 subEntity[lastIndex++] = 2;
238 subEntity[lastIndex++] = 1;
241 for (
unsigned e = 0; e < sizes_[1]-2; ++e)
243 subEntity[lastIndex++] = 0;
244 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
245 subEntity[lastIndex++] = 0;
246 subEntity[lastIndex++] = 1;
250 subEntity[lastIndex++] = 2;
251 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
252 subEntity[lastIndex++] = 3;
254 subEntity[lastIndex++] = 3;
256 assert(
size()==lastIndex);
261 void init(
const std::array<unsigned,dim>& sizes)
268 std::vector<unsigned int> codim(li_.size());
270 for (std::size_t i=0; i<codim.size(); i++)
275 std::array<unsigned int,dim> mIdx = multiindex(i);
276 for (
int j=0; j<dim; j++)
277 if (mIdx[j]==0 or mIdx[j]==sizes[j]-1)
286 std::vector<unsigned int> index(
size());
288 for (std::size_t i=0; i<index.size(); i++)
292 std::array<unsigned int,dim> mIdx = multiindex(i);
294 for (
int j=dim-1; j>=0; j--)
295 if (mIdx[j]>0 and mIdx[j]<sizes[j]-1)
296 index[i] = (sizes[j]-1)*index[i] + (mIdx[j]-1);
300 std::vector<unsigned int> subEntity(li_.size());
302 if (subEntity.size() > 0)
308 }
else if (dim==2 and sizes_[0]>1 and sizes_[1]>1) {
315 for (
size_t i=0; i<li_.size(); i++)
316 li_[i] =
LocalKey(subEntity[i], codim[i], index[i]);
322 return std::accumulate(sizes_.begin(), sizes_.end(), 1, std::multiplies<unsigned int>());
334 std::array<unsigned, dim> sizes_;
336 std::vector<LocalKey> li_;
343template<
int dim,
class LB>
348 template<
typename F,
typename C>
365template<
class GV,
class R>
368 typedef typename GV::ctype D;
369 enum {dim = GV::dimension};
382 : preBasis_(preBasis),
383 localBasis_(preBasis,*this)
389 : preBasis_(other.preBasis_),
390 localBasis_(preBasis_,*this)
399 void bind(
const std::array<unsigned,dim>& elementIdx)
402 for (
size_t i=0; i<elementIdx.size(); i++)
404 currentKnotSpan_[i] = 0;
407 while (preBasis_.knotVectors_[i][currentKnotSpan_[i]+1] < preBasis_.knotVectors_[i][currentKnotSpan_[i]]+1e-8)
408 currentKnotSpan_[i]++;
410 for (
size_t j=0; j<elementIdx[i]; j++)
412 currentKnotSpan_[i]++;
415 while (preBasis_.knotVectors_[i][currentKnotSpan_[i]+1] < preBasis_.knotVectors_[i][currentKnotSpan_[i]]+1e-8)
416 currentKnotSpan_[i]++;
420 localBasis_.offset_[i] = preBasis_.knotVectors_[i][currentKnotSpan_[i]];
421 localBasis_.scaling_[i][i] = preBasis_.knotVectors_[i][currentKnotSpan_[i]+1] - preBasis_.knotVectors_[i][currentKnotSpan_[i]];
425 std::array<unsigned int, dim> sizes;
426 for (
size_t i=0; i<dim; i++)
428 localCoefficients_.init(sizes);
440 return localCoefficients_;
446 return localInterpolation_;
453 for (
int i=0; i<dim; i++)
470 const auto& order = preBasis_.order_;
471 unsigned int r = order[i]+1;
472 if (currentKnotSpan_[i]<order[i])
473 r -= (order[i] - currentKnotSpan_[i]);
474 if ( order[i] > (preBasis_.knotVectors_[i].size() - currentKnotSpan_[i] - 2) )
475 r -= order[i] - (preBasis_.knotVectors_[i].size() - currentKnotSpan_[i] - 2);
486 std::array<unsigned,dim> currentKnotSpan_;
508 static const int dim = GV::dimension;
511 class MultiDigitCounter
518 MultiDigitCounter(
const std::array<unsigned int,dim>& limits)
521 std::fill(counter_.begin(), counter_.end(), 0);
525 MultiDigitCounter& operator++()
527 for (
int i=0; i<dim; i++)
532 if (counter_[i] < limits_[i])
541 const unsigned int& operator[](
int i)
const
547 unsigned int cycle()
const
550 for (
int i=0; i<dim; i++)
558 const std::array<unsigned int,dim> limits_;
561 std::array<unsigned int,dim> counter_;
569 using size_type = std::size_t;
571 using Node = BSplineNode<GV>;
595 const std::vector<double>& knotVector,
597 bool makeOpen =
true)
607 for (
int i=0; i<dim; i++)
612 for (
unsigned int j=0; j<order; j++)
618 for (
unsigned int j=0; j<order; j++)
649 const std::array<unsigned int,dim>& elements,
651 bool makeOpen =
true)
659 for (
int i=0; i<dim; i++)
664 for (
unsigned int j=0; j<order; j++)
668 for (
size_t j=0; j<elements[i]+1; j++)
672 for (
unsigned int j=0; j<order; j++)
706 size_type result = 1;
707 for (
int i=0; i<dim; i++)
713 template<
typename It>
718 std::array<unsigned int, dim> localSizes;
719 for (
int i=0; i<dim; i++)
720 localSizes[i] = node.finiteElement().size(i);
721 for (size_type i = 0, end = node.size() ; i < end ; ++i, ++it)
723 std::array<unsigned int,dim> localIJK =
getIJK(i, localSizes);
725 const auto currentKnotSpan = node.finiteElement().currentKnotSpan_;
726 const auto order =
order_;
728 std::array<unsigned int,dim> globalIJK;
729 for (
int i=0; i<dim; i++)
730 globalIJK[i] =
std::max((
int)currentKnotSpan[i] - (
int)order[i], 0) + localIJK[i];
733 size_type globalIdx = globalIJK[dim-1];
735 for (
int i=dim-2; i>=0; i--)
736 globalIdx = globalIdx *
size(i) + globalIJK[i];
746 unsigned int result = 1;
747 for (
size_t i=0; i<dim; i++)
753 unsigned int size (
size_t d)
const
764 const std::array<unsigned,dim>& currentKnotSpan)
const
767 std::array<std::vector<R>, dim> oneDValues;
769 for (
size_t i=0; i<dim; i++)
772 std::array<unsigned int, dim> limits;
773 for (
int i=0; i<dim; i++)
774 limits[i] = oneDValues[i].
size();
776 MultiDigitCounter ijkCounter(limits);
778 out.resize(ijkCounter.cycle());
780 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
783 for (
size_t j=0; j<dim; j++)
784 out[i] *= oneDValues[j][ijkCounter[j]];
795 const std::array<unsigned,dim>& currentKnotSpan)
const
798 std::array<unsigned int, dim> limits;
799 for (
int i=0; i<dim; i++)
802 if (currentKnotSpan[i]<
order_[i])
803 limits[i] -= (
order_[i] - currentKnotSpan[i]);
809 std::array<unsigned int, dim> offset;
810 for (
int i=0; i<dim; i++)
814 std::array<std::vector<R>, dim> oneDValues;
817 std::array<std::vector<R>, dim> lowOrderOneDValues;
819 std::array<DynamicMatrix<R>, dim> values;
821 for (
size_t i=0; i<dim; i++)
825 for (
size_t j=0; j<oneDValues[i].size(); j++)
826 oneDValues[i][j] = values[i][
order_[i]][j];
831 for (
size_t j=0; j<lowOrderOneDValues[i].size(); j++)
832 lowOrderOneDValues[i][j] = values[i][
order_[i]-1][j];
838 std::array<std::vector<R>, dim> oneDDerivatives;
839 for (
size_t i=0; i<dim; i++)
841 oneDDerivatives[i].resize(limits[i]);
844 std::fill(oneDDerivatives[i].begin(), oneDDerivatives[i].end(), R(0.0));
847 for (
size_t j=offset[i]; j<offset[i]+limits[i]; j++)
852 if (std::isnan(derivativeAddend1))
853 derivativeAddend1 = 0;
854 if (std::isnan(derivativeAddend2))
855 derivativeAddend2 = 0;
856 oneDDerivatives[i][j-offset[i]] =
order_[i] * ( derivativeAddend1 - derivativeAddend2 );
863 std::array<std::vector<R>, dim> oneDValuesShort;
865 for (
int i=0; i<dim; i++)
867 oneDValuesShort[i].resize(limits[i]);
869 for (
size_t j=0; j<limits[i]; j++)
870 oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
876 MultiDigitCounter ijkCounter(limits);
878 out.resize(ijkCounter.cycle());
881 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
882 for (
int j=0; j<dim; j++)
885 for (
int k=0; k<dim; k++)
886 out[i][0][j] *= (j==k) ? oneDDerivatives[k][ijkCounter[k]]
887 : oneDValuesShort[k][ijkCounter[k]];
893 template <
size_type k>
894 void evaluate(
const typename std::array<int,k>& directions,
897 const std::array<unsigned,dim>& currentKnotSpan)
const
899 if (k != 1 && k != 2)
903 std::array<std::vector<R>, dim> oneDValues;
904 std::array<std::vector<R>, dim> oneDDerivatives;
905 std::array<std::vector<R>, dim> oneDSecondDerivatives;
909 for (
size_t i=0; i<dim; i++)
910 evaluateAll(in[i], oneDValues[i],
true, oneDDerivatives[i],
false, oneDSecondDerivatives[i],
knotVectors_[i],
order_[i], currentKnotSpan[i]);
912 for (
size_t i=0; i<dim; i++)
916 std::array<unsigned int, dim> offset;
917 for (
int i=0; i<dim; i++)
921 std::array<unsigned int, dim> limits;
922 for (
int i=0; i<dim; i++)
927 if (currentKnotSpan[i]<
order_[i])
928 limits[i] -= (
order_[i] - currentKnotSpan[i]);
935 std::array<std::vector<R>, dim> oneDValuesShort;
937 for (
int i=0; i<dim; i++)
939 oneDValuesShort[i].resize(limits[i]);
941 for (
size_t j=0; j<limits[i]; j++)
942 oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
946 MultiDigitCounter ijkCounter(limits);
948 out.resize(ijkCounter.cycle());
953 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
956 for (
int l=0; l<dim; l++)
957 out[i][0] *= (directions[0]==l) ? oneDDerivatives[l][ijkCounter[l]]
958 : oneDValuesShort[l][ijkCounter[l]];
965 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
968 for (
int j=0; j<dim; j++)
970 if (directions[0] != directions[1])
971 if (directions[0] == j || directions[1] == j)
972 out[i][0] *= oneDDerivatives[j][ijkCounter[j]];
974 out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
976 if (directions[0] == j)
977 out[i][0] *= oneDSecondDerivatives[j][ijkCounter[j]];
979 out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
990 static std::array<unsigned int,dim>
getIJK(
typename GridView::IndexSet::IndexType idx, std::array<unsigned int,dim> elements)
992 std::array<unsigned,dim> result;
993 for (
int i=0; i<dim; i++)
995 result[i] = idx%elements[i];
1010 const std::vector<R>& knotVector,
1012 unsigned int currentKnotSpan)
1014 std::size_t outSize = order+1;
1015 if (currentKnotSpan<order)
1016 outSize -= (order - currentKnotSpan);
1017 if ( order > (knotVector.size() - currentKnotSpan - 2) )
1018 outSize -= order - (knotVector.size() - currentKnotSpan - 2);
1019 out.resize(outSize);
1030 for (
size_t i=0; i<knotVector.size()-1; i++)
1031 N[0][i] = (i == currentKnotSpan);
1033 for (
size_t r=1; r<=order; r++)
1034 for (
size_t i=0; i<knotVector.size()-r-1; i++)
1036 R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
1037 ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
1039 R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
1040 ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
1042 N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
1049 for (
size_t i=0; i<out.size(); i++) {
1050 out[i] = N[order][
std::max((
int)(currentKnotSpan - order),0) + i];
1068 const std::vector<R>& knotVector,
1070 unsigned int currentKnotSpan)
1075 N.
resize(order+1, knotVector.size()-1);
1083 for (
size_t i=0; i<knotVector.size()-1; i++)
1084 N[0][i] = (i == currentKnotSpan);
1086 for (
size_t r=1; r<=order; r++)
1087 for (
size_t i=0; i<knotVector.size()-r-1; i++)
1089 R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
1090 ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
1092 R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
1093 ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
1095 N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
1110 std::vector<R>& out,
1112 bool evaluateHessian, std::vector<R>& outHess,
1113 const std::vector<R>& knotVector,
1115 unsigned int currentKnotSpan)
1120 if (currentKnotSpan<order)
1121 limit -= (order - currentKnotSpan);
1122 if ( order > (knotVector.size() - currentKnotSpan - 2) )
1123 limit -= order - (knotVector.size() - currentKnotSpan - 2);
1126 unsigned int offset;
1127 offset =
std::max((
int)(currentKnotSpan - order),0);
1134 out.resize(knotVector.size()-order-1);
1135 for (
size_t j=0; j<out.size(); j++)
1136 out[j] = values[order][j];
1139 std::vector<R> lowOrderOneDValues;
1143 lowOrderOneDValues.
resize(knotVector.size()-(order-1)-1);
1144 for (
size_t j=0; j<lowOrderOneDValues.size(); j++)
1145 lowOrderOneDValues[j] = values[order-1][j];
1149 std::vector<R> lowOrderTwoDValues;
1151 if (order>1 && evaluateHessian)
1153 lowOrderTwoDValues.
resize(knotVector.size()-(order-2)-1);
1154 for (
size_t j=0; j<lowOrderTwoDValues.size(); j++)
1155 lowOrderTwoDValues[j] = values[order-2][j];
1164 std::fill(outJac.begin(), outJac.end(), R(0.0));
1167 for (
size_t j=offset; j<offset+limit; j++)
1169 R derivativeAddend1 = lowOrderOneDValues[j] / (knotVector[j+order]-knotVector[j]);
1170 R derivativeAddend2 = lowOrderOneDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]);
1172 if (std::isnan(derivativeAddend1))
1173 derivativeAddend1 = 0;
1174 if (std::isnan(derivativeAddend2))
1175 derivativeAddend2 = 0;
1176 outJac[j-offset] = order * ( derivativeAddend1 - derivativeAddend2 );
1182 if (evaluateHessian)
1184 outHess.resize(limit);
1187 std::fill(outHess.begin(), outHess.end(), R(0.0));
1190 for (
size_t j=offset; j<offset+limit; j++)
1192 assert(j+2 < lowOrderTwoDValues.size());
1193 R derivativeAddend1 = lowOrderTwoDValues[j] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order-1]-knotVector[j]);
1194 R derivativeAddend2 = lowOrderTwoDValues[j+1] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order]-knotVector[j+1]);
1195 R derivativeAddend3 = lowOrderTwoDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+order]-knotVector[j+1]);
1196 R derivativeAddend4 = lowOrderTwoDValues[j+2] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+1+order]-knotVector[j+2]);
1199 if (std::isnan(derivativeAddend1))
1200 derivativeAddend1 = 0;
1201 if (std::isnan(derivativeAddend2))
1202 derivativeAddend2 = 0;
1203 if (std::isnan(derivativeAddend3))
1204 derivativeAddend3 = 0;
1205 if (std::isnan(derivativeAddend4))
1206 derivativeAddend4 = 0;
1207 outHess[j-offset] = order * (order-1) * ( derivativeAddend1 - derivativeAddend2 -derivativeAddend3 + derivativeAddend4 );
1228template<
typename GV>
1230 public LeafBasisNode
1232 static const int dim = GV::dimension;
1236 using size_type = std::size_t;
1241 preBasis_(preBasis),
1242 finiteElement_(*preBasis)
1246 const Element& element()
const
1255 const FiniteElement& finiteElement()
const
1257 return finiteElement_;
1261 void bind(
const Element& e)
1264 auto elementIndex = preBasis_->gridView().indexSet().index(e);
1265 finiteElement_.bind(preBasis_->getIJK(elementIndex,preBasis_->elements_));
1266 this->setSize(finiteElement_.size());
1271 const BSplinePreBasis<GV>* preBasis_;
1273 FiniteElement finiteElement_;
1279namespace BasisFactory {
1287inline auto bSpline(
const std::vector<double>& knotVector,
1289 bool makeOpen =
true)
1291 return [&knotVector, order, makeOpen](
const auto& gridView) {
1292 return BSplinePreBasis<std::decay_t<
decltype(gridView)>>(gridView, knotVector, order, makeOpen);
1308template<
typename GV>
Construct a matrix with a dynamic size.
Definition: dynmatrix.hh:61
void resize(size_type r, size_type c, value_type v=value_type())
resize matrix to r × c
Definition: dynmatrix.hh:106
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:92
LocalBasis class in the sense of dune-localfunctions, presenting the restriction of a B-spline patch ...
Definition: bsplinebasis.hh:51
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
unsigned int order() const
Polynomial order of the shape functions.
Definition: bsplinebasis.hh:145
std::size_t size() const
Return the number of basis functions on the current knot span.
Definition: bsplinebasis.hh:152
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:102
void evaluateFunction(const FieldVector< D, dim > &in, std::vector< FieldVector< R, 1 > > &out) const
Evaluate all shape functions.
Definition: bsplinebasis.hh:75
void evaluateJacobian(const FieldVector< D, dim > &in, std::vector< FieldMatrix< D, 1, dim > > &out) const
Evaluate Jacobian of all shape functions.
Definition: bsplinebasis.hh:87
BSplineLocalBasis(const BSplinePreBasis< GV > &preBasis, const BSplineLocalFiniteElement< GV, R > &lFE)
Constructor with a given B-spline patch.
Definition: bsplinebasis.hh:66
Attaches a shape function to an entity.
Definition: bsplinebasis.hh:183
const LocalKey & localKey(std::size_t i) const
get i'th index
Definition: bsplinebasis.hh:326
std::size_t size() const
number of coefficients
Definition: bsplinebasis.hh:320
LocalFiniteElement in the sense of dune-localfunctions, for the B-spline basis on tensor-product grid...
Definition: bsplinebasis.hh:367
BSplineLocalFiniteElement(const BSplineLocalFiniteElement &other)
Copy constructor.
Definition: bsplinebasis.hh:388
const BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > & localInterpolation() const
Hand out a LocalInterpolation object.
Definition: bsplinebasis.hh:444
LocalFiniteElementTraits< BSplineLocalBasis< GV, R >, BSplineLocalCoefficients< dim >, BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > > Traits
Export various types related to this LocalFiniteElement.
Definition: bsplinebasis.hh:377
BSplineLocalFiniteElement(const BSplinePreBasis< GV > &preBasis)
Constructor with a given B-spline basis.
Definition: bsplinebasis.hh:381
const BSplineLocalCoefficients< dim > & localCoefficients() const
Hand out a LocalCoefficients object.
Definition: bsplinebasis.hh:438
void bind(const std::array< unsigned, dim > &elementIdx)
Bind LocalFiniteElement to a specific knot span of the spline patch.
Definition: bsplinebasis.hh:399
GeometryType type() const
Return the reference element that the local finite element is defined on (here, a hypercube)
Definition: bsplinebasis.hh:460
unsigned size() const
Number of shape functions in this finite element.
Definition: bsplinebasis.hh:450
unsigned int size(int i) const
Number of degrees of freedom for one coordinate direction.
Definition: bsplinebasis.hh:468
const BSplineLocalBasis< GV, R > & localBasis() const
Hand out a LocalBasis object.
Definition: bsplinebasis.hh:432
Local interpolation in the sense of dune-localfunctions, for the B-spline basis on tensor-product gri...
Definition: bsplinebasis.hh:345
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: bsplinebasis.hh:349
Pre-basis for B-spline basis.
Definition: bsplinebasis.hh:505
std::array< unsigned, dim > elements_
Number of grid elements in the different coordinate directions.
Definition: bsplinebasis.hh:1221
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:1066
void evaluateFunction(const FieldVector< typename GV::ctype, dim > &in, std::vector< FieldVector< R, 1 > > &out, const std::array< unsigned, dim > ¤tKnotSpan) const
Evaluate all B-spline basis functions at a given point.
Definition: bsplinebasis.hh:762
std::array< unsigned int, dim > order_
Order of the B-spline for each space dimension.
Definition: bsplinebasis.hh:1215
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:1109
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:1009
unsigned int size(size_t d) const
Number of shape functions in one direction.
Definition: bsplinebasis.hh:753
GV GridView
The grid view that the FE space is defined on.
Definition: bsplinebasis.hh:568
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:714
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 > ¤tKnotSpan) const
Evaluate Derivatives of all B-spline basis functions.
Definition: bsplinebasis.hh:894
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: bsplinebasis.hh:690
unsigned int dimension() const
Total number of B-spline basis functions.
Definition: bsplinebasis.hh:744
void initializeIndices()
Initialize the global indices.
Definition: bsplinebasis.hh:680
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: bsplinebasis.hh:684
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:990
Node makeNode() const
Create tree node.
Definition: bsplinebasis.hh:698
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:594
void evaluateJacobian(const FieldVector< typename GV::ctype, dim > &in, std::vector< FieldMatrix< R, 1, dim > > &out, const std::array< unsigned, dim > ¤tKnotSpan) const
Evaluate Jacobian of all B-spline basis functions.
Definition: bsplinebasis.hh:793
std::array< std::vector< double >, dim > knotVectors_
The knot vectors, one for each space dimension.
Definition: bsplinebasis.hh:1218
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition: bsplinebasis.hh:704
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:646
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:50
A generic MixIn class for PreBasis.
Definition: leafprebasismixin.hh:36
size_type size() const
Get the total dimension of the space spanned by this basis.
Definition: leafprebasismixin.hh:60
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
Describe position of one degree of freedom.
Definition: localkey.hh:24
Default exception for dummy implementations.
Definition: exceptions.hh:355
Default exception class for range errors.
Definition: exceptions.hh:346
This file implements a quadratic diagonal matrix of fixed size.
This file implements a dense matrix with dynamic numbers of rows and columns.
void umv(const X &x, Y &y) const
y += A x
Definition: diagonalmatrix.hh:296
#define DUNE_THROW(E,...)
Definition: exceptions.hh:312
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:1287
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:462
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:485
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:280
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integer_sequence< T, II..., T(IN)> push_back(std::integer_sequence< T, II... >, std::integral_constant< T, IN >={})
Append an index IN to the back of the sequence.
Definition: integersequence.hh:69
Static tag representing a codimension.
Definition: dimension.hh:24
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:35
D DomainType
domain type
Definition: localbasis.hh:43
traits helper struct
Definition: localfiniteelementtraits.hh:13
A unique label for each type of element that can occur in a grid.