DUNE PDELab (git)

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
4// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
5// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
6
7#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
9
14#include <array>
15#include <numeric>
16
19
20#include <dune/localfunctions/common/localbasis.hh>
22#include <dune/localfunctions/common/localkey.hh>
23#include <dune/localfunctions/common/localfiniteelementtraits.hh>
24#include <dune/geometry/type.hh>
25#include <dune/functions/functionspacebases/nodes.hh>
26#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
27#include <dune/functions/functionspacebases/leafprebasismixin.hh>
28
29namespace Dune
30{
31namespace Functions {
32
33// A maze of dependencies between the different parts of this. We need a few forward declarations
34template<typename GV, typename R>
35class BSplineLocalFiniteElement;
36
37template<typename GV>
38class BSplinePreBasis;
39
40
49template<class GV, class R>
51{
52 friend class BSplineLocalFiniteElement<GV,R>;
53
54 typedef typename GV::ctype D;
55 enum {dim = GV::dimension};
56public:
57
61
68 : preBasis_(preBasis),
69 lFE_(lFE)
70 {}
71
76 std::vector<FieldVector<R,1> >& out) const
77 {
78 FieldVector<D,dim> globalIn = offset_;
79 scaling_.umv(in,globalIn);
80
81 preBasis_.evaluateFunction(globalIn, out, lFE_.currentKnotSpan_);
82 }
83
88 std::vector<FieldMatrix<D,1,dim> >& out) const
89 {
90 FieldVector<D,dim> globalIn = offset_;
91 scaling_.umv(in,globalIn);
92
93 preBasis_.evaluateJacobian(globalIn, out, lFE_.currentKnotSpan_);
94
95 for (size_t i=0; i<out.size(); i++)
96 for (int j=0; j<dim; j++)
97 out[i][0][j] *= scaling_[j][j];
98 }
99
101 template<size_t k>
102 inline void evaluate (const typename std::array<int,k>& directions,
103 const typename Traits::DomainType& in,
104 std::vector<typename Traits::RangeType>& out) const
105 {
106 switch(k)
107 {
108 case 0:
109 evaluateFunction(in, out);
110 break;
111 case 1:
112 {
113 FieldVector<D,dim> globalIn = offset_;
114 scaling_.umv(in,globalIn);
115
116 preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
117
118 for (size_t i=0; i<out.size(); i++)
119 out[i][0] *= scaling_[directions[0]][directions[0]];
120 break;
121 }
122 case 2:
123 {
124 FieldVector<D,dim> globalIn = offset_;
125 scaling_.umv(in,globalIn);
126
127 preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
128
129 for (size_t i=0; i<out.size(); i++)
130 out[i][0] *= scaling_[directions[0]][directions[0]]*scaling_[directions[1]][directions[1]];
131 break;
132 }
133 default:
134 DUNE_THROW(NotImplemented, "B-Spline derivatives of order " << k << " not implemented yet!");
135 }
136 }
137
145 unsigned int order () const
146 {
147 return *std::max_element(preBasis_.order_.begin(), preBasis_.order_.end());
148 }
149
152 std::size_t size() const
153 {
154 return lFE_.size();
155 }
156
157private:
158 const BSplinePreBasis<GV>& preBasis_;
159
161
162 // Coordinates in a single knot span differ from coordinates on the B-spline patch
163 // by an affine transformation. This transformation is stored in offset_ and scaling_.
164 FieldVector<D,dim> offset_;
165 DiagonalMatrix<D,dim> scaling_;
166};
167
181template<int dim>
183{
184 // Return i as a d-digit number in the (k+1)-nary system
185 std::array<unsigned int,dim> multiindex (unsigned int i) const
186 {
187 std::array<unsigned int,dim> alpha;
188 for (int j=0; j<dim; j++)
189 {
190 alpha[j] = i % sizes_[j];
191 i = i/sizes_[j];
192 }
193 return alpha;
194 }
195
197 void setup1d(std::vector<unsigned int>& subEntity)
198 {
199 if (sizes_[0]==1)
200 {
201 subEntity[0] = 0;
202 return;
203 }
204
205 /* edge and vertex numbering
206 0----0----1
207 */
208 unsigned lastIndex=0;
209 subEntity[lastIndex++] = 0; // corner 0
210 for (unsigned i = 0; i < sizes_[0] - 2; ++i)
211 subEntity[lastIndex++] = 0; // inner dofs of element (0)
212
213 subEntity[lastIndex++] = 1; // corner 1
214
215 assert(size()==lastIndex);
216 }
217
218 void setup2d(std::vector<unsigned int>& subEntity)
219 {
220 unsigned lastIndex=0;
221
222 // LocalKey: entity number , entity codim, dof indices within each entity
223 /* edge and vertex numbering
224 2----3----3
225 | |
226 | |
227 0 1
228 | |
229 | |
230 0----2----1
231 */
232
233 // lower edge (2)
234 subEntity[lastIndex++] = 0; // corner 0
235 for (unsigned i = 0; i < sizes_[0]-2; ++i)
236 subEntity[lastIndex++] = 2; // inner dofs of lower edge (2)
237
238 subEntity[lastIndex++] = 1; // corner 1
239
240 // iterate from bottom to top over inner edge dofs
241 for (unsigned e = 0; e < sizes_[1]-2; ++e)
242 {
243 subEntity[lastIndex++] = 0; // left edge (0)
244 for (unsigned i = 0; i < sizes_[0]-2; ++i)
245 subEntity[lastIndex++] = 0; // face dofs
246 subEntity[lastIndex++] = 1; // right edge (1)
247 }
248
249 // upper edge (3)
250 subEntity[lastIndex++] = 2; // corner 2
251 for (unsigned i = 0; i < sizes_[0]-2; ++i)
252 subEntity[lastIndex++] = 3; // inner dofs of upper edge (3)
253
254 subEntity[lastIndex++] = 3; // corner 3
255
256 assert(size()==lastIndex);
257 }
258
259
260public:
261 void init(const std::array<unsigned,dim>& sizes)
262 {
263 sizes_ = sizes;
264
265 li_.resize(size());
266
267 // Set up array of codimension-per-dof-number
268 std::vector<unsigned int> codim(li_.size());
269
270 for (std::size_t i=0; i<codim.size(); i++)
271 {
272 codim[i] = 0;
273 // Codimension gets increased by 1 for each coordinate direction
274 // where dof is on boundary
275 std::array<unsigned int,dim> mIdx = multiindex(i);
276 for (int j=0; j<dim; j++)
277 if (mIdx[j]==0 or mIdx[j]==sizes[j]-1)
278 codim[i]++;
279 }
280
281 // Set up index vector (the index of the dof in the set of dofs of a given subentity)
282 // Algorithm: the 'index' has the same ordering as the dof number 'i'.
283 // To make it consecutive we interpret 'i' in the (k+1)-adic system, omit all digits
284 // that correspond to axes where the dof is on the element boundary, and transform the
285 // rest to the (k-1)-adic system.
286 std::vector<unsigned int> index(size());
287
288 for (std::size_t i=0; i<index.size(); i++)
289 {
290 index[i] = 0;
291
292 std::array<unsigned int,dim> mIdx = multiindex(i);
293
294 for (int j=dim-1; j>=0; j--)
295 if (mIdx[j]>0 and mIdx[j]<sizes[j]-1)
296 index[i] = (sizes[j]-1)*index[i] + (mIdx[j]-1);
297 }
298
299 // Set up entity and dof numbers for each (supported) dimension separately
300 std::vector<unsigned int> subEntity(li_.size());
301
302 if (subEntity.size() > 0)
303 {
304 if (dim==1) {
305
306 setup1d(subEntity);
307
308 } else if (dim==2 and sizes_[0]>1 and sizes_[1]>1) {
309
310 setup2d(subEntity);
311
312 }
313 }
314
315 for (size_t i=0; i<li_.size(); i++)
316 li_[i] = LocalKey(subEntity[i], codim[i], index[i]);
317 }
318
320 std::size_t size () const
321 {
322 return std::accumulate(sizes_.begin(), sizes_.end(), 1, std::multiplies<unsigned int>());
323 }
324
326 const LocalKey& localKey (std::size_t i) const
327 {
328 return li_[i];
329 }
330
331private:
332
333 // Number of shape functions on this element per coordinate direction
334 std::array<unsigned, dim> sizes_;
335
336 std::vector<LocalKey> li_;
337};
338
343template<int dim, class LB>
345{
346public:
348 template<typename F, typename C>
349 void interpolate (const F& f, std::vector<C>& out) const
350 {
351 DUNE_THROW(NotImplemented, "BSplineLocalInterpolation::interpolate");
352 }
353
354};
355
365template<class GV, class R>
367{
368 typedef typename GV::ctype D;
369 enum {dim = GV::dimension};
370 friend class BSplineLocalBasis<GV,R>;
371public:
372
378
382 : preBasis_(preBasis),
383 localBasis_(preBasis,*this)
384 {}
385
389 : preBasis_(other.preBasis_),
390 localBasis_(preBasis_,*this)
391 {}
392
399 void bind(const std::array<unsigned,dim>& elementIdx)
400 {
401 /* \todo In the long run we need to precompute a table for this */
402 for (size_t i=0; i<elementIdx.size(); i++)
403 {
404 currentKnotSpan_[i] = 0;
405
406 // Skip over degenerate knot spans
407 while (preBasis_.knotVectors_[i][currentKnotSpan_[i]+1] < preBasis_.knotVectors_[i][currentKnotSpan_[i]]+1e-8)
408 currentKnotSpan_[i]++;
409
410 for (size_t j=0; j<elementIdx[i]; j++)
411 {
412 currentKnotSpan_[i]++;
413
414 // Skip over degenerate knot spans
415 while (preBasis_.knotVectors_[i][currentKnotSpan_[i]+1] < preBasis_.knotVectors_[i][currentKnotSpan_[i]]+1e-8)
416 currentKnotSpan_[i]++;
417 }
418
419 // Compute the geometric transformation from knotspan-local to global coordinates
420 localBasis_.offset_[i] = preBasis_.knotVectors_[i][currentKnotSpan_[i]];
421 localBasis_.scaling_[i][i] = preBasis_.knotVectors_[i][currentKnotSpan_[i]+1] - preBasis_.knotVectors_[i][currentKnotSpan_[i]];
422 }
423
424 // Set up the LocalCoefficients object
425 std::array<unsigned int, dim> sizes;
426 for (size_t i=0; i<dim; i++)
427 sizes[i] = size(i);
428 localCoefficients_.init(sizes);
429 }
430
433 {
434 return localBasis_;
435 }
436
439 {
440 return localCoefficients_;
441 }
442
445 {
446 return localInterpolation_;
447 }
448
450 unsigned size () const
451 {
452 std::size_t r = 1;
453 for (int i=0; i<dim; i++)
454 r *= size(i);
455 return r;
456 }
457
461 {
462 return GeometryTypes::cube(dim);
463 }
464
465//private:
466
468 unsigned int size(int i) const
469 {
470 const auto& order = preBasis_.order_;
471 unsigned int r = order[i]+1; // The 'normal' value
472 if (currentKnotSpan_[i]<order[i]) // Less near the left end of the knot vector
473 r -= (order[i] - currentKnotSpan_[i]);
474 if ( order[i] > (preBasis_.knotVectors_[i].size() - currentKnotSpan_[i] - 2) )
475 r -= order[i] - (preBasis_.knotVectors_[i].size() - currentKnotSpan_[i] - 2);
476 return r;
477 }
478
479 const BSplinePreBasis<GV>& preBasis_;
480
481 BSplineLocalBasis<GV,R> localBasis_;
482 BSplineLocalCoefficients<dim> localCoefficients_;
484
485 // The knot span we are bound to
486 std::array<unsigned,dim> currentKnotSpan_;
487};
488
489
490template<typename GV>
491class BSplineNode;
492
502template<typename GV>
504 public LeafPreBasisMixin< BSplinePreBasis<GV> >
505{
507
508 static const int dim = GV::dimension;
509
511 class MultiDigitCounter
512 {
513 public:
514
518 MultiDigitCounter(const std::array<unsigned int,dim>& limits)
519 : limits_(limits)
520 {
521 std::fill(counter_.begin(), counter_.end(), 0);
522 }
523
525 MultiDigitCounter& operator++()
526 {
527 for (int i=0; i<dim; i++)
528 {
529 ++counter_[i];
530
531 // no overflow?
532 if (counter_[i] < limits_[i])
533 break;
534
535 counter_[i] = 0;
536 }
537 return *this;
538 }
539
541 const unsigned int& operator[](int i) const
542 {
543 return counter_[i];
544 }
545
547 unsigned int cycle() const
548 {
549 unsigned int r = 1;
550 for (int i=0; i<dim; i++)
551 r *= limits_[i];
552 return r;
553 }
554
555 private:
556
558 const std::array<unsigned int,dim> limits_;
559
561 std::array<unsigned int,dim> counter_;
562
563 };
564
565public:
566
568 using GridView = GV;
569 using size_type = std::size_t;
570
571 using Node = BSplineNode<GV>;
572
573 // Type used for function values
574 using R = double;
575
595 const std::vector<double>& knotVector,
596 unsigned int order,
597 bool makeOpen = true)
598 : gridView_(gridView)
599 {
600 // \todo Detection of duplicate knots
601 std::fill(elements_.begin(), elements_.end(), knotVector.size()-1);
602
603 // Mediocre sanity check: we don't know the number of grid elements in each direction.
604 // but at least we know the total number of elements.
605 assert( std::accumulate(elements_.begin(), elements_.end(), 1, std::multiplies<unsigned>()) == gridView_.size(0) );
606
607 for (int i=0; i<dim; i++)
608 {
609 // Prepend the correct number of additional knots to open the knot vector
611 if (makeOpen)
612 for (unsigned int j=0; j<order; j++)
613 knotVectors_[i].push_back(knotVector[0]);
614
615 knotVectors_[i].insert(knotVectors_[i].end(), knotVector.begin(), knotVector.end());
616
617 if (makeOpen)
618 for (unsigned int j=0; j<order; j++)
619 knotVectors_[i].push_back(knotVector.back());
620 }
621
622 std::fill(order_.begin(), order_.end(), order);
623 }
624
647 const FieldVector<double,dim>& lowerLeft,
648 const FieldVector<double,dim>& upperRight,
649 const std::array<unsigned int,dim>& elements,
650 unsigned int order,
651 bool makeOpen = true)
652 : elements_(elements),
653 gridView_(gridView)
654 {
655 // Mediocre sanity check: we don't know the number of grid elements in each direction.
656 // but at least we know the total number of elements.
657 assert( std::accumulate(elements_.begin(), elements_.end(), 1, std::multiplies<unsigned>()) == gridView_.size(0) );
658
659 for (int i=0; i<dim; i++)
660 {
661 // Prepend the correct number of additional knots to open the knot vector
663 if (makeOpen)
664 for (unsigned int j=0; j<order; j++)
665 knotVectors_[i].push_back(lowerLeft[i]);
666
667 // Construct the actual knot vector
668 for (size_t j=0; j<elements[i]+1; j++)
669 knotVectors_[i].push_back(lowerLeft[i] + j*(upperRight[i]-lowerLeft[i]) / elements[i]);
670
671 if (makeOpen)
672 for (unsigned int j=0; j<order; j++)
673 knotVectors_[i].push_back(upperRight[i]);
674 }
675
676 std::fill(order_.begin(), order_.end(), order);
677 }
678
681 {}
682
684 const GridView& gridView() const
685 {
686 return gridView_;
687 }
688
690 void update(const GridView& gv)
691 {
692 gridView_ = gv;
693 }
694
698 Node makeNode() const
699 {
700 return Node{this};
701 }
702
704 size_type maxNodeSize() const
705 {
706 size_type result = 1;
707 for (int i=0; i<dim; i++)
708 result *= order_[i]+1;
709 return result;
710 }
711
713 template<typename It>
714 It indices(const Node& node, It it) const
715 {
716 // Local degrees of freedom are arranged in a lattice.
717 // We need the lattice dimensions to be able to compute lattice coordinates from a local index
718 std::array<unsigned int, dim> localSizes;
719 for (int i=0; i<dim; i++)
720 localSizes[i] = node.finiteElement().size(i);
721 for (size_type i = 0, end = node.size() ; i < end ; ++i, ++it)
722 {
723 std::array<unsigned int,dim> localIJK = getIJK(i, localSizes);
724
725 const auto currentKnotSpan = node.finiteElement().currentKnotSpan_;
726 const auto order = order_;
727
728 std::array<unsigned int,dim> globalIJK;
729 for (int i=0; i<dim; i++)
730 globalIJK[i] = std::max((int)currentKnotSpan[i] - (int)order[i], 0) + localIJK[i]; // needs to be a signed type!
731
732 // Make one global flat index from the globalIJK tuple
733 size_type globalIdx = globalIJK[dim-1];
734
735 for (int i=dim-2; i>=0; i--)
736 globalIdx = globalIdx * size(i) + globalIJK[i];
737
738 *it = {{globalIdx}};
739 }
740 return it;
741 }
742
744 unsigned int dimension () const
745 {
746 unsigned int result = 1;
747 for (size_t i=0; i<dim; i++)
748 result *= size(i);
749 return result;
750 }
751
753 unsigned int size (size_t d) const
754 {
755 return knotVectors_[d].size() - order_[d] - 1;
756 }
757
758 using Base::size;
759
763 std::vector<FieldVector<R,1> >& out,
764 const std::array<unsigned,dim>& currentKnotSpan) const
765 {
766 // Evaluate
767 std::array<std::vector<R>, dim> oneDValues;
768
769 for (size_t i=0; i<dim; i++)
770 evaluateFunction(in[i], oneDValues[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
771
772 std::array<unsigned int, dim> limits;
773 for (int i=0; i<dim; i++)
774 limits[i] = oneDValues[i].size();
775
776 MultiDigitCounter ijkCounter(limits);
777
778 out.resize(ijkCounter.cycle());
779
780 for (size_t i=0; i<out.size(); i++, ++ijkCounter)
781 {
782 out[i] = R(1.0);
783 for (size_t j=0; j<dim; j++)
784 out[i] *= oneDValues[j][ijkCounter[j]];
785 }
786 }
787
794 std::vector<FieldMatrix<R,1,dim> >& out,
795 const std::array<unsigned,dim>& currentKnotSpan) const
796 {
797 // How many shape functions to we have in each coordinate direction?
798 std::array<unsigned int, dim> limits;
799 for (int i=0; i<dim; i++)
800 {
801 limits[i] = order_[i]+1; // The 'standard' value away from the boundaries of the knot vector
802 if (currentKnotSpan[i]<order_[i])
803 limits[i] -= (order_[i] - currentKnotSpan[i]);
804 if ( order_[i] > (knotVectors_[i].size() - currentKnotSpan[i] - 2) )
805 limits[i] -= order_[i] - (knotVectors_[i].size() - currentKnotSpan[i] - 2);
806 }
807
808 // The lowest knot spans that we need values from
809 std::array<unsigned int, dim> offset;
810 for (int i=0; i<dim; i++)
811 offset[i] = std::max((int)(currentKnotSpan[i] - order_[i]),0);
812
813 // Evaluate 1d function values (needed for the product rule)
814 std::array<std::vector<R>, dim> oneDValues;
815
816 // Evaluate 1d function values of one order lower (needed for the derivative formula)
817 std::array<std::vector<R>, dim> lowOrderOneDValues;
818
819 std::array<DynamicMatrix<R>, dim> values;
820
821 for (size_t i=0; i<dim; i++)
822 {
823 evaluateFunctionFull(in[i], values[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
824 oneDValues[i].resize(knotVectors_[i].size()-order_[i]-1);
825 for (size_t j=0; j<oneDValues[i].size(); j++)
826 oneDValues[i][j] = values[i][order_[i]][j];
827
828 if (order_[i]!=0)
829 {
830 lowOrderOneDValues[i].resize(knotVectors_[i].size()-(order_[i]-1)-1);
831 for (size_t j=0; j<lowOrderOneDValues[i].size(); j++)
832 lowOrderOneDValues[i][j] = values[i][order_[i]-1][j];
833 }
834 }
835
836
837 // Evaluate 1d function derivatives
838 std::array<std::vector<R>, dim> oneDDerivatives;
839 for (size_t i=0; i<dim; i++)
840 {
841 oneDDerivatives[i].resize(limits[i]);
842
843 if (order_[i]==0) // order-zero functions are piecewise constant, hence all derivatives are zero
844 std::fill(oneDDerivatives[i].begin(), oneDDerivatives[i].end(), R(0.0));
845 else
846 {
847 for (size_t j=offset[i]; j<offset[i]+limits[i]; j++)
848 {
849 R derivativeAddend1 = lowOrderOneDValues[i][j] / (knotVectors_[i][j+order_[i]]-knotVectors_[i][j]);
850 R derivativeAddend2 = lowOrderOneDValues[i][j+1] / (knotVectors_[i][j+order_[i]+1]-knotVectors_[i][j+1]);
851 // The two previous terms may evaluate as 0/0. This is to be interpreted as 0.
852 if (std::isnan(derivativeAddend1))
853 derivativeAddend1 = 0;
854 if (std::isnan(derivativeAddend2))
855 derivativeAddend2 = 0;
856 oneDDerivatives[i][j-offset[i]] = order_[i] * ( derivativeAddend1 - derivativeAddend2 );
857 }
858 }
859 }
860
861 // Working towards computing only the parts that we really need:
862 // Let's copy them out into a separate array
863 std::array<std::vector<R>, dim> oneDValuesShort;
864
865 for (int i=0; i<dim; i++)
866 {
867 oneDValuesShort[i].resize(limits[i]);
868
869 for (size_t j=0; j<limits[i]; j++)
870 oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
871 }
872
873
874
875 // Set up a multi-index to go from consecutive indices to integer coordinates
876 MultiDigitCounter ijkCounter(limits);
877
878 out.resize(ijkCounter.cycle());
879
880 // Complete Jacobian is given by the product rule
881 for (size_t i=0; i<out.size(); i++, ++ijkCounter)
882 for (int j=0; j<dim; j++)
883 {
884 out[i][0][j] = 1.0;
885 for (int k=0; k<dim; k++)
886 out[i][0][j] *= (j==k) ? oneDDerivatives[k][ijkCounter[k]]
887 : oneDValuesShort[k][ijkCounter[k]];
888 }
889
890 }
891
893 template <size_type k>
894 void evaluate(const typename std::array<int,k>& directions,
896 std::vector<FieldVector<R,1> >& out,
897 const std::array<unsigned,dim>& currentKnotSpan) const
898 {
899 if (k != 1 && k != 2)
900 DUNE_THROW(RangeError, "Differentiation order greater than 2 is not supported!");
901
902 // Evaluate 1d function values (needed for the product rule)
903 std::array<std::vector<R>, dim> oneDValues;
904 std::array<std::vector<R>, dim> oneDDerivatives;
905 std::array<std::vector<R>, dim> oneDSecondDerivatives;
906
907 // Evaluate 1d function derivatives
908 if (k==1)
909 for (size_t i=0; i<dim; i++)
910 evaluateAll(in[i], oneDValues[i], true, oneDDerivatives[i], false, oneDSecondDerivatives[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
911 else
912 for (size_t i=0; i<dim; i++)
913 evaluateAll(in[i], oneDValues[i], true, oneDDerivatives[i], true, oneDSecondDerivatives[i], knotVectors_[i], order_[i], currentKnotSpan[i]);
914
915 // The lowest knot spans that we need values from
916 std::array<unsigned int, dim> offset;
917 for (int i=0; i<dim; i++)
918 offset[i] = std::max((int)(currentKnotSpan[i] - order_[i]),0);
919
920 // Set up a multi-index to go from consecutive indices to integer coordinates
921 std::array<unsigned int, dim> limits;
922 for (int i=0; i<dim; i++)
923 {
924 // In a proper implementation, the following line would do
925 //limits[i] = oneDValues[i].size();
926 limits[i] = order_[i]+1; // The 'standard' value away from the boundaries of the knot vector
927 if (currentKnotSpan[i]<order_[i])
928 limits[i] -= (order_[i] - currentKnotSpan[i]);
929 if ( order_[i] > (knotVectors_[i].size() - currentKnotSpan[i] - 2) )
930 limits[i] -= order_[i] - (knotVectors_[i].size() - currentKnotSpan[i] - 2);
931 }
932
933 // Working towards computing only the parts that we really need:
934 // Let's copy them out into a separate array
935 std::array<std::vector<R>, dim> oneDValuesShort;
936
937 for (int i=0; i<dim; i++)
938 {
939 oneDValuesShort[i].resize(limits[i]);
940
941 for (size_t j=0; j<limits[i]; j++)
942 oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
943 }
944
945
946 MultiDigitCounter ijkCounter(limits);
947
948 out.resize(ijkCounter.cycle());
949
950 if (k == 1)
951 {
952 // Complete Jacobian is given by the product rule
953 for (size_t i=0; i<out.size(); i++, ++ijkCounter)
954 {
955 out[i][0] = 1.0;
956 for (int l=0; l<dim; l++)
957 out[i][0] *= (directions[0]==l) ? oneDDerivatives[l][ijkCounter[l]]
958 : oneDValuesShort[l][ijkCounter[l]];
959 }
960 }
961
962 if (k == 2)
963 {
964 // Complete derivation by deriving the tensor product
965 for (size_t i=0; i<out.size(); i++, ++ijkCounter)
966 {
967 out[i][0] = 1.0;
968 for (int j=0; j<dim; j++)
969 {
970 if (directions[0] != directions[1]) //derivation in two different variables
971 if (directions[0] == j || directions[1] == j) //the spline has to be derived (once) in this direction
972 out[i][0] *= oneDDerivatives[j][ijkCounter[j]];
973 else //no derivation in this direction
974 out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
975 else //spline is derived two times in the same direction
976 if (directions[0] == j) //the spline is derived two times in this direction
977 out[i][0] *= oneDSecondDerivatives[j][ijkCounter[j]];
978 else //no derivation in this direction
979 out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
980 }
981 }
982 }
983 }
984
985
990 static std::array<unsigned int,dim> getIJK(typename GridView::IndexSet::IndexType idx, std::array<unsigned int,dim> elements)
991 {
992 std::array<unsigned,dim> result;
993 for (int i=0; i<dim; i++)
994 {
995 result[i] = idx%elements[i];
996 idx /= elements[i];
997 }
998 return result;
999 }
1000
1009 static void evaluateFunction (const typename GV::ctype& in, std::vector<R>& out,
1010 const std::vector<R>& knotVector,
1011 unsigned int order,
1012 unsigned int currentKnotSpan)
1013 {
1014 std::size_t outSize = order+1; // The 'standard' value away from the boundaries of the knot vector
1015 if (currentKnotSpan<order) // Less near the left end of the knot vector
1016 outSize -= (order - currentKnotSpan);
1017 if ( order > (knotVector.size() - currentKnotSpan - 2) )
1018 outSize -= order - (knotVector.size() - currentKnotSpan - 2);
1019 out.resize(outSize);
1020
1021 // It's not really a matrix that is needed here, a plain 2d array would do
1022 DynamicMatrix<R> N(order+1, knotVector.size()-1);
1023
1024 // The text books on splines use the following geometric condition here to fill the array N
1025 // (see for example Cottrell, Hughes, Bazilevs, Formula (2.1). However, this condition
1026 // only works if splines are never evaluated exactly on the knots.
1027 //
1028 // for (size_t i=0; i<knotVector.size()-1; i++)
1029 // N[0][i] = (knotVector[i] <= in) and (in < knotVector[i+1]);
1030 for (size_t i=0; i<knotVector.size()-1; i++)
1031 N[0][i] = (i == currentKnotSpan);
1032
1033 for (size_t r=1; r<=order; r++)
1034 for (size_t i=0; i<knotVector.size()-r-1; i++)
1035 {
1036 R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
1037 ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
1038 : 0;
1039 R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
1040 ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
1041 : 0;
1042 N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
1043 }
1044
1049 for (size_t i=0; i<out.size(); i++) {
1050 out[i] = N[order][std::max((int)(currentKnotSpan - order),0) + i];
1051 }
1052 }
1053
1066 static void evaluateFunctionFull(const typename GV::ctype& in,
1067 DynamicMatrix<R>& out,
1068 const std::vector<R>& knotVector,
1069 unsigned int order,
1070 unsigned int currentKnotSpan)
1071 {
1072 // It's not really a matrix that is needed here, a plain 2d array would do
1073 DynamicMatrix<R>& N = out;
1074
1075 N.resize(order+1, knotVector.size()-1);
1076
1077 // The text books on splines use the following geometric condition here to fill the array N
1078 // (see for example Cottrell, Hughes, Bazilevs, Formula (2.1). However, this condition
1079 // only works if splines are never evaluated exactly on the knots.
1080 //
1081 // for (size_t i=0; i<knotVector.size()-1; i++)
1082 // N[0][i] = (knotVector[i] <= in) and (in < knotVector[i+1]);
1083 for (size_t i=0; i<knotVector.size()-1; i++)
1084 N[0][i] = (i == currentKnotSpan);
1085
1086 for (size_t r=1; r<=order; r++)
1087 for (size_t i=0; i<knotVector.size()-r-1; i++)
1088 {
1089 R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
1090 ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
1091 : 0;
1092 R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
1093 ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
1094 : 0;
1095 N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
1096 }
1097 }
1098
1099
1109 static void evaluateAll(const typename GV::ctype& in,
1110 std::vector<R>& out,
1111 bool evaluateJacobian, std::vector<R>& outJac,
1112 bool evaluateHessian, std::vector<R>& outHess,
1113 const std::vector<R>& knotVector,
1114 unsigned int order,
1115 unsigned int currentKnotSpan)
1116 {
1117 // How many shape functions to we have in each coordinate direction?
1118 unsigned int limit;
1119 limit = order+1; // The 'standard' value away from the boundaries of the knot vector
1120 if (currentKnotSpan<order)
1121 limit -= (order - currentKnotSpan);
1122 if ( order > (knotVector.size() - currentKnotSpan - 2) )
1123 limit -= order - (knotVector.size() - currentKnotSpan - 2);
1124
1125 // The lowest knot spans that we need values from
1126 unsigned int offset;
1127 offset = std::max((int)(currentKnotSpan - order),0);
1128
1129 // Evaluate 1d function values (needed for the product rule)
1130 DynamicMatrix<R> values;
1131
1132 evaluateFunctionFull(in, values, knotVector, order, currentKnotSpan);
1133
1134 out.resize(knotVector.size()-order-1);
1135 for (size_t j=0; j<out.size(); j++)
1136 out[j] = values[order][j];
1137
1138 // Evaluate 1d function values of one order lower (needed for the derivative formula)
1139 std::vector<R> lowOrderOneDValues;
1140
1141 if (order!=0)
1142 {
1143 lowOrderOneDValues.resize(knotVector.size()-(order-1)-1);
1144 for (size_t j=0; j<lowOrderOneDValues.size(); j++)
1145 lowOrderOneDValues[j] = values[order-1][j];
1146 }
1147
1148 // Evaluate 1d function values of two order lower (needed for the (second) derivative formula)
1149 std::vector<R> lowOrderTwoDValues;
1150
1151 if (order>1 && evaluateHessian)
1152 {
1153 lowOrderTwoDValues.resize(knotVector.size()-(order-2)-1);
1154 for (size_t j=0; j<lowOrderTwoDValues.size(); j++)
1155 lowOrderTwoDValues[j] = values[order-2][j];
1156 }
1157
1158 // Evaluate 1d function derivatives
1159 if (evaluateJacobian)
1160 {
1161 outJac.resize(limit);
1162
1163 if (order==0) // order-zero functions are piecewise constant, hence all derivatives are zero
1164 std::fill(outJac.begin(), outJac.end(), R(0.0));
1165 else
1166 {
1167 for (size_t j=offset; j<offset+limit; j++)
1168 {
1169 R derivativeAddend1 = lowOrderOneDValues[j] / (knotVector[j+order]-knotVector[j]);
1170 R derivativeAddend2 = lowOrderOneDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]);
1171 // The two previous terms may evaluate as 0/0. This is to be interpreted as 0.
1172 if (std::isnan(derivativeAddend1))
1173 derivativeAddend1 = 0;
1174 if (std::isnan(derivativeAddend2))
1175 derivativeAddend2 = 0;
1176 outJac[j-offset] = order * ( derivativeAddend1 - derivativeAddend2 );
1177 }
1178 }
1179 }
1180
1181 // Evaluate 1d function second derivatives
1182 if (evaluateHessian)
1183 {
1184 outHess.resize(limit);
1185
1186 if (order<2) // order-zero functions are piecewise constant, hence all derivatives are zero
1187 std::fill(outHess.begin(), outHess.end(), R(0.0));
1188 else
1189 {
1190 for (size_t j=offset; j<offset+limit; j++)
1191 {
1192 assert(j+2 < lowOrderTwoDValues.size());
1193 R derivativeAddend1 = lowOrderTwoDValues[j] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order-1]-knotVector[j]);
1194 R derivativeAddend2 = lowOrderTwoDValues[j+1] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order]-knotVector[j+1]);
1195 R derivativeAddend3 = lowOrderTwoDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+order]-knotVector[j+1]);
1196 R derivativeAddend4 = lowOrderTwoDValues[j+2] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+1+order]-knotVector[j+2]);
1197 // The two previous terms may evaluate as 0/0. This is to be interpreted as 0.
1198
1199 if (std::isnan(derivativeAddend1))
1200 derivativeAddend1 = 0;
1201 if (std::isnan(derivativeAddend2))
1202 derivativeAddend2 = 0;
1203 if (std::isnan(derivativeAddend3))
1204 derivativeAddend3 = 0;
1205 if (std::isnan(derivativeAddend4))
1206 derivativeAddend4 = 0;
1207 outHess[j-offset] = order * (order-1) * ( derivativeAddend1 - derivativeAddend2 -derivativeAddend3 + derivativeAddend4 );
1208 }
1209 }
1210 }
1211 }
1212
1213
1215 std::array<unsigned int, dim> order_;
1216
1218 std::array<std::vector<double>, dim> knotVectors_;
1219
1221 std::array<unsigned,dim> elements_;
1222
1223 GridView gridView_;
1224};
1225
1226
1227
1228template<typename GV>
1229class BSplineNode :
1230 public LeafBasisNode
1231{
1232 static const int dim = GV::dimension;
1233
1234public:
1235
1236 using size_type = std::size_t;
1237 using Element = typename GV::template Codim<0>::Entity;
1238 using FiniteElement = BSplineLocalFiniteElement<GV,double>;
1239
1240 BSplineNode(const BSplinePreBasis<GV>* preBasis) :
1241 preBasis_(preBasis),
1242 finiteElement_(*preBasis)
1243 {}
1244
1246 const Element& element() const
1247 {
1248 return element_;
1249 }
1250
1255 const FiniteElement& finiteElement() const
1256 {
1257 return finiteElement_;
1258 }
1259
1261 void bind(const Element& e)
1262 {
1263 element_ = e;
1264 auto elementIndex = preBasis_->gridView().indexSet().index(e);
1265 finiteElement_.bind(preBasis_->getIJK(elementIndex,preBasis_->elements_));
1266 this->setSize(finiteElement_.size());
1267 }
1268
1269protected:
1270
1271 const BSplinePreBasis<GV>* preBasis_;
1272
1273 FiniteElement finiteElement_;
1274 Element element_;
1275};
1276
1277
1278
1279namespace BasisFactory {
1280
1287inline auto bSpline(const std::vector<double>& knotVector,
1288 unsigned int order,
1289 bool makeOpen = true)
1290{
1291 return [&knotVector, order, makeOpen](const auto& gridView) {
1292 return BSplinePreBasis<std::decay_t<decltype(gridView)>>(gridView, knotVector, order, makeOpen);
1293 };
1294}
1295
1296} // end namespace BasisFactory
1297
1298// *****************************************************************************
1299// This is the actual global basis implementation based on the reusable parts.
1300// *****************************************************************************
1301
1308template<typename GV>
1310
1311
1312} // namespace Functions
1313
1314} // namespace Dune
1315
1316#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
Construct a matrix with a dynamic size.
Definition: dynmatrix.hh:61
void resize(size_type r, size_type c, value_type v=value_type())
resize matrix to r × c
Definition: dynmatrix.hh:106
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
LocalBasis class in the sense of dune-localfunctions, presenting the restriction of a B-spline patch ...
Definition: bsplinebasis.hh:51
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:60
unsigned int order() const
Polynomial order of the shape functions.
Definition: bsplinebasis.hh:145
std::size_t size() const
Return the number of basis functions on the current knot span.
Definition: bsplinebasis.hh:152
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:102
void evaluateFunction(const FieldVector< D, dim > &in, std::vector< FieldVector< R, 1 > > &out) const
Evaluate all shape functions.
Definition: bsplinebasis.hh:75
void evaluateJacobian(const FieldVector< D, dim > &in, std::vector< FieldMatrix< D, 1, dim > > &out) const
Evaluate Jacobian of all shape functions.
Definition: bsplinebasis.hh:87
BSplineLocalBasis(const BSplinePreBasis< GV > &preBasis, const BSplineLocalFiniteElement< GV, R > &lFE)
Constructor with a given B-spline patch.
Definition: bsplinebasis.hh:66
Attaches a shape function to an entity.
Definition: bsplinebasis.hh:183
const LocalKey & localKey(std::size_t i) const
get i'th index
Definition: bsplinebasis.hh:326
std::size_t size() const
number of coefficients
Definition: bsplinebasis.hh:320
LocalFiniteElement in the sense of dune-localfunctions, for the B-spline basis on tensor-product grid...
Definition: bsplinebasis.hh:367
BSplineLocalFiniteElement(const BSplineLocalFiniteElement &other)
Copy constructor.
Definition: bsplinebasis.hh:388
const BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > & localInterpolation() const
Hand out a LocalInterpolation object.
Definition: bsplinebasis.hh:444
LocalFiniteElementTraits< BSplineLocalBasis< GV, R >, BSplineLocalCoefficients< dim >, BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > > Traits
Export various types related to this LocalFiniteElement.
Definition: bsplinebasis.hh:377
BSplineLocalFiniteElement(const BSplinePreBasis< GV > &preBasis)
Constructor with a given B-spline basis.
Definition: bsplinebasis.hh:381
const BSplineLocalCoefficients< dim > & localCoefficients() const
Hand out a LocalCoefficients object.
Definition: bsplinebasis.hh:438
void bind(const std::array< unsigned, dim > &elementIdx)
Bind LocalFiniteElement to a specific knot span of the spline patch.
Definition: bsplinebasis.hh:399
GeometryType type() const
Return the reference element that the local finite element is defined on (here, a hypercube)
Definition: bsplinebasis.hh:460
unsigned size() const
Number of shape functions in this finite element.
Definition: bsplinebasis.hh:450
unsigned int size(int i) const
Number of degrees of freedom for one coordinate direction.
Definition: bsplinebasis.hh:468
const BSplineLocalBasis< GV, R > & localBasis() const
Hand out a LocalBasis object.
Definition: bsplinebasis.hh:432
Local interpolation in the sense of dune-localfunctions, for the B-spline basis on tensor-product gri...
Definition: bsplinebasis.hh:345
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: bsplinebasis.hh:349
Pre-basis for B-spline basis.
Definition: bsplinebasis.hh:505
std::array< unsigned, dim > elements_
Number of grid elements in the different coordinate directions.
Definition: bsplinebasis.hh:1221
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:1066
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:762
std::array< unsigned int, dim > order_
Order of the B-spline for each space dimension.
Definition: bsplinebasis.hh:1215
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:1109
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:1009
unsigned int size(size_t d) const
Number of shape functions in one direction.
Definition: bsplinebasis.hh:753
GV GridView
The grid view that the FE space is defined on.
Definition: bsplinebasis.hh:568
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:714
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:894
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: bsplinebasis.hh:690
unsigned int dimension() const
Total number of B-spline basis functions.
Definition: bsplinebasis.hh:744
void initializeIndices()
Initialize the global indices.
Definition: bsplinebasis.hh:680
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: bsplinebasis.hh:684
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:990
Node makeNode() const
Create tree node.
Definition: bsplinebasis.hh:698
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:594
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:793
std::array< std::vector< double >, dim > knotVectors_
The knot vectors, one for each space dimension.
Definition: bsplinebasis.hh:1218
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition: bsplinebasis.hh:704
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:646
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:50
A generic MixIn class for PreBasis.
Definition: leafprebasismixin.hh:36
size_type size() const
Get the total dimension of the space spanned by this basis.
Definition: leafprebasismixin.hh:60
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
Describe position of one degree of freedom.
Definition: localkey.hh:24
Default exception for dummy implementations.
Definition: exceptions.hh:263
Default exception class for range errors.
Definition: exceptions.hh:254
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:296
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
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:1287
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:462
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integer_sequence< T, II..., T(IN)> push_back(std::integer_sequence< T, II... >, std::integral_constant< T, IN >={})
Append an index IN to the back of the sequence.
Definition: integersequence.hh:69
Static tag representing a codimension.
Definition: dimension.hh:24
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:35
D DomainType
domain type
Definition: localbasis.hh:43
traits helper struct
Definition: localfiniteelementtraits.hh:13
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 23, 23:29, 2024)