6#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_MORLEYBASIS_HH
7#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_MORLEYBASIS_HH
19#include <dune/grid/common/rangegenerators.hh>
21#include <dune/localfunctions/common/localbasis.hh>
22#include <dune/localfunctions/common/localfiniteelementtraits.hh>
23#include <dune/localfunctions/common/localkey.hh>
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>
49namespace Dune::Functions
52 template<
class GV,
class R>
53 struct MorleyPreBasis;
72 template<
class GV,
class R =
double>
85 template<
class D,
class R>
86 class MorleyReferenceLocalBasis
89 static constexpr int dim = 2;
90 using Traits = H2LocalBasisTraits<D, dim, FieldVector<D, dim>, R, 1,
FieldVector<R, 1>,
107 static constexpr auto getMorleyCoefficients()
111 D sqrt2 = 0.5 * 1.414213562373095;
114 {0, 0.5, 0.5, 0.5, -1, -0.5},
115 {0, 0.5, 0.5, -0.5, -1, 0.5},
117 {0, -1., 0, 1., 0, 0},
118 {0, sqrt2, sqrt2, -sqrt2, -2. * sqrt2, -sqrt2}});
121 static constexpr auto referenceBasisCoefficients = getMorleyCoefficients();
128 static constexpr unsigned int size()
135 static constexpr unsigned int order()
145 void evaluateFunction(
const typename Traits::DomainType &in,
146 std::vector<typename Traits::RangeType> &out)
const
149 auto monomialValues = monomials(in);
150 multiplyWithCoefficentMatrix(referenceBasisCoefficients, monomialValues, out);
158 void evaluateJacobian(
const typename Traits::DomainType &in,
159 std::vector<typename Traits::JacobianType> &out)
const
162 auto monomialValues =
derivative(monomials)(in);
163 multiplyWithCoefficentMatrix(referenceBasisCoefficients, monomialValues, out);
171 void evaluateHessian(
const typename Traits::DomainType &in,
172 std::vector<typename Traits::HessianType> &out)
const
176 multiplyWithCoefficentMatrix(referenceBasisCoefficients, monomialValues, out);
185 void partial(std::array<unsigned int, dim> order,
const typename Traits::DomainType &in,
186 std::vector<typename Traits::RangeType> &out)
const
191 evaluateFunction(in, out);
192 else if (totalOrder == 1)
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];
199 else if (totalOrder == 2)
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)
209 second = std::max_element(order.begin(), order.end()) - order.begin();
211 for (
auto i : Dune::range(
size()))
212 out[i] = hessianBuffer_[i][first][second];
219 mutable std::vector<typename Traits::JacobianType> jacobiansBuffer_;
220 mutable std::vector<typename Traits::HessianType> hessianBuffer_;
227 class MorleyLocalCoefficients
230 using size_type = std::size_t;
232 MorleyLocalCoefficients()
234 for (size_type i = 0; i < 3; ++i)
237 localKeys_[3 + i] =
LocalKey(i, 1, 0);
243 static constexpr size_type
size()
250 constexpr LocalKey const &localKey(size_type i)
const
252 return localKeys_[i];
256 std::array<LocalKey, 6> localKeys_;
266 class MorleyLocalInterpolation
268 using size_type = std::size_t;
269 using FunctionalDescriptor = Dune::Functions::Impl::FunctionalDescriptor<2>;
273 MorleyLocalInterpolation()
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);
285 template <
class Element>
286 void bind(
const Element &element,
const std::bitset<3> &edgeOrientation)
288 auto geometry = element.geometry();
291 auto refElement = Dune::referenceElement<double, 2>(geometry.type());
292 for (std::size_t i = 0; i < 3; ++i)
294 localVertices_[i] = refElement.position(i, 2);
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;
300 auto edge = geometry.global(refElement.position(upper, 2))
301 - geometry.global(refElement.position(lower, 2));
303 edge /= edge.two_norm() * (edgeOrientation[i] ? -1. : 1.);
305 globalNormals_[i] = {-edge[1], edge[0]};
316 template <
class F,
class C>
317 void interpolate(
const F &f, std::vector<C> &out)
const
321 for (size_type i = 0; i < 3; ++i)
323 out[i] = f(localVertices_[i]);
324 out[3 + i] = squeeze(df(localMidpoints_[i])).dot(globalNormals_[i]);
328 const FunctionalDescriptor& functionalDescriptor(size_type i)
const
330 return descriptors_[i];
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_;
341 template<
class D,
class R>
342 using MorleyLocalBasisTraits =
typename Impl::MorleyReferenceLocalBasis<D, R>::Traits;
352 template<
class D,
class R>
353 class MorleyLocalFiniteElement
354 :
public Impl::TransformedFiniteElementMixin<MorleyLocalFiniteElement<D,R>, MorleyLocalBasisTraits<D, R>>
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;
363 using size_type = std::size_t;
365 Impl::TransformedLocalBasis<MorleyLocalFiniteElement<D,R>, MorleyLocalBasisTraits<D, R>>,
366 Impl::MorleyLocalCoefficients,
367 Impl::MorleyLocalInterpolation<D>>;
374 return coefficients_;
381 return interpolation_;
393 static constexpr size_type
size()
400 template<
class Element>
401 void bind(std::bitset<3>
const& data, Element
const &e)
403 edgeOrientation_ = data;
405 fillMatrix(e.geometry());
406 interpolation_.bind(e, edgeOrientation_);
413 Impl::MorleyReferenceLocalBasis<D, R>
const& referenceLocalBasis()
const
423 template<
class InputValues,
class OutputValues>
424 void transform(InputValues
const& inValues, OutputValues& outValues)
const
426 mat_.mv(inValues, outValues);
437 template<
class Geometry>
438 void fillMatrix(
Geometry const &geometry)
440 std::array<R, 3> B_11;
441 std::array<R, 3> B_12;
442 std::array<R, 3> l_inv;
444 std::array<Dune::FieldVector<R, 2>, 3> referenceTangents;
445 std::array<Dune::FieldVector<R, 2>, 3> globalTangents;
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)
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);
459 referenceTangents[i] = edge / edge.two_norm();
461 auto globalEdge = geometry.
global(refElement.position(upper, 2))
462 - geometry.
global(refElement.position(lower, 2));
464 l_inv[i] = 1. / globalEdge.
two_norm();
465 globalTangents[i] = globalEdge * l_inv[i];
469 for (std::size_t i = 0; i < 3; ++i)
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]);
487 for (std::size_t i = 0; i < 3; ++i)
490 for (std::size_t j = 0; j < 3; ++j)
494 mat_[j][3 + i] =
sign * B_12[i] * l_inv[i];
498 mat_[3 + i][3 + i] = (edgeOrientation_[i] ? -1. : 1.) * B_11[i];
503 typename Impl::MorleyReferenceLocalBasis<D, R> basis_;
509 std::bitset<3> edgeOrientation_;
528 template<
class GV,
class R>
530 :
public LeafBasisNode
534 using size_type = std::size_t;
536 using FiniteElement =
typename Impl::MorleyLocalFiniteElement<typename GV::ctype, R>;
539 MorleyNode(
Mapper const& m, std::vector<std::bitset<3>>
const& data)
540 : mapper_(&m), data_(&data)
542 this->setSize(finiteElement_.size());
546 Element
const &element()
const
556 FiniteElement
const &finiteElement()
const
558 return finiteElement_;
562 void bind(Element
const &e)
565 finiteElement_.bind((*data_)[mapper_->index(e)], *element_);
569 unsigned int order()
const
571 return finiteElement_.localBasis().order();
575 FiniteElement finiteElement_;
576 Element
const* element_;
577 Mapper
const* mapper_;
578 std::vector<std::bitset<3>>
const* data_;
591 template<
class GV,
class R>
598 using D =
typename GV::ctype;
599 static const std::size_t dim = GV::dimension;
604 assert(gridDim == 2);
623 using Node = MorleyNode<GridView, R>;
627 :
Base(gv, morleyMapperLayout)
645 return Node{mapper_, data_};
654 unsigned short orientation = 0;
655 auto const& idSet =
gridView.grid().globalIdSet();
657 for (
const auto &element : elements(
gridView))
660 auto elementIndex = mapper_.
index(element);
664 for (std::size_t i = 0; i < element.subEntities(dim - 1); i++)
667 auto localV0 = refElement.subEntity(i, dim - 1, 0, dim);
668 auto localV1 = refElement.subEntity(i, dim - 1, 1, dim);
671 auto globalV0 = idSet.subId(element, localV0, dim);
672 auto globalV1 = idSet.subId(element, localV1, dim);
675 if ((localV0 < localV1 && globalV0 > globalV1)
676 || (localV0 > localV1 && globalV0 < globalV1))
678 orientation |= (1 << i);
681 data_[elementIndex] = orientation;
685 SubEntityMapper mapper_;
686 std::vector<std::bitset<3>> data_;
690 namespace BasisFactory
699 template<
class R =
double>
702 return [=](
auto const &gridView) {
703 return MorleyPreBasis<std::decay_t<
decltype(gridView)>, R>(gridView);
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