\n
\n-
3#ifndef DUNE_FUNCTIONS_GRIDFUNCTIONS_GRID_FUNCTION_IMP_HH \n-
4#define DUNE_FUNCTIONS_GRIDFUNCTIONS_GRID_FUNCTION_IMP_HH \n+
3#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH \n+
4#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH \n
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n+
\n+
\n+
\n+
14#include <dune/common/dynmatrix.hh> \n
\n-
19struct HasFreeLocalFunction
\n-
\n-
\n-
22 auto require(F&& f) ->
decltype(
\n-
\n-
\n-
\n-
\n+
16#include <dune/localfunctions/common/localbasis.hh> \n+
17#include <dune/common/diagonalmatrix.hh> \n+
18#include <dune/localfunctions/common/localkey.hh> \n+
19#include <dune/localfunctions/common/localfiniteelementtraits.hh> \n+
20#include <dune/geometry/type.hh> \n+
\n+
\n+
\n+
\n+
\n+
\n
\n-
\n-
\n-
\n-
\n-
\n-
33template<
class Signature,
class DerivativeInterface,
class LocalFunctionInterface,
class EntitySet>
\n-
34class GridFunctionWrapperInterface :
\n-
35 public DifferentiableFunctionWrapperInterface<Signature, DerivativeInterface>
\n-
\n-
\n-
38 virtual LocalFunctionInterface wrappedLocalFunction()
const = 0;
\n-
\n-
40 virtual const EntitySet& wrappedEntitySet()
const = 0;
\n-
\n-
\n-
\n-
\n-
45template<
class Signature,
class DerivativeInterface,
class LocalFunctionInterface,
class EntitySet,
class B>
\n-
46class GridFunctionWrapperImplementation :
\n-
47 public DifferentiableFunctionWrapperImplementation<Signature, DerivativeInterface, B>
\n-
\n-
49 using Base = DifferentiableFunctionWrapperImplementation<Signature, DerivativeInterface, B>;
\n-
\n-
\n+
\n+
29template<
typename GV,
typename R>
\n+
\n+
\n+
\n+
\n+
\n+
\n+
44template<
class GV,
class R>
\n+
\n+
\n+
\n+
\n+
\n+
49 typedef typename GV::ctype D;
\n+
50 enum {dim = GV::dimension};
\n+
\n
\n-
53 virtual LocalFunctionInterface wrappedLocalFunction()
const \n-
\n-
55 return localFunction(this->get());
\n-
\n-
\n-
58 virtual const EntitySet& wrappedEntitySet()
const \n-
\n-
60 return this->get().entitySet();
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n+
54 typedef LocalBasisTraits<D,dim,FieldVector<D,dim>,R,1,FieldVector<R,1>,
\n+
\n+
\n+
\n+
\n+
\n+
63 : preBasis_(preBasis),
\n+
\n+
\n+
\n+
\n+
\n+
\n+
71 std::vector<FieldVector<R,1> >& out)
const \n+
\n+
73 FieldVector<D,dim> globalIn = offset_;
\n+
74 scaling_.umv(in,globalIn);
\n+
\n+
76 preBasis_.evaluateFunction(globalIn, out, lFE_.currentKnotSpan_);
\n+
\n+
\n+
\n+
\n+
\n+
83 std::vector<FieldMatrix<D,1,dim> >& out)
const \n+
\n+
85 FieldVector<D,dim> globalIn = offset_;
\n+
86 scaling_.umv(in,globalIn);
\n+
\n+
88 preBasis_.evaluateJacobian(globalIn, out, lFE_.currentKnotSpan_);
\n+
\n+
90 for (
size_t i=0; i<out.size(); i++)
\n+
91 for (
int j=0; j<dim; j++)
\n+
92 out[i][0][j] *= scaling_[j][j];
\n+
\n+
\n+
\n+
\n+
\n+
97 inline void evaluate (
const typename std::array<int,k>& directions,
\n+
98 const typename Traits::DomainType& in,
\n+
99 std::vector<typename Traits::RangeType>& out)
const \n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
108 FieldVector<D,dim> globalIn = offset_;
\n+
109 scaling_.umv(in,globalIn);
\n+
\n+
111 preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
\n+
\n+
113 for (
size_t i=0; i<out.size(); i++)
\n+
114 out[i][0] *= scaling_[directions[0]][directions[0]];
\n+
\n+
\n+
\n+
\n+
119 FieldVector<D,dim> globalIn = offset_;
\n+
120 scaling_.umv(in,globalIn);
\n+
\n+
122 preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
\n+
\n+
124 for (
size_t i=0; i<out.size(); i++)
\n+
125 out[i][0] *= scaling_[directions[0]][directions[0]]*scaling_[directions[1]][directions[1]];
\n+
\n+
\n+
\n+
129 DUNE_THROW(NotImplemented,
"B-Spline derivatives of order " << k <<
" not implemented yet!");
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
142 return *std::max_element(preBasis_.order_.begin(), preBasis_.order_.end());
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
159 FieldVector<D,dim> offset_;
\n+
160 DiagonalMatrix<D,dim> scaling_;
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
180 std::array<unsigned int,dim> multiindex (
unsigned int i)
const \n+
\n+
182 std::array<unsigned int,dim> alpha;
\n+
183 for (
int j=0; j<dim; j++)
\n+
\n+
185 alpha[j] = i % sizes_[j];
\n+
\n+
\n+
\n+
\n+
\n+
192 void setup1d(std::vector<unsigned int>& subEntity)
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
203 unsigned lastIndex=0;
\n+
204 subEntity[lastIndex++] = 0;
\n+
205 for (
unsigned i = 0; i < sizes_[0] - 2; ++i)
\n+
206 subEntity[lastIndex++] = 0;
\n+
\n+
208 subEntity[lastIndex++] = 1;
\n+
\n+
210 assert(
size()==lastIndex);
\n+
\n+
\n+
213 void setup2d(std::vector<unsigned int>& subEntity)
\n+
\n+
215 unsigned lastIndex=0;
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
229 subEntity[lastIndex++] = 0;
\n+
230 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
\n+
231 subEntity[lastIndex++] = 2;
\n+
\n+
233 subEntity[lastIndex++] = 1;
\n+
\n+
\n+
236 for (
unsigned e = 0; e < sizes_[1]-2; ++e)
\n+
\n+
238 subEntity[lastIndex++] = 0;
\n+
239 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
\n+
240 subEntity[lastIndex++] = 0;
\n+
241 subEntity[lastIndex++] = 1;
\n+
\n+
\n+
\n+
245 subEntity[lastIndex++] = 2;
\n+
246 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
\n+
247 subEntity[lastIndex++] = 3;
\n+
\n+
249 subEntity[lastIndex++] = 3;
\n+
\n+
251 assert(
size()==lastIndex);
\n+
\n+
\n+
\n+
\n+
\n+
256 void init(
const std::array<unsigned,dim>& sizes)
\n+
\n+
\n+
\n+
\n+
\n+
\n+
263 std::vector<unsigned int> codim(li_.size());
\n+
\n+
265 for (std::size_t i=0; i<codim.size(); i++)
\n+
\n+
\n+
\n+
\n+
270 std::array<unsigned int,dim> mIdx = multiindex(i);
\n+
271 for (
int j=0; j<dim; j++)
\n+
272 if (mIdx[j]==0 or mIdx[j]==sizes[j]-1)
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
281 std::vector<unsigned int> index(
size());
\n+
\n+
283 for (std::size_t i=0; i<index.size(); i++)
\n+
\n+
\n+
\n+
287 std::array<unsigned int,dim> mIdx = multiindex(i);
\n+
\n+
289 for (
int j=dim-1; j>=0; j--)
\n+
290 if (mIdx[j]>0 and mIdx[j]<sizes[j]-1)
\n+
291 index[i] = (sizes[j]-1)*index[i] + (mIdx[j]-1);
\n+
\n+
\n+
\n+
295 std::vector<unsigned int> subEntity(li_.size());
\n+
\n+
297 if (subEntity.size() > 0)
\n+
\n+
\n+
\n+
\n+
\n+
303 }
else if (dim==2 and sizes_[0]>1 and sizes_[1]>1) {
\n+
\n+
\n+
\n+
\n+
\n+
\n+
310 for (
size_t i=0; i<li_.size(); i++)
\n+
311 li_[i] = LocalKey(subEntity[i], codim[i], index[i]);
\n+
\n+
\n+
\n+
\n+
\n+
\n+
317 return std::accumulate(sizes_.begin(), sizes_.end(), 1, std::multiplies<unsigned int>());
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
329 std::array<unsigned, dim> sizes_;
\n+
\n+
331 std::vector<LocalKey> li_;
\n+
\n+
\n+
\n+
338template<
int dim,
class LB>
\n+
\n+
\n+
\n+
\n+
343 template<
typename F,
typename C>
\n+
\n+
\n+
\n+
346 DUNE_THROW(NotImplemented,
"BSplineLocalInterpolation::interpolate");
\n+
\n+
\n+
\n+
\n+
\n+
\n+
360template<
class GV,
class R>
\n+
\n+
\n+
\n+
363 typedef typename GV::ctype D;
\n+
364 enum {dim = GV::dimension};
\n+
\n+
\n+
\n+
370 typedef LocalFiniteElementTraits<BSplineLocalBasis<GV,R>,
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
394 void bind(
const std::array<unsigned,dim>& elementIdx)
\n+
\n+
\n+
397 for (
size_t i=0; i<elementIdx.size(); i++)
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
405 for (
size_t j=0; j<elementIdx[i]; j++)
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
420 std::array<unsigned int, dim> sizes;
\n+
421 for (
size_t i=0; i<dim; i++)
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
448 for (
int i=0; i<dim; i++)
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
457 return GeometryTypes::cube(dim);
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
466 unsigned int r = order[i]+1;
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
500 static const int dim = GV::dimension;
\n+
\n+
503 class MultiDigitCounter
\n+
\n+
\n+
\n+
510 MultiDigitCounter(
const std::array<unsigned int,dim>& limits)
\n+
\n+
\n+
513 std::fill(counter_.begin(), counter_.end(), 0);
\n+
\n+
\n+
517 MultiDigitCounter& operator++()
\n+
\n+
519 for (
int i=0; i<dim; i++)
\n+
\n+
\n+
\n+
\n+
524 if (counter_[i] < limits_[i])
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
533 const unsigned int& operator[](
int i)
const \n+
\n+
\n+
\n+
\n+
539 unsigned int cycle()
const \n+
\n+
\n+
542 for (
int i=0; i<dim; i++)
\n+
\n+
\n+
\n+
\n+
\n+
\n+
550 const std::array<unsigned int,dim> limits_;
\n+
\n+
553 std::array<unsigned int,dim> counter_;
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
591 const std::vector<double>& knotVector,
\n+
\n+
593 bool makeOpen =
true)
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
603 for (
int i=0; i<dim; i++)
\n+
\n+
\n+
\n+
608 for (
unsigned int j=0; j<order; j++)
\n+
\n+
\n+
\n+
\n+
\n+
614 for (
unsigned int j=0; j<order; j++)
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
643 const FieldVector<double,dim>& lowerLeft,
\n+
644 const FieldVector<double,dim>& upperRight,
\n+
645 const std::array<unsigned int,dim>& elements,
\n+
\n+
647 bool makeOpen =
true)
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
655 for (
int i=0; i<dim; i++)
\n+
\n+
\n+
\n+
660 for (
unsigned int j=0; j<order; j++)
\n+
\n+
\n+
\n+
664 for (
size_t j=0; j<elements[i]+1; j++)
\n+
665 knotVectors_[i].push_back(lowerLeft[i] + j*(upperRight[i]-lowerLeft[i]) / elements[i]);
\n+
\n+
\n+
668 for (
unsigned int j=0; j<order; j++)
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
711 template<
class ST,
int i>
\n+
\n+
\n+
\n+
714 assert(prefix.size() == 0 || prefix.size() == 1);
\n+
715 return (prefix.size() == 0) ?
size() : 0;
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
728 for (
int i=0; i<dim; i++)
\n+
\n+
\n+
\n+
\n+
\n+
734 template<
typename It>
\n+
\n+
\n+
\n+
\n+
\n+
739 std::array<unsigned int, dim> localSizes;
\n+
740 for (
int i=0; i<dim; i++)
\n+
\n+
742 for (
size_type i = 0, end = node.
size() ; i < end ; ++i, ++it)
\n+
\n+
744 std::array<unsigned int,dim> localIJK =
getIJK(i, localSizes);
\n+
\n+
\n+
747 const auto order =
order_;
\n+
\n+
749 std::array<unsigned int,dim> globalIJK;
\n+
750 for (
int i=0; i<dim; i++)
\n+
751 globalIJK[i] = std::max((
int)currentKnotSpan[i] - (
int)order[i], 0) + localIJK[i];
\n+
\n+
\n+
\n+
\n+
756 for (
int i=dim-2; i>=0; i--)
\n+
757 globalIdx = globalIdx *
size(i) + globalIJK[i];
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
767 unsigned int result = 1;
\n+
768 for (
size_t i=0; i<dim; i++)
\n+
\n+
\n+
\n+
\n+
\n+
\n+
774 unsigned int size (
size_t d)
const \n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
782 std::vector<FieldVector<R,1> >& out,
\n+
783 const std::array<unsigned,dim>& currentKnotSpan)
const \n+
\n+
\n+
786 std::array<std::vector<R>, dim> oneDValues;
\n+
\n+
788 for (
size_t i=0; i<dim; i++)
\n+
\n+
\n+
791 std::array<unsigned int, dim> limits;
\n+
792 for (
int i=0; i<dim; i++)
\n+
793 limits[i] = oneDValues[i].
size();
\n+
\n+
795 MultiDigitCounter ijkCounter(limits);
\n+
\n+
797 out.resize(ijkCounter.cycle());
\n+
\n+
799 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
\n+
\n+
\n+
802 for (
size_t j=0; j<dim; j++)
\n+
803 out[i] *= oneDValues[j][ijkCounter[j]];
\n+
\n+
\n+
\n+
\n+
\n+
\n+
813 std::vector<FieldMatrix<R,1,dim> >& out,
\n+
814 const std::array<unsigned,dim>& currentKnotSpan)
const \n+
\n+
\n+
817 std::array<unsigned int, dim> limits;
\n+
818 for (
int i=0; i<dim; i++)
\n+
\n+
\n+
821 if (currentKnotSpan[i]<
order_[i])
\n+
822 limits[i] -= (
order_[i] - currentKnotSpan[i]);
\n+
\n+
\n+
\n+
\n+
\n+
828 std::array<unsigned int, dim> offset;
\n+
829 for (
int i=0; i<dim; i++)
\n+
830 offset[i] = std::max((
int)(currentKnotSpan[i] -
order_[i]),0);
\n+
\n+
\n+
833 std::array<std::vector<R>, dim> oneDValues;
\n+
\n+
\n+
836 std::array<std::vector<R>, dim> lowOrderOneDValues;
\n+
\n+
838 std::array<DynamicMatrix<R>, dim> values;
\n+
\n+
840 for (
size_t i=0; i<dim; i++)
\n+
\n+
\n+
\n+
844 for (
size_t j=0; j<oneDValues[i].size(); j++)
\n+
845 oneDValues[i][j] = values[i][
order_[i]][j];
\n+
\n+
\n+
\n+
\n+
850 for (
size_t j=0; j<lowOrderOneDValues[i].size(); j++)
\n+
851 lowOrderOneDValues[i][j] = values[i][
order_[i]-1][j];
\n+
\n+
\n+
\n+
\n+
\n+
857 std::array<std::vector<R>, dim> oneDDerivatives;
\n+
858 for (
size_t i=0; i<dim; i++)
\n+
\n+
860 oneDDerivatives[i].resize(limits[i]);
\n+
\n+
\n+
863 std::fill(oneDDerivatives[i].begin(), oneDDerivatives[i].end(),
R(0.0));
\n+
\n+
\n+
866 for (
size_t j=offset[i]; j<offset[i]+limits[i]; j++)
\n+
\n+
\n+
\n+
\n+
871 if (std::isnan(derivativeAddend1))
\n+
872 derivativeAddend1 = 0;
\n+
873 if (std::isnan(derivativeAddend2))
\n+
874 derivativeAddend2 = 0;
\n+
875 oneDDerivatives[i][j-offset[i]] =
order_[i] * ( derivativeAddend1 - derivativeAddend2 );
\n+
\n+
\n+
\n+
\n+
\n+
\n+
882 std::array<std::vector<R>, dim> oneDValuesShort;
\n+
\n+
884 for (
int i=0; i<dim; i++)
\n+
\n+
886 oneDValuesShort[i].resize(limits[i]);
\n+
\n+
888 for (
size_t j=0; j<limits[i]; j++)
\n+
889 oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
\n+
\n+
\n+
\n+
\n+
\n+
895 MultiDigitCounter ijkCounter(limits);
\n+
\n+
897 out.resize(ijkCounter.cycle());
\n+
\n+
\n+
900 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
\n+
901 for (
int j=0; j<dim; j++)
\n+
\n+
\n+
904 for (
int k=0; k<dim; k++)
\n+
905 out[i][0][j] *= (j==k) ? oneDDerivatives[k][ijkCounter[k]]
\n+
906 : oneDValuesShort[k][ijkCounter[k]];
\n+
\n+
\n+
\n+
\n+
\n+
912 template <
size_type k>
\n+
\n+
913 void evaluate(
const typename std::array<int,k>& directions,
\n+
914 const FieldVector<typename GV::ctype,dim>& in,
\n+
915 std::vector<FieldVector<R,1> >& out,
\n+
916 const std::array<unsigned,dim>& currentKnotSpan)
const \n+
\n+
918 if (k != 1 && k != 2)
\n+
919 DUNE_THROW(RangeError,
"Differentiation order greater than 2 is not supported!");
\n+
\n+
\n+
922 std::array<std::vector<R>, dim> oneDValues;
\n+
923 std::array<std::vector<R>, dim> oneDDerivatives;
\n+
924 std::array<std::vector<R>, dim> oneDSecondDerivatives;
\n+
\n+
\n+
\n+
928 for (
size_t i=0; i<dim; i++)
\n+
929 evaluateAll(in[i], oneDValues[i],
true, oneDDerivatives[i],
false, oneDSecondDerivatives[i],
knotVectors_[i],
order_[i], currentKnotSpan[i]);
\n+
\n+
931 for (
size_t i=0; i<dim; i++)
\n+
\n+
\n+
\n+
935 std::array<unsigned int, dim> offset;
\n+
936 for (
int i=0; i<dim; i++)
\n+
937 offset[i] = std::max((
int)(currentKnotSpan[i] -
order_[i]),0);
\n+
\n+
\n+
940 std::array<unsigned int, dim> limits;
\n+
941 for (
int i=0; i<dim; i++)
\n+
\n+
\n+
\n+
\n+
946 if (currentKnotSpan[i]<
order_[i])
\n+
947 limits[i] -= (
order_[i] - currentKnotSpan[i]);
\n+
\n+
\n+
\n+
\n+
\n+
\n+
954 std::array<std::vector<R>, dim> oneDValuesShort;
\n+
\n+
956 for (
int i=0; i<dim; i++)
\n+
\n+
958 oneDValuesShort[i].resize(limits[i]);
\n+
\n+
960 for (
size_t j=0; j<limits[i]; j++)
\n+
961 oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
\n+
\n+
\n+
\n+
965 MultiDigitCounter ijkCounter(limits);
\n+
\n+
967 out.resize(ijkCounter.cycle());
\n+
\n+
\n+
\n+
\n+
972 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
\n+
\n+
\n+
975 for (
int l=0; l<dim; l++)
\n+
976 out[i][0] *= (directions[0]==l) ? oneDDerivatives[l][ijkCounter[l]]
\n+
977 : oneDValuesShort[l][ijkCounter[l]];
\n+
\n+
\n+
\n+
\n+
\n+
\n+
984 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
\n+
\n+
\n+
987 for (
int j=0; j<dim; j++)
\n+
\n+
989 if (directions[0] != directions[1])
\n+
990 if (directions[0] == j || directions[1] == j)
\n+
991 out[i][0] *= oneDDerivatives[j][ijkCounter[j]];
\n+
\n+
993 out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
\n+
\n+
995 if (directions[0] == j)
\n+
996 out[i][0] *= oneDSecondDerivatives[j][ijkCounter[j]];
\n+
\n+
998 out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
1009 static std::array<unsigned int,dim>
getIJK(
typename GridView::IndexSet::IndexType idx, std::array<unsigned int,dim> elements)
\n+
\n+
1011 std::array<unsigned,dim> result;
\n+
1012 for (
int i=0; i<dim; i++)
\n+
\n+
1014 result[i] = idx%elements[i];
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
1029 const std::vector<R>& knotVector,
\n+
\n+
1031 unsigned int currentKnotSpan)
\n+
\n+
1033 std::size_t outSize = order+1;
\n+
1034 if (currentKnotSpan<order)
\n+
1035 outSize -= (order - currentKnotSpan);
\n+
1036 if ( order > (knotVector.size() - currentKnotSpan - 2) )
\n+
1037 outSize -= order - (knotVector.size() - currentKnotSpan - 2);
\n+
1038 out.resize(outSize);
\n+
\n+
\n+
1041 DynamicMatrix<R> N(order+1, knotVector.size()-1);
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
1049 for (
size_t i=0; i<knotVector.size()-1; i++)
\n+
1050 N[0][i] = (i == currentKnotSpan);
\n+
\n+
1052 for (
size_t r=1; r<=order; r++)
\n+
1053 for (
size_t i=0; i<knotVector.size()-r-1; i++)
\n+
\n+
1055 R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
\n+
1056 ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
\n+
\n+
1058 R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
\n+
1059 ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
\n+
\n+
1061 N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
\n+
\n+
\n+
1068 for (
size_t i=0; i<out.size(); i++) {
\n+
1069 out[i] = N[order][std::max((
int)(currentKnotSpan - order),0) + i];
\n+
\n+
\n+
\n+
\n+
\n+
\n+
1086 DynamicMatrix<R>& out,
\n+
1087 const std::vector<R>& knotVector,
\n+
\n+
1089 unsigned int currentKnotSpan)
\n+
\n+
\n+
1092 DynamicMatrix<R>& N = out;
\n+
\n+
1094 N.resize(order+1, knotVector.size()-1);
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
1102 for (
size_t i=0; i<knotVector.size()-1; i++)
\n+
1103 N[0][i] = (i == currentKnotSpan);
\n+
\n+
1105 for (
size_t r=1; r<=order; r++)
\n+
1106 for (
size_t i=0; i<knotVector.size()-r-1; i++)
\n+
\n+
1108 R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
\n+
1109 ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
\n+
\n+
1111 R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
\n+
1112 ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
\n+
\n+
1114 N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
1129 std::vector<R>& out,
\n+
\n+
1131 bool evaluateHessian, std::vector<R>& outHess,
\n+
1132 const std::vector<R>& knotVector,
\n+
\n+
1134 unsigned int currentKnotSpan)
\n+
\n+
\n+
\n+
\n+
1139 if (currentKnotSpan<order)
\n+
1140 limit -= (order - currentKnotSpan);
\n+
1141 if ( order > (knotVector.size() - currentKnotSpan - 2) )
\n+
1142 limit -= order - (knotVector.size() - currentKnotSpan - 2);
\n+
\n+
\n+
1145 unsigned int offset;
\n+
1146 offset = std::max((
int)(currentKnotSpan - order),0);
\n+
\n+
\n+
1149 DynamicMatrix<R> values;
\n+
\n+
\n+
\n+
1153 out.resize(knotVector.size()-order-1);
\n+
1154 for (
size_t j=0; j<out.size(); j++)
\n+
1155 out[j] = values[order][j];
\n+
\n+
\n+
1158 std::vector<R> lowOrderOneDValues;
\n+
\n+
\n+
\n+
1162 lowOrderOneDValues.resize(knotVector.size()-(order-1)-1);
\n+
1163 for (
size_t j=0; j<lowOrderOneDValues.size(); j++)
\n+
1164 lowOrderOneDValues[j] = values[order-1][j];
\n+
\n+
\n+
\n+
1168 std::vector<R> lowOrderTwoDValues;
\n+
\n+
1170 if (order>1 && evaluateHessian)
\n+
\n+
1172 lowOrderTwoDValues.resize(knotVector.size()-(order-2)-1);
\n+
1173 for (
size_t j=0; j<lowOrderTwoDValues.size(); j++)
\n+
1174 lowOrderTwoDValues[j] = values[order-2][j];
\n+
\n+
\n+
\n+
\n+
\n+
1180 outJac.resize(limit);
\n+
\n+
\n+
1183 std::fill(outJac.begin(), outJac.end(),
R(0.0));
\n+
\n+
\n+
1186 for (
size_t j=offset; j<offset+limit; j++)
\n+
\n+
1188 R derivativeAddend1 = lowOrderOneDValues[j] / (knotVector[j+order]-knotVector[j]);
\n+
1189 R derivativeAddend2 = lowOrderOneDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]);
\n+
\n+
1191 if (std::isnan(derivativeAddend1))
\n+
1192 derivativeAddend1 = 0;
\n+
1193 if (std::isnan(derivativeAddend2))
\n+
1194 derivativeAddend2 = 0;
\n+
1195 outJac[j-offset] = order * ( derivativeAddend1 - derivativeAddend2 );
\n+
\n+
\n+
\n+
\n+
\n+
1201 if (evaluateHessian)
\n+
\n+
1203 outHess.resize(limit);
\n+
\n+
\n+
1206 std::fill(outHess.begin(), outHess.end(),
R(0.0));
\n+
\n+
\n+
1209 for (
size_t j=offset; j<offset+limit; j++)
\n+
\n+
1211 assert(j+2 < lowOrderTwoDValues.size());
\n+
1212 R derivativeAddend1 = lowOrderTwoDValues[j] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order-1]-knotVector[j]);
\n+
1213 R derivativeAddend2 = lowOrderTwoDValues[j+1] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order]-knotVector[j+1]);
\n+
1214 R derivativeAddend3 = lowOrderTwoDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+order]-knotVector[j+1]);
\n+
1215 R derivativeAddend4 = lowOrderTwoDValues[j+2] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+1+order]-knotVector[j+2]);
\n+
\n+
\n+
1218 if (std::isnan(derivativeAddend1))
\n+
1219 derivativeAddend1 = 0;
\n+
1220 if (std::isnan(derivativeAddend2))
\n+
1221 derivativeAddend2 = 0;
\n+
1222 if (std::isnan(derivativeAddend3))
\n+
1223 derivativeAddend3 = 0;
\n+
1224 if (std::isnan(derivativeAddend4))
\n+
1225 derivativeAddend4 = 0;
\n+
1226 outHess[j-offset] = order * (order-1) * ( derivativeAddend1 - derivativeAddend2 -derivativeAddend3 + derivativeAddend4 );
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
1247template<
typename GV>
\n+
\n+
\n+
\n+
\n+
1251 static const int dim = GV::dimension;
\n+
\n+
\n+
\n+
\n+
1256 using Element =
typename GV::template Codim<0>::Entity;
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
1283 auto elementIndex =
preBasis_->gridView().indexSet().index(e);
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
1298namespace BasisFactory {
\n+
\n+
\n+
1306inline auto bSpline(
const std::vector<double>& knotVector,
\n+
\n+
1308 bool makeOpen =
true)
\n+
\n+
1310 return [&knotVector, order, makeOpen](
const auto& gridView) {
\n+
1311 return BSplinePreBasis<std::decay_t<
decltype(gridView)>>(gridView, knotVector, order, makeOpen);
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
1327template<
typename GV>
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
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:1306
\n
Definition polynomial.hh:10
\n+
LocalFiniteElement in the sense of dune-localfunctions, for the B-spline basis on tensor-product grid...
Definition bsplinebasis.hh:362
\n+
BSplineLocalFiniteElement(const BSplineLocalFiniteElement &other)
Copy constructor.
Definition bsplinebasis.hh:383
\n+
const BSplinePreBasis< GV > & preBasis_
Definition bsplinebasis.hh:474
\n+
const BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > & localInterpolation() const
Hand out a LocalInterpolation object.
Definition bsplinebasis.hh:439
\n+
LocalFiniteElementTraits< BSplineLocalBasis< GV, R >, BSplineLocalCoefficients< dim >, BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > > Traits
Export various types related to this LocalFiniteElement.
Definition bsplinebasis.hh:372
\n+
std::array< unsigned, dim > currentKnotSpan_
Definition bsplinebasis.hh:481
\n+
BSplineLocalFiniteElement(const BSplinePreBasis< GV > &preBasis)
Constructor with a given B-spline basis.
Definition bsplinebasis.hh:376
\n+
const BSplineLocalCoefficients< dim > & localCoefficients() const
Hand out a LocalCoefficients object.
Definition bsplinebasis.hh:433
\n+
void bind(const std::array< unsigned, dim > &elementIdx)
Bind LocalFiniteElement to a specific knot span of the spline patch.
Definition bsplinebasis.hh:394
\n+
BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > localInterpolation_
Definition bsplinebasis.hh:478
\n+
GeometryType type() const
Return the reference element that the local finite element is defined on (here, a hypercube)
Definition bsplinebasis.hh:455
\n+
unsigned size() const
Number of shape functions in this finite element.
Definition bsplinebasis.hh:445
\n+
BSplineLocalCoefficients< dim > localCoefficients_
Definition bsplinebasis.hh:477
\n+
unsigned int size(int i) const
Number of degrees of freedom for one coordinate direction.
Definition bsplinebasis.hh:463
\n+
BSplineLocalBasis< GV, R > localBasis_
Definition bsplinebasis.hh:476
\n+
const BSplineLocalBasis< GV, R > & localBasis() const
Hand out a LocalBasis object.
Definition bsplinebasis.hh:427
\n+
Pre-basis for B-spline basis.
Definition bsplinebasis.hh:499
\n+
static constexpr size_type multiIndexBufferSize
Definition bsplinebasis.hh:567
\n+
std::array< unsigned, dim > elements_
Number of grid elements in the different coordinate directions.
Definition bsplinebasis.hh:1240
\n+
GridView gridView_
Definition bsplinebasis.hh:1242
\n+
double R
Definition bsplinebasis.hh:570
\n+
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:1085
\n+
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:781
\n+
std::array< unsigned int, dim > order_
Order of the B-spline for each space dimension.
Definition bsplinebasis.hh:1234
\n+
size_type dimension() const
Get the total dimension of the space spanned by this basis.
Definition bsplinebasis.hh:719
\n+
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:1128
\n+
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:1028
\n+
unsigned int size(size_t d) const
Number of shape functions in one direction.
Definition bsplinebasis.hh:774
\n+
GV GridView
The grid view that the FE space is defined on.
Definition bsplinebasis.hh:560
\n+
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:735
\n+
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:913
\n+
std::size_t size_type
Definition bsplinebasis.hh:561
\n+
static constexpr size_type maxMultiIndexSize
Definition bsplinebasis.hh:565
\n+
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition bsplinebasis.hh:686
\n+
void initializeIndices()
Initialize the global indices.
Definition bsplinebasis.hh:676
\n+
size_type size(const Dune::ReservedVector< ST, i > &prefix) const
Return number of possible values for next position in multi index.
Definition bsplinebasis.hh:712
\n+
unsigned int size() const
Total number of B-spline basis functions.
Definition bsplinebasis.hh:765
\n+
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition bsplinebasis.hh:680
\n+
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:1009
\n+
Node makeNode() const
Create tree node.
Definition bsplinebasis.hh:694
\n+
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:590
\n+
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:812
\n+
static constexpr size_type minMultiIndexSize
Definition bsplinebasis.hh:566
\n+
std::array< std::vector< double >, dim > knotVectors_
The knot vectors, one for each space dimension.
Definition bsplinebasis.hh:1237
\n+
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition bsplinebasis.hh:725
\n+
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:642
\n+
LocalBasis class in the sense of dune-localfunctions, presenting the restriction of a B-spline patch ...
Definition bsplinebasis.hh:46
\n+
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:55
\n+
unsigned int order() const
Polynomial order of the shape functions.
Definition bsplinebasis.hh:140
\n+
std::size_t size() const
Return the number of basis functions on the current knot span.
Definition bsplinebasis.hh:147
\n+
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:97
\n+
void evaluateFunction(const FieldVector< D, dim > &in, std::vector< FieldVector< R, 1 > > &out) const
Evaluate all shape functions.
Definition bsplinebasis.hh:70
\n+
void evaluateJacobian(const FieldVector< D, dim > &in, std::vector< FieldMatrix< D, 1, dim > > &out) const
Evaluate Jacobian of all shape functions.
Definition bsplinebasis.hh:82
\n+
BSplineLocalBasis(const BSplinePreBasis< GV > &preBasis, const BSplineLocalFiniteElement< GV, R > &lFE)
Constructor with a given B-spline patch.
Definition bsplinebasis.hh:61
\n+
Attaches a shape function to an entity.
Definition bsplinebasis.hh:178
\n+
const LocalKey & localKey(std::size_t i) const
get i'th index
Definition bsplinebasis.hh:321
\n+
void init(const std::array< unsigned, dim > &sizes)
Definition bsplinebasis.hh:256
\n+
std::size_t size() const
number of coefficients
Definition bsplinebasis.hh:315
\n+
Local interpolation in the sense of dune-localfunctions, for the B-spline basis on tensor-product gri...
Definition bsplinebasis.hh:340
\n+
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition bsplinebasis.hh:344
\n+
Definition bsplinebasis.hh:1250
\n+
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition bsplinebasis.hh:1274
\n+
typename GV::template Codim< 0 >::Entity Element
Definition bsplinebasis.hh:1256
\n+
const BSplinePreBasis< GV > * preBasis_
Definition bsplinebasis.hh:1290
\n+
Element element_
Definition bsplinebasis.hh:1293
\n+
void bind(const Element &e)
Bind to element.
Definition bsplinebasis.hh:1280
\n+
BSplineNode(const BSplinePreBasis< GV > *preBasis)
Definition bsplinebasis.hh:1259
\n+
const Element & element() const
Return current element, throw if unbound.
Definition bsplinebasis.hh:1265
\n+
FiniteElement finiteElement_
Definition bsplinebasis.hh:1292
\n+
std::size_t size_type
Definition bsplinebasis.hh:1255
\n+
Global basis for given pre-basis.
Definition defaultglobalbasis.hh:46
\n+
size_type size() const
Definition nodes.hh:142
\n+
void setSize(const size_type size)
Definition nodes.hh:164
\n+
\n
![\"doxygen\"/](\"doxygen.svg\")