\n
\n-
3#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH \n-
4#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH \n+
3#ifndef DUNE_FUNCTIONS_GRIDFUNCTIONS_GRIDVIEWENTITYSET_HH \n+
4#define DUNE_FUNCTIONS_GRIDFUNCTIONS_GRIDVIEWENTITYSET_HH \n
\n-
\n-
\n+
\n+
\n+
\n+
\n+
\n+
\n
\n-
14#include <dune/common/dynmatrix.hh> \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-
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-
54 typedef LocalBasisTraits<D,dim,FieldVector<D,dim>,R,1,FieldVector<R,1>,
\n-
\n+
\n+
21template<
class GV,
int cd>
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
32 typedef typename GridView::template Codim<codim>::Entity
Element;
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
54 return gv_.contains(e);
\n+
\n+
\n
\n-
\n-
\n-
\n-
63 : preBasis_(preBasis),
\n-
\n-
\n+
\n+
\n+
\n+
60 return gv_.size(
codim);
\n+
\n+
\n+
\n+
\n+
\n+
\n+
66 return gv_.template begin<codim>();
\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+
\n+
\n+
72 return gv_.template end<codim>();
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\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+
An entity set for all entities of given codim in a grid view.
Definition gridviewentityset.hh:23
\n+
GridViewEntitySet(const GridView &gv)
Construct GridViewEntitySet for a GridView.
Definition gridviewentityset.hh:47
\n+
GV GridView
Definition gridviewentityset.hh:26
\n+
GridView::template Codim< codim >::Entity Element
Type of Elements contained in this EntitySet.
Definition gridviewentityset.hh:32
\n+
const_iterator end() const
Create an end iterator.
Definition gridviewentityset.hh:70
\n+
const GridView & gridView() const
Return the associated GridView.
Definition gridviewentityset.hh:76
\n+
Element value_type
Definition gridviewentityset.hh:38
\n+
const_iterator begin() const
Create a begin iterator.
Definition gridviewentityset.hh:64
\n+
GridView::template Codim< codim >::Iterator const_iterator
A forward iterator.
Definition gridviewentityset.hh:41
\n+
Element::Geometry::LocalCoordinate LocalCoordinate
Type of local coordinates with respect to the Element.
Definition gridviewentityset.hh:35
\n+
size_t size() const
Return number of Elements visited by an iterator.
Definition gridviewentityset.hh:58
\n+
Element::Geometry::GlobalCoordinate GlobalCoordinate
Definition gridviewentityset.hh:36
\n+
@ codim
Definition gridviewentityset.hh:28
\n+
bool contains(const Element &e) const
Return true if e is contained in the EntitySet.
Definition gridviewentityset.hh:52
\n+
const_iterator iterator
Same as const_iterator.
Definition gridviewentityset.hh:44
\n
![\"doxygen\"/](\"doxygen.svg\")
1.9.8\n \n \n", "details": [{"source1": "html2text {}", "source2": "html2text {}", "unified_diff": "@@ -1,30 +1,29 @@\n dune-functions\u00a02.9.0\n Loading...\n Searching...\n No Matches\n * _\bd_\bu_\bn_\be\n * _\bf_\bu_\bn_\bc_\bt_\bi_\bo_\bn_\bs\n- * _\ba_\bn_\ba_\bl_\by_\bt_\bi_\bc_\bf_\bu_\bn_\bc_\bt_\bi_\bo_\bn_\bs\n-_\bC_\bl_\ba_\bs_\bs_\be_\bs | _\bN_\ba_\bm_\be_\bs_\bp_\ba_\bc_\be_\bs | _\bF_\bu_\bn_\bc_\bt_\bi_\bo_\bn_\bs\n-trigonometricfunction.hh File Reference\n-#include