dune-functions  2.6-dev
bsplinebasis.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
4 #define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
5 
10 #include <array>
11 #include <numeric>
12 
14 #include <dune/common/dynmatrix.hh>
15 
16 #include <dune/localfunctions/common/localbasis.hh>
17 #include <dune/common/diagonalmatrix.hh>
18 #include <dune/localfunctions/common/localkey.hh>
19 #include <dune/localfunctions/common/localfiniteelementtraits.hh>
20 #include <dune/geometry/type.hh>
22 
23 namespace Dune
24 {
25 namespace Functions {
26 
27 // A maze of dependencies between the different parts of this. We need lots of forward declarations
28 template<typename GV, typename R>
30 
31 template<typename GV, class MI>
33 
34 
43 template<class GV, class R>
45 {
46  friend class BSplineLocalFiniteElement<GV,R>;
47 
48  typedef typename GV::ctype D;
49  enum {dim = GV::dimension};
50 public:
51 
53  typedef LocalBasisTraits<D,dim,Dune::FieldVector<D,dim>,R,1,Dune::FieldVector<R,1>,
54  Dune::FieldMatrix<R,1,dim> > Traits;
55 
62  : preBasis_(preBasis),
63  lFE_(lFE)
64  {}
65 
69  void evaluateFunction (const FieldVector<D,dim>& in,
70  std::vector<FieldVector<R,1> >& out) const
71  {
72  FieldVector<D,dim> globalIn = offset_;
73  scaling_.umv(in,globalIn);
74 
75  preBasis_.evaluateFunction(globalIn, out, lFE_.currentKnotSpan_);
76  }
77 
81  void evaluateJacobian (const FieldVector<D,dim>& in,
82  std::vector<FieldMatrix<D,1,dim> >& out) const
83  {
84  FieldVector<D,dim> globalIn = offset_;
85  scaling_.umv(in,globalIn);
86 
87  preBasis_.evaluateJacobian(globalIn, out, lFE_.currentKnotSpan_);
88 
89  for (size_t i=0; i<out.size(); i++)
90  for (int j=0; j<dim; j++)
91  out[i][0][j] *= scaling_[j][j];
92  }
93 
95  template<size_t k>
96  inline void evaluate (const typename std::array<int,k>& directions,
97  const typename Traits::DomainType& in,
98  std::vector<typename Traits::RangeType>& out) const
99  {
100  switch(k)
101  {
102  case 0:
103  evaluateFunction(in, out);
104  break;
105  case 1:
106  {
107  FieldVector<D,dim> globalIn = offset_;
108  scaling_.umv(in,globalIn);
109 
110  preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
111 
112  for (size_t i=0; i<out.size(); i++)
113  out[i][0] *= scaling_[directions[0]][directions[0]];
114  break;
115  }
116  case 2:
117  {
118  FieldVector<D,dim> globalIn = offset_;
119  scaling_.umv(in,globalIn);
120 
121  preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
122 
123  for (size_t i=0; i<out.size(); i++)
124  out[i][0] *= scaling_[directions[0]][directions[0]]*scaling_[directions[1]][directions[1]];
125  break;
126  }
127  default:
128  DUNE_THROW(NotImplemented, "B-Spline derivatives of order " << k << " not implemented yet!");
129  }
130  }
131 
139  unsigned int order () const
140  {
141  return *std::max_element(preBasis_.order_.begin(), preBasis_.order_.end());
142  }
143 
146  std::size_t size() const
147  {
148  return lFE_.size();
149  }
150 
151 private:
153 
155 
156  // Coordinates in a single knot span differ from coordinates on the B-spline patch
157  // by an affine transformation. This transformation is stored in offset_ and scaling_.
158  FieldVector<D,dim> offset_;
159  DiagonalMatrix<D,dim> scaling_;
160 };
161 
175 template<int dim>
177 {
178  // Return i as a d-digit number in the (k+1)-nary system
179  std::array<unsigned int,dim> multiindex (unsigned int i) const
180  {
181  std::array<unsigned int,dim> alpha;
182  for (int j=0; j<dim; j++)
183  {
184  alpha[j] = i % sizes_[j];
185  i = i/sizes_[j];
186  }
187  return alpha;
188  }
189 
191  void setup1d(std::vector<unsigned int>& subEntity)
192  {
193  if (sizes_[0]==1)
194  {
195  subEntity[0] = 0;
196  return;
197  }
198 
199  /* edge and vertex numbering
200  0----0----1
201  */
202  unsigned lastIndex=0;
203  subEntity[lastIndex++] = 0; // corner 0
204  for (unsigned i = 0; i < sizes_[0] - 2; ++i)
205  subEntity[lastIndex++] = 0; // inner dofs of element (0)
206 
207  subEntity[lastIndex++] = 1; // corner 1
208 
209  assert(size()==lastIndex);
210  }
211 
212  void setup2d(std::vector<unsigned int>& subEntity)
213  {
214  unsigned lastIndex=0;
215 
216  // LocalKey: entity number , entity codim, dof indices within each entity
217  /* edge and vertex numbering
218  2----3----3
219  | |
220  | |
221  0 1
222  | |
223  | |
224  0----2----1
225  */
226 
227  // lower edge (2)
228  subEntity[lastIndex++] = 0; // corner 0
229  for (unsigned i = 0; i < sizes_[0]-2; ++i)
230  subEntity[lastIndex++] = 2; // inner dofs of lower edge (2)
231 
232  subEntity[lastIndex++] = 1; // corner 1
233 
234  // iterate from bottom to top over inner edge dofs
235  for (unsigned e = 0; e < sizes_[1]-2; ++e)
236  {
237  subEntity[lastIndex++] = 0; // left edge (0)
238  for (unsigned i = 0; i < sizes_[0]-2; ++i)
239  subEntity[lastIndex++] = 0; // face dofs
240  subEntity[lastIndex++] = 1; // right edge (1)
241  }
242 
243  // upper edge (3)
244  subEntity[lastIndex++] = 2; // corner 2
245  for (unsigned i = 0; i < sizes_[0]-2; ++i)
246  subEntity[lastIndex++] = 3; // inner dofs of upper edge (3)
247 
248  subEntity[lastIndex++] = 3; // corner 3
249 
250  assert(size()==lastIndex);
251  }
252 
253 
254 public:
255  void init(const std::array<unsigned,dim>& sizes)
256  {
257  sizes_ = sizes;
258 
259  li_.resize(size());
260 
261  // Set up array of codimension-per-dof-number
262  std::vector<unsigned int> codim(li_.size());
263 
264  for (std::size_t i=0; i<codim.size(); i++)
265  {
266  codim[i] = 0;
267  // Codimension gets increased by 1 for each coordinate direction
268  // where dof is on boundary
269  std::array<unsigned int,dim> mIdx = multiindex(i);
270  for (int j=0; j<dim; j++)
271  if (mIdx[j]==0 or mIdx[j]==sizes[j]-1)
272  codim[i]++;
273  }
274 
275  // Set up index vector (the index of the dof in the set of dofs of a given subentity)
276  // Algorithm: the 'index' has the same ordering as the dof number 'i'.
277  // To make it consecutive we interpret 'i' in the (k+1)-adic system, omit all digits
278  // that correspond to axes where the dof is on the element boundary, and transform the
279  // rest to the (k-1)-adic system.
280  std::vector<unsigned int> index(size());
281 
282  for (std::size_t i=0; i<index.size(); i++)
283  {
284  index[i] = 0;
285 
286  std::array<unsigned int,dim> mIdx = multiindex(i);
287 
288  for (int j=dim-1; j>=0; j--)
289  if (mIdx[j]>0 and mIdx[j]<sizes[j]-1)
290  index[i] = (sizes[j]-1)*index[i] + (mIdx[j]-1);
291  }
292 
293  // Set up entity and dof numbers for each (supported) dimension separately
294  std::vector<unsigned int> subEntity(li_.size());
295 
296  if (subEntity.size() > 0)
297  {
298  if (dim==1) {
299 
300  setup1d(subEntity);
301 
302  } else if (dim==2 and sizes_[0]>1 and sizes_[1]>1) {
303 
304  setup2d(subEntity);
305 
306  }
307  }
308 
309  for (size_t i=0; i<li_.size(); i++)
310  li_[i] = LocalKey(subEntity[i], codim[i], index[i]);
311  }
312 
314  std::size_t size () const
315  {
316  return std::accumulate(sizes_.begin(), sizes_.end(), 1, std::multiplies<unsigned int>());
317  }
318 
320  const LocalKey& localKey (std::size_t i) const
321  {
322  return li_[i];
323  }
324 
325 private:
326 
327  // Number of shape functions on this element per coordinate direction
328  std::array<unsigned, dim> sizes_;
329 
330  std::vector<LocalKey> li_;
331 };
332 
337 template<int dim, class LB>
339 {
340 public:
342  template<typename F, typename C>
343  void interpolate (const F& f, std::vector<C>& out) const
344  {
345  DUNE_THROW(NotImplemented, "BSplineLocalInterpolation::interpolate");
346  }
347 
348 };
349 
360 template<class GV, class R>
362 {
363  typedef typename GV::ctype D;
364  enum {dim = GV::dimension};
365  friend class BSplineLocalBasis<GV,R>;
366 public:
367 
370  typedef LocalFiniteElementTraits<BSplineLocalBasis<GV,R>,
373 
377  : preBasis_(preBasis),
378  localBasis_(preBasis,*this)
379  {}
380 
387  void bind(const std::array<unsigned,dim>& elementIdx)
388  {
389  /* \todo In the long run we need to precompute a table for this */
390  for (size_t i=0; i<elementIdx.size(); i++)
391  {
392  currentKnotSpan_[i] = 0;
393 
394  // Skip over degenerate knot spans
395  while (preBasis_.knotVectors_[i][currentKnotSpan_[i]+1] < preBasis_.knotVectors_[i][currentKnotSpan_[i]]+1e-8)
396  currentKnotSpan_[i]++;
397 
398  for (size_t j=0; j<elementIdx[i]; j++)
399  {
400  currentKnotSpan_[i]++;
401 
402  // Skip over degenerate knot spans
403  while (preBasis_.knotVectors_[i][currentKnotSpan_[i]+1] < preBasis_.knotVectors_[i][currentKnotSpan_[i]]+1e-8)
404  currentKnotSpan_[i]++;
405  }
406 
407  // Compute the geometric transformation from knotspan-local to global coordinates
408  localBasis_.offset_[i] = preBasis_.knotVectors_[i][currentKnotSpan_[i]];
409  localBasis_.scaling_[i][i] = preBasis_.knotVectors_[i][currentKnotSpan_[i]+1] - preBasis_.knotVectors_[i][currentKnotSpan_[i]];
410  }
411 
412  // Set up the LocalCoefficients object
413  std::array<unsigned int, dim> sizes;
414  for (size_t i=0; i<dim; i++)
415  sizes[i] = size(i);
416  localCoefficients_.init(sizes);
417  }
418 
421  {
422  return localBasis_;
423  }
424 
426  const BSplineLocalCoefficients<dim>& localCoefficients() const
427  {
428  return localCoefficients_;
429  }
430 
433  {
434  return localInterpolation_;
435  }
436 
438  unsigned size () const
439  {
440  std::size_t r = 1;
441  for (int i=0; i<dim; i++)
442  r *= size(i);
443  return r;
444  }
445 
448  GeometryType type () const
449  {
450  return GeometryType(GeometryType::cube,dim);
451  }
452 
453 //private:
454 
456  unsigned int size(int i) const
457  {
458  const auto& order = preBasis_.order_;
459  unsigned int r = order[i]+1; // The 'normal' value
460  if (currentKnotSpan_[i]<order[i]) // Less near the left end of the knot vector
461  r -= (order[i] - currentKnotSpan_[i]);
462  if ( order[i] > (preBasis_.knotVectors_[i].size() - currentKnotSpan_[i] - 2) )
463  r -= order[i] - (preBasis_.knotVectors_[i].size() - currentKnotSpan_[i] - 2);
464  return r;
465  }
466 
468 
470  BSplineLocalCoefficients<dim> localCoefficients_;
472 
473  // The knot span we are bound to
474  std::array<unsigned,dim> currentKnotSpan_;
475 };
476 
477 
478 template<typename GV, typename TP>
480 
481 template<typename GV, class MI, class TP>
483 
494 template<typename GV, class MI>
495 class BSplinePreBasis
496 {
497  static const int dim = GV::dimension;
498 
500  class MultiDigitCounter
501  {
502  public:
503 
507  MultiDigitCounter(const std::array<unsigned int,dim>& limits)
508  : limits_(limits)
509  {
510  std::fill(counter_.begin(), counter_.end(), 0);
511  }
512 
514  MultiDigitCounter& operator++()
515  {
516  for (int i=0; i<dim; i++)
517  {
518  ++counter_[i];
519 
520  // no overflow?
521  if (counter_[i] < limits_[i])
522  break;
523 
524  counter_[i] = 0;
525  }
526  return *this;
527  }
528 
530  const unsigned int& operator[](int i) const
531  {
532  return counter_[i];
533  }
534 
536  unsigned int cycle() const
537  {
538  unsigned int r = 1;
539  for (int i=0; i<dim; i++)
540  r *= limits_[i];
541  return r;
542  }
543 
544  private:
545 
547  const std::array<unsigned int,dim> limits_;
548 
550  std::array<unsigned int,dim> counter_;
551 
552  };
553 
554 public:
555 
557  using GridView = GV;
558  using size_type = std::size_t;
559 
560  template<class TP>
562 
563  template<class TP>
565 
567  using MultiIndex = MI;
568 
569  using SizePrefix = Dune::ReservedVector<size_type, 1>;
570 
571  // Type used for function values
572  using R = double;
573 
592  BSplinePreBasis(const GridView& gridView,
593  const std::vector<double>& knotVector,
594  unsigned int order,
595  bool makeOpen = true)
596  : gridView_(gridView)
597  {
598  // \todo Detection of duplicate knots
599  std::fill(elements_.begin(), elements_.end(), knotVector.size()-1);
600 
601  // Mediocre sanity check: we don't know the number of grid elements in each direction.
602  // but at least we know the total number of elements.
603  assert( std::accumulate(elements_.begin(), elements_.end(), 1, std::multiplies<unsigned>()) == gridView_.size(0) );
604 
605  for (int i=0; i<dim; i++)
606  {
607  // Prepend the correct number of additional knots to open the knot vector
609  if (makeOpen)
610  for (unsigned int j=0; j<order; j++)
611  knotVectors_[i].push_back(knotVector[0]);
612 
613  knotVectors_[i].insert(knotVectors_[i].end(), knotVector.begin(), knotVector.end());
614 
615  if (makeOpen)
616  for (unsigned int j=0; j<order; j++)
617  knotVectors_[i].push_back(knotVector.back());
618  }
619 
620  std::fill(order_.begin(), order_.end(), order);
621  }
622 
644  BSplinePreBasis(const GridView& gridView,
645  const FieldVector<double,dim>& lowerLeft,
646  const FieldVector<double,dim>& upperRight,
647  const std::array<unsigned int,dim>& elements,
648  unsigned int order,
649  bool makeOpen = true)
650  : elements_(elements),
651  gridView_(gridView)
652  {
653  // Mediocre sanity check: we don't know the number of grid elements in each direction.
654  // but at least we know the total number of elements.
655  assert( std::accumulate(elements_.begin(), elements_.end(), 1, std::multiplies<unsigned>()) == gridView_.size(0) );
656 
657  for (int i=0; i<dim; i++)
658  {
659  // Prepend the correct number of additional knots to open the knot vector
661  if (makeOpen)
662  for (unsigned int j=0; j<order; j++)
663  knotVectors_[i].push_back(lowerLeft[i]);
664 
665  // Construct the actual knot vector
666  for (size_t j=0; j<elements[i]+1; j++)
667  knotVectors_[i].push_back(lowerLeft[i] + j*(upperRight[i]-lowerLeft[i]) / elements[i]);
668 
669  if (makeOpen)
670  for (unsigned int j=0; j<order; j++)
671  knotVectors_[i].push_back(upperRight[i]);
672  }
673 
674  std::fill(order_.begin(), order_.end(), order);
675  }
676 
679  {}
680 
682  const GridView& gridView() const
683  {
684  return gridView_;
685  }
686 
688  void update(const GridView& gv)
689  {
690  gridView_ = gv;
691  }
692 
703  template<class TP>
704  Node<TP> node(const TP& tp) const
705  {
706  return Node<TP>{tp,this};
707  }
708 
718  template<class TP>
720  {
721  return IndexSet<TP>{*this};
722  }
723 
725  size_type size(const SizePrefix prefix) const
726  {
727  assert(prefix.size() == 0 || prefix.size() == 1);
728  return (prefix.size() == 0) ? size() : 0;
729  }
730 
733  {
734  return size();
735  }
736 
739  {
740  size_type result = 1;
741  for (int i=0; i<dim; i++)
742  result *= order_[i]+1;
743  return result;
744  }
745 
747  unsigned int size () const
748  {
749  unsigned int result = 1;
750  for (size_t i=0; i<dim; i++)
751  result *= size(i);
752  return result;
753  }
754 
756  unsigned int size (size_t d) const
757  {
758  return knotVectors_[d].size() - order_[d] - 1;
759  }
760 
763  void evaluateFunction (const FieldVector<typename GV::ctype,dim>& in,
764  std::vector<FieldVector<R,1> >& out,
765  const std::array<unsigned,dim>& currentKnotSpan) const
766  {
767  // Evaluate
768  std::array<std::vector<R>, dim> oneDValues;
769 
770  for (size_t i=0; i<dim; i++)
771  evaluateFunction(in[i], oneDValues[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
772 
773  std::array<unsigned int, dim> limits;
774  for (int i=0; i<dim; i++)
775  limits[i] = oneDValues[i].size();
776 
777  MultiDigitCounter ijkCounter(limits);
778 
779  out.resize(ijkCounter.cycle());
780 
781  for (size_t i=0; i<out.size(); i++, ++ijkCounter)
782  {
783  out[i] = R(1.0);
784  for (size_t j=0; j<dim; j++)
785  out[i] *= oneDValues[j][ijkCounter[j]];
786  }
787  }
788 
794  void evaluateJacobian (const FieldVector<typename GV::ctype,dim>& in,
795  std::vector<FieldMatrix<R,1,dim> >& out,
796  const std::array<unsigned,dim>& currentKnotSpan) const
797  {
798  // How many shape functions to we have in each coordinate direction?
799  std::array<unsigned int, dim> limits;
800  for (int i=0; i<dim; i++)
801  {
802  limits[i] = order_[i]+1; // The 'standard' value away from the boundaries of the knot vector
803  if (currentKnotSpan[i]<order_[i])
804  limits[i] -= (order_[i] - currentKnotSpan[i]);
805  if ( order_[i] > (knotVectors_[i].size() - currentKnotSpan[i] - 2) )
806  limits[i] -= order_[i] - (knotVectors_[i].size() - currentKnotSpan[i] - 2);
807  }
808 
809  // The lowest knot spans that we need values from
810  std::array<unsigned int, dim> offset;
811  for (int i=0; i<dim; i++)
812  offset[i] = std::max((int)(currentKnotSpan[i] - order_[i]),0);
813 
814  // Evaluate 1d function values (needed for the product rule)
815  std::array<std::vector<R>, dim> oneDValues;
816 
817  // Evaluate 1d function values of one order lower (needed for the derivative formula)
818  std::array<std::vector<R>, dim> lowOrderOneDValues;
819 
820  std::array<DynamicMatrix<R>, dim> values;
821 
822  for (size_t i=0; i<dim; i++)
823  {
824  evaluateFunctionFull(in[i], values[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
825  oneDValues[i].resize(knotVectors_[i].size()-order_[i]-1);
826  for (size_t j=0; j<oneDValues[i].size(); j++)
827  oneDValues[i][j] = values[i][order_[i]][j];
828 
829  if (order_[i]!=0)
830  {
831  lowOrderOneDValues[i].resize(knotVectors_[i].size()-(order_[i]-1)-1);
832  for (size_t j=0; j<lowOrderOneDValues[i].size(); j++)
833  lowOrderOneDValues[i][j] = values[i][order_[i]-1][j];
834  }
835  }
836 
837 
838  // Evaluate 1d function derivatives
839  std::array<std::vector<R>, dim> oneDDerivatives;
840  for (size_t i=0; i<dim; i++)
841  {
842  oneDDerivatives[i].resize(limits[i]);
843 
844  if (order_[i]==0) // order-zero functions are piecewise constant, hence all derivatives are zero
845  std::fill(oneDDerivatives[i].begin(), oneDDerivatives[i].end(), R(0.0));
846  else
847  {
848  for (size_t j=offset[i]; j<offset[i]+limits[i]; j++)
849  {
850  R derivativeAddend1 = lowOrderOneDValues[i][j] / (knotVectors_[i][j+order_[i]]-knotVectors_[i][j]);
851  R derivativeAddend2 = lowOrderOneDValues[i][j+1] / (knotVectors_[i][j+order_[i]+1]-knotVectors_[i][j+1]);
852  // The two previous terms may evaluate as 0/0. This is to be interpreted as 0.
853  if (std::isnan(derivativeAddend1))
854  derivativeAddend1 = 0;
855  if (std::isnan(derivativeAddend2))
856  derivativeAddend2 = 0;
857  oneDDerivatives[i][j-offset[i]] = order_[i] * ( derivativeAddend1 - derivativeAddend2 );
858  }
859  }
860  }
861 
862  // Working towards computing only the parts that we really need:
863  // Let's copy them out into a separate array
864  std::array<std::vector<R>, dim> oneDValuesShort;
865 
866  for (int i=0; i<dim; i++)
867  {
868  oneDValuesShort[i].resize(limits[i]);
869 
870  for (size_t j=0; j<limits[i]; j++)
871  oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
872  }
873 
874 
875 
876  // Set up a multi-index to go from consecutive indices to integer coordinates
877  MultiDigitCounter ijkCounter(limits);
878 
879  out.resize(ijkCounter.cycle());
880 
881  // Complete Jacobian is given by the product rule
882  for (size_t i=0; i<out.size(); i++, ++ijkCounter)
883  for (int j=0; j<dim; j++)
884  {
885  out[i][0][j] = 1.0;
886  for (int k=0; k<dim; k++)
887  out[i][0][j] *= (j==k) ? oneDDerivatives[k][ijkCounter[k]]
888  : oneDValuesShort[k][ijkCounter[k]];
889  }
890 
891  }
892 
894  template <size_type k>
895  void evaluate(const typename std::array<int,k>& directions,
896  const FieldVector<typename GV::ctype,dim>& in,
897  std::vector<FieldVector<R,1> >& out,
898  const std::array<unsigned,dim>& currentKnotSpan) const
899  {
900  if (k != 1 && k != 2)
901  DUNE_THROW(RangeError, "Differentiation order greater than 2 is not supported!");
902 
903  // Evaluate 1d function values (needed for the product rule)
904  std::array<std::vector<R>, dim> oneDValues;
905  std::array<std::vector<R>, dim> oneDDerivatives;
906  std::array<std::vector<R>, dim> oneDSecondDerivatives;
907 
908  // Evaluate 1d function derivatives
909  if (k==1)
910  for (size_t i=0; i<dim; i++)
911  evaluateAll(in[i], oneDValues[i], true, oneDDerivatives[i], false, oneDSecondDerivatives[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
912  else
913  for (size_t i=0; i<dim; i++)
914  evaluateAll(in[i], oneDValues[i], true, oneDDerivatives[i], true, oneDSecondDerivatives[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
915 
916  // The lowest knot spans that we need values from
917  std::array<unsigned int, dim> offset;
918  for (int i=0; i<dim; i++)
919  offset[i] = std::max((int)(currentKnotSpan[i] - order_[i]),0);
920 
921  // Set up a multi-index to go from consecutive indices to integer coordinates
922  std::array<unsigned int, dim> limits;
923  for (int i=0; i<dim; i++)
924  {
925  // In a proper implementation, the following line would do
926  //limits[i] = oneDValues[i].size();
927  limits[i] = order_[i]+1; // The 'standard' value away from the boundaries of the knot vector
928  if (currentKnotSpan[i]<order_[i])
929  limits[i] -= (order_[i] - currentKnotSpan[i]);
930  if ( order_[i] > (knotVectors_[i].size() - currentKnotSpan[i] - 2) )
931  limits[i] -= order_[i] - (knotVectors_[i].size() - currentKnotSpan[i] - 2);
932  }
933 
934  // Working towards computing only the parts that we really need:
935  // Let's copy them out into a separate array
936  std::array<std::vector<R>, dim> oneDValuesShort;
937 
938  for (int i=0; i<dim; i++)
939  {
940  oneDValuesShort[i].resize(limits[i]);
941 
942  for (size_t j=0; j<limits[i]; j++)
943  oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
944  }
945 
946 
947  MultiDigitCounter ijkCounter(limits);
948 
949  out.resize(ijkCounter.cycle());
950 
951  if (k == 1)
952  {
953  // Complete Jacobian is given by the product rule
954  for (size_t i=0; i<out.size(); i++, ++ijkCounter)
955  {
956  out[i][0] = 1.0;
957  for (int l=0; l<dim; l++)
958  out[i][0] *= (directions[0]==l) ? oneDDerivatives[l][ijkCounter[l]]
959  : oneDValuesShort[l][ijkCounter[l]];
960  }
961  }
962 
963  if (k == 2)
964  {
965  // Complete derivation by deriving the tensor product
966  for (size_t i=0; i<out.size(); i++, ++ijkCounter)
967  {
968  out[i][0] = 1.0;
969  for (int j=0; j<dim; j++)
970  {
971  if (directions[0] != directions[1]) //derivation in two different variables
972  if (directions[0] == j || directions[1] == j) //the spline has to be derived (once) in this direction
973  out[i][0] *= oneDDerivatives[j][ijkCounter[j]];
974  else //no derivation in this direction
975  out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
976  else //spline is derived two times in the same direction
977  if (directions[0] == j) //the spline is derived two times in this direction
978  out[i][0] *= oneDSecondDerivatives[j][ijkCounter[j]];
979  else //no derivation in this direction
980  out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
981  }
982  }
983  }
984  }
985 
986 
991  static std::array<unsigned int,dim> getIJK(typename GridView::IndexSet::IndexType idx, std::array<unsigned int,dim> elements)
992  {
993  std::array<unsigned,dim> result;
994  for (int i=0; i<dim; i++)
995  {
996  result[i] = idx%elements[i];
997  idx /= elements[i];
998  }
999  return result;
1000  }
1001 
1010  static void evaluateFunction (const typename GV::ctype& in, std::vector<R>& out,
1011  const std::vector<R>& knotVector,
1012  unsigned int order,
1013  unsigned int currentKnotSpan)
1014  {
1015  std::size_t outSize = order+1; // The 'standard' value away from the boundaries of the knot vector
1016  if (currentKnotSpan<order) // Less near the left end of the knot vector
1017  outSize -= (order - currentKnotSpan);
1018  if ( order > (knotVector.size() - currentKnotSpan - 2) )
1019  outSize -= order - (knotVector.size() - currentKnotSpan - 2);
1020  out.resize(outSize);
1021 
1022  // It's not really a matrix that is needed here, a plain 2d array would do
1023  DynamicMatrix<R> N(order+1, knotVector.size()-1);
1024 
1025  // The text books on splines use the following geometric condition here to fill the array N
1026  // (see for example Cottrell, Hughes, Bazilevs, Formula (2.1). However, this condition
1027  // only works if splines are never evaluated exactly on the knots.
1028  //
1029  // for (size_t i=0; i<knotVector.size()-1; i++)
1030  // N[0][i] = (knotVector[i] <= in) and (in < knotVector[i+1]);
1031  for (size_t i=0; i<knotVector.size()-1; i++)
1032  N[0][i] = (i == currentKnotSpan);
1033 
1034  for (size_t r=1; r<=order; r++)
1035  for (size_t i=0; i<knotVector.size()-r-1; i++)
1036  {
1037  R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
1038  ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
1039  : 0;
1040  R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
1041  ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
1042  : 0;
1043  N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
1044  }
1045 
1050  for (size_t i=0; i<out.size(); i++) {
1051  out[i] = N[order][std::max((int)(currentKnotSpan - order),0) + i];
1052  }
1053  }
1054 
1067  static void evaluateFunctionFull(const typename GV::ctype& in,
1068  DynamicMatrix<R>& out,
1069  const std::vector<R>& knotVector,
1070  unsigned int order,
1071  unsigned int currentKnotSpan)
1072  {
1073  // It's not really a matrix that is needed here, a plain 2d array would do
1074  DynamicMatrix<R>& N = out;
1075 
1076  N.resize(order+1, knotVector.size()-1);
1077 
1078  // The text books on splines use the following geometric condition here to fill the array N
1079  // (see for example Cottrell, Hughes, Bazilevs, Formula (2.1). However, this condition
1080  // only works if splines are never evaluated exactly on the knots.
1081  //
1082  // for (size_t i=0; i<knotVector.size()-1; i++)
1083  // N[0][i] = (knotVector[i] <= in) and (in < knotVector[i+1]);
1084  for (size_t i=0; i<knotVector.size()-1; i++)
1085  N[0][i] = (i == currentKnotSpan);
1086 
1087  for (size_t r=1; r<=order; r++)
1088  for (size_t i=0; i<knotVector.size()-r-1; i++)
1089  {
1090  R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
1091  ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
1092  : 0;
1093  R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
1094  ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
1095  : 0;
1096  N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
1097  }
1098  }
1099 
1100 
1110  static void evaluateAll(const typename GV::ctype& in,
1111  std::vector<R>& out,
1112  bool evaluateJacobian, std::vector<R>& outJac,
1113  bool evaluateHessian, std::vector<R>& outHess,
1114  const std::vector<R>& knotVector,
1115  unsigned int order,
1116  unsigned int currentKnotSpan)
1117  {
1118  // How many shape functions to we have in each coordinate direction?
1119  unsigned int limit;
1120  limit = order+1; // The 'standard' value away from the boundaries of the knot vector
1121  if (currentKnotSpan<order)
1122  limit -= (order - currentKnotSpan);
1123  if ( order > (knotVector.size() - currentKnotSpan - 2) )
1124  limit -= order - (knotVector.size() - currentKnotSpan - 2);
1125 
1126  // The lowest knot spans that we need values from
1127  unsigned int offset;
1128  offset = std::max((int)(currentKnotSpan - order),0);
1129 
1130  // Evaluate 1d function values (needed for the product rule)
1131  DynamicMatrix<R> values;
1132 
1133  evaluateFunctionFull(in, values, knotVector, order, currentKnotSpan);
1134 
1135  out.resize(knotVector.size()-order-1);
1136  for (size_t j=0; j<out.size(); j++)
1137  out[j] = values[order][j];
1138 
1139  // Evaluate 1d function values of one order lower (needed for the derivative formula)
1140  std::vector<R> lowOrderOneDValues;
1141 
1142  if (order!=0)
1143  {
1144  lowOrderOneDValues.resize(knotVector.size()-(order-1)-1);
1145  for (size_t j=0; j<lowOrderOneDValues.size(); j++)
1146  lowOrderOneDValues[j] = values[order-1][j];
1147  }
1148 
1149  // Evaluate 1d function values of two order lower (needed for the (second) derivative formula)
1150  std::vector<R> lowOrderTwoDValues;
1151 
1152  if (order>1 && evaluateHessian)
1153  {
1154  lowOrderTwoDValues.resize(knotVector.size()-(order-2)-1);
1155  for (size_t j=0; j<lowOrderTwoDValues.size(); j++)
1156  lowOrderTwoDValues[j] = values[order-2][j];
1157  }
1158 
1159  // Evaluate 1d function derivatives
1160  if (evaluateJacobian)
1161  {
1162  outJac.resize(limit);
1163 
1164  if (order==0) // order-zero functions are piecewise constant, hence all derivatives are zero
1165  std::fill(outJac.begin(), outJac.end(), R(0.0));
1166  else
1167  {
1168  for (size_t j=offset; j<offset+limit; j++)
1169  {
1170  R derivativeAddend1 = lowOrderOneDValues[j] / (knotVector[j+order]-knotVector[j]);
1171  R derivativeAddend2 = lowOrderOneDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]);
1172  // The two previous terms may evaluate as 0/0. This is to be interpreted as 0.
1173  if (std::isnan(derivativeAddend1))
1174  derivativeAddend1 = 0;
1175  if (std::isnan(derivativeAddend2))
1176  derivativeAddend2 = 0;
1177  outJac[j-offset] = order * ( derivativeAddend1 - derivativeAddend2 );
1178  }
1179  }
1180  }
1181 
1182  // Evaluate 1d function second derivatives
1183  if (evaluateHessian)
1184  {
1185  outHess.resize(limit);
1186 
1187  if (order<2) // order-zero functions are piecewise constant, hence all derivatives are zero
1188  std::fill(outHess.begin(), outHess.end(), R(0.0));
1189  else
1190  {
1191  for (size_t j=offset; j<offset+limit; j++)
1192  {
1193  assert(j+2 < lowOrderTwoDValues.size());
1194  R derivativeAddend1 = lowOrderTwoDValues[j] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order-1]-knotVector[j]);
1195  R derivativeAddend2 = lowOrderTwoDValues[j+1] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order]-knotVector[j+1]);
1196  R derivativeAddend3 = lowOrderTwoDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+order]-knotVector[j+1]);
1197  R derivativeAddend4 = lowOrderTwoDValues[j+2] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+1+order]-knotVector[j+2]);
1198  // The two previous terms may evaluate as 0/0. This is to be interpreted as 0.
1199 
1200  if (std::isnan(derivativeAddend1))
1201  derivativeAddend1 = 0;
1202  if (std::isnan(derivativeAddend2))
1203  derivativeAddend2 = 0;
1204  if (std::isnan(derivativeAddend3))
1205  derivativeAddend3 = 0;
1206  if (std::isnan(derivativeAddend4))
1207  derivativeAddend4 = 0;
1208  outHess[j-offset] = order * (order-1) * ( derivativeAddend1 - derivativeAddend2 -derivativeAddend3 + derivativeAddend4 );
1209  }
1210  }
1211  }
1212  }
1213 
1214 
1216  std::array<unsigned int, dim> order_;
1217 
1219  std::array<std::vector<double>, dim> knotVectors_;
1220 
1222  std::array<unsigned,dim> elements_;
1223 
1225 };
1226 
1227 
1228 
1229 template<typename GV, typename TP>
1230 class BSplineNode :
1231  public LeafBasisNode<std::size_t, TP>
1232 {
1233  static const int dim = GV::dimension;
1234 
1235  using Base = LeafBasisNode<std::size_t,TP>;
1236 
1237 public:
1238 
1239  using size_type = std::size_t;
1240  using TreePath = TP;
1241  using Element = typename GV::template Codim<0>::Entity;
1243 
1244  BSplineNode(const TreePath& treePath, const BSplinePreBasis<GV, FlatMultiIndex<std::size_t>>* preBasis) :
1245  Base(treePath),
1246  preBasis_(preBasis),
1247  finiteElement_(*preBasis)
1248  {}
1249 
1251  const Element& element() const
1252  {
1253  return element_;
1254  }
1255 
1261  {
1262  return finiteElement_;
1263  }
1264 
1266  void bind(const Element& e)
1267  {
1268  element_ = e;
1269  auto elementIndex = preBasis_->gridView().indexSet().index(e);
1270  finiteElement_.bind(preBasis_->getIJK(elementIndex,preBasis_->elements_));
1271  this->setSize(finiteElement_.size());
1272  }
1273 
1274 protected:
1275 
1277 
1280 };
1281 
1282 
1283 
1284 template<typename GV, class MI, class TP>
1285 class BSplineNodeIndexSet
1286 {
1287  enum {dim = GV::dimension};
1288 
1289 public:
1290 
1291  using size_type = std::size_t;
1292 
1294  using MultiIndex = MI;
1295 
1297 
1298  using Node = typename PreBasis::template Node<TP>;
1299 
1300  BSplineNodeIndexSet(const PreBasis& preBasis) :
1301  preBasis_(&preBasis)
1302  {}
1303 
1309  void bind(const Node& node)
1310  {
1311  node_ = &node;
1312  // Local degrees of freedom are arranged in a lattice.
1313  // We need the lattice dimensions to be able to compute lattice coordinates from a local index
1314  for (int i=0; i<dim; i++)
1315  localSizes_[i] = node_->finiteElement().size(i);
1316  }
1317 
1320  void unbind()
1321  {
1322  node_ = nullptr;
1323  }
1324 
1327  size_type size() const
1328  {
1329  return node_->finiteElement().size();
1330  }
1331 
1333  template<typename It>
1334  It indices(It it) const
1335  {
1336  for (size_type i = 0, end = size() ; i < end ; ++i, ++it)
1337  {
1338  std::array<unsigned int,dim> localIJK = preBasis_->getIJK(i, localSizes_);
1339 
1340  const auto currentKnotSpan = node_->finiteElement().currentKnotSpan_;
1341  const auto order = preBasis_->order_;
1342 
1343  std::array<unsigned int,dim> globalIJK;
1344  for (int i=0; i<dim; i++)
1345  globalIJK[i] = std::max((int)currentKnotSpan[i] - (int)order[i], 0) + localIJK[i]; // needs to be a signed type!
1346 
1347  // Make one global flat index from the globalIJK tuple
1348  size_type globalIdx = globalIJK[dim-1];
1349 
1350  for (int i=dim-2; i>=0; i--)
1351  globalIdx = globalIdx * preBasis_->size(i) + globalIJK[i];
1352 
1353  *it = {{globalIdx}};
1354  }
1355  return it;
1356  }
1357 
1358 protected:
1360 
1361  const Node* node_;
1362 
1363  std::array<unsigned int, dim> localSizes_;
1364 };
1365 
1366 
1367 
1368 // *****************************************************************************
1369 // This is the actual global basis implementation based on the reusable parts.
1370 // *****************************************************************************
1371 
1378 template<typename GV>
1380 
1381 
1382 } // namespace Functions
1383 
1384 } // namespace Dune
1385 
1386 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
const PreBasis * preBasis_
Definition: bsplinebasis.hh:1359
std::array< unsigned int, dim > localSizes_
Definition: bsplinebasis.hh:1363
typename PreBasis::template Node< TP > Node
Definition: bsplinebasis.hh:1298
void bind(const Element &e)
Bind to element.
Definition: bsplinebasis.hh:1266
typename GV::template Codim< 0 >::Entity Element
Definition: bsplinebasis.hh:1241
const size_type offset_
Definition: nodes.hh:40
IndexSet< TP > indexSet() const
Create tree node index set with given root tree path.
Definition: bsplinebasis.hh:719
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: bsplinebasis.hh:1294
std::size_t size() const
number of coefficients
Definition: bsplinebasis.hh:314
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: bsplinebasis.hh:682
Node< TP > node(const TP &tp) const
Create tree node with given root tree path.
Definition: bsplinebasis.hh:704
LocalBasis class in the sense of dune-localfunctions, presenting the restriction of a B-spline patch ...
Definition: bsplinebasis.hh:44
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:592
size_type size(const SizePrefix prefix) const
Return number of possible values for next position in multi index.
Definition: bsplinebasis.hh:725
void bind(const std::array< unsigned, dim > &elementIdx)
Bind LocalFiniteElement to a specific knot span of the spline patch.
Definition: bsplinebasis.hh:387
const BSplineLocalCoefficients< dim > & localCoefficients() const
Hand out a LocalCoefficients object.
Definition: bsplinebasis.hh:426
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:644
std::size_t size_type
Definition: bsplinebasis.hh:1239
BSplineNodeIndexSet(const PreBasis &preBasis)
Definition: bsplinebasis.hh:1300
std::array< unsigned, dim > elements_
Number of grid elements in the different coordinate directions.
Definition: bsplinebasis.hh:1222
LocalBasisTraits< D, dim, Dune::FieldVector< D, dim >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, dim > > Traits
export type traits for function signature
Definition: bsplinebasis.hh:54
void init(const std::array< unsigned, dim > &sizes)
Definition: bsplinebasis.hh:255
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition: bsplinebasis.hh:738
FiniteElement finiteElement_
Definition: bsplinebasis.hh:1278
GeometryType type() const
Return the reference element that the local finite element is defined on (here, a hypercube) ...
Definition: bsplinebasis.hh:448
unsigned int size() const
Total number of B-spline basis functions.
Definition: bsplinebasis.hh:747
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:96
unsigned int size(size_t d) const
Number of shape functions in one direction.
Definition: bsplinebasis.hh:756
Local interpolation in the sense of dune-localfunctions, for the B-spline basis on tensor-product gri...
Definition: bsplinebasis.hh:338
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:42
void unbind()
Unbind the view.
Definition: bsplinebasis.hh:1320
Pre-basis for B-spline basis.
Definition: bsplinebasis.hh:32
void evaluateJacobian(const FieldVector< typename GV::ctype, dim > &in, std::vector< FieldMatrix< R, 1, dim > > &out, const std::array< unsigned, dim > &currentKnotSpan) const
Evaluate Jacobian of all B-spline basis functions.
Definition: bsplinebasis.hh:794
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: bsplinebasis.hh:688
const BSplinePreBasis< GV, FlatMultiIndex< std::size_t > > * preBasis_
Definition: bsplinebasis.hh:1276
Definition: nodes.hh:189
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: bsplinebasis.hh:343
std::array< unsigned int, dim > order_
Order of the B-spline for each space dimension.
Definition: bsplinebasis.hh:1216
void bind(const Node &node)
Bind the view to a grid element.
Definition: bsplinebasis.hh:1309
size_type size() const
Size of subtree rooted in this node (element-local)
Definition: bsplinebasis.hh:1327
A multi-index class with only one level.
Definition: flatmultiindex.hh:31
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:991
Definition: bsplinebasis.hh:482
std::array< std::vector< double >, dim > knotVectors_
The knot vectors, one for each space dimension.
Definition: bsplinebasis.hh:1219
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:1110
Attaches a shape function to an entity.
Definition: bsplinebasis.hh:176
LocalFiniteElement in the sense of dune-localfunctions, for the B-spline basis on tensor-product grid...
Definition: bsplinebasis.hh:29
const Node * node_
Definition: bsplinebasis.hh:1361
unsigned int order() const
Polynomial order of the shape functions.
Definition: bsplinebasis.hh:139
BSplineLocalBasis< GV, R > localBasis_
Definition: bsplinebasis.hh:469
void evaluateFunction(const FieldVector< typename GV::ctype, dim > &in, std::vector< FieldVector< R, 1 > > &out, const std::array< unsigned, dim > &currentKnotSpan) const
Evaluate all B-spline basis functions at a given point.
Definition: bsplinebasis.hh:763
const BSplinePreBasis< GV, FlatMultiIndex< std::size_t > > & preBasis_
Definition: bsplinebasis.hh:467
unsigned size() const
Number of shape functions in this finite element.
Definition: bsplinebasis.hh:438
const BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > & localInterpolation() const
Hand out a LocalInterpolation object.
Definition: bsplinebasis.hh:432
void evaluateJacobian(const FieldVector< D, dim > &in, std::vector< FieldMatrix< D, 1, dim > > &out) const
Evaluate Jacobian of all shape functions.
Definition: bsplinebasis.hh:81
GV GridView
The grid view that the FE space is defined on.
Definition: bsplinebasis.hh:557
Definition: bsplinebasis.hh:479
void evaluateFunction(const FieldVector< D, dim > &in, std::vector< FieldVector< R, 1 > > &out) const
Evaluate all shape functions.
Definition: bsplinebasis.hh:69
std::size_t size_type
Definition: bsplinebasis.hh:1291
std::array< unsigned, dim > currentKnotSpan_
Definition: bsplinebasis.hh:474
GridView gridView_
Definition: bsplinebasis.hh:1224
void initializeIndices()
Initialize the global indices.
Definition: bsplinebasis.hh:678
Definition: polynomial.hh:7
const LocalKey & localKey(std::size_t i) const
get i&#39;th index
Definition: bsplinebasis.hh:320
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:1067
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition: bsplinebasis.hh:1260
BSplineLocalBasis(const BSplinePreBasis< GV, FlatMultiIndex< std::size_t >> &preBasis, const BSplineLocalFiniteElement< GV, R > &lFE)
Constructor with a given B-spline patch.
Definition: bsplinebasis.hh:60
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:1010
const BSplineLocalBasis< GV, R > & localBasis() const
Hand out a LocalBasis object.
Definition: bsplinebasis.hh:420
Dune::ReservedVector< size_type, 1 > SizePrefix
Definition: bsplinebasis.hh:569
std::size_t size() const
Return the number of basis functions on the current knot span.
Definition: bsplinebasis.hh:146
BSplineLocalFiniteElement(const BSplinePreBasis< GV, FlatMultiIndex< std::size_t >> &preBasis)
Constructor with a given B-spline basis.
Definition: bsplinebasis.hh:376
BSplineLocalCoefficients< dim > localCoefficients_
Definition: bsplinebasis.hh:470
BSplineNode(const TreePath &treePath, const BSplinePreBasis< GV, FlatMultiIndex< std::size_t >> *preBasis)
Definition: bsplinebasis.hh:1244
Element element_
Definition: bsplinebasis.hh:1279
TP TreePath
Definition: bsplinebasis.hh:1240
It indices(It it) const
Maps from subtree index set [0..size-1] to a globally unique multi index in global basis...
Definition: bsplinebasis.hh:1334
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 > &currentKnotSpan) const
Evaluate Derivatives of all B-spline basis functions.
Definition: bsplinebasis.hh:895
size_type dimension() const
Get the total dimension of the space spanned by this basis.
Definition: bsplinebasis.hh:732
LocalFiniteElementTraits< BSplineLocalBasis< GV, R >, BSplineLocalCoefficients< dim >, BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > > Traits
Export various types related to this LocalFiniteElement.
Definition: bsplinebasis.hh:372
unsigned int size(int i) const
Number of degrees of freedom for one coordinate direction.
Definition: bsplinebasis.hh:456
BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > localInterpolation_
Definition: bsplinebasis.hh:471
const Element & element() const
Return current element, throw if unbound.
Definition: bsplinebasis.hh:1251