dune-functions  2.7.0
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>
24 
25 namespace Dune
26 {
27 namespace Functions {
28 
29 // A maze of dependencies between the different parts of this. We need a few forward declarations
30 template<typename GV, typename R, typename MI>
32 
33 template<typename GV, class MI>
35 
36 
45 template<class GV, class R, class MI>
47 {
48  friend class BSplineLocalFiniteElement<GV,R,MI>;
49 
50  typedef typename GV::ctype D;
51  enum {dim = GV::dimension};
52 public:
53 
55  typedef LocalBasisTraits<D,dim,FieldVector<D,dim>,R,1,FieldVector<R,1>,
56  FieldMatrix<R,1,dim> > Traits;
57 
64  : preBasis_(preBasis),
65  lFE_(lFE)
66  {}
67 
71  void evaluateFunction (const FieldVector<D,dim>& in,
72  std::vector<FieldVector<R,1> >& out) const
73  {
74  FieldVector<D,dim> globalIn = offset_;
75  scaling_.umv(in,globalIn);
76 
77  preBasis_.evaluateFunction(globalIn, out, lFE_.currentKnotSpan_);
78  }
79 
83  void evaluateJacobian (const FieldVector<D,dim>& in,
84  std::vector<FieldMatrix<D,1,dim> >& out) const
85  {
86  FieldVector<D,dim> globalIn = offset_;
87  scaling_.umv(in,globalIn);
88 
89  preBasis_.evaluateJacobian(globalIn, out, lFE_.currentKnotSpan_);
90 
91  for (size_t i=0; i<out.size(); i++)
92  for (int j=0; j<dim; j++)
93  out[i][0][j] *= scaling_[j][j];
94  }
95 
97  template<size_t k>
98  inline void evaluate (const typename std::array<int,k>& directions,
99  const typename Traits::DomainType& in,
100  std::vector<typename Traits::RangeType>& out) const
101  {
102  switch(k)
103  {
104  case 0:
105  evaluateFunction(in, out);
106  break;
107  case 1:
108  {
109  FieldVector<D,dim> globalIn = offset_;
110  scaling_.umv(in,globalIn);
111 
112  preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
113 
114  for (size_t i=0; i<out.size(); i++)
115  out[i][0] *= scaling_[directions[0]][directions[0]];
116  break;
117  }
118  case 2:
119  {
120  FieldVector<D,dim> globalIn = offset_;
121  scaling_.umv(in,globalIn);
122 
123  preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
124 
125  for (size_t i=0; i<out.size(); i++)
126  out[i][0] *= scaling_[directions[0]][directions[0]]*scaling_[directions[1]][directions[1]];
127  break;
128  }
129  default:
130  DUNE_THROW(NotImplemented, "B-Spline derivatives of order " << k << " not implemented yet!");
131  }
132  }
133 
141  unsigned int order () const
142  {
143  return *std::max_element(preBasis_.order_.begin(), preBasis_.order_.end());
144  }
145 
148  std::size_t size() const
149  {
150  return lFE_.size();
151  }
152 
153 private:
154  const BSplinePreBasis<GV,MI>& preBasis_;
155 
157 
158  // Coordinates in a single knot span differ from coordinates on the B-spline patch
159  // by an affine transformation. This transformation is stored in offset_ and scaling_.
160  FieldVector<D,dim> offset_;
161  DiagonalMatrix<D,dim> scaling_;
162 };
163 
177 template<int dim>
179 {
180  // Return i as a d-digit number in the (k+1)-nary system
181  std::array<unsigned int,dim> multiindex (unsigned int i) const
182  {
183  std::array<unsigned int,dim> alpha;
184  for (int j=0; j<dim; j++)
185  {
186  alpha[j] = i % sizes_[j];
187  i = i/sizes_[j];
188  }
189  return alpha;
190  }
191 
193  void setup1d(std::vector<unsigned int>& subEntity)
194  {
195  if (sizes_[0]==1)
196  {
197  subEntity[0] = 0;
198  return;
199  }
200 
201  /* edge and vertex numbering
202  0----0----1
203  */
204  unsigned lastIndex=0;
205  subEntity[lastIndex++] = 0; // corner 0
206  for (unsigned i = 0; i < sizes_[0] - 2; ++i)
207  subEntity[lastIndex++] = 0; // inner dofs of element (0)
208 
209  subEntity[lastIndex++] = 1; // corner 1
210 
211  assert(size()==lastIndex);
212  }
213 
214  void setup2d(std::vector<unsigned int>& subEntity)
215  {
216  unsigned lastIndex=0;
217 
218  // LocalKey: entity number , entity codim, dof indices within each entity
219  /* edge and vertex numbering
220  2----3----3
221  | |
222  | |
223  0 1
224  | |
225  | |
226  0----2----1
227  */
228 
229  // lower edge (2)
230  subEntity[lastIndex++] = 0; // corner 0
231  for (unsigned i = 0; i < sizes_[0]-2; ++i)
232  subEntity[lastIndex++] = 2; // inner dofs of lower edge (2)
233 
234  subEntity[lastIndex++] = 1; // corner 1
235 
236  // iterate from bottom to top over inner edge dofs
237  for (unsigned e = 0; e < sizes_[1]-2; ++e)
238  {
239  subEntity[lastIndex++] = 0; // left edge (0)
240  for (unsigned i = 0; i < sizes_[0]-2; ++i)
241  subEntity[lastIndex++] = 0; // face dofs
242  subEntity[lastIndex++] = 1; // right edge (1)
243  }
244 
245  // upper edge (3)
246  subEntity[lastIndex++] = 2; // corner 2
247  for (unsigned i = 0; i < sizes_[0]-2; ++i)
248  subEntity[lastIndex++] = 3; // inner dofs of upper edge (3)
249 
250  subEntity[lastIndex++] = 3; // corner 3
251 
252  assert(size()==lastIndex);
253  }
254 
255 
256 public:
257  void init(const std::array<unsigned,dim>& sizes)
258  {
259  sizes_ = sizes;
260 
261  li_.resize(size());
262 
263  // Set up array of codimension-per-dof-number
264  std::vector<unsigned int> codim(li_.size());
265 
266  for (std::size_t i=0; i<codim.size(); i++)
267  {
268  codim[i] = 0;
269  // Codimension gets increased by 1 for each coordinate direction
270  // where dof is on boundary
271  std::array<unsigned int,dim> mIdx = multiindex(i);
272  for (int j=0; j<dim; j++)
273  if (mIdx[j]==0 or mIdx[j]==sizes[j]-1)
274  codim[i]++;
275  }
276 
277  // Set up index vector (the index of the dof in the set of dofs of a given subentity)
278  // Algorithm: the 'index' has the same ordering as the dof number 'i'.
279  // To make it consecutive we interpret 'i' in the (k+1)-adic system, omit all digits
280  // that correspond to axes where the dof is on the element boundary, and transform the
281  // rest to the (k-1)-adic system.
282  std::vector<unsigned int> index(size());
283 
284  for (std::size_t i=0; i<index.size(); i++)
285  {
286  index[i] = 0;
287 
288  std::array<unsigned int,dim> mIdx = multiindex(i);
289 
290  for (int j=dim-1; j>=0; j--)
291  if (mIdx[j]>0 and mIdx[j]<sizes[j]-1)
292  index[i] = (sizes[j]-1)*index[i] + (mIdx[j]-1);
293  }
294 
295  // Set up entity and dof numbers for each (supported) dimension separately
296  std::vector<unsigned int> subEntity(li_.size());
297 
298  if (subEntity.size() > 0)
299  {
300  if (dim==1) {
301 
302  setup1d(subEntity);
303 
304  } else if (dim==2 and sizes_[0]>1 and sizes_[1]>1) {
305 
306  setup2d(subEntity);
307 
308  }
309  }
310 
311  for (size_t i=0; i<li_.size(); i++)
312  li_[i] = LocalKey(subEntity[i], codim[i], index[i]);
313  }
314 
316  std::size_t size () const
317  {
318  return std::accumulate(sizes_.begin(), sizes_.end(), 1, std::multiplies<unsigned int>());
319  }
320 
322  const LocalKey& localKey (std::size_t i) const
323  {
324  return li_[i];
325  }
326 
327 private:
328 
329  // Number of shape functions on this element per coordinate direction
330  std::array<unsigned, dim> sizes_;
331 
332  std::vector<LocalKey> li_;
333 };
334 
339 template<int dim, class LB>
341 {
342 public:
344  template<typename F, typename C>
345  void interpolate (const F& f, std::vector<C>& out) const
346  {
347  DUNE_THROW(NotImplemented, "BSplineLocalInterpolation::interpolate");
348  }
349 
350 };
351 
362 template<class GV, class R, class MI>
363 class BSplineLocalFiniteElement
364 {
365  typedef typename GV::ctype D;
366  enum {dim = GV::dimension};
367  friend class BSplineLocalBasis<GV,R,MI>;
368 public:
369 
372  typedef LocalFiniteElementTraits<BSplineLocalBasis<GV,R,MI>,
375 
379  : preBasis_(preBasis),
380  localBasis_(preBasis,*this)
381  {}
382 
389  void bind(const std::array<unsigned,dim>& elementIdx)
390  {
391  /* \todo In the long run we need to precompute a table for this */
392  for (size_t i=0; i<elementIdx.size(); i++)
393  {
394  currentKnotSpan_[i] = 0;
395 
396  // Skip over degenerate knot spans
397  while (preBasis_.knotVectors_[i][currentKnotSpan_[i]+1] < preBasis_.knotVectors_[i][currentKnotSpan_[i]]+1e-8)
398  currentKnotSpan_[i]++;
399 
400  for (size_t j=0; j<elementIdx[i]; j++)
401  {
402  currentKnotSpan_[i]++;
403 
404  // Skip over degenerate knot spans
405  while (preBasis_.knotVectors_[i][currentKnotSpan_[i]+1] < preBasis_.knotVectors_[i][currentKnotSpan_[i]]+1e-8)
406  currentKnotSpan_[i]++;
407  }
408 
409  // Compute the geometric transformation from knotspan-local to global coordinates
410  localBasis_.offset_[i] = preBasis_.knotVectors_[i][currentKnotSpan_[i]];
411  localBasis_.scaling_[i][i] = preBasis_.knotVectors_[i][currentKnotSpan_[i]+1] - preBasis_.knotVectors_[i][currentKnotSpan_[i]];
412  }
413 
414  // Set up the LocalCoefficients object
415  std::array<unsigned int, dim> sizes;
416  for (size_t i=0; i<dim; i++)
417  sizes[i] = size(i);
418  localCoefficients_.init(sizes);
419  }
420 
423  {
424  return localBasis_;
425  }
426 
429  {
430  return localCoefficients_;
431  }
432 
435  {
436  return localInterpolation_;
437  }
438 
440  unsigned size () const
441  {
442  std::size_t r = 1;
443  for (int i=0; i<dim; i++)
444  r *= size(i);
445  return r;
446  }
447 
450  GeometryType type () const
451  {
452  return GeometryTypes::cube(dim);
453  }
454 
455 //private:
456 
458  unsigned int size(int i) const
459  {
460  const auto& order = preBasis_.order_;
461  unsigned int r = order[i]+1; // The 'normal' value
462  if (currentKnotSpan_[i]<order[i]) // Less near the left end of the knot vector
463  r -= (order[i] - currentKnotSpan_[i]);
464  if ( order[i] > (preBasis_.knotVectors_[i].size() - currentKnotSpan_[i] - 2) )
465  r -= order[i] - (preBasis_.knotVectors_[i].size() - currentKnotSpan_[i] - 2);
466  return r;
467  }
468 
470 
474 
475  // The knot span we are bound to
476  std::array<unsigned,dim> currentKnotSpan_;
477 };
478 
479 
480 template<typename GV, typename MI>
482 
483 template<typename GV, class MI>
485 
496 template<typename GV, class MI>
497 class BSplinePreBasis
498 {
499  static const int dim = GV::dimension;
500 
502  class MultiDigitCounter
503  {
504  public:
505 
509  MultiDigitCounter(const std::array<unsigned int,dim>& limits)
510  : limits_(limits)
511  {
512  std::fill(counter_.begin(), counter_.end(), 0);
513  }
514 
516  MultiDigitCounter& operator++()
517  {
518  for (int i=0; i<dim; i++)
519  {
520  ++counter_[i];
521 
522  // no overflow?
523  if (counter_[i] < limits_[i])
524  break;
525 
526  counter_[i] = 0;
527  }
528  return *this;
529  }
530 
532  const unsigned int& operator[](int i) const
533  {
534  return counter_[i];
535  }
536 
538  unsigned int cycle() const
539  {
540  unsigned int r = 1;
541  for (int i=0; i<dim; i++)
542  r *= limits_[i];
543  return r;
544  }
545 
546  private:
547 
549  const std::array<unsigned int,dim> limits_;
550 
552  std::array<unsigned int,dim> counter_;
553 
554  };
555 
556 public:
557 
559  using GridView = GV;
560  using size_type = std::size_t;
561 
563 
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 
593  const std::vector<double>& knotVector,
594  unsigned int order,
595  bool makeOpen = true)
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 
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),
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 
696  Node makeNode() const
697  {
698  return Node{this};
699  }
700 
708  {
709  return IndexSet{*this};
710  }
711 
713  size_type size(const SizePrefix prefix) const
714  {
715  assert(prefix.size() == 0 || prefix.size() == 1);
716  return (prefix.size() == 0) ? size() : 0;
717  }
718 
721  {
722  return size();
723  }
724 
727  {
728  size_type result = 1;
729  for (int i=0; i<dim; i++)
730  result *= order_[i]+1;
731  return result;
732  }
733 
735  unsigned int size () const
736  {
737  unsigned int result = 1;
738  for (size_t i=0; i<dim; i++)
739  result *= size(i);
740  return result;
741  }
742 
744  unsigned int size (size_t d) const
745  {
746  return knotVectors_[d].size() - order_[d] - 1;
747  }
748 
751  void evaluateFunction (const FieldVector<typename GV::ctype,dim>& in,
752  std::vector<FieldVector<R,1> >& out,
753  const std::array<unsigned,dim>& currentKnotSpan) const
754  {
755  // Evaluate
756  std::array<std::vector<R>, dim> oneDValues;
757 
758  for (size_t i=0; i<dim; i++)
759  evaluateFunction(in[i], oneDValues[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
760 
761  std::array<unsigned int, dim> limits;
762  for (int i=0; i<dim; i++)
763  limits[i] = oneDValues[i].size();
764 
765  MultiDigitCounter ijkCounter(limits);
766 
767  out.resize(ijkCounter.cycle());
768 
769  for (size_t i=0; i<out.size(); i++, ++ijkCounter)
770  {
771  out[i] = R(1.0);
772  for (size_t j=0; j<dim; j++)
773  out[i] *= oneDValues[j][ijkCounter[j]];
774  }
775  }
776 
782  void evaluateJacobian (const FieldVector<typename GV::ctype,dim>& in,
783  std::vector<FieldMatrix<R,1,dim> >& out,
784  const std::array<unsigned,dim>& currentKnotSpan) const
785  {
786  // How many shape functions to we have in each coordinate direction?
787  std::array<unsigned int, dim> limits;
788  for (int i=0; i<dim; i++)
789  {
790  limits[i] = order_[i]+1; // The 'standard' value away from the boundaries of the knot vector
791  if (currentKnotSpan[i]<order_[i])
792  limits[i] -= (order_[i] - currentKnotSpan[i]);
793  if ( order_[i] > (knotVectors_[i].size() - currentKnotSpan[i] - 2) )
794  limits[i] -= order_[i] - (knotVectors_[i].size() - currentKnotSpan[i] - 2);
795  }
796 
797  // The lowest knot spans that we need values from
798  std::array<unsigned int, dim> offset;
799  for (int i=0; i<dim; i++)
800  offset[i] = std::max((int)(currentKnotSpan[i] - order_[i]),0);
801 
802  // Evaluate 1d function values (needed for the product rule)
803  std::array<std::vector<R>, dim> oneDValues;
804 
805  // Evaluate 1d function values of one order lower (needed for the derivative formula)
806  std::array<std::vector<R>, dim> lowOrderOneDValues;
807 
808  std::array<DynamicMatrix<R>, dim> values;
809 
810  for (size_t i=0; i<dim; i++)
811  {
812  evaluateFunctionFull(in[i], values[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
813  oneDValues[i].resize(knotVectors_[i].size()-order_[i]-1);
814  for (size_t j=0; j<oneDValues[i].size(); j++)
815  oneDValues[i][j] = values[i][order_[i]][j];
816 
817  if (order_[i]!=0)
818  {
819  lowOrderOneDValues[i].resize(knotVectors_[i].size()-(order_[i]-1)-1);
820  for (size_t j=0; j<lowOrderOneDValues[i].size(); j++)
821  lowOrderOneDValues[i][j] = values[i][order_[i]-1][j];
822  }
823  }
824 
825 
826  // Evaluate 1d function derivatives
827  std::array<std::vector<R>, dim> oneDDerivatives;
828  for (size_t i=0; i<dim; i++)
829  {
830  oneDDerivatives[i].resize(limits[i]);
831 
832  if (order_[i]==0) // order-zero functions are piecewise constant, hence all derivatives are zero
833  std::fill(oneDDerivatives[i].begin(), oneDDerivatives[i].end(), R(0.0));
834  else
835  {
836  for (size_t j=offset[i]; j<offset[i]+limits[i]; j++)
837  {
838  R derivativeAddend1 = lowOrderOneDValues[i][j] / (knotVectors_[i][j+order_[i]]-knotVectors_[i][j]);
839  R derivativeAddend2 = lowOrderOneDValues[i][j+1] / (knotVectors_[i][j+order_[i]+1]-knotVectors_[i][j+1]);
840  // The two previous terms may evaluate as 0/0. This is to be interpreted as 0.
841  if (std::isnan(derivativeAddend1))
842  derivativeAddend1 = 0;
843  if (std::isnan(derivativeAddend2))
844  derivativeAddend2 = 0;
845  oneDDerivatives[i][j-offset[i]] = order_[i] * ( derivativeAddend1 - derivativeAddend2 );
846  }
847  }
848  }
849 
850  // Working towards computing only the parts that we really need:
851  // Let's copy them out into a separate array
852  std::array<std::vector<R>, dim> oneDValuesShort;
853 
854  for (int i=0; i<dim; i++)
855  {
856  oneDValuesShort[i].resize(limits[i]);
857 
858  for (size_t j=0; j<limits[i]; j++)
859  oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
860  }
861 
862 
863 
864  // Set up a multi-index to go from consecutive indices to integer coordinates
865  MultiDigitCounter ijkCounter(limits);
866 
867  out.resize(ijkCounter.cycle());
868 
869  // Complete Jacobian is given by the product rule
870  for (size_t i=0; i<out.size(); i++, ++ijkCounter)
871  for (int j=0; j<dim; j++)
872  {
873  out[i][0][j] = 1.0;
874  for (int k=0; k<dim; k++)
875  out[i][0][j] *= (j==k) ? oneDDerivatives[k][ijkCounter[k]]
876  : oneDValuesShort[k][ijkCounter[k]];
877  }
878 
879  }
880 
882  template <size_type k>
883  void evaluate(const typename std::array<int,k>& directions,
884  const FieldVector<typename GV::ctype,dim>& in,
885  std::vector<FieldVector<R,1> >& out,
886  const std::array<unsigned,dim>& currentKnotSpan) const
887  {
888  if (k != 1 && k != 2)
889  DUNE_THROW(RangeError, "Differentiation order greater than 2 is not supported!");
890 
891  // Evaluate 1d function values (needed for the product rule)
892  std::array<std::vector<R>, dim> oneDValues;
893  std::array<std::vector<R>, dim> oneDDerivatives;
894  std::array<std::vector<R>, dim> oneDSecondDerivatives;
895 
896  // Evaluate 1d function derivatives
897  if (k==1)
898  for (size_t i=0; i<dim; i++)
899  evaluateAll(in[i], oneDValues[i], true, oneDDerivatives[i], false, oneDSecondDerivatives[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
900  else
901  for (size_t i=0; i<dim; i++)
902  evaluateAll(in[i], oneDValues[i], true, oneDDerivatives[i], true, oneDSecondDerivatives[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
903 
904  // The lowest knot spans that we need values from
905  std::array<unsigned int, dim> offset;
906  for (int i=0; i<dim; i++)
907  offset[i] = std::max((int)(currentKnotSpan[i] - order_[i]),0);
908 
909  // Set up a multi-index to go from consecutive indices to integer coordinates
910  std::array<unsigned int, dim> limits;
911  for (int i=0; i<dim; i++)
912  {
913  // In a proper implementation, the following line would do
914  //limits[i] = oneDValues[i].size();
915  limits[i] = order_[i]+1; // The 'standard' value away from the boundaries of the knot vector
916  if (currentKnotSpan[i]<order_[i])
917  limits[i] -= (order_[i] - currentKnotSpan[i]);
918  if ( order_[i] > (knotVectors_[i].size() - currentKnotSpan[i] - 2) )
919  limits[i] -= order_[i] - (knotVectors_[i].size() - currentKnotSpan[i] - 2);
920  }
921 
922  // Working towards computing only the parts that we really need:
923  // Let's copy them out into a separate array
924  std::array<std::vector<R>, dim> oneDValuesShort;
925 
926  for (int i=0; i<dim; i++)
927  {
928  oneDValuesShort[i].resize(limits[i]);
929 
930  for (size_t j=0; j<limits[i]; j++)
931  oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
932  }
933 
934 
935  MultiDigitCounter ijkCounter(limits);
936 
937  out.resize(ijkCounter.cycle());
938 
939  if (k == 1)
940  {
941  // Complete Jacobian is given by the product rule
942  for (size_t i=0; i<out.size(); i++, ++ijkCounter)
943  {
944  out[i][0] = 1.0;
945  for (int l=0; l<dim; l++)
946  out[i][0] *= (directions[0]==l) ? oneDDerivatives[l][ijkCounter[l]]
947  : oneDValuesShort[l][ijkCounter[l]];
948  }
949  }
950 
951  if (k == 2)
952  {
953  // Complete derivation by deriving the tensor product
954  for (size_t i=0; i<out.size(); i++, ++ijkCounter)
955  {
956  out[i][0] = 1.0;
957  for (int j=0; j<dim; j++)
958  {
959  if (directions[0] != directions[1]) //derivation in two different variables
960  if (directions[0] == j || directions[1] == j) //the spline has to be derived (once) in this direction
961  out[i][0] *= oneDDerivatives[j][ijkCounter[j]];
962  else //no derivation in this direction
963  out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
964  else //spline is derived two times in the same direction
965  if (directions[0] == j) //the spline is derived two times in this direction
966  out[i][0] *= oneDSecondDerivatives[j][ijkCounter[j]];
967  else //no derivation in this direction
968  out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
969  }
970  }
971  }
972  }
973 
974 
979  static std::array<unsigned int,dim> getIJK(typename GridView::IndexSet::IndexType idx, std::array<unsigned int,dim> elements)
980  {
981  std::array<unsigned,dim> result;
982  for (int i=0; i<dim; i++)
983  {
984  result[i] = idx%elements[i];
985  idx /= elements[i];
986  }
987  return result;
988  }
989 
998  static void evaluateFunction (const typename GV::ctype& in, std::vector<R>& out,
999  const std::vector<R>& knotVector,
1000  unsigned int order,
1001  unsigned int currentKnotSpan)
1002  {
1003  std::size_t outSize = order+1; // The 'standard' value away from the boundaries of the knot vector
1004  if (currentKnotSpan<order) // Less near the left end of the knot vector
1005  outSize -= (order - currentKnotSpan);
1006  if ( order > (knotVector.size() - currentKnotSpan - 2) )
1007  outSize -= order - (knotVector.size() - currentKnotSpan - 2);
1008  out.resize(outSize);
1009 
1010  // It's not really a matrix that is needed here, a plain 2d array would do
1011  DynamicMatrix<R> N(order+1, knotVector.size()-1);
1012 
1013  // The text books on splines use the following geometric condition here to fill the array N
1014  // (see for example Cottrell, Hughes, Bazilevs, Formula (2.1). However, this condition
1015  // only works if splines are never evaluated exactly on the knots.
1016  //
1017  // for (size_t i=0; i<knotVector.size()-1; i++)
1018  // N[0][i] = (knotVector[i] <= in) and (in < knotVector[i+1]);
1019  for (size_t i=0; i<knotVector.size()-1; i++)
1020  N[0][i] = (i == currentKnotSpan);
1021 
1022  for (size_t r=1; r<=order; r++)
1023  for (size_t i=0; i<knotVector.size()-r-1; i++)
1024  {
1025  R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
1026  ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
1027  : 0;
1028  R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
1029  ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
1030  : 0;
1031  N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
1032  }
1033 
1038  for (size_t i=0; i<out.size(); i++) {
1039  out[i] = N[order][std::max((int)(currentKnotSpan - order),0) + i];
1040  }
1041  }
1042 
1055  static void evaluateFunctionFull(const typename GV::ctype& in,
1056  DynamicMatrix<R>& out,
1057  const std::vector<R>& knotVector,
1058  unsigned int order,
1059  unsigned int currentKnotSpan)
1060  {
1061  // It's not really a matrix that is needed here, a plain 2d array would do
1062  DynamicMatrix<R>& N = out;
1063 
1064  N.resize(order+1, knotVector.size()-1);
1065 
1066  // The text books on splines use the following geometric condition here to fill the array N
1067  // (see for example Cottrell, Hughes, Bazilevs, Formula (2.1). However, this condition
1068  // only works if splines are never evaluated exactly on the knots.
1069  //
1070  // for (size_t i=0; i<knotVector.size()-1; i++)
1071  // N[0][i] = (knotVector[i] <= in) and (in < knotVector[i+1]);
1072  for (size_t i=0; i<knotVector.size()-1; i++)
1073  N[0][i] = (i == currentKnotSpan);
1074 
1075  for (size_t r=1; r<=order; r++)
1076  for (size_t i=0; i<knotVector.size()-r-1; i++)
1077  {
1078  R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
1079  ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
1080  : 0;
1081  R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
1082  ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
1083  : 0;
1084  N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
1085  }
1086  }
1087 
1088 
1098  static void evaluateAll(const typename GV::ctype& in,
1099  std::vector<R>& out,
1100  bool evaluateJacobian, std::vector<R>& outJac,
1101  bool evaluateHessian, std::vector<R>& outHess,
1102  const std::vector<R>& knotVector,
1103  unsigned int order,
1104  unsigned int currentKnotSpan)
1105  {
1106  // How many shape functions to we have in each coordinate direction?
1107  unsigned int limit;
1108  limit = order+1; // The 'standard' value away from the boundaries of the knot vector
1109  if (currentKnotSpan<order)
1110  limit -= (order - currentKnotSpan);
1111  if ( order > (knotVector.size() - currentKnotSpan - 2) )
1112  limit -= order - (knotVector.size() - currentKnotSpan - 2);
1113 
1114  // The lowest knot spans that we need values from
1115  unsigned int offset;
1116  offset = std::max((int)(currentKnotSpan - order),0);
1117 
1118  // Evaluate 1d function values (needed for the product rule)
1119  DynamicMatrix<R> values;
1120 
1121  evaluateFunctionFull(in, values, knotVector, order, currentKnotSpan);
1122 
1123  out.resize(knotVector.size()-order-1);
1124  for (size_t j=0; j<out.size(); j++)
1125  out[j] = values[order][j];
1126 
1127  // Evaluate 1d function values of one order lower (needed for the derivative formula)
1128  std::vector<R> lowOrderOneDValues;
1129 
1130  if (order!=0)
1131  {
1132  lowOrderOneDValues.resize(knotVector.size()-(order-1)-1);
1133  for (size_t j=0; j<lowOrderOneDValues.size(); j++)
1134  lowOrderOneDValues[j] = values[order-1][j];
1135  }
1136 
1137  // Evaluate 1d function values of two order lower (needed for the (second) derivative formula)
1138  std::vector<R> lowOrderTwoDValues;
1139 
1140  if (order>1 && evaluateHessian)
1141  {
1142  lowOrderTwoDValues.resize(knotVector.size()-(order-2)-1);
1143  for (size_t j=0; j<lowOrderTwoDValues.size(); j++)
1144  lowOrderTwoDValues[j] = values[order-2][j];
1145  }
1146 
1147  // Evaluate 1d function derivatives
1148  if (evaluateJacobian)
1149  {
1150  outJac.resize(limit);
1151 
1152  if (order==0) // order-zero functions are piecewise constant, hence all derivatives are zero
1153  std::fill(outJac.begin(), outJac.end(), R(0.0));
1154  else
1155  {
1156  for (size_t j=offset; j<offset+limit; j++)
1157  {
1158  R derivativeAddend1 = lowOrderOneDValues[j] / (knotVector[j+order]-knotVector[j]);
1159  R derivativeAddend2 = lowOrderOneDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]);
1160  // The two previous terms may evaluate as 0/0. This is to be interpreted as 0.
1161  if (std::isnan(derivativeAddend1))
1162  derivativeAddend1 = 0;
1163  if (std::isnan(derivativeAddend2))
1164  derivativeAddend2 = 0;
1165  outJac[j-offset] = order * ( derivativeAddend1 - derivativeAddend2 );
1166  }
1167  }
1168  }
1169 
1170  // Evaluate 1d function second derivatives
1171  if (evaluateHessian)
1172  {
1173  outHess.resize(limit);
1174 
1175  if (order<2) // order-zero functions are piecewise constant, hence all derivatives are zero
1176  std::fill(outHess.begin(), outHess.end(), R(0.0));
1177  else
1178  {
1179  for (size_t j=offset; j<offset+limit; j++)
1180  {
1181  assert(j+2 < lowOrderTwoDValues.size());
1182  R derivativeAddend1 = lowOrderTwoDValues[j] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order-1]-knotVector[j]);
1183  R derivativeAddend2 = lowOrderTwoDValues[j+1] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order]-knotVector[j+1]);
1184  R derivativeAddend3 = lowOrderTwoDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+order]-knotVector[j+1]);
1185  R derivativeAddend4 = lowOrderTwoDValues[j+2] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+1+order]-knotVector[j+2]);
1186  // The two previous terms may evaluate as 0/0. This is to be interpreted as 0.
1187 
1188  if (std::isnan(derivativeAddend1))
1189  derivativeAddend1 = 0;
1190  if (std::isnan(derivativeAddend2))
1191  derivativeAddend2 = 0;
1192  if (std::isnan(derivativeAddend3))
1193  derivativeAddend3 = 0;
1194  if (std::isnan(derivativeAddend4))
1195  derivativeAddend4 = 0;
1196  outHess[j-offset] = order * (order-1) * ( derivativeAddend1 - derivativeAddend2 -derivativeAddend3 + derivativeAddend4 );
1197  }
1198  }
1199  }
1200  }
1201 
1202 
1204  std::array<unsigned int, dim> order_;
1205 
1207  std::array<std::vector<double>, dim> knotVectors_;
1208 
1210  std::array<unsigned,dim> elements_;
1211 
1213 };
1214 
1215 
1216 
1217 template<typename GV, typename MI>
1218 class BSplineNode :
1219  public LeafBasisNode
1220 {
1221  static const int dim = GV::dimension;
1222 
1223 public:
1224 
1225  using size_type = std::size_t;
1226  using Element = typename GV::template Codim<0>::Entity;
1228 
1230  preBasis_(preBasis),
1231  finiteElement_(*preBasis)
1232  {}
1233 
1235  const Element& element() const
1236  {
1237  return element_;
1238  }
1239 
1245  {
1246  return finiteElement_;
1247  }
1248 
1250  void bind(const Element& e)
1251  {
1252  element_ = e;
1253  auto elementIndex = preBasis_->gridView().indexSet().index(e);
1254  finiteElement_.bind(preBasis_->getIJK(elementIndex,preBasis_->elements_));
1255  this->setSize(finiteElement_.size());
1256  }
1257 
1258 protected:
1259 
1261 
1264 };
1265 
1266 
1267 
1268 template<typename GV, class MI>
1269 class BSplineNodeIndexSet
1270 {
1271  enum {dim = GV::dimension};
1272 
1273 public:
1274 
1275  using size_type = std::size_t;
1276 
1278  using MultiIndex = MI;
1279 
1281 
1283 
1284  BSplineNodeIndexSet(const PreBasis& preBasis) :
1285  preBasis_(&preBasis)
1286  {}
1287 
1293  void bind(const Node& node)
1294  {
1295  node_ = &node;
1296  // Local degrees of freedom are arranged in a lattice.
1297  // We need the lattice dimensions to be able to compute lattice coordinates from a local index
1298  for (int i=0; i<dim; i++)
1299  localSizes_[i] = node_->finiteElement().size(i);
1300  }
1301 
1304  void unbind()
1305  {
1306  node_ = nullptr;
1307  }
1308 
1311  size_type size() const
1312  {
1313  return node_->finiteElement().size();
1314  }
1315 
1317  template<typename It>
1318  It indices(It it) const
1319  {
1320  for (size_type i = 0, end = size() ; i < end ; ++i, ++it)
1321  {
1322  std::array<unsigned int,dim> localIJK = preBasis_->getIJK(i, localSizes_);
1323 
1324  const auto currentKnotSpan = node_->finiteElement().currentKnotSpan_;
1325  const auto order = preBasis_->order_;
1326 
1327  std::array<unsigned int,dim> globalIJK;
1328  for (int i=0; i<dim; i++)
1329  globalIJK[i] = std::max((int)currentKnotSpan[i] - (int)order[i], 0) + localIJK[i]; // needs to be a signed type!
1330 
1331  // Make one global flat index from the globalIJK tuple
1332  size_type globalIdx = globalIJK[dim-1];
1333 
1334  for (int i=dim-2; i>=0; i--)
1335  globalIdx = globalIdx * preBasis_->size(i) + globalIJK[i];
1336 
1337  *it = {{globalIdx}};
1338  }
1339  return it;
1340  }
1341 
1342 protected:
1344 
1345  const Node* node_;
1346 
1347  std::array<unsigned int, dim> localSizes_;
1348 };
1349 
1350 
1351 
1352 namespace BasisFactory {
1353 
1354 namespace Imp {
1355 
1356 class BSplinePreBasisFactory
1357 {
1358 public:
1359  static const std::size_t requiredMultiIndexSize=1;
1360 
1361  BSplinePreBasisFactory(const std::vector<double>& knotVector,
1362  unsigned int order,
1363  bool makeOpen = true)
1364  : knotVector_(knotVector),
1365  order_(order),
1366  makeOpen_(makeOpen)
1367  {}
1368 
1369  template<class MultiIndex, class GridView>
1370  auto makePreBasis(const GridView& gridView) const
1371  {
1372  return BSplinePreBasis<GridView, MultiIndex>(gridView, knotVector_, order_, makeOpen_);
1373  }
1374 
1375 private:
1376  const std::vector<double>& knotVector_;
1377  unsigned int order_;
1378  bool makeOpen_;
1379 };
1380 
1381 } // end namespace BasisFactory::Imp
1382 
1389 auto bSpline(const std::vector<double>& knotVector,
1390  unsigned int order,
1391  bool makeOpen = true)
1392 {
1393  return Imp::BSplinePreBasisFactory(knotVector, order, makeOpen);
1394 }
1395 
1396 } // end namespace BasisFactory
1397 
1398 // *****************************************************************************
1399 // This is the actual global basis implementation based on the reusable parts.
1400 // *****************************************************************************
1401 
1408 template<typename GV>
1410 
1411 
1412 } // namespace Functions
1413 
1414 } // namespace Dune
1415 
1416 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
Dune::Functions::BSplineNodeIndexSet::unbind
void unbind()
Unbind the view.
Definition: bsplinebasis.hh:1304
Dune::Functions::BSplineNodeIndexSet::localSizes_
std::array< unsigned int, dim > localSizes_
Definition: bsplinebasis.hh:1347
Dune::Functions::BSplineNode::size_type
std::size_t size_type
Definition: bsplinebasis.hh:1225
Dune::Functions::BSplineLocalFiniteElement::localBasis_
BSplineLocalBasis< GV, R, MI > localBasis_
Definition: bsplinebasis.hh:471
Dune::Functions::BSplinePreBasis::size
size_type size(const SizePrefix prefix) const
Return number of possible values for next position in multi index.
Definition: bsplinebasis.hh:713
Dune::Functions::BSplinePreBasis::GridView
GV GridView
The grid view that the FE space is defined on.
Definition: bsplinebasis.hh:559
Dune::Functions::BSplinePreBasis::BSplinePreBasis
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
Dune::Functions::BSplinePreBasis::size
unsigned int size(size_t d) const
Number of shape functions in one direction.
Definition: bsplinebasis.hh:744
Dune::Functions::BSplineLocalFiniteElement::size
unsigned int size(int i) const
Number of degrees of freedom for one coordinate direction.
Definition: bsplinebasis.hh:458
Dune::Functions::BSplinePreBasis::makeNode
Node makeNode() const
Create tree node.
Definition: bsplinebasis.hh:696
Dune::Functions::BSplinePreBasis::R
double R
Definition: bsplinebasis.hh:572
Dune::Functions::BSplinePreBasis::evaluateFunctionFull
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:1055
Dune::Functions::BSplinePreBasis::MultiIndex
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: bsplinebasis.hh:567
Dune::Functions::BSplineNodeIndexSet
Definition: bsplinebasis.hh:484
Dune::Functions::BSplinePreBasis::evaluate
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:883
Dune::Functions::BSplineLocalBasis
LocalBasis class in the sense of dune-localfunctions, presenting the restriction of a B-spline patch ...
Definition: bsplinebasis.hh:46
Dune::Functions::BSplineLocalFiniteElement::type
GeometryType type() const
Return the reference element that the local finite element is defined on (here, a hypercube)
Definition: bsplinebasis.hh:450
Dune::Functions::BSplinePreBasis::maxNodeSize
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition: bsplinebasis.hh:726
Dune::Functions::BSplineLocalFiniteElement::BSplineLocalFiniteElement
BSplineLocalFiniteElement(const BSplinePreBasis< GV, MI > &preBasis)
Constructor with a given B-spline basis.
Definition: bsplinebasis.hh:378
Dune::Functions::BSplinePreBasis::size
unsigned int size() const
Total number of B-spline basis functions.
Definition: bsplinebasis.hh:735
Dune::Functions::BSplineLocalFiniteElement::localCoefficients_
BSplineLocalCoefficients< dim > localCoefficients_
Definition: bsplinebasis.hh:472
Dune::Functions::BSplineNode
Definition: bsplinebasis.hh:481
Dune::Functions::BSplineNodeIndexSet::size
size_type size() const
Size of subtree rooted in this node (element-local)
Definition: bsplinebasis.hh:1311
Dune::Functions::BSplineNode::preBasis_
const BSplinePreBasis< GV, MI > * preBasis_
Definition: bsplinebasis.hh:1260
Dune::Functions::BSplineNodeIndexSet::preBasis_
const PreBasis * preBasis_
Definition: bsplinebasis.hh:1343
Dune::Functions::BSplineLocalFiniteElement::bind
void bind(const std::array< unsigned, dim > &elementIdx)
Bind LocalFiniteElement to a specific knot span of the spline patch.
Definition: bsplinebasis.hh:389
Dune::Functions::BSplineLocalFiniteElement::localInterpolation_
BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R, MI > > localInterpolation_
Definition: bsplinebasis.hh:473
Dune::Functions::BSplineLocalBasis::Traits
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:56
Dune::Functions::BSplinePreBasis::gridView
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: bsplinebasis.hh:682
Dune::Functions::BSplineNodeIndexSet::bind
void bind(const Node &node)
Bind the view to a grid element.
Definition: bsplinebasis.hh:1293
Dune::Functions::BSplineLocalInterpolation::interpolate
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: bsplinebasis.hh:345
Dune::Functions::BasisFactory::bSpline
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:1389
Dune::Functions::BSplineNodeIndexSet::BSplineNodeIndexSet
BSplineNodeIndexSet(const PreBasis &preBasis)
Definition: bsplinebasis.hh:1284
Dune::Functions::BSplinePreBasis::evaluateAll
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:1098
Dune::Functions::BSplineLocalBasis::evaluateJacobian
void evaluateJacobian(const FieldVector< D, dim > &in, std::vector< FieldMatrix< D, 1, dim > > &out) const
Evaluate Jacobian of all shape functions.
Definition: bsplinebasis.hh:83
Dune::Functions::BSplinePreBasis::size_type
std::size_t size_type
Definition: bsplinebasis.hh:560
Dune::Functions::BSplineNode::element_
Element element_
Definition: bsplinebasis.hh:1263
Dune::Functions::BSplinePreBasis
Pre-basis for B-spline basis.
Definition: bsplinebasis.hh:34
Dune::Functions::DefaultGlobalBasis
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:46
Dune::Functions::BSplineNodeIndexSet::indices
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:1318
Dune::Functions::BSplineLocalFiniteElement::currentKnotSpan_
std::array< unsigned, dim > currentKnotSpan_
Definition: bsplinebasis.hh:476
Dune::Functions::BSplinePreBasis::update
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: bsplinebasis.hh:688
Dune::Functions::BSplinePreBasis::BSplinePreBasis
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
Dune::Functions::BSplinePreBasis::elements_
std::array< unsigned, dim > elements_
Number of grid elements in the different coordinate directions.
Definition: bsplinebasis.hh:1210
Dune::Functions::BSplinePreBasis::gridView_
GridView gridView_
Definition: bsplinebasis.hh:1212
Dune::Functions::BSplineLocalFiniteElement::size
unsigned size() const
Number of shape functions in this finite element.
Definition: bsplinebasis.hh:440
Dune::Functions::BSplineLocalBasis::BSplineLocalBasis
BSplineLocalBasis(const BSplinePreBasis< GV, MI > &preBasis, const BSplineLocalFiniteElement< GV, R, MI > &lFE)
Constructor with a given B-spline patch.
Definition: bsplinebasis.hh:62
Dune::Functions::BSplineNodeIndexSet::MultiIndex
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: bsplinebasis.hh:1278
Dune::Functions::BSplineNode::element
const Element & element() const
Return current element, throw if unbound.
Definition: bsplinebasis.hh:1235
Dune::Functions::BSplineLocalCoefficients
Attaches a shape function to an entity.
Definition: bsplinebasis.hh:178
Dune::Functions::BSplineNodeIndexSet::size_type
std::size_t size_type
Definition: bsplinebasis.hh:1275
Dune::Functions::BSplineLocalFiniteElement::Traits
LocalFiniteElementTraits< BSplineLocalBasis< GV, R, MI >, BSplineLocalCoefficients< dim >, BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R, MI > > > Traits
Export various types related to this LocalFiniteElement.
Definition: bsplinebasis.hh:374
Dune::Functions::BSplineLocalFiniteElement
LocalFiniteElement in the sense of dune-localfunctions, for the B-spline basis on tensor-product grid...
Definition: bsplinebasis.hh:31
Dune::Functions::BSplinePreBasis::order_
std::array< unsigned int, dim > order_
Order of the B-spline for each space dimension.
Definition: bsplinebasis.hh:1204
Dune::Functions::BSplineLocalFiniteElement::preBasis_
const BSplinePreBasis< GV, MI > & preBasis_
Definition: bsplinebasis.hh:469
nodes.hh
Dune::Functions::BSplineLocalFiniteElement::localCoefficients
const BSplineLocalCoefficients< dim > & localCoefficients() const
Hand out a LocalCoefficients object.
Definition: bsplinebasis.hh:428
Dune::Functions::BSplineNode::bind
void bind(const Element &e)
Bind to element.
Definition: bsplinebasis.hh:1250
defaultglobalbasis.hh
Dune
Definition: polynomial.hh:10
Dune::Functions::BSplinePreBasis::evaluateJacobian
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:782
Dune::Functions::BSplineLocalCoefficients::localKey
const LocalKey & localKey(std::size_t i) const
get i'th index
Definition: bsplinebasis.hh:322
Dune::Functions::BSplinePreBasis::evaluateFunction
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:998
Dune::Functions::BSplineNode::BSplineNode
BSplineNode(const BSplinePreBasis< GV, MI > *preBasis)
Definition: bsplinebasis.hh:1229
Dune::Functions::BSplineLocalBasis::evaluate
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:98
Dune::Functions::BSplineLocalFiniteElement::localBasis
const BSplineLocalBasis< GV, R, MI > & localBasis() const
Hand out a LocalBasis object.
Definition: bsplinebasis.hh:422
Dune::Functions::BSplineLocalBasis::order
unsigned int order() const
Polynomial order of the shape functions.
Definition: bsplinebasis.hh:141
Dune::Functions::BSplinePreBasis::getIJK
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:979
Dune::Functions::BSplinePreBasis::knotVectors_
std::array< std::vector< double >, dim > knotVectors_
The knot vectors, one for each space dimension.
Definition: bsplinebasis.hh:1207
Dune::Functions::BSplinePreBasis::initializeIndices
void initializeIndices()
Initialize the global indices.
Definition: bsplinebasis.hh:678
Dune::Functions::BSplineNode::finiteElement
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition: bsplinebasis.hh:1244
Dune::Functions::BSplinePreBasis::makeIndexSet
IndexSet makeIndexSet() const
Create tree node index set.
Definition: bsplinebasis.hh:707
Dune::Functions::BSplineLocalCoefficients::init
void init(const std::array< unsigned, dim > &sizes)
Definition: bsplinebasis.hh:257
flatmultiindex.hh
Dune::Functions::BSplineLocalCoefficients::size
std::size_t size() const
number of coefficients
Definition: bsplinebasis.hh:316
Dune::Functions::BSplineLocalFiniteElement::localInterpolation
const BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R, MI > > & localInterpolation() const
Hand out a LocalInterpolation object.
Definition: bsplinebasis.hh:434
Dune::Functions::BSplineLocalBasis::size
std::size_t size() const
Return the number of basis functions on the current knot span.
Definition: bsplinebasis.hh:148
Dune::Functions::LeafBasisNode
Definition: nodes.hh:179
Dune::Functions::BSplineLocalBasis::evaluateFunction
void evaluateFunction(const FieldVector< D, dim > &in, std::vector< FieldVector< R, 1 > > &out) const
Evaluate all shape functions.
Definition: bsplinebasis.hh:71
Dune::Functions::BSplineLocalInterpolation
Local interpolation in the sense of dune-localfunctions, for the B-spline basis on tensor-product gri...
Definition: bsplinebasis.hh:340
Dune::Functions::BSplinePreBasis::dimension
size_type dimension() const
Get the total dimension of the space spanned by this basis.
Definition: bsplinebasis.hh:720
Dune::Functions::BSplinePreBasis::evaluateFunction
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:751
Dune::Functions::BSplineNode::finiteElement_
FiniteElement finiteElement_
Definition: bsplinebasis.hh:1262
Dune::Functions::BSplineNode::Element
typename GV::template Codim< 0 >::Entity Element
Definition: bsplinebasis.hh:1226
Dune::Functions::BSplineNodeIndexSet::node_
const Node * node_
Definition: bsplinebasis.hh:1345
Dune::Functions::BSplinePreBasis::SizePrefix
Dune::ReservedVector< size_type, 1 > SizePrefix
Definition: bsplinebasis.hh:569