7#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_CUBICHERMITEBASIS_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_CUBICHERMITEBASIS_HH
15#include <dune/grid/common/rangegenerators.hh>
17#include <dune/localfunctions/common/localfiniteelementvariant.hh>
19#include <dune/localfunctions/crouzeixraviart.hh>
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>
38template<
class DF,
class RF,
unsigned int dim,
bool reduced>
39class CubicHermiteLocalBasis
42 static constexpr auto makeReferenceBasisCoefficients() {
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}
63 if constexpr ((dim==2) and (reduced))
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},
91 static constexpr auto referenceBasisCoefficients = makeReferenceBasisCoefficients();
99 template<
class LambdaRefValues,
class Entry>
100 void transformToElementBasis(
const LambdaRefValues& refValues, std::vector<Entry>& out)
const
102 if constexpr (dim==1)
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];;
111 if constexpr (dim==2)
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];
131 using LocalSubEntityMeshSize = std::vector<double>;
139 using OrderArray = std::array<unsigned int, dim>;
141 static constexpr unsigned int size()
143 return decltype(referenceBasisCoefficients)::rows;
146 inline void evaluateFunction(
const Domain& x, std::vector<Range>& values)
const
148 static constexpr auto monomials = MonomialSet<RF, dim, 3>{};
149 auto monomialValues = monomials(x);
151 referenceBasisCoefficients.mv(monomialValues, referenceValues);
152 transformToElementBasis(referenceValues, values);
155 inline void evaluateJacobian(
const Domain& x, std::vector<Jacobian>& jacobians)
const
157 static constexpr auto monomials = MonomialSet<RF, dim, 3>{};
158 auto monomialJacobians =
derivative(monomials)(x);
160 referenceBasisCoefficients.mv(monomialJacobians, referenceJacobians);
161 transformToElementBasis(referenceJacobians, jacobians);
164 void partial(
const OrderArray& order,
const Domain& x, std::vector<Range>& out)
const
168 evaluateFunction(x, out);
170 DUNE_THROW(RangeError,
"partial() not implemented for given order");
173 unsigned int order()
const
178 template<
class Element>
179 void bind(
const Element& element,
const LocalSubEntityMeshSize& localSubEntityMeshSize) {
180 localSubEntityMeshSize_ = &localSubEntityMeshSize;
182 elementJacobian_ = element.geometry().jacobian(center);
186 ElementJacobian elementJacobian_;
187 const LocalSubEntityMeshSize* localSubEntityMeshSize_;
192template<
unsigned int dim,
bool reduced>
193struct CubicHermiteLocalCoefficients
196 static constexpr auto makeLocalKeys() {
197 if constexpr (dim==1)
204 if constexpr ((dim==2) and (not reduced))
217 if constexpr ((dim==2) and (reduced))
231 using LocalKeys = std::decay_t<
decltype(makeLocalKeys())>;
235 std::size_t
size()
const
237 return localKeys_.size();
240 const LocalKey& localKey(std::size_t i)
const
242 assert( i < localKeys_.size() );
243 return localKeys_[i];
246 LocalKeys localKeys_ = makeLocalKeys();
251template<
class DF,
class RF,
unsigned int dim,
bool reduced>
252class CubicHermiteLocalInterpolation
255 using LocalSubEntityMeshSize = std::vector<double>;
259 template<
class Element>
260 void bind(
const Element& element,
const LocalSubEntityMeshSize& localSubEntityMeshSize) {
261 localSubEntityMeshSize_ = &localSubEntityMeshSize;
264 template<
class F,
class C>
265 void interpolate(
const F& f, std::vector<C>& out)
const
269 if constexpr (dim==1)
273 out[1] = df(0) * (*localSubEntityMeshSize_)[1];
275 out[3] = df(1) * (*localSubEntityMeshSize_)[3];
277 if constexpr (dim==2)
279 if constexpr (not reduced)
282 out[9] = f(Domain({1.0/3.0,1.0/3.0}));
284 if constexpr (reduced)
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];
301 const LocalSubEntityMeshSize* localSubEntityMeshSize_;
313template<
class DF,
class RF,
unsigned int dim,
bool reduced>
314class CubicHermiteLocalFiniteElement
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>;
323 using Traits = LocalFiniteElementTraits<LocalBasis, LocalCoefficients, LocalInterpolation>;
325 const LocalBasis& localBasis()
const {
329 const LocalCoefficients& localCoefficients()
const {
330 return localCoefficients_;
333 const LocalInterpolation& localInterpolation()
const {
334 return localInterpolation_;
337 unsigned int size()
const {
338 return localBasis_.size();
345 template<
class Element,
class Mapper,
class MeshSizeContainer>
346 void bind(
const Element& element,
const Mapper& mapper,
const MeshSizeContainer& subEntityMeshSize) {
348 localSubEntityMeshSize_.resize(localCoefficients_.size());
349 for(
auto i : Dune::range(
size()))
351 auto localKey = localCoefficients_.localKey(i);
352 auto subEntityIndex = mapper.subIndex(element, localKey.subEntity(), localKey.codim());
353 localSubEntityMeshSize_[i] = subEntityMeshSize[subEntityIndex];
356 localBasis_.bind(element, localSubEntityMeshSize_);
357 localInterpolation_.bind(element, localSubEntityMeshSize_);
358 type_ = element.type();
362 LocalSubEntityMeshSize localSubEntityMeshSize_;
363 LocalBasis localBasis_;
364 LocalCoefficients localCoefficients_;
365 LocalInterpolation localInterpolation_;
387template<
typename GV,
bool reduced>
388class CubicHermiteNode :
391 static const int gridDim = GV::dimension;
393 using CubeFiniteElement = Impl::CubicHermiteLocalFiniteElement<typename GV::ctype, double, gridDim, reduced>;
395 using Base = LeafBasisNode;
399 using size_type = std::size_t;
400 using Element =
typename GV::template Codim<0>::Entity;
401 using FiniteElement = CubeFiniteElement;
404 CubicHermiteNode(
const SubEntityMapper& subEntityMapper,
const std::vector<double>& subEntityMeshSize) :
407 subEntityMapper_(&subEntityMapper),
408 subEntityMeshSize_(&subEntityMeshSize)
412 const Element& element()
const
421 const FiniteElement& finiteElement()
const
423 return finiteElement_;
427 void bind(
const Element& e)
430 finiteElement_.bind(*element_, *subEntityMapper_, *subEntityMeshSize_);
431 this->setSize(finiteElement_.size());
436 FiniteElement finiteElement_;
437 const Element* element_;
438 const SubEntityMapper* subEntityMapper_;
439 const std::vector<double>* subEntityMeshSize_;
451template<
typename GV,
bool reduced = false>
458 static constexpr auto cubicHermiteMapperLayout(
Dune::GeometryType type,
int gridDim) {
468 return (cubicHermiteMapperLayout(type, gridDim) > 0);
474 using Node = CubicHermiteNode<GridView, reduced>;
477 Base(gv, cubicHermiteMapperLayout),
478 subEntityMapper_(gv, subEntityLayout)
480 static_assert(
GridView::dimension<=2,
"CubicHermitePreBasis is only implemented for dim<=2");
481 subEntityMeshSize_ = Impl::computeAverageSubEntityMeshSize(subEntityMapper_);
487 subEntityMapper_.
update(gv);
488 subEntityMeshSize_ = Impl::computeAverageSubEntityMeshSize(subEntityMapper_);
493 return Node{subEntityMapper_, subEntityMeshSize_};
497 std::vector<double> subEntityMeshSize_;
503namespace BasisFactory {
510template<
class Dummy=
void>
513 return [](
const auto& gridView) {
523template<
class Dummy=
void>
526 return [](
const auto& gridView) {
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[].