DUNE PDELab (git)

cubichermitebasis.hh
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_CUBICHERMITEBASIS_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_CUBICHERMITEBASIS_HH
9
12
15#include <dune/grid/common/rangegenerators.hh>
16
17#include <dune/localfunctions/common/localfiniteelementvariant.hh>
19#include <dune/localfunctions/crouzeixraviart.hh>
20
21#include <dune/functions/analyticfunctions/monomialset.hh>
22#include <dune/functions/functionspacebases/nodes.hh>
23#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
24#include <dune/functions/functionspacebases/leafprebasismappermixin.hh>
25#include <dune/functions/common/mapperutilities.hh>
26
27namespace Dune {
28namespace Functions {
29
30namespace Impl {
31
32// *****************************************************************************
33// * CubicHermiteLocalFiniteElement
34// *****************************************************************************
35
36
37
38template<class DF, class RF, unsigned int dim, bool reduced>
39class CubicHermiteLocalBasis
40{
41
42 static constexpr auto makeReferenceBasisCoefficients() {
43 if constexpr (dim==1)
45 { 1, 0, -3, 2},
46 { 0, 1, -2, 1},
47 { 0, 0, 3, -2},
48 { 0, 0, -1, 1}
49 };
50 if constexpr ((dim==2) and (not reduced))
52 { 1, 0, 0, -3, -13, -3, 2, 13, 13, 2},
53 { 0, 1, 0, -2, -3, 0, 1, 3, 2, 0},
54 { 0, 0, 1, 0, -3, -2, 0, 2, 3, 1},
55 { 0, 0, 0, 3, -7, 0, -2, 7, 7, 0},
56 { 0, 0, 0, -1, 2, 0, 1, -2, -2, 0},
57 { 0, 0, 0, 0, -1, 0, 0, 2, 1, 0},
58 { 0, 0, 0, 0, -7, 3, 0, 7, 7, -2},
59 { 0, 0, 0, 0, -1, 0, 0, 1, 2, 0},
60 { 0, 0, 0, 0, 2, -1, 0, -2, -2, 1},
61 { 0, 0, 0, 0, 27, 0, 0, -27, -27, 0}
62 };
63 if constexpr ((dim==2) and (reduced))
64 {
65 auto w = std::array{1./3, 1./18, 1./18, 1./3, -1./9, 1./18, 1./3, 1./18, -1./9};
67 { 1, 0, 0, -3, -13 + w[0]*27, -3, 2, 13 - w[0]*27, 13 - w[0]*27, 2},
68 { 0, 1, 0, -2, -3 + w[1]*27, 0, 1, 3 - w[1]*27, 2 - w[1]*27, 0},
69 { 0, 0, 1, 0, -3 + w[2]*27, -2, 0, 2 - w[2]*27, 3 - w[2]*27, 1},
70 { 0, 0, 0, 3, -7 + w[3]*27, 0, -2, 7 - w[3]*27, 7 - w[3]*27, 0},
71 { 0, 0, 0, -1, 2 + w[4]*27, 0, 1, -2 - w[4]*27, -2 - w[4]*27, 0},
72 { 0, 0, 0, 0, -1 + w[5]*27, 0, 0, 2 - w[5]*27, 1 - w[5]*27, 0},
73 { 0, 0, 0, 0, -7 + w[6]*27, 3, 0, 7 - w[6]*27, 7 - w[6]*27, -2},
74 { 0, 0, 0, 0, -1 + w[7]*27, 0, 0, 1 - w[7]*27, 2 - w[7]*27, 0},
75 { 0, 0, 0, 0, 2 + w[8]*27, -1, 0, -2 - w[8]*27, -2 - w[8]*27, 1},
76 };
77 }
78 }
79
80 // These are the coefficients of the cubic Hermite basis functions
81 // on the reference element wrt. the monomials. These have been computed
82 // by solving with the corresponiding Vandermonde-matrix for the reference
83 // element in advance.
84 // The basis functions can be evaluated by first evaluating the monomials
85 // and then transforming their values with these coefficients using
86 //
87 // referenceBasisCoefficients.mv(evaluateMonomialValues(x), values);
88 //
89 // Surprisingly, storing them as static constexpr member is slightly faster
90 // than using the static constexpr function directly.
91 static constexpr auto referenceBasisCoefficients = makeReferenceBasisCoefficients();
92
93 // This transforms the function or derivative values from the basis
94 // functions on the reference element to those on the grid element.
95 // Since Hermite elements do not form an affine family, the transformation
96 // of derivative DOFs involves the Jacobian of the grid element transformation.
97 // To avoid blowup of condition numbers due to h-dependend basis functions,
98 // an h-dependent rescaling is applied.
99 template<class LambdaRefValues, class Entry>
100 void transformToElementBasis(const LambdaRefValues& refValues, std::vector<Entry>& out) const
101 {
102 if constexpr (dim==1)
103 {
104 const auto& J = elementJacobian_;
105 out.resize(refValues.size());
106 out[0][0] = refValues[0];
107 out[1][0] = J*refValues[1] / (*localSubEntityMeshSize_)[1];
108 out[2][0] = refValues[2];
109 out[3][0] = J*refValues[3] / (*localSubEntityMeshSize_)[1];;
110 }
111 if constexpr (dim==2)
112 {
113 const auto& J = elementJacobian_;
114 out.resize(refValues.size());
115 out[0][0] = refValues[0];
116 out[1][0] = (J[0][0]*refValues[1] + J[0][1]*refValues[2]) / (*localSubEntityMeshSize_)[1];
117 out[2][0] = (J[1][0]*refValues[1] + J[1][1]*refValues[2]) / (*localSubEntityMeshSize_)[2];
118 out[3][0] = refValues[3];
119 out[4][0] = (J[0][0]*refValues[4] + J[0][1]*refValues[5]) / (*localSubEntityMeshSize_)[4];
120 out[5][0] = (J[1][0]*refValues[4] + J[1][1]*refValues[5]) / (*localSubEntityMeshSize_)[5];
121 out[6][0] = refValues[6];
122 out[7][0] = (J[0][0]*refValues[7] + J[0][1]*refValues[8]) / (*localSubEntityMeshSize_)[7];
123 out[8][0] = (J[1][0]*refValues[7] + J[1][1]*refValues[8]) / (*localSubEntityMeshSize_)[8];
124 if constexpr (not reduced)
125 out[9][0] = refValues[9];
126 }
127 }
128
129 using ElementJacobian = Dune::FieldMatrix<DF, dim,dim>;
130
131 using LocalSubEntityMeshSize = std::vector<double>;
132
133public:
134
135 using Domain = Dune::FieldVector<DF, dim>;
136 using Range = Dune::FieldVector<RF, 1>;
137 using Jacobian = Dune::FieldMatrix<RF, 1, dim>;
139 using OrderArray = std::array<unsigned int, dim>;
140
141 static constexpr unsigned int size()
142 {
143 return decltype(referenceBasisCoefficients)::rows;
144 }
145
146 inline void evaluateFunction(const Domain& x, std::vector<Range>& values) const
147 {
148 static constexpr auto monomials = MonomialSet<RF, dim, 3>{};
149 auto monomialValues = monomials(x);
150 auto referenceValues = Dune::FieldVector<RF, size()>{};
151 referenceBasisCoefficients.mv(monomialValues, referenceValues);
152 transformToElementBasis(referenceValues, values);
153 }
154
155 inline void evaluateJacobian(const Domain& x, std::vector<Jacobian>& jacobians) const
156 {
157 static constexpr auto monomials = MonomialSet<RF, dim, 3>{};
158 auto monomialJacobians = derivative(monomials)(x);
159 auto referenceJacobians = Dune::FieldMatrix<RF, size(), dim>{};
160 referenceBasisCoefficients.mv(monomialJacobians, referenceJacobians);
161 transformToElementBasis(referenceJacobians, jacobians);
162 }
163
164 void partial(const OrderArray& order, const Domain& x, std::vector<Range>& out) const
165 {
166 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
167 if (totalOrder == 0)
168 evaluateFunction(x, out);
169 else
170 DUNE_THROW(RangeError, "partial() not implemented for given order");
171 }
172
173 unsigned int order() const
174 {
175 return 3;
176 }
177
178 template<class Element>
179 void bind(const Element& element, const LocalSubEntityMeshSize& localSubEntityMeshSize) {
180 localSubEntityMeshSize_ = &localSubEntityMeshSize;
181 auto center = Dune::ReferenceElements<DF, dim>::simplex().position(0, 0);
182 elementJacobian_ = element.geometry().jacobian(center);
183 }
184
185private:
186 ElementJacobian elementJacobian_;
187 const LocalSubEntityMeshSize* localSubEntityMeshSize_;
188};
189
190
191
192template<unsigned int dim, bool reduced>
193struct CubicHermiteLocalCoefficients
194{
195
196 static constexpr auto makeLocalKeys() {
197 if constexpr (dim==1)
198 return std::array{
199 LocalKey(0, 1, 0),
200 LocalKey(0, 1, 1),
201 LocalKey(1, 1, 0),
202 LocalKey(1, 1, 1)
203 };
204 if constexpr ((dim==2) and (not reduced))
205 return std::array{
206 LocalKey(0, 2, 0),
207 LocalKey(0, 2, 1),
208 LocalKey(0, 2, 2),
209 LocalKey(1, 2, 0),
210 LocalKey(1, 2, 1),
211 LocalKey(1, 2, 2),
212 LocalKey(2, 2, 0),
213 LocalKey(2, 2, 1),
214 LocalKey(2, 2, 2),
215 LocalKey(0, 0, 0)
216 };
217 if constexpr ((dim==2) and (reduced))
218 return std::array{
219 LocalKey(0, 2, 0),
220 LocalKey(0, 2, 1),
221 LocalKey(0, 2, 2),
222 LocalKey(1, 2, 0),
223 LocalKey(1, 2, 1),
224 LocalKey(1, 2, 2),
225 LocalKey(2, 2, 0),
226 LocalKey(2, 2, 1),
227 LocalKey(2, 2, 2)
228 };
229 }
230
231 using LocalKeys = std::decay_t<decltype(makeLocalKeys())>;
232
233public:
234
235 std::size_t size() const
236 {
237 return localKeys_.size();
238 }
239
240 const LocalKey& localKey(std::size_t i) const
241 {
242 assert( i < localKeys_.size() );
243 return localKeys_[i];
244 }
245private:
246 LocalKeys localKeys_ = makeLocalKeys();
247};
248
249
250
251template<class DF, class RF, unsigned int dim, bool reduced>
252class CubicHermiteLocalInterpolation
253{
254 using ElementJacobianInverse = Dune::FieldMatrix<DF, dim,dim>;
255 using LocalSubEntityMeshSize = std::vector<double>;
256
257public:
258
259 template<class Element>
260 void bind(const Element& element, const LocalSubEntityMeshSize& localSubEntityMeshSize) {
261 localSubEntityMeshSize_ = &localSubEntityMeshSize;
262 }
263
264 template<class F, class C>
265 void interpolate(const F& f, std::vector<C>& out) const
266 {
267 using Domain = Dune::FieldVector<DF, dim>;
268 auto&& df = derivative(f);
269 if constexpr (dim==1)
270 {
271 out.resize(4);
272 out[0] = f(0);
273 out[1] = df(0) * (*localSubEntityMeshSize_)[1];
274 out[2] = f(1);
275 out[3] = df(1) * (*localSubEntityMeshSize_)[3];
276 }
277 if constexpr (dim==2)
278 {
279 if constexpr (not reduced)
280 {
281 out.resize(10);
282 out[9] = f(Domain({1.0/3.0,1.0/3.0}));
283 }
284 if constexpr (reduced)
285 out.resize(9);
286 auto J0 = df(Domain({0,0}));
287 auto J1 = df(Domain({1,0}));
288 auto J2 = df(Domain({0,1}));
289 out[0] = f(Domain({0,0}));
290 out[1] = J0[0][0] * (*localSubEntityMeshSize_)[1];
291 out[2] = J0[0][1] * (*localSubEntityMeshSize_)[2];
292 out[3] = f(Domain({1,0}));
293 out[4] = J1[0][0] * (*localSubEntityMeshSize_)[4];
294 out[5] = J1[0][1] * (*localSubEntityMeshSize_)[5];
295 out[6] = f(Domain({0,1}));
296 out[7] = J2[0][0] * (*localSubEntityMeshSize_)[7];
297 out[8] = J2[0][1] * (*localSubEntityMeshSize_)[8];
298 }
299 }
300private:
301 const LocalSubEntityMeshSize* localSubEntityMeshSize_;
302};
303
304
305
313template<class DF, class RF, unsigned int dim, bool reduced>
314class CubicHermiteLocalFiniteElement
315{
316 using LocalBasis = CubicHermiteLocalBasis<DF, RF, dim, reduced>;
317 using LocalCoefficients = CubicHermiteLocalCoefficients<dim, reduced>;
318 using LocalInterpolation = CubicHermiteLocalInterpolation<DF, RF, dim, reduced>;
319 using LocalSubEntityMeshSize = std::vector<double>;
320
321public:
322
323 using Traits = LocalFiniteElementTraits<LocalBasis, LocalCoefficients, LocalInterpolation>;
324
325 const LocalBasis& localBasis() const {
326 return localBasis_;
327 }
328
329 const LocalCoefficients& localCoefficients() const {
330 return localCoefficients_;
331 }
332
333 const LocalInterpolation& localInterpolation() const {
334 return localInterpolation_;
335 }
336
337 unsigned int size() const {
338 return localBasis_.size();
339 }
340
341 GeometryType type() const {
342 return type_;
343 }
344
345 template<class Element, class Mapper, class MeshSizeContainer>
346 void bind(const Element& element, const Mapper& mapper, const MeshSizeContainer& subEntityMeshSize) {
347 // Update local mesh size cache
348 localSubEntityMeshSize_.resize(localCoefficients_.size());
349 for(auto i : Dune::range(size()))
350 {
351 auto localKey = localCoefficients_.localKey(i);
352 auto subEntityIndex = mapper.subIndex(element, localKey.subEntity(), localKey.codim());
353 localSubEntityMeshSize_[i] = subEntityMeshSize[subEntityIndex];
354 }
355
356 localBasis_.bind(element, localSubEntityMeshSize_);
357 localInterpolation_.bind(element, localSubEntityMeshSize_);
358 type_ = element.type();
359 }
360
361private:
362 LocalSubEntityMeshSize localSubEntityMeshSize_;
363 LocalBasis localBasis_;
364 LocalCoefficients localCoefficients_;
365 LocalInterpolation localInterpolation_;
366 GeometryType type_;
367};
368
369
370
371
372} // namespace Impl in Dune::Functions::
373
374
375
376// *****************************************************************************
377// This is the reusable part of the basis. It contains
378//
379// CubicHermitePreBasis
380// CubicHermiteNode
381//
382// The pre-basis allows to create the others and is the owner of possible shared
383// state. These components do _not_ depend on the global basis and local view
384// and can be used without a global basis.
385// *****************************************************************************
386
387template<typename GV, bool reduced>
388class CubicHermiteNode :
389 public LeafBasisNode
390{
391 static const int gridDim = GV::dimension;
392
393 using CubeFiniteElement = Impl::CubicHermiteLocalFiniteElement<typename GV::ctype, double, gridDim, reduced>;
394 using SubEntityMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GV>;
395 using Base = LeafBasisNode;
396
397public:
398
399 using size_type = std::size_t;
400 using Element = typename GV::template Codim<0>::Entity;
401 using FiniteElement = CubeFiniteElement;
402 using Base::size;
403
404 CubicHermiteNode(const SubEntityMapper& subEntityMapper, const std::vector<double>& subEntityMeshSize) :
405 finiteElement_(),
406 element_(nullptr),
407 subEntityMapper_(&subEntityMapper),
408 subEntityMeshSize_(&subEntityMeshSize)
409 {}
410
412 const Element& element() const
413 {
414 return *element_;
415 }
416
421 const FiniteElement& finiteElement() const
422 {
423 return finiteElement_;
424 }
425
427 void bind(const Element& e)
428 {
429 element_ = &e;
430 finiteElement_.bind(*element_, *subEntityMapper_, *subEntityMeshSize_);
431 this->setSize(finiteElement_.size());
432 }
433
434protected:
435
436 FiniteElement finiteElement_;
437 const Element* element_;
438 const SubEntityMapper* subEntityMapper_;
439 const std::vector<double>* subEntityMeshSize_;
440};
441
442
443
451template<typename GV, bool reduced = false>
453{
455 using Base::mapper_;
457
458 static constexpr auto cubicHermiteMapperLayout(Dune::GeometryType type, int gridDim) {
459 if (type.isVertex())
460 return 1 + gridDim;
461 if ((type.isTriangle()) and (not reduced))
462 return 1;
463 else
464 return 0;
465 }
466
467 static constexpr auto subEntityLayout(Dune::GeometryType type, int gridDim) {
468 return (cubicHermiteMapperLayout(type, gridDim) > 0);
469 }
470
471public:
472
473 using GridView = typename Base::GridView;
474 using Node = CubicHermiteNode<GridView, reduced>;
475
476 CubicHermitePreBasis(const GridView& gv) :
477 Base(gv, cubicHermiteMapperLayout),
478 subEntityMapper_(gv, subEntityLayout)
479 {
480 static_assert(GridView::dimension<=2, "CubicHermitePreBasis is only implemented for dim<=2");
481 subEntityMeshSize_ = Impl::computeAverageSubEntityMeshSize(subEntityMapper_);
482 }
483
484 void update(const GridView& gv)
485 {
486 Base::update(gv);
487 subEntityMapper_.update(gv);
488 subEntityMeshSize_ = Impl::computeAverageSubEntityMeshSize(subEntityMapper_);
489 }
490
491 Node makeNode() const
492 {
493 return Node{subEntityMapper_, subEntityMeshSize_};
494 }
495
496private:
497 std::vector<double> subEntityMeshSize_;
498 SubEntityMapper subEntityMapper_;
499};
500
501
502
503namespace BasisFactory {
504
510template<class Dummy=void>
512{
513 return [](const auto& gridView) {
514 return CubicHermitePreBasis<std::decay_t<decltype(gridView)>>(gridView);
515 };
516}
517
523template<class Dummy=void>
525{
526 return [](const auto& gridView) {
527 return CubicHermitePreBasis<std::decay_t<decltype(gridView)>, true>(gridView);
528 };
529}
530
531} // end namespace BasisFactory
532
533
534
535
542template<typename GV>
543using CubicHermiteBasis = DefaultGlobalBasis<CubicHermitePreBasis<GV>>;
544
545
546
553template<typename GV>
555
556
557
558} // end namespace Functions
559} // end namespace Dune
560
561#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_CUBICHERMITEBASIS_HH
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:50
A generic MixIn class for PreBasis with flat indices computed from a mapper.
Definition: leafprebasismappermixin.hh:62
void update(const GridView &gv)
Update the stored GridView.
Definition: leafprebasismappermixin.hh:101
GV GridView
Type of the associated GridView.
Definition: leafprebasismappermixin.hh:68
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
constexpr bool isVertex() const
Return true if entity is a vertex.
Definition: type.hh:279
constexpr bool isTriangle() const
Return true if entity is a triangle.
Definition: type.hh:289
Implementation class for a multiple codim and multiple geometry type mapper.
Definition: mcmgmapper.hh:129
void update(const GV &gridView)
Recalculates indices after grid adaptation.
Definition: mcmgmapper.hh:308
A set of traits classes to store static information about grid implementation.
A few common exception classes.
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:132
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
TrigonometricFunction< K, -cosFactor, sinFactor > derivative(const TrigonometricFunction< K, sinFactor, cosFactor > &f)
Obtain derivative of TrigonometricFunction function.
Definition: trigonometricfunction.hh:43
auto cubicHermite()
Create a pre-basis factory that can create a CubicHermite pre-basis.
Definition: cubichermitebasis.hh:511
DefaultGlobalBasis< CubicHermitePreBasis< GV, R, reduced > > CubicHermiteBasis
Nodal basis of a scalar cubic Hermite finite element space.
Definition: cubichermitebasis.hh:77
auto reducedCubicHermite()
Create a pre-basis factory that can create a CubicHermite pre-basis.
Definition: cubichermitebasis.hh:524
static constexpr int dimension
The dimension of the grid.
Definition: gridview.hh:134
void interpolate(const F &f, const GFS &gfs, XG &xg)
interpolation from a given grid function
Definition: interpolate.hh:177
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
Mapper for multiple codim and multiple geometry types.
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
Convenience header that includes all available Rannacher-Turek LocalFiniteElements.
A pre-basis for a Hermitebasis.
Definition: cubichermitebasis.hh:453
CubicHermiteNode< GridView, R, reduced > Node
Template mapping root tree path to type of created tree node.
Definition: cubichermitebasis.hh:785
Node makeNode() const
Create tree node.
Definition: cubichermitebasis.hh:810
GV GridView
The grid view that the FE basis is defined on.
Definition: cubichermitebasis.hh:779
void update(GridView const &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: cubichermitebasis.hh:800
CubicHermitePreBasis(const GV &gv)
Constructor for a given grid view object.
Definition: cubichermitebasis.hh:790
static const ReferenceElement & simplex()
get simplex reference elements
Definition: referenceelements.hh:162
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:35
Provides the TupleVector class that augments std::tuple by operator[].
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)