dune-functions  2.5-dev
bsplinebasis.hh
Go to the documentation of this file.
1 #ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
2 #define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
3 
8 #include <numeric>
9 
11 #include <dune/common/dynmatrix.hh>
12 
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>
19 
20 namespace Dune
21 {
22 namespace Functions {
23 
24 // A maze of dependencies between the different parts of this. We need lots of forward declarations
25 template<typename GV, typename R>
27 
28 template<typename GV, class MI>
30 
31 
40 template<class GV, class R>
42 {
43  friend class BSplineLocalFiniteElement<GV,R>;
44 
45  typedef typename GV::ctype D;
46  enum {dim = GV::dimension};
47 public:
48 
50  typedef LocalBasisTraits<D,dim,Dune::FieldVector<D,dim>,R,1,Dune::FieldVector<R,1>,
51  Dune::FieldMatrix<R,1,dim>, 2> Traits;
52 
59  : nodeFactory_(nodeFactory),
60  lFE_(lFE)
61  {}
62 
66  void evaluateFunction (const FieldVector<D,dim>& in,
67  std::vector<FieldVector<R,1> >& out) const
68  {
69  FieldVector<D,dim> globalIn = offset_;
70  scaling_.umv(in,globalIn);
71 
72  nodeFactory_.evaluateFunction(globalIn, out, lFE_.currentKnotSpan_);
73  }
74 
78  void evaluateJacobian (const FieldVector<D,dim>& in,
79  std::vector<FieldMatrix<D,1,dim> >& out) const
80  {
81  FieldVector<D,dim> globalIn = offset_;
82  scaling_.umv(in,globalIn);
83 
84  nodeFactory_.evaluateJacobian(globalIn, out, lFE_.currentKnotSpan_);
85 
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];
89  }
90 
92  template<size_t k>
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
96  {
97  switch(k)
98  {
99  case 0:
100  evaluateFunction(in, out);
101  break;
102  case 1:
103  {
104  FieldVector<D,dim> globalIn = offset_;
105  scaling_.umv(in,globalIn);
106 
107  nodeFactory_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
108 
109  for (size_t i=0; i<out.size(); i++)
110  out[i][0] *= scaling_[directions[0]][directions[0]];
111  break;
112  }
113  case 2:
114  {
115  FieldVector<D,dim> globalIn = offset_;
116  scaling_.umv(in,globalIn);
117 
118  nodeFactory_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
119 
120  for (size_t i=0; i<out.size(); i++)
121  out[i][0] *= scaling_[directions[0]][directions[0]]*scaling_[directions[1]][directions[1]];
122  break;
123  }
124  default:
125  DUNE_THROW(NotImplemented, "B-Spline derivatives of order " << k << " not implemented yet!");
126  }
127  }
128 
136  unsigned int order () const
137  {
138  return *std::max_element(nodeFactory_.order_.begin(), nodeFactory_.order_.end());
139  }
140 
143  std::size_t size() const
144  {
145  return lFE_.size();
146  }
147 
148 private:
150 
152 
153  // Coordinates in a single knot span differ from coordinates on the B-spline patch
154  // by an affine transformation. This transformation is stored in offset_ and scaling_.
155  FieldVector<D,dim> offset_;
156  DiagonalMatrix<D,dim> scaling_;
157 };
158 
172 template<int dim>
174 {
175  // Return i as a d-digit number in the (k+1)-nary system
176  std::array<unsigned int,dim> multiindex (unsigned int i) const
177  {
178  std::array<unsigned int,dim> alpha;
179  for (int j=0; j<dim; j++)
180  {
181  alpha[j] = i % sizes_[j];
182  i = i/sizes_[j];
183  }
184  return alpha;
185  }
186 
188  void setup1d(std::vector<unsigned int>& subEntity)
189  {
190  if (sizes_[0]==1)
191  {
192  subEntity[0] = 0;
193  return;
194  }
195 
196  /* edge and vertex numbering
197  0----0----1
198  */
199  unsigned lastIndex=0;
200  subEntity[lastIndex++] = 0; // corner 0
201  for (unsigned i = 0; i < sizes_[0] - 2; ++i)
202  subEntity[lastIndex++] = 0; // inner dofs of element (0)
203 
204  subEntity[lastIndex++] = 1; // corner 1
205 
206  assert(size()==lastIndex);
207  }
208 
209  void setup2d(std::vector<unsigned int>& subEntity)
210  {
211  unsigned lastIndex=0;
212 
213  // LocalKey: entity number , entity codim, dof indices within each entity
214  /* edge and vertex numbering
215  2----3----3
216  | |
217  | |
218  0 1
219  | |
220  | |
221  0----2----1
222  */
223 
224  // lower edge (2)
225  subEntity[lastIndex++] = 0; // corner 0
226  for (unsigned i = 0; i < sizes_[0]-2; ++i)
227  subEntity[lastIndex++] = 2; // inner dofs of lower edge (2)
228 
229  subEntity[lastIndex++] = 1; // corner 1
230 
231  // iterate from bottom to top over inner edge dofs
232  for (unsigned e = 0; e < sizes_[1]-2; ++e)
233  {
234  subEntity[lastIndex++] = 0; // left edge (0)
235  for (unsigned i = 0; i < sizes_[0]-2; ++i)
236  subEntity[lastIndex++] = 0; // face dofs
237  subEntity[lastIndex++] = 1; // right edge (1)
238  }
239 
240  // upper edge (3)
241  subEntity[lastIndex++] = 2; // corner 2
242  for (unsigned i = 0; i < sizes_[0]-2; ++i)
243  subEntity[lastIndex++] = 3; // inner dofs of upper edge (3)
244 
245  subEntity[lastIndex++] = 3; // corner 3
246 
247  assert(size()==lastIndex);
248  }
249 
250 
251 public:
252  void init(const std::array<unsigned,dim>& sizes)
253  {
254  sizes_ = sizes;
255 
256  li_.resize(size());
257 
258  // Set up array of codimension-per-dof-number
259  std::vector<unsigned int> codim(li_.size());
260 
261  for (std::size_t i=0; i<codim.size(); i++)
262  {
263  codim[i] = 0;
264  // Codimension gets increased by 1 for each coordinate direction
265  // where dof is on boundary
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)
269  codim[i]++;
270  }
271 
272  // Set up index vector (the index of the dof in the set of dofs of a given subentity)
273  // Algorithm: the 'index' has the same ordering as the dof number 'i'.
274  // To make it consecutive we interpret 'i' in the (k+1)-adic system, omit all digits
275  // that correspond to axes where the dof is on the element boundary, and transform the
276  // rest to the (k-1)-adic system.
277  std::vector<unsigned int> index(size());
278 
279  for (std::size_t i=0; i<index.size(); i++)
280  {
281  index[i] = 0;
282 
283  std::array<unsigned int,dim> mIdx = multiindex(i);
284 
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);
288  }
289 
290  // Set up entity and dof numbers for each (supported) dimension separately
291  std::vector<unsigned int> subEntity(li_.size());
292 
293  if (subEntity.size() > 0)
294  {
295  if (dim==1) {
296 
297  setup1d(subEntity);
298 
299  } else if (dim==2 and sizes_[0]>1 and sizes_[1]>1) {
300 
301  setup2d(subEntity);
302 
303  }
304  }
305 
306  for (size_t i=0; i<li_.size(); i++)
307  li_[i] = LocalKey(subEntity[i], codim[i], index[i]);
308  }
309 
311  std::size_t size () const
312  {
313  return std::accumulate(sizes_.begin(), sizes_.end(), 1, std::multiplies<unsigned int>());
314  }
315 
317  const LocalKey& localKey (std::size_t i) const
318  {
319  return li_[i];
320  }
321 
322 private:
323 
324  // Number of shape functions on this element per coordinate direction
325  std::array<unsigned, dim> sizes_;
326 
327  std::vector<LocalKey> li_;
328 };
329 
334 template<int dim, class LB>
336 {
337 public:
339  template<typename F, typename C>
340  void interpolate (const F& f, std::vector<C>& out) const
341  {
342  DUNE_THROW(NotImplemented, "BSplineLocalInterpolation::interpolate");
343  }
344 
345 };
346 
357 template<class GV, class R>
359 {
360  typedef typename GV::ctype D;
361  enum {dim = GV::dimension};
362  friend class BSplineLocalBasis<GV,R>;
363 public:
364 
367  typedef LocalFiniteElementTraits<BSplineLocalBasis<GV,R>,
370 
374  : nodeFactory_(nodeFactory),
375  localBasis_(nodeFactory,*this)
376  {}
377 
384  void bind(const std::array<uint,dim>& elementIdx)
385  {
386  /* \todo In the long run we need to precompute a table for this */
387  for (size_t i=0; i<elementIdx.size(); i++)
388  {
389  currentKnotSpan_[i] = 0;
390 
391  // Skip over degenerate knot spans
392  while (nodeFactory_.knotVectors_[i][currentKnotSpan_[i]+1] < nodeFactory_.knotVectors_[i][currentKnotSpan_[i]]+1e-8)
393  currentKnotSpan_[i]++;
394 
395  for (size_t j=0; j<elementIdx[i]; j++)
396  {
397  currentKnotSpan_[i]++;
398 
399  // Skip over degenerate knot spans
400  while (nodeFactory_.knotVectors_[i][currentKnotSpan_[i]+1] < nodeFactory_.knotVectors_[i][currentKnotSpan_[i]]+1e-8)
401  currentKnotSpan_[i]++;
402  }
403 
404  // Compute the geometric transformation from knotspan-local to global coordinates
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]];
407  }
408 
409  // Set up the LocalCoefficients object
410  std::array<unsigned int, dim> sizes;
411  for (size_t i=0; i<dim; i++)
412  sizes[i] = size(i);
413  localCoefficients_.init(sizes);
414  }
415 
418  {
419  return localBasis_;
420  }
421 
423  const BSplineLocalCoefficients<dim>& localCoefficients() const
424  {
425  return localCoefficients_;
426  }
427 
430  {
431  return localInterpolation_;
432  }
433 
435  uint size () const
436  {
437  std::size_t r = 1;
438  for (int i=0; i<dim; i++)
439  r *= size(i);
440  return r;
441  }
442 
445  GeometryType type () const
446  {
447  return GeometryType(GeometryType::cube,dim);
448  }
449 
450 //private:
451 
453  unsigned int size(int i) const
454  {
455  const auto& order = nodeFactory_.order_;
456  unsigned int r = order[i]+1; // The 'normal' value
457  if (currentKnotSpan_[i]<order[i]) // Less near the left end of the knot vector
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);
461  return r;
462  }
463 
465 
467  BSplineLocalCoefficients<dim> localCoefficients_;
469 
470  // The knot span we are bound to
471  std::array<uint,dim> currentKnotSpan_;
472 };
473 
474 
475 template<typename GV, typename TP>
477 
478 template<typename GV, class MI, class TP>
480 
491 template<typename GV, class MI>
492 class BSplineNodeFactory
493 {
494  static const int dim = GV::dimension;
495 
497  class MultiDigitCounter
498  {
499  public:
500 
504  MultiDigitCounter(const std::array<unsigned int,dim>& limits)
505  : limits_(limits)
506  {
507  std::fill(counter_.begin(), counter_.end(), 0);
508  }
509 
511  MultiDigitCounter& operator++()
512  {
513  for (int i=0; i<dim; i++)
514  {
515  ++counter_[i];
516 
517  // no overflow?
518  if (counter_[i] < limits_[i])
519  break;
520 
521  counter_[i] = 0;
522  }
523  return *this;
524  }
525 
527  const unsigned int& operator[](int i) const
528  {
529  return counter_[i];
530  }
531 
533  unsigned int cycle() const
534  {
535  unsigned int r = 1;
536  for (int i=0; i<dim; i++)
537  r *= limits_[i];
538  return r;
539  }
540 
541  private:
542 
544  const std::array<unsigned int,dim> limits_;
545 
547  std::array<unsigned int,dim> counter_;
548 
549  };
550 
551 public:
552 
554  using GridView = GV;
555  using size_type = std::size_t;
556 
557  template<class TP>
559 
560  template<class TP>
562 
564  using MultiIndex = MI;
565 
566  using SizePrefix = Dune::ReservedVector<size_type, 1>;
567 
568  // Type used for function values
569  using R = double;
570 
589  BSplineNodeFactory(const GridView& gridView,
590  const std::vector<double>& knotVector,
591  unsigned int order,
592  bool makeOpen = true)
593  : gridView_(gridView)
594  {
595  // \todo Detection of duplicate knots
596  std::fill(elements_.begin(), elements_.end(), knotVector.size()-1);
597 
598  // Mediocre sanity check: we don't know the number of grid elements in each direction.
599  // but at least we know the total number of elements.
600  assert( std::accumulate(elements_.begin(), elements_.end(), 1, std::multiplies<uint>()) == gridView_.size(0) );
601 
602  for (int i=0; i<dim; i++)
603  {
604  // Prepend the correct number of additional knots to open the knot vector
606  if (makeOpen)
607  for (unsigned int j=0; j<order; j++)
608  knotVectors_[i].push_back(knotVector[0]);
609 
610  knotVectors_[i].insert(knotVectors_[i].end(), knotVector.begin(), knotVector.end());
611 
612  if (makeOpen)
613  for (unsigned int j=0; j<order; j++)
614  knotVectors_[i].push_back(knotVector.back());
615  }
616 
617  std::fill(order_.begin(), order_.end(), order);
618  }
619 
641  BSplineNodeFactory(const GridView& gridView,
642  const FieldVector<double,dim>& lowerLeft,
643  const FieldVector<double,dim>& upperRight,
644  const array<unsigned int,dim>& elements,
645  unsigned int order,
646  bool makeOpen = true)
647  : elements_(elements),
648  gridView_(gridView)
649  {
650  // Mediocre sanity check: we don't know the number of grid elements in each direction.
651  // but at least we know the total number of elements.
652  assert( std::accumulate(elements_.begin(), elements_.end(), 1, std::multiplies<uint>()) == gridView_.size(0) );
653 
654  for (int i=0; i<dim; i++)
655  {
656  // Prepend the correct number of additional knots to open the knot vector
658  if (makeOpen)
659  for (unsigned int j=0; j<order; j++)
660  knotVectors_[i].push_back(lowerLeft[i]);
661 
662  // Construct the actual knot vector
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]);
665 
666  if (makeOpen)
667  for (unsigned int j=0; j<order; j++)
668  knotVectors_[i].push_back(upperRight[i]);
669  }
670 
671  std::fill(order_.begin(), order_.end(), order);
672  }
673 
676  {}
677 
679  const GridView& gridView() const
680  {
681  return gridView_;
682  }
683 
685  void update(const GridView& gv)
686  {
687  gridView_ = gv;
688  }
689 
700  template<class TP>
701  Node<TP> node(const TP& tp) const
702  {
703  return Node<TP>{tp,this};
704  }
705 
715  template<class TP>
717  {
718  return IndexSet<TP>{*this};
719  }
720 
722  size_type size(const SizePrefix prefix) const
723  {
724  if (prefix.size() == 0)
725  return size();
726  assert(false);
727  }
728 
731  {
732  return size();
733  }
734 
737  {
738  size_type result = 1;
739  for (int i=0; i<dim; i++)
740  result *= order_[i]+1;
741  return result;
742  }
743 
745  unsigned int size () const
746  {
747  unsigned int result = 1;
748  for (size_t i=0; i<dim; i++)
749  result *= size(i);
750  return result;
751  }
752 
754  unsigned int size (size_t d) const
755  {
756  return knotVectors_[d].size() - order_[d] - 1;
757  }
758 
761  void evaluateFunction (const FieldVector<typename GV::ctype,dim>& in,
762  std::vector<FieldVector<R,1> >& out,
763  const std::array<uint,dim>& currentKnotSpan) const
764  {
765  // Evaluate
766  Dune::array<std::vector<R>, dim> oneDValues;
767 
768  for (size_t i=0; i<dim; i++)
769  evaluateFunction(in[i], oneDValues[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
770 
771  std::array<unsigned int, dim> limits;
772  for (int i=0; i<dim; i++)
773  limits[i] = oneDValues[i].size();
774 
775  MultiDigitCounter ijkCounter(limits);
776 
777  out.resize(ijkCounter.cycle());
778 
779  for (size_t i=0; i<out.size(); i++, ++ijkCounter)
780  {
781  out[i] = R(1.0);
782  for (size_t j=0; j<dim; j++)
783  out[i] *= oneDValues[j][ijkCounter[j]];
784  }
785  }
786 
792  void evaluateJacobian (const FieldVector<typename GV::ctype,dim>& in,
793  std::vector<FieldMatrix<R,1,dim> >& out,
794  const std::array<uint,dim>& currentKnotSpan) const
795  {
796  // How many shape functions to we have in each coordinate direction?
797  std::array<unsigned int, dim> limits;
798  for (int i=0; i<dim; i++)
799  {
800  limits[i] = order_[i]+1; // The 'standard' value away from the boundaries of the knot vector
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);
805  }
806 
807  // The lowest knot spans that we need values from
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);
811 
812  // Evaluate 1d function values (needed for the product rule)
813  Dune::array<std::vector<R>, dim> oneDValues;
814 
815  // Evaluate 1d function values of one order lower (needed for the derivative formula)
816  Dune::array<std::vector<R>, dim> lowOrderOneDValues;
817 
818  Dune::array<DynamicMatrix<R>, dim> values;
819 
820  for (size_t i=0; i<dim; i++)
821  {
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];
826 
827  if (order_[i]!=0)
828  {
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];
832  }
833  }
834 
835 
836  // Evaluate 1d function derivatives
837  Dune::array<std::vector<R>, dim> oneDDerivatives;
838  for (size_t i=0; i<dim; i++)
839  {
840  oneDDerivatives[i].resize(limits[i]);
841 
842  if (order_[i]==0) // order-zero functions are piecewise constant, hence all derivatives are zero
843  std::fill(oneDDerivatives[i].begin(), oneDDerivatives[i].end(), R(0.0));
844  else
845  {
846  for (size_t j=offset[i]; j<offset[i]+limits[i]; j++)
847  {
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]);
850  // The two previous terms may evaluate as 0/0. This is to be interpreted as 0.
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 );
856  }
857  }
858  }
859 
860  // Working towards computing only the parts that we really need:
861  // Let's copy them out into a separate array
862  Dune::array<std::vector<R>, dim> oneDValuesShort;
863 
864  for (int i=0; i<dim; i++)
865  {
866  oneDValuesShort[i].resize(limits[i]);
867 
868  for (size_t j=0; j<limits[i]; j++)
869  oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
870  }
871 
872 
873 
874  // Set up a multi-index to go from consecutive indices to integer coordinates
875  MultiDigitCounter ijkCounter(limits);
876 
877  out.resize(ijkCounter.cycle());
878 
879  // Complete Jacobian is given by the product rule
880  for (size_t i=0; i<out.size(); i++, ++ijkCounter)
881  for (int j=0; j<dim; j++)
882  {
883  out[i][0][j] = 1.0;
884  for (int k=0; k<dim; k++)
885  out[i][0][j] *= (j==k) ? oneDDerivatives[k][ijkCounter[k]]
886  : oneDValuesShort[k][ijkCounter[k]];
887  }
888 
889  }
890 
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
897  {
898  if (k != 1 && k != 2)
899  DUNE_THROW(RangeError, "Differentiation order greater than 2 is not supported!");
900 
901  // Evaluate 1d function values (needed for the product rule)
902  std::array<std::vector<R>, dim> oneDValues;
903  std::array<std::vector<R>, dim> oneDDerivatives;
904  std::array<std::vector<R>, dim> oneDSecondDerivatives;
905 
906  // Evaluate 1d function derivatives
907  if (k==1)
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]);
910  else
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]);
913 
914  // The lowest knot spans that we need values from
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);
918 
919  // Set up a multi-index to go from consecutive indices to integer coordinates
920  std::array<unsigned int, dim> limits;
921  for (int i=0; i<dim; i++)
922  {
923  // In a proper implementation, the following line would do
924  //limits[i] = oneDValues[i].size();
925  limits[i] = order_[i]+1; // The 'standard' value away from the boundaries of the knot vector
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);
930  }
931 
932  // Working towards computing only the parts that we really need:
933  // Let's copy them out into a separate array
934  std::array<std::vector<R>, dim> oneDValuesShort;
935 
936  for (int i=0; i<dim; i++)
937  {
938  oneDValuesShort[i].resize(limits[i]);
939 
940  for (size_t j=0; j<limits[i]; j++)
941  oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
942  }
943 
944 
945  MultiDigitCounter ijkCounter(limits);
946 
947  out.resize(ijkCounter.cycle());
948 
949  if (k == 1)
950  {
951  // Complete Jacobian is given by the product rule
952  for (size_t i=0; i<out.size(); i++, ++ijkCounter)
953  {
954  out[i][0] = 1.0;
955  for (int l=0; l<dim; l++)
956  out[i][0] *= (directions[0]==l) ? oneDDerivatives[l][ijkCounter[l]]
957  : oneDValuesShort[l][ijkCounter[l]];
958  }
959  }
960 
961  if (k == 2)
962  {
963  // Complete derivation by deriving the tensor product
964  for (size_t i=0; i<out.size(); i++, ++ijkCounter)
965  {
966  out[i][0] = 1.0;
967  for (int j=0; j<dim; j++)
968  {
969  if (directions[0] != directions[1]) //derivation in two different variables
970  if (directions[0] == j || directions[1] == j) //the spline has to be derived (once) in this direction
971  out[i][0] *= oneDDerivatives[j][ijkCounter[j]];
972  else //no derivation in this direction
973  out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
974  else //spline is derived two times in the same direction
975  if (directions[0] == j) //the spline is derived two times in this direction
976  out[i][0] *= oneDSecondDerivatives[j][ijkCounter[j]];
977  else //no derivation in this direction
978  out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
979  }
980  }
981  }
982  }
983 
984 
989  static std::array<unsigned int,dim> getIJK(typename GridView::IndexSet::IndexType idx, std::array<unsigned int,dim> elements)
990  {
991  std::array<uint,dim> result;
992  for (int i=0; i<dim; i++)
993  {
994  result[i] = idx%elements[i];
995  idx /= elements[i];
996  }
997  return result;
998  }
999 
1008  static void evaluateFunction (const typename GV::ctype& in, std::vector<R>& out,
1009  const std::vector<R>& knotVector,
1010  unsigned int order,
1011  unsigned int currentKnotSpan)
1012  {
1013  std::size_t outSize = order+1; // The 'standard' value away from the boundaries of the knot vector
1014  if (currentKnotSpan<order) // Less near the left end of the knot vector
1015  outSize -= (order - currentKnotSpan);
1016  if ( order > (knotVector.size() - currentKnotSpan - 2) )
1017  outSize -= order - (knotVector.size() - currentKnotSpan - 2);
1018  out.resize(outSize);
1019 
1020  // It's not really a matrix that is needed here, a plain 2d array would do
1021  DynamicMatrix<R> N(order+1, knotVector.size()-1);
1022 
1023  // The text books on splines use the following geometric condition here to fill the array N
1024  // (see for example Cottrell, Hughes, Bazilevs, Formula (2.1). However, this condition
1025  // only works if splines are never evaluated exactly on the knots.
1026  //
1027  // for (size_t i=0; i<knotVector.size()-1; i++)
1028  // N[0][i] = (knotVector[i] <= in) and (in < knotVector[i+1]);
1029  for (size_t i=0; i<knotVector.size()-1; i++)
1030  N[0][i] = (i == currentKnotSpan);
1031 
1032  for (size_t r=1; r<=order; r++)
1033  for (size_t i=0; i<knotVector.size()-r-1; i++)
1034  {
1035  R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
1036  ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
1037  : 0;
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])
1040  : 0;
1041  N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
1042  }
1043 
1048  for (size_t i=0; i<out.size(); i++) {
1049  out[i] = N[order][std::max((int)(currentKnotSpan - order),0) + i];
1050  }
1051  }
1052 
1065  static void evaluateFunctionFull(const typename GV::ctype& in,
1066  DynamicMatrix<R>& out,
1067  const std::vector<R>& knotVector,
1068  unsigned int order,
1069  unsigned int currentKnotSpan)
1070  {
1071  // It's not really a matrix that is needed here, a plain 2d array would do
1072  DynamicMatrix<R>& N = out;
1073 
1074  N.resize(order+1, knotVector.size()-1);
1075 
1076  // The text books on splines use the following geometric condition here to fill the array N
1077  // (see for example Cottrell, Hughes, Bazilevs, Formula (2.1). However, this condition
1078  // only works if splines are never evaluated exactly on the knots.
1079  //
1080  // for (size_t i=0; i<knotVector.size()-1; i++)
1081  // N[0][i] = (knotVector[i] <= in) and (in < knotVector[i+1]);
1082  for (size_t i=0; i<knotVector.size()-1; i++)
1083  N[0][i] = (i == currentKnotSpan);
1084 
1085  for (size_t r=1; r<=order; r++)
1086  for (size_t i=0; i<knotVector.size()-r-1; i++)
1087  {
1088  R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
1089  ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
1090  : 0;
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])
1093  : 0;
1094  N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
1095  }
1096  }
1097 
1098 
1108  static void evaluateAll(const typename GV::ctype& in,
1109  std::vector<R>& out,
1110  bool evaluateJacobian, std::vector<R>& outJac,
1111  bool evaluateHessian, std::vector<R>& outHess,
1112  const std::vector<R>& knotVector,
1113  unsigned int order,
1114  unsigned int currentKnotSpan)
1115  {
1116  // How many shape functions to we have in each coordinate direction?
1117  unsigned int limit;
1118  limit = order+1; // The 'standard' value away from the boundaries of the knot vector
1119  if (currentKnotSpan<order)
1120  limit -= (order - currentKnotSpan);
1121  if ( order > (knotVector.size() - currentKnotSpan - 2) )
1122  limit -= order - (knotVector.size() - currentKnotSpan - 2);
1123 
1124  // The lowest knot spans that we need values from
1125  unsigned int offset;
1126  offset = std::max((int)(currentKnotSpan - order),0);
1127 
1128  // Evaluate 1d function values (needed for the product rule)
1129  DynamicMatrix<R> values;
1130 
1131  evaluateFunctionFull(in, values, knotVector, order, currentKnotSpan);
1132 
1133  out.resize(knotVector.size()-order-1);
1134  for (size_t j=0; j<out.size(); j++)
1135  out[j] = values[order][j];
1136 
1137  // Evaluate 1d function values of one order lower (needed for the derivative formula)
1138  std::vector<R> lowOrderOneDValues;
1139 
1140  if (order!=0)
1141  {
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];
1145  }
1146 
1147  // Evaluate 1d function values of two order lower (needed for the (second) derivative formula)
1148  std::vector<R> lowOrderTwoDValues;
1149 
1150  if (order>1 && evaluateHessian)
1151  {
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];
1155  }
1156 
1157  // Evaluate 1d function derivatives
1158  if (evaluateJacobian)
1159  {
1160  outJac.resize(limit);
1161 
1162  if (order==0) // order-zero functions are piecewise constant, hence all derivatives are zero
1163  std::fill(outJac.begin(), outJac.end(), R(0.0));
1164  else
1165  {
1166  for (size_t j=offset; j<offset+limit; j++)
1167  {
1168  R derivativeAddend1 = lowOrderOneDValues[j] / (knotVector[j+order]-knotVector[j]);
1169  R derivativeAddend2 = lowOrderOneDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]);
1170  // The two previous terms may evaluate as 0/0. This is to be interpreted as 0.
1171  if (std::isnan(derivativeAddend1))
1172  derivativeAddend1 = 0;
1173  if (std::isnan(derivativeAddend2))
1174  derivativeAddend2 = 0;
1175  outJac[j-offset] = order * ( derivativeAddend1 - derivativeAddend2 );
1176  }
1177  }
1178  }
1179 
1180  // Evaluate 1d function second derivatives
1181  if (evaluateHessian)
1182  {
1183  outHess.resize(limit);
1184 
1185  if (order<2) // order-zero functions are piecewise constant, hence all derivatives are zero
1186  std::fill(outHess.begin(), outHess.end(), R(0.0));
1187  else
1188  {
1189  for (size_t j=offset; j<offset+limit; j++)
1190  {
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]);
1196  // The two previous terms may evaluate as 0/0. This is to be interpreted as 0.
1197 
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 );
1207  }
1208  }
1209  }
1210  }
1211 
1212 
1214  array<unsigned int, dim> order_;
1215 
1217  array<std::vector<double>, dim> knotVectors_;
1218 
1220  std::array<uint,dim> elements_;
1221 
1223 };
1224 
1225 
1226 
1227 template<typename GV, typename TP>
1228 class BSplineNode :
1229  public LeafBasisNode<std::size_t, TP>
1230 {
1231  static const int dim = GV::dimension;
1232 
1233  using Base = LeafBasisNode<std::size_t,TP>;
1234 
1235 public:
1236 
1237  using size_type = std::size_t;
1238  using TreePath = TP;
1239  using Element = typename GV::template Codim<0>::Entity;
1241 
1242  BSplineNode(const TreePath& treePath, const BSplineNodeFactory<GV, FlatMultiIndex<std::size_t>>* nodeFactory) :
1243  Base(treePath),
1244  nodeFactory_(nodeFactory),
1245  finiteElement_(*nodeFactory)
1246  {}
1247 
1249  const Element& element() const
1250  {
1251  return element_;
1252  }
1253 
1259  {
1260  return finiteElement_;
1261  }
1262 
1264  void bind(const Element& e)
1265  {
1266  element_ = e;
1267  auto elementIndex = nodeFactory_->gridView().indexSet().index(e);
1268  finiteElement_.bind(nodeFactory_->getIJK(elementIndex,nodeFactory_->elements_));
1269  this->setSize(finiteElement_.size());
1270  }
1271 
1272 protected:
1273 
1275 
1278 };
1279 
1280 
1281 
1282 template<typename GV, class MI, class TP>
1283 class BSplineNodeIndexSet
1284 {
1285  enum {dim = GV::dimension};
1286 
1287 public:
1288 
1289  using size_type = std::size_t;
1290 
1292  using MultiIndex = MI;
1293 
1295 
1296  using Node = typename NodeFactory::template Node<TP>;
1297 
1298  BSplineNodeIndexSet(const NodeFactory& nodeFactory) :
1299  nodeFactory_(&nodeFactory)
1300  {}
1301 
1307  void bind(const Node& node)
1308  {
1309  node_ = &node;
1310  // Local degrees of freedom are arranged in a lattice.
1311  // We need the lattice dimensions to be able to compute lattice coordinates from a local index
1312  for (int i=0; i<dim; i++)
1313  localSizes_[i] = node_->finiteElement().size(i);
1314  }
1315 
1318  void unbind()
1319  {
1320  node_ = nullptr;
1321  }
1322 
1325  size_type size() const
1326  {
1327  return node_->finiteElement().size();
1328  }
1329 
1332  {
1333  std::array<unsigned int,dim> localIJK = nodeFactory_->getIJK(i, localSizes_);
1334 
1335  const auto currentKnotSpan = node_->finiteElement().currentKnotSpan_;
1336  const auto order = nodeFactory_->order_;
1337 
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]; // needs to be a signed type!
1341 
1342  // Make one global flat index from the globalIJK tuple
1343  size_type globalIdx = globalIJK[dim-1];
1344 
1345  for (int i=dim-2; i>=0; i--)
1346  globalIdx = globalIdx * nodeFactory_->size(i) + globalIJK[i];
1347 
1348  return { globalIdx };
1349  }
1350 
1351 protected:
1353 
1354  const Node* node_;
1355 
1356  std::array<unsigned int, dim> localSizes_;
1357 };
1358 
1359 
1360 
1361 // *****************************************************************************
1362 // This is the actual global basis implementation based on the reusable parts.
1363 // *****************************************************************************
1364 
1371 template<typename GV>
1373 
1374 
1375 } // namespace Functions
1376 
1377 } // namespace Dune
1378 
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
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&#39;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
Definition: nodes.hh:189
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 > &currentKnotSpan) 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 > &currentKnotSpan) 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 > &currentKnotSpan) 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