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>
28 template<
typename GV,
typename R>
31 template<
typename GV,
class MI>
43 template<
class GV,
class R>
48 typedef typename GV::ctype D;
49 enum {dim = GV::dimension};
53 typedef LocalBasisTraits<D,dim,Dune::FieldVector<D,dim>,R,1,Dune::FieldVector<R,1>,
62 : preBasis_(preBasis),
70 std::vector<FieldVector<R,1> >& out)
const
72 FieldVector<D,dim> globalIn =
offset_;
73 scaling_.umv(in,globalIn);
75 preBasis_.evaluateFunction(globalIn, out, lFE_.currentKnotSpan_);
82 std::vector<FieldMatrix<D,1,dim> >& out)
const
84 FieldVector<D,dim> globalIn =
offset_;
85 scaling_.umv(in,globalIn);
87 preBasis_.evaluateJacobian(globalIn, out, lFE_.currentKnotSpan_);
89 for (
size_t i=0; i<out.size(); i++)
90 for (
int j=0; j<dim; j++)
91 out[i][0][j] *= scaling_[j][j];
96 inline void evaluate (
const typename std::array<int,k>& directions,
97 const typename Traits::DomainType& in,
98 std::vector<typename Traits::RangeType>& out)
const
107 FieldVector<D,dim> globalIn =
offset_;
108 scaling_.umv(in,globalIn);
110 preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
112 for (
size_t i=0; i<out.size(); i++)
113 out[i][0] *= scaling_[directions[0]][directions[0]];
118 FieldVector<D,dim> globalIn =
offset_;
119 scaling_.umv(in,globalIn);
121 preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
123 for (
size_t i=0; i<out.size(); i++)
124 out[i][0] *= scaling_[directions[0]][directions[0]]*scaling_[directions[1]][directions[1]];
128 DUNE_THROW(NotImplemented,
"B-Spline derivatives of order " << k <<
" not implemented yet!");
141 return *std::max_element(preBasis_.order_.begin(), preBasis_.order_.end());
159 DiagonalMatrix<D,dim> scaling_;
179 std::array<unsigned int,dim> multiindex (
unsigned int i)
const
181 std::array<unsigned int,dim> alpha;
182 for (
int j=0; j<dim; j++)
184 alpha[j] = i % sizes_[j];
191 void setup1d(std::vector<unsigned int>& subEntity)
202 unsigned lastIndex=0;
203 subEntity[lastIndex++] = 0;
204 for (
unsigned i = 0; i < sizes_[0] - 2; ++i)
205 subEntity[lastIndex++] = 0;
207 subEntity[lastIndex++] = 1;
209 assert(
size()==lastIndex);
212 void setup2d(std::vector<unsigned int>& subEntity)
214 unsigned lastIndex=0;
228 subEntity[lastIndex++] = 0;
229 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
230 subEntity[lastIndex++] = 2;
232 subEntity[lastIndex++] = 1;
235 for (
unsigned e = 0; e < sizes_[1]-2; ++e)
237 subEntity[lastIndex++] = 0;
238 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
239 subEntity[lastIndex++] = 0;
240 subEntity[lastIndex++] = 1;
244 subEntity[lastIndex++] = 2;
245 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
246 subEntity[lastIndex++] = 3;
248 subEntity[lastIndex++] = 3;
250 assert(
size()==lastIndex);
255 void init(
const std::array<unsigned,dim>& sizes)
262 std::vector<unsigned int> codim(li_.size());
264 for (std::size_t i=0; i<codim.size(); i++)
269 std::array<unsigned int,dim> mIdx = multiindex(i);
270 for (
int j=0; j<dim; j++)
271 if (mIdx[j]==0 or mIdx[j]==sizes[j]-1)
280 std::vector<unsigned int> index(
size());
282 for (std::size_t i=0; i<index.size(); i++)
286 std::array<unsigned int,dim> mIdx = multiindex(i);
288 for (
int j=dim-1; j>=0; j--)
289 if (mIdx[j]>0 and mIdx[j]<sizes[j]-1)
290 index[i] = (sizes[j]-1)*index[i] + (mIdx[j]-1);
294 std::vector<unsigned int> subEntity(li_.size());
296 if (subEntity.size() > 0)
302 }
else if (dim==2 and sizes_[0]>1 and sizes_[1]>1) {
309 for (
size_t i=0; i<li_.size(); i++)
310 li_[i] = LocalKey(subEntity[i], codim[i], index[i]);
316 return std::accumulate(sizes_.begin(), sizes_.end(), 1, std::multiplies<unsigned int>());
328 std::array<unsigned, dim> sizes_;
330 std::vector<LocalKey> li_;
337 template<
int dim,
class LB>
342 template<
typename F,
typename C>
345 DUNE_THROW(NotImplemented,
"BSplineLocalInterpolation::interpolate");
360 template<
class GV,
class R>
361 class BSplineLocalFiniteElement
363 typedef typename GV::ctype D;
364 enum {dim = GV::dimension};
370 typedef LocalFiniteElementTraits<BSplineLocalBasis<GV,R>,
387 void bind(
const std::array<unsigned,dim>& elementIdx)
390 for (
size_t i=0; i<elementIdx.size(); i++)
398 for (
size_t j=0; j<elementIdx[i]; j++)
413 std::array<unsigned int, dim> sizes;
414 for (
size_t i=0; i<dim; i++)
441 for (
int i=0; i<dim; i++)
450 return GeometryType(GeometryType::cube,dim);
459 unsigned int r = order[i]+1;
478 template<
typename GV,
typename TP>
481 template<
typename GV,
class MI,
class TP>
494 template<
typename GV,
class MI>
497 static const int dim = GV::dimension;
500 class MultiDigitCounter
507 MultiDigitCounter(
const std::array<unsigned int,dim>& limits)
510 std::fill(counter_.begin(), counter_.end(), 0);
514 MultiDigitCounter& operator++()
516 for (
int i=0; i<dim; i++)
521 if (counter_[i] < limits_[i])
530 const unsigned int& operator[](
int i)
const
536 unsigned int cycle()
const
539 for (
int i=0; i<dim; i++)
547 const std::array<unsigned int,dim> limits_;
550 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++)
727 assert(prefix.size() == 0 || prefix.size() == 1);
728 return (prefix.size() == 0) ?
size() : 0;
741 for (
int i=0; i<dim; i++)
749 unsigned int result = 1;
750 for (
size_t i=0; i<dim; i++)
756 unsigned int size (
size_t d)
const
764 std::vector<FieldVector<R,1> >& out,
765 const std::array<unsigned,dim>& currentKnotSpan)
const
768 std::array<std::vector<R>, dim> oneDValues;
770 for (
size_t i=0; i<dim; i++)
773 std::array<unsigned int, dim> limits;
774 for (
int i=0; i<dim; i++)
775 limits[i] = oneDValues[i].
size();
777 MultiDigitCounter ijkCounter(limits);
779 out.resize(ijkCounter.cycle());
781 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
784 for (
size_t j=0; j<dim; j++)
785 out[i] *= oneDValues[j][ijkCounter[j]];
795 std::vector<FieldMatrix<R,1,dim> >& out,
796 const std::array<unsigned,dim>& currentKnotSpan)
const
799 std::array<unsigned int, dim> limits;
800 for (
int i=0; i<dim; i++)
803 if (currentKnotSpan[i]<
order_[i])
804 limits[i] -= (
order_[i] - currentKnotSpan[i]);
810 std::array<unsigned int, dim> offset;
811 for (
int i=0; i<dim; i++)
812 offset[i] = std::max((
int)(currentKnotSpan[i] -
order_[i]),0);
815 std::array<std::vector<R>, dim> oneDValues;
818 std::array<std::vector<R>, dim> lowOrderOneDValues;
820 std::array<DynamicMatrix<R>, dim> values;
822 for (
size_t i=0; i<dim; i++)
826 for (
size_t j=0; j<oneDValues[i].size(); j++)
827 oneDValues[i][j] = values[i][
order_[i]][j];
832 for (
size_t j=0; j<lowOrderOneDValues[i].size(); j++)
833 lowOrderOneDValues[i][j] = values[i][
order_[i]-1][j];
839 std::array<std::vector<R>, dim> oneDDerivatives;
840 for (
size_t i=0; i<dim; i++)
842 oneDDerivatives[i].resize(limits[i]);
845 std::fill(oneDDerivatives[i].begin(), oneDDerivatives[i].end(),
R(0.0));
848 for (
size_t j=offset[i]; j<offset[i]+limits[i]; j++)
853 if (std::isnan(derivativeAddend1))
854 derivativeAddend1 = 0;
855 if (std::isnan(derivativeAddend2))
856 derivativeAddend2 = 0;
857 oneDDerivatives[i][j-offset[i]] =
order_[i] * ( derivativeAddend1 - derivativeAddend2 );
864 std::array<std::vector<R>, dim> oneDValuesShort;
866 for (
int i=0; i<dim; i++)
868 oneDValuesShort[i].resize(limits[i]);
870 for (
size_t j=0; j<limits[i]; j++)
871 oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
877 MultiDigitCounter ijkCounter(limits);
879 out.resize(ijkCounter.cycle());
882 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
883 for (
int j=0; j<dim; j++)
886 for (
int k=0; k<dim; k++)
887 out[i][0][j] *= (j==k) ? oneDDerivatives[k][ijkCounter[k]]
888 : oneDValuesShort[k][ijkCounter[k]];
894 template <
size_type k>
895 void evaluate(
const typename std::array<int,k>& directions,
896 const FieldVector<typename GV::ctype,dim>& in,
897 std::vector<FieldVector<R,1> >& out,
898 const std::array<unsigned,dim>& currentKnotSpan)
const
900 if (k != 1 && k != 2)
901 DUNE_THROW(RangeError,
"Differentiation order greater than 2 is not supported!");
904 std::array<std::vector<R>, dim> oneDValues;
905 std::array<std::vector<R>, dim> oneDDerivatives;
906 std::array<std::vector<R>, dim> oneDSecondDerivatives;
910 for (
size_t i=0; i<dim; i++)
911 evaluateAll(in[i], oneDValues[i],
true, oneDDerivatives[i],
false, oneDSecondDerivatives[i],
knotVectors_[i],
order_[i], currentKnotSpan[i]);
913 for (
size_t i=0; i<dim; i++)
917 std::array<unsigned int, dim> offset;
918 for (
int i=0; i<dim; i++)
919 offset[i] = std::max((
int)(currentKnotSpan[i] -
order_[i]),0);
922 std::array<unsigned int, dim> limits;
923 for (
int i=0; i<dim; i++)
928 if (currentKnotSpan[i]<
order_[i])
929 limits[i] -= (
order_[i] - currentKnotSpan[i]);
936 std::array<std::vector<R>, dim> oneDValuesShort;
938 for (
int i=0; i<dim; i++)
940 oneDValuesShort[i].resize(limits[i]);
942 for (
size_t j=0; j<limits[i]; j++)
943 oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
947 MultiDigitCounter ijkCounter(limits);
949 out.resize(ijkCounter.cycle());
954 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
957 for (
int l=0; l<dim; l++)
958 out[i][0] *= (directions[0]==l) ? oneDDerivatives[l][ijkCounter[l]]
959 : oneDValuesShort[l][ijkCounter[l]];
966 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
969 for (
int j=0; j<dim; j++)
971 if (directions[0] != directions[1])
972 if (directions[0] == j || directions[1] == j)
973 out[i][0] *= oneDDerivatives[j][ijkCounter[j]];
975 out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
977 if (directions[0] == j)
978 out[i][0] *= oneDSecondDerivatives[j][ijkCounter[j]];
980 out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
991 static std::array<unsigned int,dim>
getIJK(
typename GridView::IndexSet::IndexType idx, std::array<unsigned int,dim> elements)
993 std::array<unsigned,dim> result;
994 for (
int i=0; i<dim; i++)
996 result[i] = idx%elements[i];
1011 const std::vector<R>& knotVector,
1013 unsigned int currentKnotSpan)
1015 std::size_t outSize = order+1;
1016 if (currentKnotSpan<order)
1017 outSize -= (order - currentKnotSpan);
1018 if ( order > (knotVector.size() - currentKnotSpan - 2) )
1019 outSize -= order - (knotVector.size() - currentKnotSpan - 2);
1020 out.resize(outSize);
1023 DynamicMatrix<R> N(order+1, knotVector.size()-1);
1031 for (
size_t i=0; i<knotVector.size()-1; i++)
1032 N[0][i] = (i == currentKnotSpan);
1034 for (
size_t r=1; r<=order; r++)
1035 for (
size_t i=0; i<knotVector.size()-r-1; i++)
1037 R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
1038 ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
1040 R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
1041 ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
1043 N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
1050 for (
size_t i=0; i<out.size(); i++) {
1051 out[i] = N[order][std::max((
int)(currentKnotSpan - order),0) + i];
1068 DynamicMatrix<R>& out,
1069 const std::vector<R>& knotVector,
1071 unsigned int currentKnotSpan)
1074 DynamicMatrix<R>& N = out;
1076 N.resize(order+1, knotVector.size()-1);
1084 for (
size_t i=0; i<knotVector.size()-1; i++)
1085 N[0][i] = (i == currentKnotSpan);
1087 for (
size_t r=1; r<=order; r++)
1088 for (
size_t i=0; i<knotVector.size()-r-1; i++)
1090 R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
1091 ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
1093 R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
1094 ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
1096 N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
1111 std::vector<R>& out,
1113 bool evaluateHessian, std::vector<R>& outHess,
1114 const std::vector<R>& knotVector,
1116 unsigned int currentKnotSpan)
1121 if (currentKnotSpan<order)
1122 limit -= (order - currentKnotSpan);
1123 if ( order > (knotVector.size() - currentKnotSpan - 2) )
1124 limit -= order - (knotVector.size() - currentKnotSpan - 2);
1127 unsigned int offset;
1128 offset = std::max((
int)(currentKnotSpan - order),0);
1131 DynamicMatrix<R> values;
1135 out.resize(knotVector.size()-order-1);
1136 for (
size_t j=0; j<out.size(); j++)
1137 out[j] = values[order][j];
1140 std::vector<R> lowOrderOneDValues;
1144 lowOrderOneDValues.resize(knotVector.size()-(order-1)-1);
1145 for (
size_t j=0; j<lowOrderOneDValues.size(); j++)
1146 lowOrderOneDValues[j] = values[order-1][j];
1150 std::vector<R> lowOrderTwoDValues;
1152 if (order>1 && evaluateHessian)
1154 lowOrderTwoDValues.resize(knotVector.size()-(order-2)-1);
1155 for (
size_t j=0; j<lowOrderTwoDValues.size(); j++)
1156 lowOrderTwoDValues[j] = values[order-2][j];
1162 outJac.resize(limit);
1165 std::fill(outJac.begin(), outJac.end(),
R(0.0));
1168 for (
size_t j=offset; j<offset+limit; j++)
1170 R derivativeAddend1 = lowOrderOneDValues[j] / (knotVector[j+order]-knotVector[j]);
1171 R derivativeAddend2 = lowOrderOneDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]);
1173 if (std::isnan(derivativeAddend1))
1174 derivativeAddend1 = 0;
1175 if (std::isnan(derivativeAddend2))
1176 derivativeAddend2 = 0;
1177 outJac[j-offset] = order * ( derivativeAddend1 - derivativeAddend2 );
1183 if (evaluateHessian)
1185 outHess.resize(limit);
1188 std::fill(outHess.begin(), outHess.end(),
R(0.0));
1191 for (
size_t j=offset; j<offset+limit; j++)
1193 assert(j+2 < lowOrderTwoDValues.size());
1194 R derivativeAddend1 = lowOrderTwoDValues[j] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order-1]-knotVector[j]);
1195 R derivativeAddend2 = lowOrderTwoDValues[j+1] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order]-knotVector[j+1]);
1196 R derivativeAddend3 = lowOrderTwoDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+order]-knotVector[j+1]);
1197 R derivativeAddend4 = lowOrderTwoDValues[j+2] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+1+order]-knotVector[j+2]);
1200 if (std::isnan(derivativeAddend1))
1201 derivativeAddend1 = 0;
1202 if (std::isnan(derivativeAddend2))
1203 derivativeAddend2 = 0;
1204 if (std::isnan(derivativeAddend3))
1205 derivativeAddend3 = 0;
1206 if (std::isnan(derivativeAddend4))
1207 derivativeAddend4 = 0;
1208 outHess[j-offset] = order * (order-1) * ( derivativeAddend1 - derivativeAddend2 -derivativeAddend3 + derivativeAddend4 );
1229 template<
typename GV,
typename TP>
1233 static const int dim = GV::dimension;
1241 using Element =
typename GV::template Codim<0>::Entity;
1269 auto elementIndex =
preBasis_->gridView().indexSet().index(e);
1284 template<
typename GV,
class MI,
class TP>
1287 enum {dim = GV::dimension};
1314 for (
int i=0; i<dim; i++)
1329 return node_->finiteElement().size();
1333 template<
typename It>
1340 const auto currentKnotSpan =
node_->finiteElement().currentKnotSpan_;
1343 std::array<unsigned int,dim> globalIJK;
1344 for (
int i=0; i<dim; i++)
1345 globalIJK[i] = std::max((
int)currentKnotSpan[i] - (
int)order[i], 0) + localIJK[i];
1350 for (
int i=dim-2; i>=0; i--)
1353 *it = {{globalIdx}};
1378 template<
typename GV>
1386 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH