1 #ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH 2 #define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH 11 #include <dune/common/dynmatrix.hh> 13 #include <dune/localfunctions/common/localbasis.hh> 14 #include <dune/common/diagonalmatrix.hh> 15 #include <dune/localfunctions/common/localkey.hh> 16 #include <dune/localfunctions/common/localfiniteelementtraits.hh> 17 #include <dune/geometry/type.hh> 25 template<
typename GV,
typename R>
28 template<
typename GV,
class MI>
40 template<
class GV,
class R>
45 typedef typename GV::ctype D;
46 enum {dim = GV::dimension};
50 typedef LocalBasisTraits<D,dim,Dune::FieldVector<D,dim>,R,1,Dune::FieldVector<R,1>,
51 Dune::FieldMatrix<R,1,dim>, 2>
Traits;
59 : nodeFactory_(nodeFactory),
67 std::vector<FieldVector<R,1> >& out)
const 69 FieldVector<D,dim> globalIn =
offset_;
70 scaling_.umv(in,globalIn);
72 nodeFactory_.evaluateFunction(globalIn, out, lFE_.currentKnotSpan_);
79 std::vector<FieldMatrix<D,1,dim> >& out)
const 81 FieldVector<D,dim> globalIn =
offset_;
82 scaling_.umv(in,globalIn);
84 nodeFactory_.evaluateJacobian(globalIn, out, lFE_.currentKnotSpan_);
86 for (
size_t i=0; i<out.size(); i++)
87 for (
int j=0; j<dim; j++)
88 out[i][0][j] *= scaling_[j][j];
93 inline void evaluate (
const typename Dune::array<int,k>& directions,
94 const typename Traits::DomainType& in,
95 std::vector<typename Traits::RangeType>& out)
const 104 FieldVector<D,dim> globalIn =
offset_;
105 scaling_.umv(in,globalIn);
107 nodeFactory_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
109 for (
size_t i=0; i<out.size(); i++)
110 out[i][0] *= scaling_[directions[0]][directions[0]];
115 FieldVector<D,dim> globalIn =
offset_;
116 scaling_.umv(in,globalIn);
118 nodeFactory_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
120 for (
size_t i=0; i<out.size(); i++)
121 out[i][0] *= scaling_[directions[0]][directions[0]]*scaling_[directions[1]][directions[1]];
125 DUNE_THROW(NotImplemented,
"B-Spline derivatives of order " << k <<
" not implemented yet!");
138 return *std::max_element(nodeFactory_.order_.begin(), nodeFactory_.order_.end());
156 DiagonalMatrix<D,dim> scaling_;
176 std::array<unsigned int,dim> multiindex (
unsigned int i)
const 178 std::array<unsigned int,dim> alpha;
179 for (
int j=0; j<dim; j++)
181 alpha[j] = i % sizes_[j];
188 void setup1d(std::vector<unsigned int>& subEntity)
199 unsigned lastIndex=0;
200 subEntity[lastIndex++] = 0;
201 for (
unsigned i = 0; i < sizes_[0] - 2; ++i)
202 subEntity[lastIndex++] = 0;
204 subEntity[lastIndex++] = 1;
206 assert(
size()==lastIndex);
209 void setup2d(std::vector<unsigned int>& subEntity)
211 unsigned lastIndex=0;
225 subEntity[lastIndex++] = 0;
226 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
227 subEntity[lastIndex++] = 2;
229 subEntity[lastIndex++] = 1;
232 for (
unsigned e = 0; e < sizes_[1]-2; ++e)
234 subEntity[lastIndex++] = 0;
235 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
236 subEntity[lastIndex++] = 0;
237 subEntity[lastIndex++] = 1;
241 subEntity[lastIndex++] = 2;
242 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
243 subEntity[lastIndex++] = 3;
245 subEntity[lastIndex++] = 3;
247 assert(
size()==lastIndex);
252 void init(
const std::array<unsigned,dim>& sizes)
259 std::vector<unsigned int> codim(li_.size());
261 for (std::size_t i=0; i<codim.size(); i++)
266 std::array<unsigned int,dim> mIdx = multiindex(i);
267 for (
int j=0; j<dim; j++)
268 if (mIdx[j]==0 or mIdx[j]==sizes[j]-1)
277 std::vector<unsigned int> index(
size());
279 for (std::size_t i=0; i<index.size(); i++)
283 std::array<unsigned int,dim> mIdx = multiindex(i);
285 for (
int j=dim-1; j>=0; j--)
286 if (mIdx[j]>0 and mIdx[j]<sizes[j]-1)
287 index[i] = (sizes[j]-1)*index[i] + (mIdx[j]-1);
291 std::vector<unsigned int> subEntity(li_.size());
293 if (subEntity.size() > 0)
299 }
else if (dim==2 and sizes_[0]>1 and sizes_[1]>1) {
306 for (
size_t i=0; i<li_.size(); i++)
307 li_[i] = LocalKey(subEntity[i], codim[i], index[i]);
313 return std::accumulate(sizes_.begin(), sizes_.end(), 1, std::multiplies<unsigned int>());
325 std::array<unsigned, dim> sizes_;
327 std::vector<LocalKey> li_;
334 template<
int dim,
class LB>
339 template<
typename F,
typename C>
342 DUNE_THROW(NotImplemented,
"BSplineLocalInterpolation::interpolate");
357 template<
class GV,
class R>
360 typedef typename GV::ctype D;
361 enum {dim = GV::dimension};
367 typedef LocalFiniteElementTraits<BSplineLocalBasis<GV,R>,
374 : nodeFactory_(nodeFactory),
375 localBasis_(nodeFactory,*this)
384 void bind(
const std::array<uint,dim>& elementIdx)
387 for (
size_t i=0; i<elementIdx.size(); i++)
389 currentKnotSpan_[i] = 0;
392 while (nodeFactory_.knotVectors_[i][currentKnotSpan_[i]+1] < nodeFactory_.knotVectors_[i][currentKnotSpan_[i]]+1e-8)
393 currentKnotSpan_[i]++;
395 for (
size_t j=0; j<elementIdx[i]; j++)
397 currentKnotSpan_[i]++;
400 while (nodeFactory_.knotVectors_[i][currentKnotSpan_[i]+1] < nodeFactory_.knotVectors_[i][currentKnotSpan_[i]]+1e-8)
401 currentKnotSpan_[i]++;
405 localBasis_.offset_[i] = nodeFactory_.knotVectors_[i][currentKnotSpan_[i]];
406 localBasis_.scaling_[i][i] = nodeFactory_.knotVectors_[i][currentKnotSpan_[i]+1] - nodeFactory_.knotVectors_[i][currentKnotSpan_[i]];
410 std::array<unsigned int, dim> sizes;
411 for (
size_t i=0; i<dim; i++)
413 localCoefficients_.init(sizes);
425 return localCoefficients_;
431 return localInterpolation_;
438 for (
int i=0; i<dim; i++)
447 return GeometryType(GeometryType::cube,dim);
455 const auto&
order = nodeFactory_.order_;
456 unsigned int r =
order[i]+1;
457 if (currentKnotSpan_[i]<
order[i])
458 r -= (
order[i] - currentKnotSpan_[i]);
459 if (
order[i] > (nodeFactory_.knotVectors_[i].size() - currentKnotSpan_[i] - 2) )
460 r -=
order[i] - (nodeFactory_.knotVectors_[i].size() - currentKnotSpan_[i] - 2);
475 template<
typename GV,
typename TP>
478 template<
typename GV,
class MI,
class TP>
491 template<
typename GV,
class MI>
494 static const int dim = GV::dimension;
497 class MultiDigitCounter
504 MultiDigitCounter(
const std::array<unsigned int,dim>& limits)
507 std::fill(counter_.begin(), counter_.end(), 0);
511 MultiDigitCounter& operator++()
513 for (
int i=0; i<dim; i++)
518 if (counter_[i] < limits_[i])
527 const unsigned int& operator[](
int i)
const 533 unsigned int cycle()
const 536 for (
int i=0; i<dim; i++)
544 const std::array<unsigned int,dim> limits_;
547 std::array<unsigned int,dim> counter_;
590 const std::vector<double>& knotVector,
592 bool makeOpen =
true)
593 : gridView_(gridView)
596 std::fill(elements_.begin(), elements_.end(), knotVector.size()-1);
600 assert( std::accumulate(elements_.begin(), elements_.end(), 1, std::multiplies<uint>()) == gridView_.size(0) );
602 for (
int i=0; i<dim; i++)
607 for (
unsigned int j=0; j<
order; j++)
608 knotVectors_[i].push_back(knotVector[0]);
610 knotVectors_[i].insert(knotVectors_[i].end(), knotVector.begin(), knotVector.end());
613 for (
unsigned int j=0; j<
order; j++)
614 knotVectors_[i].push_back(knotVector.back());
617 std::fill(order_.begin(), order_.end(),
order);
642 const FieldVector<double,dim>& lowerLeft,
643 const FieldVector<double,dim>& upperRight,
644 const array<unsigned int,dim>& elements,
646 bool makeOpen =
true)
647 : elements_(elements),
652 assert( std::accumulate(elements_.begin(), elements_.end(), 1, std::multiplies<uint>()) == gridView_.size(0) );
654 for (
int i=0; i<dim; i++)
659 for (
unsigned int j=0; j<
order; j++)
660 knotVectors_[i].push_back(lowerLeft[i]);
663 for (
size_t j=0; j<elements[i]+1; j++)
664 knotVectors_[i].push_back(lowerLeft[i] + j*(upperRight[i]-lowerLeft[i]) / elements[i]);
667 for (
unsigned int j=0; j<
order; j++)
668 knotVectors_[i].push_back(upperRight[i]);
671 std::fill(order_.begin(), order_.end(),
order);
724 if (prefix.size() == 0)
739 for (
int i=0; i<dim; i++)
740 result *= order_[i]+1;
747 unsigned int result = 1;
748 for (
size_t i=0; i<dim; i++)
754 unsigned int size (
size_t d)
const 756 return knotVectors_[d].size() - order_[d] - 1;
762 std::vector<FieldVector<R,1> >& out,
763 const std::array<uint,dim>& currentKnotSpan)
const 766 Dune::array<std::vector<R>, dim> oneDValues;
768 for (
size_t i=0; i<dim; i++)
769 evaluateFunction(in[i], oneDValues[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
771 std::array<unsigned int, dim> limits;
772 for (
int i=0; i<dim; i++)
773 limits[i] = oneDValues[i].
size();
775 MultiDigitCounter ijkCounter(limits);
777 out.resize(ijkCounter.cycle());
779 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
782 for (
size_t j=0; j<dim; j++)
783 out[i] *= oneDValues[j][ijkCounter[j]];
793 std::vector<FieldMatrix<R,1,dim> >& out,
794 const std::array<uint,dim>& currentKnotSpan)
const 797 std::array<unsigned int, dim> limits;
798 for (
int i=0; i<dim; i++)
800 limits[i] = order_[i]+1;
801 if (currentKnotSpan[i]<order_[i])
802 limits[i] -= (order_[i] - currentKnotSpan[i]);
803 if ( order_[i] > (knotVectors_[i].
size() - currentKnotSpan[i] - 2) )
804 limits[i] -= order_[i] - (knotVectors_[i].
size() - currentKnotSpan[i] - 2);
808 std::array<unsigned int, dim> offset;
809 for (
int i=0; i<dim; i++)
810 offset[i] = std::max((
int)(currentKnotSpan[i] - order_[i]),0);
813 Dune::array<std::vector<R>, dim> oneDValues;
816 Dune::array<std::vector<R>, dim> lowOrderOneDValues;
818 Dune::array<DynamicMatrix<R>, dim> values;
820 for (
size_t i=0; i<dim; i++)
822 evaluateFunctionFull(in[i], values[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
823 oneDValues[i].resize(knotVectors_[i].
size()-order_[i]-1);
824 for (
size_t j=0; j<oneDValues[i].size(); j++)
825 oneDValues[i][j] = values[i][order_[i]][j];
829 lowOrderOneDValues[i].resize(knotVectors_[i].
size()-(order_[i]-1)-1);
830 for (
size_t j=0; j<lowOrderOneDValues[i].size(); j++)
831 lowOrderOneDValues[i][j] = values[i][order_[i]-1][j];
837 Dune::array<std::vector<R>, dim> oneDDerivatives;
838 for (
size_t i=0; i<dim; i++)
840 oneDDerivatives[i].resize(limits[i]);
843 std::fill(oneDDerivatives[i].begin(), oneDDerivatives[i].end(),
R(0.0));
846 for (
size_t j=offset[i]; j<offset[i]+limits[i]; j++)
848 R derivativeAddend1 = lowOrderOneDValues[i][j] / (knotVectors_[i][j+order_[i]]-knotVectors_[i][j]);
849 R derivativeAddend2 = lowOrderOneDValues[i][j+1] / (knotVectors_[i][j+order_[i]+1]-knotVectors_[i][j+1]);
851 if (std::isnan(derivativeAddend1))
852 derivativeAddend1 = 0;
853 if (std::isnan(derivativeAddend2))
854 derivativeAddend2 = 0;
855 oneDDerivatives[i][j-offset[i]] = order_[i] * ( derivativeAddend1 - derivativeAddend2 );
862 Dune::array<std::vector<R>, dim> oneDValuesShort;
864 for (
int i=0; i<dim; i++)
866 oneDValuesShort[i].resize(limits[i]);
868 for (
size_t j=0; j<limits[i]; j++)
869 oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
875 MultiDigitCounter ijkCounter(limits);
877 out.resize(ijkCounter.cycle());
880 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
881 for (
int j=0; j<dim; j++)
884 for (
int k=0; k<dim; k++)
885 out[i][0][j] *= (j==k) ? oneDDerivatives[k][ijkCounter[k]]
886 : oneDValuesShort[k][ijkCounter[k]];
892 template <
size_type k>
893 void evaluate(
const typename std::array<int,k>& directions,
894 const FieldVector<typename GV::ctype,dim>& in,
895 std::vector<FieldVector<R,1> >& out,
896 const std::array<uint,dim>& currentKnotSpan)
const 898 if (k != 1 && k != 2)
899 DUNE_THROW(RangeError,
"Differentiation order greater than 2 is not supported!");
902 std::array<std::vector<R>, dim> oneDValues;
903 std::array<std::vector<R>, dim> oneDDerivatives;
904 std::array<std::vector<R>, dim> oneDSecondDerivatives;
908 for (
size_t i=0; i<dim; i++)
909 evaluateAll(in[i], oneDValues[i],
true, oneDDerivatives[i],
false, oneDSecondDerivatives[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
911 for (
size_t i=0; i<dim; i++)
912 evaluateAll(in[i], oneDValues[i],
true, oneDDerivatives[i],
true, oneDSecondDerivatives[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
915 std::array<unsigned int, dim> offset;
916 for (
int i=0; i<dim; i++)
917 offset[i] = std::max((
int)(currentKnotSpan[i] - order_[i]),0);
920 std::array<unsigned int, dim> limits;
921 for (
int i=0; i<dim; i++)
925 limits[i] = order_[i]+1;
926 if (currentKnotSpan[i]<order_[i])
927 limits[i] -= (order_[i] - currentKnotSpan[i]);
928 if ( order_[i] > (knotVectors_[i].
size() - currentKnotSpan[i] - 2) )
929 limits[i] -= order_[i] - (knotVectors_[i].
size() - currentKnotSpan[i] - 2);
934 std::array<std::vector<R>, dim> oneDValuesShort;
936 for (
int i=0; i<dim; i++)
938 oneDValuesShort[i].resize(limits[i]);
940 for (
size_t j=0; j<limits[i]; j++)
941 oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
945 MultiDigitCounter ijkCounter(limits);
947 out.resize(ijkCounter.cycle());
952 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
955 for (
int l=0; l<dim; l++)
956 out[i][0] *= (directions[0]==l) ? oneDDerivatives[l][ijkCounter[l]]
957 : oneDValuesShort[l][ijkCounter[l]];
964 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
967 for (
int j=0; j<dim; j++)
969 if (directions[0] != directions[1])
970 if (directions[0] == j || directions[1] == j)
971 out[i][0] *= oneDDerivatives[j][ijkCounter[j]];
973 out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
975 if (directions[0] == j)
976 out[i][0] *= oneDSecondDerivatives[j][ijkCounter[j]];
978 out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
989 static std::array<unsigned int,dim>
getIJK(
typename GridView::IndexSet::IndexType idx, std::array<unsigned int,dim> elements)
991 std::array<uint,dim> result;
992 for (
int i=0; i<dim; i++)
994 result[i] = idx%elements[i];
1009 const std::vector<R>& knotVector,
1011 unsigned int currentKnotSpan)
1013 std::size_t outSize = order+1;
1014 if (currentKnotSpan<order)
1015 outSize -= (order - currentKnotSpan);
1016 if ( order > (knotVector.size() - currentKnotSpan - 2) )
1017 outSize -= order - (knotVector.size() - currentKnotSpan - 2);
1018 out.resize(outSize);
1021 DynamicMatrix<R> N(order+1, knotVector.size()-1);
1029 for (
size_t i=0; i<knotVector.size()-1; i++)
1030 N[0][i] = (i == currentKnotSpan);
1032 for (
size_t r=1; r<=
order; r++)
1033 for (
size_t i=0; i<knotVector.size()-r-1; i++)
1035 R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
1036 ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
1038 R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
1039 ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
1041 N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
1048 for (
size_t i=0; i<out.size(); i++) {
1049 out[i] = N[
order][std::max((
int)(currentKnotSpan - order),0) + i];
1066 DynamicMatrix<R>& out,
1067 const std::vector<R>& knotVector,
1069 unsigned int currentKnotSpan)
1072 DynamicMatrix<R>& N = out;
1074 N.resize(order+1, knotVector.size()-1);
1082 for (
size_t i=0; i<knotVector.size()-1; i++)
1083 N[0][i] = (i == currentKnotSpan);
1085 for (
size_t r=1; r<=
order; r++)
1086 for (
size_t i=0; i<knotVector.size()-r-1; i++)
1088 R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
1089 ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
1091 R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
1092 ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
1094 N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
1109 std::vector<R>& out,
1111 bool evaluateHessian, std::vector<R>& outHess,
1112 const std::vector<R>& knotVector,
1114 unsigned int currentKnotSpan)
1119 if (currentKnotSpan<order)
1120 limit -= (order - currentKnotSpan);
1121 if ( order > (knotVector.size() - currentKnotSpan - 2) )
1122 limit -= order - (knotVector.size() - currentKnotSpan - 2);
1125 unsigned int offset;
1126 offset = std::max((
int)(currentKnotSpan - order),0);
1129 DynamicMatrix<R> values;
1131 evaluateFunctionFull(in, values, knotVector, order, currentKnotSpan);
1133 out.resize(knotVector.size()-order-1);
1134 for (
size_t j=0; j<out.size(); j++)
1135 out[j] = values[order][j];
1138 std::vector<R> lowOrderOneDValues;
1142 lowOrderOneDValues.resize(knotVector.size()-(order-1)-1);
1143 for (
size_t j=0; j<lowOrderOneDValues.size(); j++)
1144 lowOrderOneDValues[j] = values[order-1][j];
1148 std::vector<R> lowOrderTwoDValues;
1150 if (order>1 && evaluateHessian)
1152 lowOrderTwoDValues.resize(knotVector.size()-(order-2)-1);
1153 for (
size_t j=0; j<lowOrderTwoDValues.size(); j++)
1154 lowOrderTwoDValues[j] = values[order-2][j];
1158 if (evaluateJacobian)
1160 outJac.resize(limit);
1163 std::fill(outJac.begin(), outJac.end(),
R(0.0));
1166 for (
size_t j=offset; j<offset+limit; j++)
1168 R derivativeAddend1 = lowOrderOneDValues[j] / (knotVector[j+
order]-knotVector[j]);
1169 R derivativeAddend2 = lowOrderOneDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]);
1171 if (std::isnan(derivativeAddend1))
1172 derivativeAddend1 = 0;
1173 if (std::isnan(derivativeAddend2))
1174 derivativeAddend2 = 0;
1175 outJac[j-offset] = order * ( derivativeAddend1 - derivativeAddend2 );
1181 if (evaluateHessian)
1183 outHess.resize(limit);
1186 std::fill(outHess.begin(), outHess.end(),
R(0.0));
1189 for (
size_t j=offset; j<offset+limit; j++)
1191 assert(j+2 < lowOrderTwoDValues.size());
1192 R derivativeAddend1 = lowOrderTwoDValues[j] / (knotVector[j+
order]-knotVector[j]) / (knotVector[j+order-1]-knotVector[j]);
1193 R derivativeAddend2 = lowOrderTwoDValues[j+1] / (knotVector[j+
order]-knotVector[j]) / (knotVector[j+order]-knotVector[j+1]);
1194 R derivativeAddend3 = lowOrderTwoDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+order]-knotVector[j+1]);
1195 R derivativeAddend4 = lowOrderTwoDValues[j+2] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+1+order]-knotVector[j+2]);
1198 if (std::isnan(derivativeAddend1))
1199 derivativeAddend1 = 0;
1200 if (std::isnan(derivativeAddend2))
1201 derivativeAddend2 = 0;
1202 if (std::isnan(derivativeAddend3))
1203 derivativeAddend3 = 0;
1204 if (std::isnan(derivativeAddend4))
1205 derivativeAddend4 = 0;
1206 outHess[j-offset] = order * (order-1) * ( derivativeAddend1 - derivativeAddend2 -derivativeAddend3 + derivativeAddend4 );
1227 template<
typename GV,
typename TP>
1231 static const int dim = GV::dimension;
1239 using Element =
typename GV::template Codim<0>::Entity;
1244 nodeFactory_(nodeFactory),
1245 finiteElement_(*nodeFactory)
1260 return finiteElement_;
1267 auto elementIndex = nodeFactory_->gridView().indexSet().index(e);
1268 finiteElement_.bind(nodeFactory_->getIJK(elementIndex,nodeFactory_->elements_));
1269 this->setSize(finiteElement_.size());
1282 template<
typename GV,
class MI,
class TP>
1285 enum {dim = GV::dimension};
1299 nodeFactory_(&nodeFactory)
1312 for (
int i=0; i<dim; i++)
1313 localSizes_[i] = node_->finiteElement().size(i);
1327 return node_->finiteElement().size();
1333 std::array<unsigned int,dim> localIJK = nodeFactory_->getIJK(i, localSizes_);
1335 const auto currentKnotSpan = node_->finiteElement().currentKnotSpan_;
1336 const auto order = nodeFactory_->order_;
1338 std::array<unsigned int,dim> globalIJK;
1339 for (
int i=0; i<dim; i++)
1340 globalIJK[i] = std::max((
int)currentKnotSpan[i] - (
int)
order[i], 0) + localIJK[i];
1345 for (
int i=dim-2; i>=0; i--)
1346 globalIdx = globalIdx * nodeFactory_->size(i) + globalIJK[i];
1348 return { globalIdx };
1371 template<
typename GV>
1379 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH GridView gridView_
Definition: bsplinebasis.hh:1222
GeometryType type() const
Return the reference element that the local finite element is defined on (here, a hypercube) ...
Definition: bsplinebasis.hh:445
Local interpolation in the sense of dune-localfunctions, for the B-spline basis on tensor-product gri...
Definition: bsplinebasis.hh:335
MultiIndex index(size_type i) const
Maps from subtree index set [0..size-1] to a globally unique multi index in global basis...
Definition: bsplinebasis.hh:1331
array< std::vector< double >, dim > knotVectors_
The knot vectors, one for each space dimension.
Definition: bsplinebasis.hh:1217
const size_type offset_
Definition: nodes.hh:40
const NodeFactory * nodeFactory_
Definition: bsplinebasis.hh:1352
IndexSet< TP > indexSet() const
Create tree node index set with given root tree path.
Definition: bsplinebasis.hh:716
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: bsplinebasis.hh:679
double R
Definition: bsplinebasis.hh:569
Dune::ReservedVector< size_type, 1 > SizePrefix
Definition: bsplinebasis.hh:566
const Node * node_
Definition: bsplinebasis.hh:1354
const LocalKey & localKey(std::size_t i) const
get i'th index
Definition: bsplinebasis.hh:317
BSplineLocalBasis< GV, R > localBasis_
Definition: bsplinebasis.hh:466
uint size() const
Number of shape functions in this finite element.
Definition: bsplinebasis.hh:435
Node factory for B-spline basis.
Definition: bsplinebasis.hh:29
BSplineLocalBasis(const BSplineNodeFactory< GV, FlatMultiIndex< std::size_t >> &nodeFactory, const BSplineLocalFiniteElement< GV, R > &lFE)
Constructor with a given B-spline patch.
Definition: bsplinebasis.hh:57
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition: bsplinebasis.hh:736
void evaluateJacobian(const FieldVector< D, dim > &in, std::vector< FieldMatrix< D, 1, dim > > &out) const
Evaluate Jacobian of all shape functions.
Definition: bsplinebasis.hh:78
LocalFiniteElementTraits< BSplineLocalBasis< GV, R >, BSplineLocalCoefficients< dim >, BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > > Traits
Export various types related to this LocalFiniteElement.
Definition: bsplinebasis.hh:369
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:1065
std::size_t size() const
Return the number of basis functions on the current knot span.
Definition: bsplinebasis.hh:143
Element element_
Definition: bsplinebasis.hh:1277
GV GridView
The grid view that the FE space is defined on.
Definition: bsplinebasis.hh:554
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: bsplinebasis.hh:1292
Definition: bsplinebasis.hh:479
BSplineLocalFiniteElement(const BSplineNodeFactory< GV, FlatMultiIndex< std::size_t >> &nodeFactory)
Constructor with a given B-spline basis.
Definition: bsplinebasis.hh:373
BSplineNodeIndexSet(const NodeFactory &nodeFactory)
Definition: bsplinebasis.hh:1298
FiniteElement finiteElement_
Definition: bsplinebasis.hh:1276
Global basis for given node factory.
Definition: defaultglobalbasis.hh:42
Definition: polynomial.hh:7
Attaches a shape function to an entity.
Definition: bsplinebasis.hh:173
LocalFiniteElement in the sense of dune-localfunctions, for the B-spline basis on tensor-product grid...
Definition: bsplinebasis.hh:26
unsigned int order() const
Polynomial order of the shape functions.
Definition: bsplinebasis.hh:136
const BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > & localInterpolation() const
Hand out a LocalInterpolation object.
Definition: bsplinebasis.hh:429
const BSplineLocalCoefficients< dim > & localCoefficients() const
Hand out a LocalCoefficients object.
Definition: bsplinebasis.hh:423
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:989
array< unsigned int, dim > order_
Order of the B-spline for each space dimension.
Definition: bsplinebasis.hh:1214
BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > localInterpolation_
Definition: bsplinebasis.hh:468
std::array< unsigned int, dim > localSizes_
Definition: bsplinebasis.hh:1356
Definition: bsplinebasis.hh:476
BSplineLocalCoefficients< dim > localCoefficients_
Definition: bsplinebasis.hh:467
std::size_t size_type
Definition: bsplinebasis.hh:1289
BSplineNode(const TreePath &treePath, const BSplineNodeFactory< GV, FlatMultiIndex< std::size_t >> *nodeFactory)
Definition: bsplinebasis.hh:1242
unsigned int size(int i) const
Number of degrees of freedom for one coordinate direction.
Definition: bsplinebasis.hh:453
void unbind()
Unbind the view.
Definition: bsplinebasis.hh:1318
void bind(const Element &e)
Bind to element.
Definition: bsplinebasis.hh:1264
void evaluate(const typename Dune::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:93
std::size_t size() const
number of coefficients
Definition: bsplinebasis.hh:311
void init(const std::array< unsigned, dim > &sizes)
Definition: bsplinebasis.hh:252
typename NodeFactory::template Node< TP > Node
Definition: bsplinebasis.hh:1296
size_type size() const
Size of subtree rooted in this node (element-local)
Definition: bsplinebasis.hh:1325
const BSplineLocalBasis< GV, R > & localBasis() const
Hand out a LocalBasis object.
Definition: bsplinebasis.hh:417
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: bsplinebasis.hh:340
size_type size(const SizePrefix prefix) const
Return number of possible values for next position in multi index.
Definition: bsplinebasis.hh:722
std::array< uint, dim > elements_
Number of grid elements in the different coordinate directions.
Definition: bsplinebasis.hh:1220
Node< TP > node(const TP &tp) const
Create tree node with given root tree path.
Definition: bsplinebasis.hh:701
void bind(const std::array< uint, dim > &elementIdx)
Bind LocalFiniteElement to a specific knot span of the spline patch.
Definition: bsplinebasis.hh:384
BSplineNodeFactory(const GridView &gridView, const FieldVector< double, dim > &lowerLeft, const FieldVector< double, dim > &upperRight, const 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:641
const BSplineNodeFactory< GV, FlatMultiIndex< std::size_t > > & nodeFactory_
Definition: bsplinebasis.hh:464
std::size_t size_type
Definition: bsplinebasis.hh:1237
void bind(const Node &node)
Bind the view to a grid element.
Definition: bsplinebasis.hh:1307
const Element & element() const
Return current element, throw if unbound.
Definition: bsplinebasis.hh:1249
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:1008
BSplineNodeFactory(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:589
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< uint, dim > ¤tKnotSpan) const
Evaluate Derivatives of all B-spline basis functions.
Definition: bsplinebasis.hh:893
unsigned int size() const
Total number of B-spline basis functions.
Definition: bsplinebasis.hh:745
const BSplineNodeFactory< GV, FlatMultiIndex< std::size_t > > * nodeFactory_
Definition: bsplinebasis.hh:1274
unsigned int size(size_t d) const
Number of shape functions in one direction.
Definition: bsplinebasis.hh:754
std::array< uint, dim > currentKnotSpan_
Definition: bsplinebasis.hh:471
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: bsplinebasis.hh:685
void evaluateFunction(const FieldVector< D, dim > &in, std::vector< FieldVector< R, 1 > > &out) const
Evaluate all shape functions.
Definition: bsplinebasis.hh:66
LocalBasis class in the sense of dune-localfunctions, presenting the restriction of a B-spline patch ...
Definition: bsplinebasis.hh:41
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:1108
typename GV::template Codim< 0 >::Entity Element
Definition: bsplinebasis.hh:1239
A multi index class with only one level.
Definition: flatmultiindex.hh:23
void evaluateJacobian(const FieldVector< typename GV::ctype, dim > &in, std::vector< FieldMatrix< R, 1, dim > > &out, const std::array< uint, dim > ¤tKnotSpan) const
Evaluate Jacobian of all B-spline basis functions.
Definition: bsplinebasis.hh:792
void initializeIndices()
Initialize the global indices.
Definition: bsplinebasis.hh:675
size_type dimension() const
Get the total dimension of the space spanned by this basis.
Definition: bsplinebasis.hh:730
TP TreePath
Definition: bsplinebasis.hh:1238
void evaluateFunction(const FieldVector< typename GV::ctype, dim > &in, std::vector< FieldVector< R, 1 > > &out, const std::array< uint, dim > ¤tKnotSpan) const
Evaluate all B-spline basis functions at a given point.
Definition: bsplinebasis.hh:761
LocalBasisTraits< D, dim, Dune::FieldVector< D, dim >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, dim >, 2 > Traits
export type traits for function signature
Definition: bsplinebasis.hh:51
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition: bsplinebasis.hh:1258
std::size_t size_type
Definition: bsplinebasis.hh:555