3 #ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
4 #define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
14 #include <dune/common/dynmatrix.hh>
16 #include <dune/localfunctions/common/localbasis.hh>
17 #include <dune/common/diagonalmatrix.hh>
18 #include <dune/localfunctions/common/localkey.hh>
19 #include <dune/localfunctions/common/localfiniteelementtraits.hh>
20 #include <dune/geometry/type.hh>
30 template<
typename GV,
typename R,
typename MI>
33 template<
typename GV,
class MI>
45 template<
class GV,
class R,
class MI>
50 typedef typename GV::ctype D;
51 enum {dim = GV::dimension};
55 typedef LocalBasisTraits<D,dim,FieldVector<D,dim>,R,1,FieldVector<R,1>,
64 : preBasis_(preBasis),
72 std::vector<FieldVector<R,1> >& out)
const
74 FieldVector<D,dim> globalIn = offset_;
75 scaling_.umv(in,globalIn);
77 preBasis_.evaluateFunction(globalIn, out, lFE_.currentKnotSpan_);
84 std::vector<FieldMatrix<D,1,dim> >& out)
const
86 FieldVector<D,dim> globalIn = offset_;
87 scaling_.umv(in,globalIn);
89 preBasis_.evaluateJacobian(globalIn, out, lFE_.currentKnotSpan_);
91 for (
size_t i=0; i<out.size(); i++)
92 for (
int j=0; j<dim; j++)
93 out[i][0][j] *= scaling_[j][j];
98 inline void evaluate (
const typename std::array<int,k>& directions,
99 const typename Traits::DomainType& in,
100 std::vector<typename Traits::RangeType>& out)
const
109 FieldVector<D,dim> globalIn = offset_;
110 scaling_.umv(in,globalIn);
112 preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
114 for (
size_t i=0; i<out.size(); i++)
115 out[i][0] *= scaling_[directions[0]][directions[0]];
120 FieldVector<D,dim> globalIn = offset_;
121 scaling_.umv(in,globalIn);
123 preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
125 for (
size_t i=0; i<out.size(); i++)
126 out[i][0] *= scaling_[directions[0]][directions[0]]*scaling_[directions[1]][directions[1]];
130 DUNE_THROW(NotImplemented,
"B-Spline derivatives of order " << k <<
" not implemented yet!");
143 return *std::max_element(preBasis_.order_.begin(), preBasis_.order_.end());
160 FieldVector<D,dim> offset_;
161 DiagonalMatrix<D,dim> scaling_;
181 std::array<unsigned int,dim> multiindex (
unsigned int i)
const
183 std::array<unsigned int,dim> alpha;
184 for (
int j=0; j<dim; j++)
186 alpha[j] = i % sizes_[j];
193 void setup1d(std::vector<unsigned int>& subEntity)
204 unsigned lastIndex=0;
205 subEntity[lastIndex++] = 0;
206 for (
unsigned i = 0; i < sizes_[0] - 2; ++i)
207 subEntity[lastIndex++] = 0;
209 subEntity[lastIndex++] = 1;
211 assert(
size()==lastIndex);
214 void setup2d(std::vector<unsigned int>& subEntity)
216 unsigned lastIndex=0;
230 subEntity[lastIndex++] = 0;
231 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
232 subEntity[lastIndex++] = 2;
234 subEntity[lastIndex++] = 1;
237 for (
unsigned e = 0; e < sizes_[1]-2; ++e)
239 subEntity[lastIndex++] = 0;
240 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
241 subEntity[lastIndex++] = 0;
242 subEntity[lastIndex++] = 1;
246 subEntity[lastIndex++] = 2;
247 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
248 subEntity[lastIndex++] = 3;
250 subEntity[lastIndex++] = 3;
252 assert(
size()==lastIndex);
257 void init(
const std::array<unsigned,dim>& sizes)
264 std::vector<unsigned int> codim(li_.size());
266 for (std::size_t i=0; i<codim.size(); i++)
271 std::array<unsigned int,dim> mIdx = multiindex(i);
272 for (
int j=0; j<dim; j++)
273 if (mIdx[j]==0 or mIdx[j]==sizes[j]-1)
282 std::vector<unsigned int> index(
size());
284 for (std::size_t i=0; i<index.size(); i++)
288 std::array<unsigned int,dim> mIdx = multiindex(i);
290 for (
int j=dim-1; j>=0; j--)
291 if (mIdx[j]>0 and mIdx[j]<sizes[j]-1)
292 index[i] = (sizes[j]-1)*index[i] + (mIdx[j]-1);
296 std::vector<unsigned int> subEntity(li_.size());
298 if (subEntity.size() > 0)
304 }
else if (dim==2 and sizes_[0]>1 and sizes_[1]>1) {
311 for (
size_t i=0; i<li_.size(); i++)
312 li_[i] = LocalKey(subEntity[i], codim[i], index[i]);
318 return std::accumulate(sizes_.begin(), sizes_.end(), 1, std::multiplies<unsigned int>());
330 std::array<unsigned, dim> sizes_;
332 std::vector<LocalKey> li_;
339 template<
int dim,
class LB>
344 template<
typename F,
typename C>
347 DUNE_THROW(NotImplemented,
"BSplineLocalInterpolation::interpolate");
362 template<
class GV,
class R,
class MI>
365 typedef typename GV::ctype D;
366 enum {dim = GV::dimension};
372 typedef LocalFiniteElementTraits<BSplineLocalBasis<GV,R,MI>,
389 void bind(
const std::array<unsigned,dim>& elementIdx)
392 for (
size_t i=0; i<elementIdx.size(); i++)
400 for (
size_t j=0; j<elementIdx[i]; j++)
415 std::array<unsigned int, dim> sizes;
416 for (
size_t i=0; i<dim; i++)
443 for (
int i=0; i<dim; i++)
452 return GeometryTypes::cube(dim);
461 unsigned int r = order[i]+1;
480 template<
typename GV,
typename MI>
483 template<
typename GV,
class MI>
496 template<
typename GV,
class MI>
499 static const int dim = GV::dimension;
502 class MultiDigitCounter
509 MultiDigitCounter(
const std::array<unsigned int,dim>& limits)
512 std::fill(counter_.begin(), counter_.end(), 0);
516 MultiDigitCounter& operator++()
518 for (
int i=0; i<dim; i++)
523 if (counter_[i] < limits_[i])
532 const unsigned int& operator[](
int i)
const
538 unsigned int cycle()
const
541 for (
int i=0; i<dim; i++)
549 const std::array<unsigned int,dim> limits_;
552 std::array<unsigned int,dim> counter_;
593 const std::vector<double>& knotVector,
595 bool makeOpen =
true)
605 for (
int i=0; i<dim; i++)
610 for (
unsigned int j=0; j<order; j++)
616 for (
unsigned int j=0; j<order; j++)
645 const FieldVector<double,dim>& lowerLeft,
646 const FieldVector<double,dim>& upperRight,
647 const std::array<unsigned int,dim>& elements,
649 bool makeOpen =
true)
657 for (
int i=0; i<dim; i++)
662 for (
unsigned int j=0; j<order; j++)
666 for (
size_t j=0; j<elements[i]+1; j++)
667 knotVectors_[i].push_back(lowerLeft[i] + j*(upperRight[i]-lowerLeft[i]) / elements[i]);
670 for (
unsigned int j=0; j<order; j++)
715 assert(prefix.size() == 0 || prefix.size() == 1);
716 return (prefix.size() == 0) ?
size() : 0;
729 for (
int i=0; i<dim; i++)
737 unsigned int result = 1;
738 for (
size_t i=0; i<dim; i++)
744 unsigned int size (
size_t d)
const
752 std::vector<FieldVector<R,1> >& out,
753 const std::array<unsigned,dim>& currentKnotSpan)
const
756 std::array<std::vector<R>, dim> oneDValues;
758 for (
size_t i=0; i<dim; i++)
761 std::array<unsigned int, dim> limits;
762 for (
int i=0; i<dim; i++)
763 limits[i] = oneDValues[i].
size();
765 MultiDigitCounter ijkCounter(limits);
767 out.resize(ijkCounter.cycle());
769 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
772 for (
size_t j=0; j<dim; j++)
773 out[i] *= oneDValues[j][ijkCounter[j]];
783 std::vector<FieldMatrix<R,1,dim> >& out,
784 const std::array<unsigned,dim>& currentKnotSpan)
const
787 std::array<unsigned int, dim> limits;
788 for (
int i=0; i<dim; i++)
791 if (currentKnotSpan[i]<
order_[i])
792 limits[i] -= (
order_[i] - currentKnotSpan[i]);
798 std::array<unsigned int, dim> offset;
799 for (
int i=0; i<dim; i++)
800 offset[i] = std::max((
int)(currentKnotSpan[i] -
order_[i]),0);
803 std::array<std::vector<R>, dim> oneDValues;
806 std::array<std::vector<R>, dim> lowOrderOneDValues;
808 std::array<DynamicMatrix<R>, dim> values;
810 for (
size_t i=0; i<dim; i++)
814 for (
size_t j=0; j<oneDValues[i].size(); j++)
815 oneDValues[i][j] = values[i][
order_[i]][j];
820 for (
size_t j=0; j<lowOrderOneDValues[i].size(); j++)
821 lowOrderOneDValues[i][j] = values[i][
order_[i]-1][j];
827 std::array<std::vector<R>, dim> oneDDerivatives;
828 for (
size_t i=0; i<dim; i++)
830 oneDDerivatives[i].resize(limits[i]);
833 std::fill(oneDDerivatives[i].begin(), oneDDerivatives[i].end(),
R(0.0));
836 for (
size_t j=offset[i]; j<offset[i]+limits[i]; j++)
841 if (std::isnan(derivativeAddend1))
842 derivativeAddend1 = 0;
843 if (std::isnan(derivativeAddend2))
844 derivativeAddend2 = 0;
845 oneDDerivatives[i][j-offset[i]] =
order_[i] * ( derivativeAddend1 - derivativeAddend2 );
852 std::array<std::vector<R>, dim> oneDValuesShort;
854 for (
int i=0; i<dim; i++)
856 oneDValuesShort[i].resize(limits[i]);
858 for (
size_t j=0; j<limits[i]; j++)
859 oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
865 MultiDigitCounter ijkCounter(limits);
867 out.resize(ijkCounter.cycle());
870 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
871 for (
int j=0; j<dim; j++)
874 for (
int k=0; k<dim; k++)
875 out[i][0][j] *= (j==k) ? oneDDerivatives[k][ijkCounter[k]]
876 : oneDValuesShort[k][ijkCounter[k]];
882 template <
size_type k>
883 void evaluate(
const typename std::array<int,k>& directions,
884 const FieldVector<typename GV::ctype,dim>& in,
885 std::vector<FieldVector<R,1> >& out,
886 const std::array<unsigned,dim>& currentKnotSpan)
const
888 if (k != 1 && k != 2)
889 DUNE_THROW(RangeError,
"Differentiation order greater than 2 is not supported!");
892 std::array<std::vector<R>, dim> oneDValues;
893 std::array<std::vector<R>, dim> oneDDerivatives;
894 std::array<std::vector<R>, dim> oneDSecondDerivatives;
898 for (
size_t i=0; i<dim; i++)
899 evaluateAll(in[i], oneDValues[i],
true, oneDDerivatives[i],
false, oneDSecondDerivatives[i],
knotVectors_[i],
order_[i], currentKnotSpan[i]);
901 for (
size_t i=0; i<dim; i++)
905 std::array<unsigned int, dim> offset;
906 for (
int i=0; i<dim; i++)
907 offset[i] = std::max((
int)(currentKnotSpan[i] -
order_[i]),0);
910 std::array<unsigned int, dim> limits;
911 for (
int i=0; i<dim; i++)
916 if (currentKnotSpan[i]<
order_[i])
917 limits[i] -= (
order_[i] - currentKnotSpan[i]);
924 std::array<std::vector<R>, dim> oneDValuesShort;
926 for (
int i=0; i<dim; i++)
928 oneDValuesShort[i].resize(limits[i]);
930 for (
size_t j=0; j<limits[i]; j++)
931 oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
935 MultiDigitCounter ijkCounter(limits);
937 out.resize(ijkCounter.cycle());
942 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
945 for (
int l=0; l<dim; l++)
946 out[i][0] *= (directions[0]==l) ? oneDDerivatives[l][ijkCounter[l]]
947 : oneDValuesShort[l][ijkCounter[l]];
954 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
957 for (
int j=0; j<dim; j++)
959 if (directions[0] != directions[1])
960 if (directions[0] == j || directions[1] == j)
961 out[i][0] *= oneDDerivatives[j][ijkCounter[j]];
963 out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
965 if (directions[0] == j)
966 out[i][0] *= oneDSecondDerivatives[j][ijkCounter[j]];
968 out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
979 static std::array<unsigned int,dim>
getIJK(
typename GridView::IndexSet::IndexType idx, std::array<unsigned int,dim> elements)
981 std::array<unsigned,dim> result;
982 for (
int i=0; i<dim; i++)
984 result[i] = idx%elements[i];
999 const std::vector<R>& knotVector,
1001 unsigned int currentKnotSpan)
1003 std::size_t outSize = order+1;
1004 if (currentKnotSpan<order)
1005 outSize -= (order - currentKnotSpan);
1006 if ( order > (knotVector.size() - currentKnotSpan - 2) )
1007 outSize -= order - (knotVector.size() - currentKnotSpan - 2);
1008 out.resize(outSize);
1011 DynamicMatrix<R> N(order+1, knotVector.size()-1);
1019 for (
size_t i=0; i<knotVector.size()-1; i++)
1020 N[0][i] = (i == currentKnotSpan);
1022 for (
size_t r=1; r<=order; r++)
1023 for (
size_t i=0; i<knotVector.size()-r-1; i++)
1025 R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
1026 ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
1028 R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
1029 ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
1031 N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
1038 for (
size_t i=0; i<out.size(); i++) {
1039 out[i] = N[order][std::max((
int)(currentKnotSpan - order),0) + i];
1056 DynamicMatrix<R>& out,
1057 const std::vector<R>& knotVector,
1059 unsigned int currentKnotSpan)
1062 DynamicMatrix<R>& N = out;
1064 N.resize(order+1, knotVector.size()-1);
1072 for (
size_t i=0; i<knotVector.size()-1; i++)
1073 N[0][i] = (i == currentKnotSpan);
1075 for (
size_t r=1; r<=order; r++)
1076 for (
size_t i=0; i<knotVector.size()-r-1; i++)
1078 R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
1079 ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
1081 R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
1082 ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
1084 N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
1099 std::vector<R>& out,
1101 bool evaluateHessian, std::vector<R>& outHess,
1102 const std::vector<R>& knotVector,
1104 unsigned int currentKnotSpan)
1109 if (currentKnotSpan<order)
1110 limit -= (order - currentKnotSpan);
1111 if ( order > (knotVector.size() - currentKnotSpan - 2) )
1112 limit -= order - (knotVector.size() - currentKnotSpan - 2);
1115 unsigned int offset;
1116 offset = std::max((
int)(currentKnotSpan - order),0);
1119 DynamicMatrix<R> values;
1123 out.resize(knotVector.size()-order-1);
1124 for (
size_t j=0; j<out.size(); j++)
1125 out[j] = values[order][j];
1128 std::vector<R> lowOrderOneDValues;
1132 lowOrderOneDValues.resize(knotVector.size()-(order-1)-1);
1133 for (
size_t j=0; j<lowOrderOneDValues.size(); j++)
1134 lowOrderOneDValues[j] = values[order-1][j];
1138 std::vector<R> lowOrderTwoDValues;
1140 if (order>1 && evaluateHessian)
1142 lowOrderTwoDValues.resize(knotVector.size()-(order-2)-1);
1143 for (
size_t j=0; j<lowOrderTwoDValues.size(); j++)
1144 lowOrderTwoDValues[j] = values[order-2][j];
1150 outJac.resize(limit);
1153 std::fill(outJac.begin(), outJac.end(),
R(0.0));
1156 for (
size_t j=offset; j<offset+limit; j++)
1158 R derivativeAddend1 = lowOrderOneDValues[j] / (knotVector[j+order]-knotVector[j]);
1159 R derivativeAddend2 = lowOrderOneDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]);
1161 if (std::isnan(derivativeAddend1))
1162 derivativeAddend1 = 0;
1163 if (std::isnan(derivativeAddend2))
1164 derivativeAddend2 = 0;
1165 outJac[j-offset] = order * ( derivativeAddend1 - derivativeAddend2 );
1171 if (evaluateHessian)
1173 outHess.resize(limit);
1176 std::fill(outHess.begin(), outHess.end(),
R(0.0));
1179 for (
size_t j=offset; j<offset+limit; j++)
1181 assert(j+2 < lowOrderTwoDValues.size());
1182 R derivativeAddend1 = lowOrderTwoDValues[j] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order-1]-knotVector[j]);
1183 R derivativeAddend2 = lowOrderTwoDValues[j+1] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order]-knotVector[j+1]);
1184 R derivativeAddend3 = lowOrderTwoDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+order]-knotVector[j+1]);
1185 R derivativeAddend4 = lowOrderTwoDValues[j+2] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+1+order]-knotVector[j+2]);
1188 if (std::isnan(derivativeAddend1))
1189 derivativeAddend1 = 0;
1190 if (std::isnan(derivativeAddend2))
1191 derivativeAddend2 = 0;
1192 if (std::isnan(derivativeAddend3))
1193 derivativeAddend3 = 0;
1194 if (std::isnan(derivativeAddend4))
1195 derivativeAddend4 = 0;
1196 outHess[j-offset] = order * (order-1) * ( derivativeAddend1 - derivativeAddend2 -derivativeAddend3 + derivativeAddend4 );
1217 template<
typename GV,
typename MI>
1221 static const int dim = GV::dimension;
1226 using Element =
typename GV::template Codim<0>::Entity;
1253 auto elementIndex =
preBasis_->gridView().indexSet().index(e);
1268 template<
typename GV,
class MI>
1271 enum {dim = GV::dimension};
1298 for (
int i=0; i<dim; i++)
1317 template<
typename It>
1327 std::array<unsigned int,dim> globalIJK;
1328 for (
int i=0; i<dim; i++)
1329 globalIJK[i] = std::max((
int)currentKnotSpan[i] - (
int)order[i], 0) + localIJK[i];
1334 for (
int i=dim-2; i>=0; i--)
1337 *it = {{globalIdx}};
1352 namespace BasisFactory {
1356 class BSplinePreBasisFactory
1359 static const std::size_t requiredMultiIndexSize=1;
1361 BSplinePreBasisFactory(
const std::vector<double>& knotVector,
1363 bool makeOpen =
true)
1364 : knotVector_(knotVector),
1369 template<
class MultiIndex,
class Gr
idView>
1370 auto makePreBasis(
const GridView& gridView)
const
1376 const std::vector<double>& knotVector_;
1377 unsigned int order_;
1389 auto bSpline(
const std::vector<double>& knotVector,
1391 bool makeOpen =
true)
1393 return Imp::BSplinePreBasisFactory(knotVector, order, makeOpen);
1408 template<
typename GV>
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:1389
Definition: polynomial.hh:10
LocalFiniteElement in the sense of dune-localfunctions, for the B-spline basis on tensor-product grid...
Definition: bsplinebasis.hh:364
const BSplinePreBasis< GV, MI > & preBasis_
Definition: bsplinebasis.hh:469
unsigned size() const
Number of shape functions in this finite element.
Definition: bsplinebasis.hh:440
std::array< unsigned, dim > currentKnotSpan_
Definition: bsplinebasis.hh:476
GeometryType type() const
Return the reference element that the local finite element is defined on (here, a hypercube)
Definition: bsplinebasis.hh:450
BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R, MI > > localInterpolation_
Definition: bsplinebasis.hh:473
void bind(const std::array< unsigned, dim > &elementIdx)
Bind LocalFiniteElement to a specific knot span of the spline patch.
Definition: bsplinebasis.hh:389
LocalFiniteElementTraits< BSplineLocalBasis< GV, R, MI >, BSplineLocalCoefficients< dim >, BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R, MI > > > Traits
Export various types related to this LocalFiniteElement.
Definition: bsplinebasis.hh:374
const BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R, MI > > & localInterpolation() const
Hand out a LocalInterpolation object.
Definition: bsplinebasis.hh:434
const BSplineLocalBasis< GV, R, MI > & localBasis() const
Hand out a LocalBasis object.
Definition: bsplinebasis.hh:422
BSplineLocalBasis< GV, R, MI > localBasis_
Definition: bsplinebasis.hh:471
BSplineLocalFiniteElement(const BSplinePreBasis< GV, MI > &preBasis)
Constructor with a given B-spline basis.
Definition: bsplinebasis.hh:378
unsigned int size(int i) const
Number of degrees of freedom for one coordinate direction.
Definition: bsplinebasis.hh:458
const BSplineLocalCoefficients< dim > & localCoefficients() const
Hand out a LocalCoefficients object.
Definition: bsplinebasis.hh:428
BSplineLocalCoefficients< dim > localCoefficients_
Definition: bsplinebasis.hh:472
Pre-basis for B-spline basis.
Definition: bsplinebasis.hh:498
double R
Definition: bsplinebasis.hh:572
GV GridView
The grid view that the FE space is defined on.
Definition: bsplinebasis.hh:559
unsigned int size(size_t d) const
Number of shape functions in one direction.
Definition: bsplinebasis.hh:744
size_type dimension() const
Get the total dimension of the space spanned by this basis.
Definition: bsplinebasis.hh:720
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: bsplinebasis.hh:688
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:782
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:998
std::size_t size_type
Definition: bsplinebasis.hh:560
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:592
unsigned int size() const
Total number of B-spline basis functions.
Definition: bsplinebasis.hh:735
size_type size(const SizePrefix prefix) const
Return number of possible values for next position in multi index.
Definition: bsplinebasis.hh:713
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: bsplinebasis.hh:567
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition: bsplinebasis.hh:726
void initializeIndices()
Initialize the global indices.
Definition: bsplinebasis.hh:678
std::array< unsigned, dim > elements_
Number of grid elements in the different coordinate directions.
Definition: bsplinebasis.hh:1210
GridView gridView_
Definition: bsplinebasis.hh:1212
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:1098
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: bsplinebasis.hh:682
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:751
Dune::ReservedVector< size_type, 1 > SizePrefix
Definition: bsplinebasis.hh:569
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:644
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:1055
IndexSet makeIndexSet() const
Create tree node index set.
Definition: bsplinebasis.hh:707
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:883
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:979
Node makeNode() const
Create tree node.
Definition: bsplinebasis.hh:696
std::array< std::vector< double >, dim > knotVectors_
The knot vectors, one for each space dimension.
Definition: bsplinebasis.hh:1207
std::array< unsigned int, dim > order_
Order of the B-spline for each space dimension.
Definition: bsplinebasis.hh:1204
LocalBasis class in the sense of dune-localfunctions, presenting the restriction of a B-spline patch ...
Definition: bsplinebasis.hh:47
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:56
BSplineLocalBasis(const BSplinePreBasis< GV, MI > &preBasis, const BSplineLocalFiniteElement< GV, R, MI > &lFE)
Constructor with a given B-spline patch.
Definition: bsplinebasis.hh:62
void evaluateFunction(const FieldVector< D, dim > &in, std::vector< FieldVector< R, 1 > > &out) const
Evaluate all shape functions.
Definition: bsplinebasis.hh:71
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:98
unsigned int order() const
Polynomial order of the shape functions.
Definition: bsplinebasis.hh:141
void evaluateJacobian(const FieldVector< D, dim > &in, std::vector< FieldMatrix< D, 1, dim > > &out) const
Evaluate Jacobian of all shape functions.
Definition: bsplinebasis.hh:83
std::size_t size() const
Return the number of basis functions on the current knot span.
Definition: bsplinebasis.hh:148
Attaches a shape function to an entity.
Definition: bsplinebasis.hh:179
const LocalKey & localKey(std::size_t i) const
get i'th index
Definition: bsplinebasis.hh:322
void init(const std::array< unsigned, dim > &sizes)
Definition: bsplinebasis.hh:257
std::size_t size() const
number of coefficients
Definition: bsplinebasis.hh:316
Local interpolation in the sense of dune-localfunctions, for the B-spline basis on tensor-product gri...
Definition: bsplinebasis.hh:341
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: bsplinebasis.hh:345
Definition: bsplinebasis.hh:1220
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition: bsplinebasis.hh:1244
BSplineNode(const BSplinePreBasis< GV, MI > *preBasis)
Definition: bsplinebasis.hh:1229
void bind(const Element &e)
Bind to element.
Definition: bsplinebasis.hh:1250
const Element & element() const
Return current element, throw if unbound.
Definition: bsplinebasis.hh:1235
std::size_t size_type
Definition: bsplinebasis.hh:1225
const BSplinePreBasis< GV, MI > * preBasis_
Definition: bsplinebasis.hh:1260
FiniteElement finiteElement_
Definition: bsplinebasis.hh:1262
Element element_
Definition: bsplinebasis.hh:1263
typename GV::template Codim< 0 >::Entity Element
Definition: bsplinebasis.hh:1226
Definition: bsplinebasis.hh:1270
std::size_t size_type
Definition: bsplinebasis.hh:1275
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: bsplinebasis.hh:1278
BSplineNodeIndexSet(const PreBasis &preBasis)
Definition: bsplinebasis.hh:1284
void bind(const Node &node)
Bind the view to a grid element.
Definition: bsplinebasis.hh:1293
It indices(It it) const
Maps from subtree index set [0..size-1] to a globally unique multi index in global basis.
Definition: bsplinebasis.hh:1318
const Node * node_
Definition: bsplinebasis.hh:1345
std::array< unsigned int, dim > localSizes_
Definition: bsplinebasis.hh:1347
void unbind()
Unbind the view.
Definition: bsplinebasis.hh:1304
size_type size() const
Size of subtree rooted in this node (element-local)
Definition: bsplinebasis.hh:1311
const PreBasis * preBasis_
Definition: bsplinebasis.hh:1343
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:47
void setSize(const size_type size)
Definition: nodes.hh:160