Dune Core Modules (2.8.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
15
16#include <dune/localfunctions/common/localbasis.hh>
18#include <dune/localfunctions/common/localkey.hh>
19#include <dune/localfunctions/common/localfiniteelementtraits.hh>
20#include <dune/geometry/type.hh>
21#include <dune/functions/functionspacebases/nodes.hh>
22#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
23#include <dune/functions/functionspacebases/flatmultiindex.hh>
24
25namespace Dune
26{
27namespace Functions {
28
29// A maze of dependencies between the different parts of this. We need a few forward declarations
30template<typename GV, typename R, typename MI>
31class BSplineLocalFiniteElement;
32
33template<typename GV, class MI>
34class BSplinePreBasis;
35
36
45template<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};
52public:
53
57
64 : preBasis_(preBasis),
65 lFE_(lFE)
66 {}
67
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
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
153private:
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
177template<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
256public:
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
327private:
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
339template<int dim, class LB>
341{
342public:
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
362template<class GV, class R, class MI>
364{
365 typedef typename GV::ctype D;
366 enum {dim = GV::dimension};
367 friend class BSplineLocalBasis<GV,R,MI>;
368public:
369
375
379 : preBasis_(preBasis),
380 localBasis_(preBasis,*this)
381 {}
382
386 : preBasis_(other.preBasis_),
387 localBasis_(preBasis_,*this)
388 {}
389
396 void bind(const std::array<unsigned,dim>& elementIdx)
397 {
398 /* \todo In the long run we need to precompute a table for this */
399 for (size_t i=0; i<elementIdx.size(); i++)
400 {
401 currentKnotSpan_[i] = 0;
402
403 // Skip over degenerate knot spans
404 while (preBasis_.knotVectors_[i][currentKnotSpan_[i]+1] < preBasis_.knotVectors_[i][currentKnotSpan_[i]]+1e-8)
405 currentKnotSpan_[i]++;
406
407 for (size_t j=0; j<elementIdx[i]; j++)
408 {
409 currentKnotSpan_[i]++;
410
411 // Skip over degenerate knot spans
412 while (preBasis_.knotVectors_[i][currentKnotSpan_[i]+1] < preBasis_.knotVectors_[i][currentKnotSpan_[i]]+1e-8)
413 currentKnotSpan_[i]++;
414 }
415
416 // Compute the geometric transformation from knotspan-local to global coordinates
417 localBasis_.offset_[i] = preBasis_.knotVectors_[i][currentKnotSpan_[i]];
418 localBasis_.scaling_[i][i] = preBasis_.knotVectors_[i][currentKnotSpan_[i]+1] - preBasis_.knotVectors_[i][currentKnotSpan_[i]];
419 }
420
421 // Set up the LocalCoefficients object
422 std::array<unsigned int, dim> sizes;
423 for (size_t i=0; i<dim; i++)
424 sizes[i] = size(i);
425 localCoefficients_.init(sizes);
426 }
427
430 {
431 return localBasis_;
432 }
433
436 {
437 return localCoefficients_;
438 }
439
442 {
443 return localInterpolation_;
444 }
445
447 unsigned size () const
448 {
449 std::size_t r = 1;
450 for (int i=0; i<dim; i++)
451 r *= size(i);
452 return r;
453 }
454
458 {
459 return GeometryTypes::cube(dim);
460 }
461
462//private:
463
465 unsigned int size(int i) const
466 {
467 const auto& order = preBasis_.order_;
468 unsigned int r = order[i]+1; // The 'normal' value
469 if (currentKnotSpan_[i]<order[i]) // Less near the left end of the knot vector
470 r -= (order[i] - currentKnotSpan_[i]);
471 if ( order[i] > (preBasis_.knotVectors_[i].size() - currentKnotSpan_[i] - 2) )
472 r -= order[i] - (preBasis_.knotVectors_[i].size() - currentKnotSpan_[i] - 2);
473 return r;
474 }
475
476 const BSplinePreBasis<GV,MI>& preBasis_;
477
478 BSplineLocalBasis<GV,R,MI> localBasis_;
479 BSplineLocalCoefficients<dim> localCoefficients_;
481
482 // The knot span we are bound to
483 std::array<unsigned,dim> currentKnotSpan_;
484};
485
486
487template<typename GV, typename MI>
488class BSplineNode;
489
500template<typename GV, class MI>
502{
503 static const int dim = GV::dimension;
504
506 class MultiDigitCounter
507 {
508 public:
509
513 MultiDigitCounter(const std::array<unsigned int,dim>& limits)
514 : limits_(limits)
515 {
516 std::fill(counter_.begin(), counter_.end(), 0);
517 }
518
520 MultiDigitCounter& operator++()
521 {
522 for (int i=0; i<dim; i++)
523 {
524 ++counter_[i];
525
526 // no overflow?
527 if (counter_[i] < limits_[i])
528 break;
529
530 counter_[i] = 0;
531 }
532 return *this;
533 }
534
536 const unsigned int& operator[](int i) const
537 {
538 return counter_[i];
539 }
540
542 unsigned int cycle() const
543 {
544 unsigned int r = 1;
545 for (int i=0; i<dim; i++)
546 r *= limits_[i];
547 return r;
548 }
549
550 private:
551
553 const std::array<unsigned int,dim> limits_;
554
556 std::array<unsigned int,dim> counter_;
557
558 };
559
560public:
561
563 using GridView = GV;
564 using size_type = std::size_t;
565
566 using Node = BSplineNode<GV, MI>;
567
569 using IndexSet = Impl::DefaultNodeIndexSet<BSplinePreBasis>;
570
572 using MultiIndex = MI;
573
575
576 // Type used for function values
577 using R = double;
578
598 const std::vector<double>& knotVector,
599 unsigned int order,
600 bool makeOpen = true)
601 : gridView_(gridView)
602 {
603 // \todo Detection of duplicate knots
604 std::fill(elements_.begin(), elements_.end(), knotVector.size()-1);
605
606 // Mediocre sanity check: we don't know the number of grid elements in each direction.
607 // but at least we know the total number of elements.
608 assert( std::accumulate(elements_.begin(), elements_.end(), 1, std::multiplies<unsigned>()) == gridView_.size(0) );
609
610 for (int i=0; i<dim; i++)
611 {
612 // Prepend the correct number of additional knots to open the knot vector
614 if (makeOpen)
615 for (unsigned int j=0; j<order; j++)
616 knotVectors_[i].push_back(knotVector[0]);
617
618 knotVectors_[i].insert(knotVectors_[i].end(), knotVector.begin(), knotVector.end());
619
620 if (makeOpen)
621 for (unsigned int j=0; j<order; j++)
622 knotVectors_[i].push_back(knotVector.back());
623 }
624
625 std::fill(order_.begin(), order_.end(), order);
626 }
627
650 const FieldVector<double,dim>& lowerLeft,
651 const FieldVector<double,dim>& upperRight,
652 const std::array<unsigned int,dim>& elements,
653 unsigned int order,
654 bool makeOpen = true)
655 : elements_(elements),
656 gridView_(gridView)
657 {
658 // Mediocre sanity check: we don't know the number of grid elements in each direction.
659 // but at least we know the total number of elements.
660 assert( std::accumulate(elements_.begin(), elements_.end(), 1, std::multiplies<unsigned>()) == gridView_.size(0) );
661
662 for (int i=0; i<dim; i++)
663 {
664 // Prepend the correct number of additional knots to open the knot vector
666 if (makeOpen)
667 for (unsigned int j=0; j<order; j++)
668 knotVectors_[i].push_back(lowerLeft[i]);
669
670 // Construct the actual knot vector
671 for (size_t j=0; j<elements[i]+1; j++)
672 knotVectors_[i].push_back(lowerLeft[i] + j*(upperRight[i]-lowerLeft[i]) / elements[i]);
673
674 if (makeOpen)
675 for (unsigned int j=0; j<order; j++)
676 knotVectors_[i].push_back(upperRight[i]);
677 }
678
679 std::fill(order_.begin(), order_.end(), order);
680 }
681
684 {}
685
687 const GridView& gridView() const
688 {
689 return gridView_;
690 }
691
693 void update(const GridView& gv)
694 {
695 gridView_ = gv;
696 }
697
701 Node makeNode() const
702 {
703 return Node{this};
704 }
705
713 [[deprecated("Warning: The IndexSet typedef and the makeIndexSet method are deprecated. "\
714 "As a replacement use the indices() method of the PreBasis directly.")]]
716 {
717 return IndexSet{*this};
718 }
719
721 size_type size(const SizePrefix prefix) const
722 {
723 assert(prefix.size() == 0 || prefix.size() == 1);
724 return (prefix.size() == 0) ? size() : 0;
725 }
726
728 size_type dimension() const
729 {
730 return size();
731 }
732
734 size_type maxNodeSize() const
735 {
736 size_type result = 1;
737 for (int i=0; i<dim; i++)
738 result *= order_[i]+1;
739 return result;
740 }
741
743 template<typename It>
744 It indices(const Node& node, It it) const
745 {
746 // Local degrees of freedom are arranged in a lattice.
747 // We need the lattice dimensions to be able to compute lattice coordinates from a local index
748 std::array<unsigned int, dim> localSizes;
749 for (int i=0; i<dim; i++)
750 localSizes[i] = node.finiteElement().size(i);
751 for (size_type i = 0, end = node.size() ; i < end ; ++i, ++it)
752 {
753 std::array<unsigned int,dim> localIJK = getIJK(i, localSizes);
754
755 const auto currentKnotSpan = node.finiteElement().currentKnotSpan_;
756 const auto order = order_;
757
758 std::array<unsigned int,dim> globalIJK;
759 for (int i=0; i<dim; i++)
760 globalIJK[i] = std::max((int)currentKnotSpan[i] - (int)order[i], 0) + localIJK[i]; // needs to be a signed type!
761
762 // Make one global flat index from the globalIJK tuple
763 size_type globalIdx = globalIJK[dim-1];
764
765 for (int i=dim-2; i>=0; i--)
766 globalIdx = globalIdx * size(i) + globalIJK[i];
767
768 *it = {{globalIdx}};
769 }
770 return it;
771 }
772
774 unsigned int size () const
775 {
776 unsigned int result = 1;
777 for (size_t i=0; i<dim; i++)
778 result *= size(i);
779 return result;
780 }
781
783 unsigned int size (size_t d) const
784 {
785 return knotVectors_[d].size() - order_[d] - 1;
786 }
787
791 std::vector<FieldVector<R,1> >& out,
792 const std::array<unsigned,dim>& currentKnotSpan) const
793 {
794 // Evaluate
795 std::array<std::vector<R>, dim> oneDValues;
796
797 for (size_t i=0; i<dim; i++)
798 evaluateFunction(in[i], oneDValues[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
799
800 std::array<unsigned int, dim> limits;
801 for (int i=0; i<dim; i++)
802 limits[i] = oneDValues[i].size();
803
804 MultiDigitCounter ijkCounter(limits);
805
806 out.resize(ijkCounter.cycle());
807
808 for (size_t i=0; i<out.size(); i++, ++ijkCounter)
809 {
810 out[i] = R(1.0);
811 for (size_t j=0; j<dim; j++)
812 out[i] *= oneDValues[j][ijkCounter[j]];
813 }
814 }
815
822 std::vector<FieldMatrix<R,1,dim> >& out,
823 const std::array<unsigned,dim>& currentKnotSpan) const
824 {
825 // How many shape functions to we have in each coordinate direction?
826 std::array<unsigned int, dim> limits;
827 for (int i=0; i<dim; i++)
828 {
829 limits[i] = order_[i]+1; // The 'standard' value away from the boundaries of the knot vector
830 if (currentKnotSpan[i]<order_[i])
831 limits[i] -= (order_[i] - currentKnotSpan[i]);
832 if ( order_[i] > (knotVectors_[i].size() - currentKnotSpan[i] - 2) )
833 limits[i] -= order_[i] - (knotVectors_[i].size() - currentKnotSpan[i] - 2);
834 }
835
836 // The lowest knot spans that we need values from
837 std::array<unsigned int, dim> offset;
838 for (int i=0; i<dim; i++)
839 offset[i] = std::max((int)(currentKnotSpan[i] - order_[i]),0);
840
841 // Evaluate 1d function values (needed for the product rule)
842 std::array<std::vector<R>, dim> oneDValues;
843
844 // Evaluate 1d function values of one order lower (needed for the derivative formula)
845 std::array<std::vector<R>, dim> lowOrderOneDValues;
846
847 std::array<DynamicMatrix<R>, dim> values;
848
849 for (size_t i=0; i<dim; i++)
850 {
851 evaluateFunctionFull(in[i], values[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
852 oneDValues[i].resize(knotVectors_[i].size()-order_[i]-1);
853 for (size_t j=0; j<oneDValues[i].size(); j++)
854 oneDValues[i][j] = values[i][order_[i]][j];
855
856 if (order_[i]!=0)
857 {
858 lowOrderOneDValues[i].resize(knotVectors_[i].size()-(order_[i]-1)-1);
859 for (size_t j=0; j<lowOrderOneDValues[i].size(); j++)
860 lowOrderOneDValues[i][j] = values[i][order_[i]-1][j];
861 }
862 }
863
864
865 // Evaluate 1d function derivatives
866 std::array<std::vector<R>, dim> oneDDerivatives;
867 for (size_t i=0; i<dim; i++)
868 {
869 oneDDerivatives[i].resize(limits[i]);
870
871 if (order_[i]==0) // order-zero functions are piecewise constant, hence all derivatives are zero
872 std::fill(oneDDerivatives[i].begin(), oneDDerivatives[i].end(), R(0.0));
873 else
874 {
875 for (size_t j=offset[i]; j<offset[i]+limits[i]; j++)
876 {
877 R derivativeAddend1 = lowOrderOneDValues[i][j] / (knotVectors_[i][j+order_[i]]-knotVectors_[i][j]);
878 R derivativeAddend2 = lowOrderOneDValues[i][j+1] / (knotVectors_[i][j+order_[i]+1]-knotVectors_[i][j+1]);
879 // The two previous terms may evaluate as 0/0. This is to be interpreted as 0.
880 if (std::isnan(derivativeAddend1))
881 derivativeAddend1 = 0;
882 if (std::isnan(derivativeAddend2))
883 derivativeAddend2 = 0;
884 oneDDerivatives[i][j-offset[i]] = order_[i] * ( derivativeAddend1 - derivativeAddend2 );
885 }
886 }
887 }
888
889 // Working towards computing only the parts that we really need:
890 // Let's copy them out into a separate array
891 std::array<std::vector<R>, dim> oneDValuesShort;
892
893 for (int i=0; i<dim; i++)
894 {
895 oneDValuesShort[i].resize(limits[i]);
896
897 for (size_t j=0; j<limits[i]; j++)
898 oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
899 }
900
901
902
903 // Set up a multi-index to go from consecutive indices to integer coordinates
904 MultiDigitCounter ijkCounter(limits);
905
906 out.resize(ijkCounter.cycle());
907
908 // Complete Jacobian is given by the product rule
909 for (size_t i=0; i<out.size(); i++, ++ijkCounter)
910 for (int j=0; j<dim; j++)
911 {
912 out[i][0][j] = 1.0;
913 for (int k=0; k<dim; k++)
914 out[i][0][j] *= (j==k) ? oneDDerivatives[k][ijkCounter[k]]
915 : oneDValuesShort[k][ijkCounter[k]];
916 }
917
918 }
919
921 template <size_type k>
922 void evaluate(const typename std::array<int,k>& directions,
924 std::vector<FieldVector<R,1> >& out,
925 const std::array<unsigned,dim>& currentKnotSpan) const
926 {
927 if (k != 1 && k != 2)
928 DUNE_THROW(RangeError, "Differentiation order greater than 2 is not supported!");
929
930 // Evaluate 1d function values (needed for the product rule)
931 std::array<std::vector<R>, dim> oneDValues;
932 std::array<std::vector<R>, dim> oneDDerivatives;
933 std::array<std::vector<R>, dim> oneDSecondDerivatives;
934
935 // Evaluate 1d function derivatives
936 if (k==1)
937 for (size_t i=0; i<dim; i++)
938 evaluateAll(in[i], oneDValues[i], true, oneDDerivatives[i], false, oneDSecondDerivatives[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
939 else
940 for (size_t i=0; i<dim; i++)
941 evaluateAll(in[i], oneDValues[i], true, oneDDerivatives[i], true, oneDSecondDerivatives[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
942
943 // The lowest knot spans that we need values from
944 std::array<unsigned int, dim> offset;
945 for (int i=0; i<dim; i++)
946 offset[i] = std::max((int)(currentKnotSpan[i] - order_[i]),0);
947
948 // Set up a multi-index to go from consecutive indices to integer coordinates
949 std::array<unsigned int, dim> limits;
950 for (int i=0; i<dim; i++)
951 {
952 // In a proper implementation, the following line would do
953 //limits[i] = oneDValues[i].size();
954 limits[i] = order_[i]+1; // The 'standard' value away from the boundaries of the knot vector
955 if (currentKnotSpan[i]<order_[i])
956 limits[i] -= (order_[i] - currentKnotSpan[i]);
957 if ( order_[i] > (knotVectors_[i].size() - currentKnotSpan[i] - 2) )
958 limits[i] -= order_[i] - (knotVectors_[i].size() - currentKnotSpan[i] - 2);
959 }
960
961 // Working towards computing only the parts that we really need:
962 // Let's copy them out into a separate array
963 std::array<std::vector<R>, dim> oneDValuesShort;
964
965 for (int i=0; i<dim; i++)
966 {
967 oneDValuesShort[i].resize(limits[i]);
968
969 for (size_t j=0; j<limits[i]; j++)
970 oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
971 }
972
973
974 MultiDigitCounter ijkCounter(limits);
975
976 out.resize(ijkCounter.cycle());
977
978 if (k == 1)
979 {
980 // Complete Jacobian is given by the product rule
981 for (size_t i=0; i<out.size(); i++, ++ijkCounter)
982 {
983 out[i][0] = 1.0;
984 for (int l=0; l<dim; l++)
985 out[i][0] *= (directions[0]==l) ? oneDDerivatives[l][ijkCounter[l]]
986 : oneDValuesShort[l][ijkCounter[l]];
987 }
988 }
989
990 if (k == 2)
991 {
992 // Complete derivation by deriving the tensor product
993 for (size_t i=0; i<out.size(); i++, ++ijkCounter)
994 {
995 out[i][0] = 1.0;
996 for (int j=0; j<dim; j++)
997 {
998 if (directions[0] != directions[1]) //derivation in two different variables
999 if (directions[0] == j || directions[1] == j) //the spline has to be derived (once) in this direction
1000 out[i][0] *= oneDDerivatives[j][ijkCounter[j]];
1001 else //no derivation in this direction
1002 out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
1003 else //spline is derived two times in the same direction
1004 if (directions[0] == j) //the spline is derived two times in this direction
1005 out[i][0] *= oneDSecondDerivatives[j][ijkCounter[j]];
1006 else //no derivation in this direction
1007 out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
1008 }
1009 }
1010 }
1011 }
1012
1013
1018 static std::array<unsigned int,dim> getIJK(typename GridView::IndexSet::IndexType idx, std::array<unsigned int,dim> elements)
1019 {
1020 std::array<unsigned,dim> result;
1021 for (int i=0; i<dim; i++)
1022 {
1023 result[i] = idx%elements[i];
1024 idx /= elements[i];
1025 }
1026 return result;
1027 }
1028
1037 static void evaluateFunction (const typename GV::ctype& in, std::vector<R>& out,
1038 const std::vector<R>& knotVector,
1039 unsigned int order,
1040 unsigned int currentKnotSpan)
1041 {
1042 std::size_t outSize = order+1; // The 'standard' value away from the boundaries of the knot vector
1043 if (currentKnotSpan<order) // Less near the left end of the knot vector
1044 outSize -= (order - currentKnotSpan);
1045 if ( order > (knotVector.size() - currentKnotSpan - 2) )
1046 outSize -= order - (knotVector.size() - currentKnotSpan - 2);
1047 out.resize(outSize);
1048
1049 // It's not really a matrix that is needed here, a plain 2d array would do
1050 DynamicMatrix<R> N(order+1, knotVector.size()-1);
1051
1052 // The text books on splines use the following geometric condition here to fill the array N
1053 // (see for example Cottrell, Hughes, Bazilevs, Formula (2.1). However, this condition
1054 // only works if splines are never evaluated exactly on the knots.
1055 //
1056 // for (size_t i=0; i<knotVector.size()-1; i++)
1057 // N[0][i] = (knotVector[i] <= in) and (in < knotVector[i+1]);
1058 for (size_t i=0; i<knotVector.size()-1; i++)
1059 N[0][i] = (i == currentKnotSpan);
1060
1061 for (size_t r=1; r<=order; r++)
1062 for (size_t i=0; i<knotVector.size()-r-1; i++)
1063 {
1064 R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
1065 ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
1066 : 0;
1067 R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
1068 ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
1069 : 0;
1070 N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
1071 }
1072
1077 for (size_t i=0; i<out.size(); i++) {
1078 out[i] = N[order][std::max((int)(currentKnotSpan - order),0) + i];
1079 }
1080 }
1081
1094 static void evaluateFunctionFull(const typename GV::ctype& in,
1095 DynamicMatrix<R>& out,
1096 const std::vector<R>& knotVector,
1097 unsigned int order,
1098 unsigned int currentKnotSpan)
1099 {
1100 // It's not really a matrix that is needed here, a plain 2d array would do
1101 DynamicMatrix<R>& N = out;
1102
1103 N.resize(order+1, knotVector.size()-1);
1104
1105 // The text books on splines use the following geometric condition here to fill the array N
1106 // (see for example Cottrell, Hughes, Bazilevs, Formula (2.1). However, this condition
1107 // only works if splines are never evaluated exactly on the knots.
1108 //
1109 // for (size_t i=0; i<knotVector.size()-1; i++)
1110 // N[0][i] = (knotVector[i] <= in) and (in < knotVector[i+1]);
1111 for (size_t i=0; i<knotVector.size()-1; i++)
1112 N[0][i] = (i == currentKnotSpan);
1113
1114 for (size_t r=1; r<=order; r++)
1115 for (size_t i=0; i<knotVector.size()-r-1; i++)
1116 {
1117 R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
1118 ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
1119 : 0;
1120 R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
1121 ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
1122 : 0;
1123 N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
1124 }
1125 }
1126
1127
1137 static void evaluateAll(const typename GV::ctype& in,
1138 std::vector<R>& out,
1139 bool evaluateJacobian, std::vector<R>& outJac,
1140 bool evaluateHessian, std::vector<R>& outHess,
1141 const std::vector<R>& knotVector,
1142 unsigned int order,
1143 unsigned int currentKnotSpan)
1144 {
1145 // How many shape functions to we have in each coordinate direction?
1146 unsigned int limit;
1147 limit = order+1; // The 'standard' value away from the boundaries of the knot vector
1148 if (currentKnotSpan<order)
1149 limit -= (order - currentKnotSpan);
1150 if ( order > (knotVector.size() - currentKnotSpan - 2) )
1151 limit -= order - (knotVector.size() - currentKnotSpan - 2);
1152
1153 // The lowest knot spans that we need values from
1154 unsigned int offset;
1155 offset = std::max((int)(currentKnotSpan - order),0);
1156
1157 // Evaluate 1d function values (needed for the product rule)
1158 DynamicMatrix<R> values;
1159
1160 evaluateFunctionFull(in, values, knotVector, order, currentKnotSpan);
1161
1162 out.resize(knotVector.size()-order-1);
1163 for (size_t j=0; j<out.size(); j++)
1164 out[j] = values[order][j];
1165
1166 // Evaluate 1d function values of one order lower (needed for the derivative formula)
1167 std::vector<R> lowOrderOneDValues;
1168
1169 if (order!=0)
1170 {
1171 lowOrderOneDValues.resize(knotVector.size()-(order-1)-1);
1172 for (size_t j=0; j<lowOrderOneDValues.size(); j++)
1173 lowOrderOneDValues[j] = values[order-1][j];
1174 }
1175
1176 // Evaluate 1d function values of two order lower (needed for the (second) derivative formula)
1177 std::vector<R> lowOrderTwoDValues;
1178
1179 if (order>1 && evaluateHessian)
1180 {
1181 lowOrderTwoDValues.resize(knotVector.size()-(order-2)-1);
1182 for (size_t j=0; j<lowOrderTwoDValues.size(); j++)
1183 lowOrderTwoDValues[j] = values[order-2][j];
1184 }
1185
1186 // Evaluate 1d function derivatives
1187 if (evaluateJacobian)
1188 {
1189 outJac.resize(limit);
1190
1191 if (order==0) // order-zero functions are piecewise constant, hence all derivatives are zero
1192 std::fill(outJac.begin(), outJac.end(), R(0.0));
1193 else
1194 {
1195 for (size_t j=offset; j<offset+limit; j++)
1196 {
1197 R derivativeAddend1 = lowOrderOneDValues[j] / (knotVector[j+order]-knotVector[j]);
1198 R derivativeAddend2 = lowOrderOneDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]);
1199 // The two previous terms may evaluate as 0/0. This is to be interpreted as 0.
1200 if (std::isnan(derivativeAddend1))
1201 derivativeAddend1 = 0;
1202 if (std::isnan(derivativeAddend2))
1203 derivativeAddend2 = 0;
1204 outJac[j-offset] = order * ( derivativeAddend1 - derivativeAddend2 );
1205 }
1206 }
1207 }
1208
1209 // Evaluate 1d function second derivatives
1210 if (evaluateHessian)
1211 {
1212 outHess.resize(limit);
1213
1214 if (order<2) // order-zero functions are piecewise constant, hence all derivatives are zero
1215 std::fill(outHess.begin(), outHess.end(), R(0.0));
1216 else
1217 {
1218 for (size_t j=offset; j<offset+limit; j++)
1219 {
1220 assert(j+2 < lowOrderTwoDValues.size());
1221 R derivativeAddend1 = lowOrderTwoDValues[j] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order-1]-knotVector[j]);
1222 R derivativeAddend2 = lowOrderTwoDValues[j+1] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order]-knotVector[j+1]);
1223 R derivativeAddend3 = lowOrderTwoDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+order]-knotVector[j+1]);
1224 R derivativeAddend4 = lowOrderTwoDValues[j+2] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+1+order]-knotVector[j+2]);
1225 // The two previous terms may evaluate as 0/0. This is to be interpreted as 0.
1226
1227 if (std::isnan(derivativeAddend1))
1228 derivativeAddend1 = 0;
1229 if (std::isnan(derivativeAddend2))
1230 derivativeAddend2 = 0;
1231 if (std::isnan(derivativeAddend3))
1232 derivativeAddend3 = 0;
1233 if (std::isnan(derivativeAddend4))
1234 derivativeAddend4 = 0;
1235 outHess[j-offset] = order * (order-1) * ( derivativeAddend1 - derivativeAddend2 -derivativeAddend3 + derivativeAddend4 );
1236 }
1237 }
1238 }
1239 }
1240
1241
1243 std::array<unsigned int, dim> order_;
1244
1246 std::array<std::vector<double>, dim> knotVectors_;
1247
1249 std::array<unsigned,dim> elements_;
1250
1251 GridView gridView_;
1252};
1253
1254
1255
1256template<typename GV, typename MI>
1257class BSplineNode :
1258 public LeafBasisNode
1259{
1260 static const int dim = GV::dimension;
1261
1262public:
1263
1264 using size_type = std::size_t;
1265 using Element = typename GV::template Codim<0>::Entity;
1266 using FiniteElement = BSplineLocalFiniteElement<GV,double,MI>;
1267
1268 BSplineNode(const BSplinePreBasis<GV, MI>* preBasis) :
1269 preBasis_(preBasis),
1270 finiteElement_(*preBasis)
1271 {}
1272
1274 const Element& element() const
1275 {
1276 return element_;
1277 }
1278
1283 const FiniteElement& finiteElement() const
1284 {
1285 return finiteElement_;
1286 }
1287
1289 void bind(const Element& e)
1290 {
1291 element_ = e;
1292 auto elementIndex = preBasis_->gridView().indexSet().index(e);
1293 finiteElement_.bind(preBasis_->getIJK(elementIndex,preBasis_->elements_));
1294 this->setSize(finiteElement_.size());
1295 }
1296
1297protected:
1298
1299 const BSplinePreBasis<GV,MI>* preBasis_;
1300
1301 FiniteElement finiteElement_;
1302 Element element_;
1303};
1304
1305
1306
1307namespace BasisFactory {
1308
1309namespace Imp {
1310
1311class BSplinePreBasisFactory
1312{
1313public:
1314 static const std::size_t requiredMultiIndexSize=1;
1315
1316 BSplinePreBasisFactory(const std::vector<double>& knotVector,
1317 unsigned int order,
1318 bool makeOpen = true)
1319 : knotVector_(knotVector),
1320 order_(order),
1321 makeOpen_(makeOpen)
1322 {}
1323
1324 template<class MultiIndex, class GridView>
1325 auto makePreBasis(const GridView& gridView) const
1326 {
1327 return BSplinePreBasis<GridView, MultiIndex>(gridView, knotVector_, order_, makeOpen_);
1328 }
1329
1330private:
1331 const std::vector<double>& knotVector_;
1332 unsigned int order_;
1333 bool makeOpen_;
1334};
1335
1336} // end namespace BasisFactory::Imp
1337
1344auto bSpline(const std::vector<double>& knotVector,
1345 unsigned int order,
1346 bool makeOpen = true)
1347{
1348 return Imp::BSplinePreBasisFactory(knotVector, order, makeOpen);
1349}
1350
1351} // end namespace BasisFactory
1352
1353// *****************************************************************************
1354// This is the actual global basis implementation based on the reusable parts.
1355// *****************************************************************************
1356
1363template<typename GV>
1365
1366
1367} // namespace Functions
1368
1369} // namespace Dune
1370
1371#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
Construct a matrix with a dynamic size.
Definition: dynmatrix.hh:59
void resize(size_type r, size_type c, value_type v=value_type())
resize matrix to r × c
Definition: dynmatrix.hh:104
A dense n x m matrix.
Definition: fmatrix.hh:69
vector space out of a tensor product of fields.
Definition: fvector.hh:95
LocalBasis class in the sense of dune-localfunctions, presenting the restriction of a B-spline patch ...
Definition: bsplinebasis.hh:47
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
BSplineLocalBasis(const BSplinePreBasis< GV, MI > &preBasis, const BSplineLocalFiniteElement< GV, R, MI > &lFE)
Constructor with a given B-spline patch.
Definition: bsplinebasis.hh:62
void evaluateFunction(const FieldVector< D, dim > &in, std::vector< FieldVector< R, 1 > > &out) const
Evaluate all shape functions.
Definition: bsplinebasis.hh:71
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
unsigned int order() const
Polynomial order of the shape functions.
Definition: bsplinebasis.hh:141
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
std::size_t size() const
Return the number of basis functions on the current knot span.
Definition: bsplinebasis.hh:148
Attaches a shape function to an entity.
Definition: bsplinebasis.hh:179
const LocalKey & localKey(std::size_t i) const
get i'th index
Definition: bsplinebasis.hh:322
std::size_t size() const
number of coefficients
Definition: bsplinebasis.hh:316
LocalFiniteElement in the sense of dune-localfunctions, for the B-spline basis on tensor-product grid...
Definition: bsplinebasis.hh:364
const BSplineLocalBasis< GV, R, MI > & localBasis() const
Hand out a LocalBasis object.
Definition: bsplinebasis.hh:429
unsigned size() const
Number of shape functions in this finite element.
Definition: bsplinebasis.hh:447
GeometryType type() const
Return the reference element that the local finite element is defined on (here, a hypercube)
Definition: bsplinebasis.hh:457
void bind(const std::array< unsigned, dim > &elementIdx)
Bind LocalFiniteElement to a specific knot span of the spline patch.
Definition: bsplinebasis.hh:396
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
const BSplineLocalCoefficients< dim > & localCoefficients() const
Hand out a LocalCoefficients object.
Definition: bsplinebasis.hh:435
const BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R, MI > > & localInterpolation() const
Hand out a LocalInterpolation object.
Definition: bsplinebasis.hh:441
BSplineLocalFiniteElement(const BSplineLocalFiniteElement &other)
Copy constructor.
Definition: bsplinebasis.hh:385
BSplineLocalFiniteElement(const BSplinePreBasis< GV, MI > &preBasis)
Constructor with a given B-spline basis.
Definition: bsplinebasis.hh:378
unsigned int size(int i) const
Number of degrees of freedom for one coordinate direction.
Definition: bsplinebasis.hh:465
Local interpolation in the sense of dune-localfunctions, for the B-spline basis on tensor-product gri...
Definition: bsplinebasis.hh:341
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: bsplinebasis.hh:345
Pre-basis for B-spline basis.
Definition: bsplinebasis.hh:502
GV GridView
The grid view that the FE space is defined on.
Definition: bsplinebasis.hh:563
unsigned int size(size_t d) const
Number of shape functions in one direction.
Definition: bsplinebasis.hh:783
size_type dimension() const
Get the total dimension of the space spanned by this basis.
Definition: bsplinebasis.hh:728
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: bsplinebasis.hh:693
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:821
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:1037
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:597
unsigned int size() const
Total number of B-spline basis functions.
Definition: bsplinebasis.hh:774
size_type size(const SizePrefix prefix) const
Return number of possible values for next position in multi index.
Definition: bsplinebasis.hh:721
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: bsplinebasis.hh:572
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition: bsplinebasis.hh:734
void initializeIndices()
Initialize the global indices.
Definition: bsplinebasis.hh:683
std::array< unsigned, dim > elements_
Number of grid elements in the different coordinate directions.
Definition: bsplinebasis.hh:1249
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: bsplinebasis.hh:687
It indices(const Node &node, It it) const
Maps from subtree index set [0..size-1] to a globally unique multi index in global basis.
Definition: bsplinebasis.hh:744
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:1137
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:790
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:649
Impl::DefaultNodeIndexSet< BSplinePreBasis > IndexSet
Type of created tree node index set.
Definition: bsplinebasis.hh:569
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:1094
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:1018
IndexSet makeIndexSet() const
Create tree node index set.
Definition: bsplinebasis.hh:715
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:922
Node makeNode() const
Create tree node.
Definition: bsplinebasis.hh:701
std::array< std::vector< double >, dim > knotVectors_
The knot vectors, one for each space dimension.
Definition: bsplinebasis.hh:1246
std::array< unsigned int, dim > order_
Order of the B-spline for each space dimension.
Definition: bsplinebasis.hh:1243
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:47
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:123
Describe position of one degree of freedom.
Definition: localkey.hh:21
Default exception for dummy implementations.
Definition: exceptions.hh:261
Default exception class for range errors.
Definition: exceptions.hh:252
A Vector class with statically reserved memory.
Definition: reservedvector.hh:43
size_type size() const
Returns number of elements in the vector.
Definition: reservedvector.hh:184
This file implements a quadratic diagonal matrix of fixed size.
This file implements a dense matrix with dynamic numbers of rows and columns.
void umv(const X &x, Y &y) const
y += A x
Definition: diagonalmatrix.hh:295
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
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:1344
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:470
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:289
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
constexpr HybridTreePath< T..., std::size_t > push_back(const HybridTreePath< T... > &tp, std::size_t i)
Appends a run time index to a HybridTreePath.
Definition: treepath.hh:278
Dune namespace.
Definition: alignedallocator.hh:11
Static tag representing a codimension.
Definition: dimension.hh:22
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:32
D DomainType
domain type
Definition: localbasis.hh:43
traits helper struct
Definition: localfiniteelementtraits.hh:11
A unique label for each type of element that can occur in a grid.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)