DUNE PDELab (git)

morleybasis.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#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_MORLEYBASIS_HH
7#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_MORLEYBASIS_HH
8
9#include <algorithm>
10#include <type_traits>
11#include <vector>
12#include <array>
13#include <bitset>
14
17
19#include <dune/grid/common/rangegenerators.hh>
20
21#include <dune/localfunctions/common/localbasis.hh>
22#include <dune/localfunctions/common/localfiniteelementtraits.hh>
23#include <dune/localfunctions/common/localkey.hh>
24
25#include <dune/functions/analyticfunctions/monomialset.hh>
27#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
28#include <dune/functions/functionspacebases/functionaldescriptor.hh>
29#include <dune/functions/functionspacebases/leafprebasismappermixin.hh>
30#include <dune/functions/functionspacebases/nodes.hh>
31#include <dune/functions/functionspacebases/transformedfiniteelementmixin.hh>
32
49namespace Dune::Functions
50{
51
52 template<class GV, class R>
53 struct MorleyPreBasis;
54
72 template<class GV, class R = double>
74
75
76 namespace Impl
77 {
78
85 template<class D, class R>
86 class MorleyReferenceLocalBasis
87 {
88 public:
89 static constexpr int dim = 2;
90 using Traits = H2LocalBasisTraits<D, dim, FieldVector<D, dim>, R, 1, FieldVector<R, 1>,
92
93 private:
94
107 static constexpr auto getMorleyCoefficients()
108 {
109 // Define std::sqrt(2.) manually in double precision,
110 // because std::sqrt is not constexpr before C++26.
111 D sqrt2 = 0.5 * 1.414213562373095;
112
113 return Dune::FieldMatrix<D, 6,6>({{1, -1, -1, 0, 2, 0},
114 {0, 0.5, 0.5, 0.5, -1, -0.5},
115 {0, 0.5, 0.5, -0.5, -1, 0.5},
116 {0, 0, 1, 0, 0, -1},
117 {0, -1., 0, 1., 0, 0},
118 {0, sqrt2, sqrt2, -sqrt2, -2. * sqrt2, -sqrt2}});
119 }
120
121 static constexpr auto referenceBasisCoefficients = getMorleyCoefficients();
123
124 public:
125
128 static constexpr unsigned int size()
129 {
130 return 6;
131 }
132
135 static constexpr unsigned int order()
136 {
137 return 2;
138 }
139
145 void evaluateFunction(const typename Traits::DomainType &in,
146 std::vector<typename Traits::RangeType> &out) const
147 {
148 out.resize(size());
149 auto monomialValues = monomials(in);
150 multiplyWithCoefficentMatrix(referenceBasisCoefficients, monomialValues, out);
151 }
152
158 void evaluateJacobian(const typename Traits::DomainType &in,
159 std::vector<typename Traits::JacobianType> &out) const
160 {
161 out.resize(size());
162 auto monomialValues = derivative(monomials)(in);
163 multiplyWithCoefficentMatrix(referenceBasisCoefficients, monomialValues, out);
164 }
165
171 void evaluateHessian(const typename Traits::DomainType &in,
172 std::vector<typename Traits::HessianType> &out) const
173 {
174 out.resize(size());
175 auto monomialValues = derivative(derivative(monomials))(in);
176 multiplyWithCoefficentMatrix(referenceBasisCoefficients, monomialValues, out);
177 }
178
185 void partial(std::array<unsigned int, dim> order, const typename Traits::DomainType &in,
186 std::vector<typename Traits::RangeType> &out) const
187 {
188 out.resize(size());
189 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
190 if (totalOrder == 0)
191 evaluateFunction(in, out);
192 else if (totalOrder == 1)
193 {
194 evaluateJacobian(in,jacobiansBuffer_);
195 std::size_t which = std::max_element(order.begin(), order.end()) - order.begin();
196 for (auto i : Dune::range(size()))
197 out[i] = jacobiansBuffer_[i][0][which];
198 }
199 else if (totalOrder == 2)
200 {
201 evaluateHessian(in, hessianBuffer_);
202 std::size_t first, second;
203 first = std::max_element(order.begin(), order.end()) - order.begin();
204 if (order[first] == 2)
205 second = first;
206 else
207 {
208 order[first] = 0;
209 second = std::max_element(order.begin(), order.end()) - order.begin();
210 }
211 for (auto i : Dune::range(size()))
212 out[i] = hessianBuffer_[i][first][second];
213 }
214 else
215 DUNE_THROW(RangeError, "partial() not implemented for given order");
216 }
217
218 private:
219 mutable std::vector<typename Traits::JacobianType> jacobiansBuffer_;
220 mutable std::vector<typename Traits::HessianType> hessianBuffer_;
221 };
222
223
227 class MorleyLocalCoefficients
228 {
229 public:
230 using size_type = std::size_t;
231
232 MorleyLocalCoefficients()
233 {
234 for (size_type i = 0; i < 3; ++i)
235 {
236 localKeys_[i] = LocalKey(i, 2, 0); // vertex dofs
237 localKeys_[3 + i] = LocalKey(i, 1, 0); // edge dofs
238 }
239 }
240
243 static constexpr size_type size()
244 {
245 return 6;
246 }
247
250 constexpr LocalKey const &localKey(size_type i) const
251 {
252 return localKeys_[i];
253 }
254
255 private:
256 std::array<LocalKey, 6> localKeys_;
257 };
258
259
265 template<class D>
266 class MorleyLocalInterpolation
267 {
268 using size_type = std::size_t;
269 using FunctionalDescriptor = Dune::Functions::Impl::FunctionalDescriptor<2>;
270
271 public:
272
273 MorleyLocalInterpolation()
274 {
275 descriptors_[0] = FunctionalDescriptor();
276 descriptors_[1] = FunctionalDescriptor();
277 descriptors_[2] = FunctionalDescriptor();
278 descriptors_[3] = FunctionalDescriptor(1);
279 descriptors_[4] = FunctionalDescriptor(1);
280 descriptors_[5] = FunctionalDescriptor(1);
281 }
282
285 template <class Element>
286 void bind(const Element &element, const std::bitset<3> &edgeOrientation)
287 {
288 auto geometry = element.geometry();
289
290 // get global Normals and midpoints
291 auto refElement = Dune::referenceElement<double, 2>(geometry.type());
292 for (std::size_t i = 0; i < 3; ++i)
293 {
294 localVertices_[i] = refElement.position(i, 2);
295
296 localMidpoints_[i] = refElement.position(i, 1);
297 std::size_t lower = (i == 2) ? 1 : 0;
298 std::size_t upper = (i == 0) ? 1 : 2;
299
300 auto edge = geometry.global(refElement.position(upper, 2))
301 - geometry.global(refElement.position(lower, 2));
302 // normalize and orient
303 edge /= edge.two_norm() * (edgeOrientation[i] ? -1. : 1.);
304 // Rotation by pi/2. Note that Kirby rotates by 3*pi/2
305 globalNormals_[i] = {-edge[1], edge[0]};
306 }
307 }
308
316 template <class F, class C>
317 void interpolate(const F &f, std::vector<C> &out) const
318 {
319 out.resize(6);
320 auto&& df = derivative(f);
321 for (size_type i = 0; i < 3; ++i)
322 {
323 out[i] = f(localVertices_[i]);
324 out[3 + i] = squeeze(df(localMidpoints_[i])).dot(globalNormals_[i]);
325 }
326 }
327
328 const FunctionalDescriptor& functionalDescriptor(size_type i) const
329 {
330 return descriptors_[i];
331 }
332
333 protected:
334 std::array<Dune::FieldVector<D, 2>, 3> globalNormals_;
335 std::array<Dune::FieldVector<D, 2>, 3> localMidpoints_;
336 std::array<Dune::FieldVector<D, 2>, 3> localVertices_;
337 std::array<FunctionalDescriptor, 6> descriptors_;
338
339 };
340
341 template<class D, class R>
342 using MorleyLocalBasisTraits = typename Impl::MorleyReferenceLocalBasis<D, R>::Traits;
343
352 template<class D, class R>
353 class MorleyLocalFiniteElement
354 : public Impl::TransformedFiniteElementMixin<MorleyLocalFiniteElement<D,R>, MorleyLocalBasisTraits<D, R>>
355 {
356 using Base = Impl::TransformedFiniteElementMixin< MorleyLocalFiniteElement<D,R>, MorleyLocalBasisTraits<D, R>>;
357 friend class Impl::TransformedLocalBasis<MorleyLocalFiniteElement<D,R>, MorleyLocalBasisTraits<D, R>>;
358 static constexpr int dim = 2;
359 public:
360
363 using size_type = std::size_t;
364 using Traits = LocalFiniteElementTraits<
365 Impl::TransformedLocalBasis<MorleyLocalFiniteElement<D,R>, MorleyLocalBasisTraits<D, R>>,
366 Impl::MorleyLocalCoefficients,
367 Impl::MorleyLocalInterpolation<D>>;
368
372 const typename Traits::LocalCoefficientsType &localCoefficients() const
373 {
374 return coefficients_;
375 }
376
379 const typename Traits::LocalInterpolationType &localInterpolation() const
380 {
381 return interpolation_;
382 }
383
386 static constexpr GeometryType type()
387 {
388 return GeometryTypes::simplex(dim);
389 }
390
393 static constexpr size_type size()
394 {
395 return 6;
396 }
397
400 template<class Element>
401 void bind(std::bitset<3> const& data, Element const &e)
402 {
403 edgeOrientation_ = data;
404
405 fillMatrix(e.geometry());
406 interpolation_.bind(e, edgeOrientation_);
407 }
408
409 protected:
410
413 Impl::MorleyReferenceLocalBasis<D, R> const& referenceLocalBasis() const
414 {
415 return basis_;
416 }
417
423 template<class InputValues, class OutputValues>
424 void transform(InputValues const& inValues, OutputValues& outValues) const
425 {
426 mat_.mv(inValues, outValues);
427 }
428
429 private:
430
437 template<class Geometry>
438 void fillMatrix(Geometry const &geometry)
439 {
440 std::array<R, 3> B_11;
441 std::array<R, 3> B_12;
442 std::array<R, 3> l_inv;
443
444 std::array<Dune::FieldVector<R, 2>, 3> referenceTangents;
445 std::array<Dune::FieldVector<R, 2>, 3> globalTangents;
446
447 // By default, edges point from the vertex with the smaller index
448 // to the vertex with the larger index.
449
450 // get local and global Tangents
451 auto refElement = Dune::referenceElement<double, 2>(geometry.type());
452 auto x = refElement.position(0,0);
453 for (std::size_t i = 0; i < 3; ++i)
454 {
455 std::size_t lower = (i == 2) ? 1 : 0;
456 std::size_t upper = (i == 0) ? 1 : 2;
457 auto edge = refElement.position(upper, 2) - refElement.position(lower, 2);
458
459 referenceTangents[i] = edge / edge.two_norm();
460
461 auto globalEdge = geometry.global(refElement.position(upper, 2))
462 - geometry.global(refElement.position(lower, 2));
463
464 l_inv[i] = 1. / globalEdge.two_norm();
465 globalTangents[i] = globalEdge * l_inv[i];
466 }
467
468 auto jacobianTransposed = geometry.jacobianTransposed(x);
469 for (std::size_t i = 0; i < 3; ++i)
470 {
471 B_11[i] = -referenceTangents[i][1]
472 *(-globalTangents[i][1] * jacobianTransposed[0][0]
473 + globalTangents[i][0] * jacobianTransposed[0][1])
474 + referenceTangents[i][0]
475 *(-globalTangents[i][1] * jacobianTransposed[1][0]
476 + globalTangents[i][0] * jacobianTransposed[1][1]);
477 B_12[i] = -referenceTangents[i][1]
478 *(globalTangents[i][0] * jacobianTransposed[0][0]
479 + globalTangents[i][1] * jacobianTransposed[0][1])
480 + referenceTangents[i][0]
481 *(globalTangents[i][0] * jacobianTransposed[1][0]
482 + globalTangents[i][1] * jacobianTransposed[1][1]);
483 }
484
485 // Actually setup matrix
486 int sign = -1;
487 for (std::size_t i = 0; i < 3; ++i)
488 {
489 mat_[i][i] = 1.;
490 for (std::size_t j = 0; j < 3; ++j)
491 {
492 if (j != (2 - i)) // dune specific edge order
493 {
494 mat_[j][3 + i] = sign * B_12[i] * l_inv[i];
495 sign *= -1;
496 }
497 }
498 mat_[3 + i][3 + i] = (edgeOrientation_[i] ? -1. : 1.) * B_11[i];
499 }
500 }
501
502 // a finite element consists of a basis, coeffiecents and an interpolation
503 typename Impl::MorleyReferenceLocalBasis<D, R> basis_;
504 typename Traits::LocalCoefficientsType coefficients_;
505 typename Traits::LocalInterpolationType interpolation_;
506 // This is the matrix M in Kirbys paper
508 // the local state, i.e. a collection of global information restricted to this element
509 std::bitset<3> edgeOrientation_;
510
511 };
512
513 } // end namespace Impl
514
515
516
517 // *****************************************************************************
518 // This is the reusable part of the basis. It contains
519 //
520 // MorleyPreBasis
521 // MorleyNode
522 //
523 // The pre-basis allows to create the others and is the owner of possible shared
524 // state. These components do _not_ depend on the global basis and local view
525 // and can be used without a global basis.
526 // *****************************************************************************
527
528 template<class GV, class R>
529 class MorleyNode
530 : public LeafBasisNode
531 {
533 public:
534 using size_type = std::size_t;
535 using Element = typename GV::template Codim<0>::Entity;
536 using FiniteElement = typename Impl::MorleyLocalFiniteElement<typename GV::ctype, R>;
537
538
539 MorleyNode(Mapper const& m, std::vector<std::bitset<3>> const& data)
540 : mapper_(&m), data_(&data)
541 {
542 this->setSize(finiteElement_.size());
543 }
544
546 Element const &element() const
547 {
548 return *element_;
549 }
550
556 FiniteElement const &finiteElement() const
557 {
558 return finiteElement_;
559 }
560
562 void bind(Element const &e)
563 {
564 element_ = &e;
565 finiteElement_.bind((*data_)[mapper_->index(e)], *element_);
566 }
567
569 unsigned int order() const
570 {
571 return finiteElement_.localBasis().order();
572 }
573
574 protected:
575 FiniteElement finiteElement_;
576 Element const* element_;
577 Mapper const* mapper_;
578 std::vector<std::bitset<3>> const* data_;
579 };
580
581
591 template<class GV, class R>
593 : public LeafPreBasisMapperMixin<GV>
594 {
597 using Element = typename GV::template Codim<0>::Entity;
598 using D = typename GV::ctype;
599 static const std::size_t dim = GV::dimension;
600
601 // helper methods to assign each subentity the number of dofs. Used by the LeafPreBasisMapperMixin.
602 static constexpr auto morleyMapperLayout(Dune::GeometryType type, int gridDim)
603 {
604 assert(gridDim == 2);
605 if (type.isVertex())
606 return 1; // one evaluation dof per vertex
607 if (type.isLine())
608 return 1;
609 if ((type.isTriangle()) )
610 return 0;
611 else
612 return 0;
613 }
614
615 public:
617 using GridView = GV;
618
620 using size_type = std::size_t;
621
623 using Node = MorleyNode<GridView, R>;
624
626 MorleyPreBasis(const GV &gv)
627 : Base(gv, morleyMapperLayout)
628 , mapper_({gv, mcmgElementLayout()})
629 {
630 updateState(gv);
631 }
632
634 void update(GridView const &gv)
635 {
636 Base::update(gv);
637 updateState(gv);
638 }
639
644 {
645 return Node{mapper_, data_};
646 }
647
648 protected:
649
650 void updateState(GridView const &gridView)
651 {
652 data_.resize(gridView.size(0));
653 // compute orientation for all elements
654 unsigned short orientation = 0;
655 auto const& idSet = gridView.grid().globalIdSet();
656
657 for (const auto &element : elements(gridView))
658 {
659 const auto &refElement = referenceElement(element);
660 auto elementIndex = mapper_.index(element);
661
662 orientation = 0;
663
664 for (std::size_t i = 0; i < element.subEntities(dim - 1); i++)
665 {
666 // Local vertex indices within the element
667 auto localV0 = refElement.subEntity(i, dim - 1, 0, dim);
668 auto localV1 = refElement.subEntity(i, dim - 1, 1, dim);
669
670 // Global vertex indices within the grid
671 auto globalV0 = idSet.subId(element, localV0, dim);
672 auto globalV1 = idSet.subId(element, localV1, dim);
673
674 // The edge is flipped if the local ordering disagrees with global ordering
675 if ((localV0 < localV1 && globalV0 > globalV1)
676 || (localV0 > localV1 && globalV0 < globalV1))
677 {
678 orientation |= (1 << i);
679 }
680 }
681 data_[elementIndex] = orientation;
682 }
683 }
684
685 SubEntityMapper mapper_;
686 std::vector<std::bitset<3>> data_;
687
688 }; // class MorleyPreBasis
689
690 namespace BasisFactory
691 {
692
699 template<class R = double>
700 auto morley()
701 {
702 return [=](auto const &gridView) {
703 return MorleyPreBasis<std::decay_t<decltype(gridView)>, R>(gridView);
704 };
705 }
706
707 } // end namespace BasisFactory
708
709} // end namespace Dune::Functions
710
711#endif
FieldTraits< value_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition: densevector.hh:641
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
const GridView & gridView() const
Export the stored GridView.
Definition: leafprebasismappermixin.hh:95
std::size_t size_type
Type used for index digits.
Definition: leafprebasismappermixin.hh:71
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
constexpr bool isLine() const
Return true if entity is a line segment.
Definition: type.hh:284
Wrapper class for geometries.
Definition: geometry.hh:71
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
Return the transposed of the Jacobian.
Definition: geometry.hh:302
GeometryType type() const
Return the type of the reference element. The type can be used to access the Dune::ReferenceElement.
Definition: geometry.hh:194
GlobalCoordinate global(const LocalCoordinate &local) const
Evaluate the map .
Definition: geometry.hh:228
Grid view abstract base class.
Definition: gridview.hh:66
Describe position of one degree of freedom.
Definition: localkey.hh:24
Mapper interface.
Definition: mapper.hh:110
Implementation class for a multiple codim and multiple geometry type mapper.
Definition: mcmgmapper.hh:129
Index index(const EntityType &e) const
Map entity to starting index in array for dof block.
Definition: mcmgmapper.hh:171
Default exception class for range errors.
Definition: exceptions.hh:254
This file provides an implementation of the cubic Hermite finite element in 1 to 3 dimensions.
A few common exception classes.
Implements a vector constructed from a given type representing a field and a compile-time given size.
#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
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
constexpr GeometryType simplex(unsigned int dim)
Returns a GeometryType representing a simplex of dimension dim.
Definition: type.hh:453
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
MCMGLayout mcmgElementLayout()
layout for elements (codim-0 entities)
Definition: mcmgmapper.hh:97
Mapper for multiple codim and multiple geometry types.
auto morley()
construct a PreBasisFactory for the Morley Finite Element
Definition: morleybasis.hh:700
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
int sign(const T &val)
Return the sign of the value.
Definition: math.hh:180
Static tag representing a codimension.
Definition: dimension.hh:24
Function, which evaluates all monomials up to degree maxDegree in a given coordinate.
Definition: monomialset.hh:64
A pre-basis for a Morleybasis.
Definition: morleybasis.hh:594
void update(GridView const &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: morleybasis.hh:634
Node makeNode() const
Create tree node.
Definition: morleybasis.hh:643
MorleyNode< GridView, R > Node
Template mapping root tree path to type of created tree node.
Definition: morleybasis.hh:623
MorleyPreBasis(const GV &gv)
Constructor for a given grid view object.
Definition: morleybasis.hh:626
traits helper struct
Definition: localfiniteelementtraits.hh:13
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:20
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:24
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)