7#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_CUBICHERMITEBASIS_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_CUBICHERMITEBASIS_HH
10#include <dune/common/exceptions.hh>
11#include <dune/common/tuplevector.hh>
13#include <dune/grid/common/capabilities.hh>
14#include <dune/grid/common/mcmgmapper.hh>
15#include <dune/grid/common/rangegenerators.hh>
17#include <dune/localfunctions/common/localfiniteelementvariant.hh>
18#include <dune/localfunctions/rannacherturek.hh>
19#include <dune/localfunctions/crouzeixraviart.hh>
21#include <dune/functions/functionspacebases/nodes.hh>
22#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
23#include <dune/functions/functionspacebases/leafprebasismappermixin.hh>
36 template<
class Gr
idView>
37 auto subIndexSet(
const Dune::MultipleCodimMultipleGeomTypeMapper<GridView>& mapper,
const typename GridView::template Codim<0>::Entity& element)
39 using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
40 using Index =
typename Mapper::Index;
41 constexpr auto dimension = GridView::dimension;
42 auto subIndices = std::vector<Index>();
43 auto referenceElement = Dune::referenceElement<double, dimension>(element.type());
44 for(
auto codim : Dune::range(dimension+1))
46 for(
auto subEntity : Dune::range(referenceElement.size(codim)))
48 std::size_t c = mapper.layout()(referenceElement.type(subEntity, codim), dimension);
51 std::size_t firstIndex = mapper.subIndex(element, subEntity, codim);
52 for(
auto j : Dune::range(firstIndex, firstIndex+c))
54 subIndices.push_back(j);
66 template<
class Mapper>
67 auto computeAverageSubEntityMeshSize(
const Mapper& mapper)
69 constexpr auto dimension = Mapper::GridView::dimension;
71 std::vector<unsigned int> adjacentElements(mapper.size(), 0);
72 std::vector<double> subEntityMeshSize(mapper.size(), 0.0);
73 for(
const auto& element : Dune::elements(mapper.gridView()))
75 auto A = element.geometry().volume();
76 for(
auto i : Impl::subIndexSet(mapper, element))
78 subEntityMeshSize[i] += A;
79 ++(adjacentElements[i]);
82 for(
auto i : Dune::range(mapper.size()))
83 subEntityMeshSize[i] = std::pow(subEntityMeshSize[i]/adjacentElements[i], 1./dimension);
84 return subEntityMeshSize;
95static constexpr auto evaluateMonomialValues(
const Dune::FieldVector<K,1>& x)
97 using Range = Dune::FieldVector<K,1>;
98 constexpr std::size_t maxOrder=3;
99 constexpr std::size_t size = (maxOrder+1);
100 auto xPowers = std::array<double,maxOrder+1>{};
102 for(
auto k: Dune::range(maxOrder))
103 xPowers[k+1] = xPowers[k]*x[0];
104 auto y = Dune::FieldVector<Range,size>{};
105 for(
auto order : Dune::range(maxOrder+1))
106 y[order] = xPowers[order];
112static constexpr auto evaluateMonomialJacobians(
const Dune::FieldVector<K,1>& x)
114 using Jacobian = Dune::FieldMatrix<K,1,1>;
115 constexpr std::size_t maxOrder=3;
116 constexpr std::size_t size = (maxOrder+1);
117 auto xPowers = std::array<double,maxOrder+1>{};
119 for(
auto k: Dune::range(maxOrder))
120 xPowers[k+1] = xPowers[k]*x[0];
121 auto y = Dune::FieldVector<Jacobian,size>{};
122 for(
auto order : Dune::range(std::size_t(1), maxOrder+1))
123 y[order][0][2] = order*xPowers[order-1];
129static constexpr auto evaluateMonomialValues(
const Dune::FieldVector<K,2>& x)
131 using Range = Dune::FieldVector<K,1>;
132 constexpr std::size_t maxOrder=3;
133 constexpr std::size_t dim=2;
134 constexpr std::size_t size = (maxOrder+1)*(maxOrder+2)/2;
135 auto xPowers = std::array<std::array<double,maxOrder+1>,dim>{};
136 for(
auto j: Dune::range(dim))
139 for(
auto k: Dune::range(maxOrder))
140 xPowers[j][k+1] = xPowers[j][k]*x[j];
142 auto y = Dune::FieldVector<Range,size>{};
144 for(
auto order : Dune::range(maxOrder+1))
146 for(
auto k : Dune::range(order+1))
148 y[index] = xPowers[0][order-k]*xPowers[1][k];
157static constexpr auto evaluateMonomialJacobians(
const Dune::FieldVector<K,2>& x)
159 using Jacobian = Dune::FieldMatrix<K,1,2>;
160 constexpr std::size_t maxOrder=3;
161 constexpr std::size_t dim=2;
162 constexpr std::size_t size = (maxOrder+1)*(maxOrder+2)/2;
163 auto xPowers = std::array<std::array<double,maxOrder+1>,dim>{};
164 for(
auto j: Dune::range(dim))
167 for(
auto k: Dune::range(maxOrder))
168 xPowers[j][k+1] = xPowers[j][k]*x[j];
170 auto y = Dune::FieldVector<Jacobian,size>{};
172 for(
auto order : Dune::range(maxOrder+1))
174 for(
auto k : Dune::range(order+1))
177 y[index][0][0] = (order-k)*xPowers[0][order-k-1]*xPowers[1][k];
179 y[index][0][1] = k*xPowers[0][order-k]*xPowers[1][k-1];
194template<
class DF,
class RF,
unsigned int dim,
bool reduced>
195class CubicHermiteLocalBasis
198 static constexpr auto makeReferenceBasisCoefficients() {
199 if constexpr (dim==1)
200 return Dune::FieldMatrix<int,4,4>{
206 if constexpr ((dim==2) and (not reduced))
207 return Dune::FieldMatrix<int,10,10>{
208 { 1, 0, 0, -3, -13, -3, 2, 13, 13, 2},
209 { 0, 1, 0, -2, -3, 0, 1, 3, 2, 0},
210 { 0, 0, 1, 0, -3, -2, 0, 2, 3, 1},
211 { 0, 0, 0, 3, -7, 0, -2, 7, 7, 0},
212 { 0, 0, 0, -1, 2, 0, 1, -2, -2, 0},
213 { 0, 0, 0, 0, -1, 0, 0, 2, 1, 0},
214 { 0, 0, 0, 0, -7, 3, 0, 7, 7, -2},
215 { 0, 0, 0, 0, -1, 0, 0, 1, 2, 0},
216 { 0, 0, 0, 0, 2, -1, 0, -2, -2, 1},
217 { 0, 0, 0, 0, 27, 0, 0, -27, -27, 0}
219 if constexpr ((dim==2) and (reduced))
221 auto w = std::array{1./3, 1./18, 1./18, 1./3, -1./9, 1./18, 1./3, 1./18, -1./9};
222 return Dune::FieldMatrix<double,9,10>{
223 { 1, 0, 0, -3, -13 + w[0]*27, -3, 2, 13 - w[0]*27, 13 - w[0]*27, 2},
224 { 0, 1, 0, -2, -3 + w[1]*27, 0, 1, 3 - w[1]*27, 2 - w[1]*27, 0},
225 { 0, 0, 1, 0, -3 + w[2]*27, -2, 0, 2 - w[2]*27, 3 - w[2]*27, 1},
226 { 0, 0, 0, 3, -7 + w[3]*27, 0, -2, 7 - w[3]*27, 7 - w[3]*27, 0},
227 { 0, 0, 0, -1, 2 + w[4]*27, 0, 1, -2 - w[4]*27, -2 - w[4]*27, 0},
228 { 0, 0, 0, 0, -1 + w[5]*27, 0, 0, 2 - w[5]*27, 1 - w[5]*27, 0},
229 { 0, 0, 0, 0, -7 + w[6]*27, 3, 0, 7 - w[6]*27, 7 - w[6]*27, -2},
230 { 0, 0, 0, 0, -1 + w[7]*27, 0, 0, 1 - w[7]*27, 2 - w[7]*27, 0},
231 { 0, 0, 0, 0, 2 + w[8]*27, -1, 0, -2 - w[8]*27, -2 - w[8]*27, 1},
247 static constexpr auto referenceBasisCoefficients = makeReferenceBasisCoefficients();
255 template<
class LambdaRefValues,
class Entry>
256 void transformToElementBasis(
const LambdaRefValues& refValues, std::vector<Entry>& out)
const
258 if constexpr (dim==1)
260 const auto& J = elementJacobian_;
261 out.resize(refValues.size());
262 out[0] = refValues[0];
263 out[1] = J*refValues[1] / (*localSubEntityMeshSize_)[1];
264 out[2] = refValues[2];
265 out[3] = J*refValues[3] / (*localSubEntityMeshSize_)[1];;
267 if constexpr (dim==2)
269 const auto& J = elementJacobian_;
270 out.resize(refValues.size());
271 out[0] = refValues[0];
272 out[1] = (J[0][0]*refValues[1] + J[0][1]*refValues[2]) / (*localSubEntityMeshSize_)[1];
273 out[2] = (J[1][0]*refValues[1] + J[1][1]*refValues[2]) / (*localSubEntityMeshSize_)[2];
274 out[3] = refValues[3];
275 out[4] = (J[0][0]*refValues[4] + J[0][1]*refValues[5]) / (*localSubEntityMeshSize_)[4];
276 out[5] = (J[1][0]*refValues[4] + J[1][1]*refValues[5]) / (*localSubEntityMeshSize_)[5];
277 out[6] = refValues[6];
278 out[7] = (J[0][0]*refValues[7] + J[0][1]*refValues[8]) / (*localSubEntityMeshSize_)[7];
279 out[8] = (J[1][0]*refValues[7] + J[1][1]*refValues[8]) / (*localSubEntityMeshSize_)[8];
280 if constexpr (not reduced)
281 out[9] = refValues[9];
285 using ElementJacobian = Dune::FieldMatrix<DF, dim,dim>;
287 using LocalSubEntityMeshSize = std::vector<double>;
291 using Domain = Dune::FieldVector<DF, dim>;
292 using Range = Dune::FieldVector<RF, 1>;
293 using Jacobian = Dune::FieldMatrix<RF, 1, dim>;
294 using Traits = Dune::LocalBasisTraits<DF, dim, Domain, RF, 1, Range, Jacobian>;
295 using OrderArray = std::array<unsigned int, dim>;
297 static constexpr unsigned int size()
299 return decltype(referenceBasisCoefficients)::rows;
302 inline void evaluateFunction(
const Domain& x, std::vector<Range>& values)
const
304 auto monomialValues = evaluateMonomialValues(x);
305 auto referenceValues = Dune::FieldVector<Range, size()>{};
306 referenceBasisCoefficients.mv(monomialValues, referenceValues);
307 transformToElementBasis(referenceValues, values);
310 inline void evaluateJacobian(
const Domain& x, std::vector<Jacobian>& jacobians)
const
312 auto monomialJacobians = evaluateMonomialJacobians(x);
313 auto referenceJacobians = Dune::FieldVector<Jacobian, size()>{};
314 referenceBasisCoefficients.mv(monomialJacobians, referenceJacobians);
315 transformToElementBasis(referenceJacobians, jacobians);
318 void partial(
const OrderArray& order,
const Domain& x, std::vector<Range>& out)
const
320 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
322 evaluateFunction(x, out);
323 DUNE_THROW(RangeError,
"partial() not implemented for given order");
326 unsigned int order()
const
331 template<
class Element>
332 void bind(
const Element& element,
const LocalSubEntityMeshSize& localSubEntityMeshSize) {
333 localSubEntityMeshSize_ = &localSubEntityMeshSize;
334 auto center = Dune::ReferenceElements<DF, dim>::simplex().position(0, 0);
335 elementJacobian_ = element.geometry().jacobian(center);
339 ElementJacobian elementJacobian_;
340 const LocalSubEntityMeshSize* localSubEntityMeshSize_;
345template<
unsigned int dim,
bool reduced>
346struct CubicHermiteLocalCoefficients
349 static constexpr auto makeLocalKeys() {
350 if constexpr (dim==1)
357 if constexpr ((dim==2) and (not reduced))
370 if constexpr ((dim==2) and (reduced))
384 using LocalKeys = std::decay_t<
decltype(makeLocalKeys())>;
388 std::size_t size()
const
390 return localKeys_.size();
393 const LocalKey& localKey(std::size_t i)
const
395 assert( i < localKeys_.size() );
396 return localKeys_[i];
399 LocalKeys localKeys_ = makeLocalKeys();
404template<
class DF,
class RF,
unsigned int dim,
bool reduced>
405class CubicHermiteLocalInterpolation
407 using ElementJacobianInverse = Dune::FieldMatrix<DF, dim,dim>;
408 using LocalSubEntityMeshSize = std::vector<double>;
412 template<
class Element>
413 void bind(
const Element& element,
const LocalSubEntityMeshSize& localSubEntityMeshSize) {
414 localSubEntityMeshSize_ = &localSubEntityMeshSize;
417 template<
class F,
class C>
418 void interpolate(
const F& f, std::vector<C>& out)
const
420 using Domain = Dune::FieldVector<DF, dim>;
422 if constexpr (dim==1)
426 out[1] = df(0) * (*localSubEntityMeshSize_)[1];
428 out[3] = df(1) * (*localSubEntityMeshSize_)[3];
430 if constexpr (dim==2)
432 if constexpr (not reduced)
435 out[9] = f(Domain({1.0/3.0,1.0/3.0}));
437 if constexpr (reduced)
439 auto J0 = df(Domain({0,0}));
440 auto J1 = df(Domain({1,0}));
441 auto J2 = df(Domain({0,1}));
442 out[0] = f(Domain({0,0}));
443 out[1] = J0[0][0] * (*localSubEntityMeshSize_)[1];
444 out[2] = J0[0][1] * (*localSubEntityMeshSize_)[2];
445 out[3] = f(Domain({1,0}));
446 out[4] = J1[0][0] * (*localSubEntityMeshSize_)[4];
447 out[5] = J1[0][1] * (*localSubEntityMeshSize_)[5];
448 out[6] = f(Domain({0,1}));
449 out[7] = J2[0][0] * (*localSubEntityMeshSize_)[7];
450 out[8] = J2[0][1] * (*localSubEntityMeshSize_)[8];
454 const LocalSubEntityMeshSize* localSubEntityMeshSize_;
466template<
class DF,
class RF,
unsigned int dim,
bool reduced>
467class CubicHermiteLocalFiniteElement
469 using LocalBasis = CubicHermiteLocalBasis<DF, RF, dim, reduced>;
470 using LocalCoefficients = CubicHermiteLocalCoefficients<dim, reduced>;
471 using LocalInterpolation = CubicHermiteLocalInterpolation<DF, RF, dim, reduced>;
472 using LocalSubEntityMeshSize = std::vector<double>;
476 using Traits = LocalFiniteElementTraits<LocalBasis, LocalCoefficients, LocalInterpolation>;
478 const LocalBasis& localBasis()
const {
482 const LocalCoefficients& localCoefficients()
const {
483 return localCoefficients_;
486 const LocalInterpolation& localInterpolation()
const {
487 return localInterpolation_;
490 unsigned int size()
const {
491 return localBasis_.size();
494 GeometryType type()
const {
498 template<
class Element,
class Mapper,
class MeshSizeContainer>
499 void bind(
const Element& element,
const Mapper& mapper,
const MeshSizeContainer& subEntityMeshSize) {
501 localSubEntityMeshSize_.resize(localCoefficients_.size());
502 for(
auto i : Dune::range(size()))
504 auto localKey = localCoefficients_.localKey(i);
505 auto subEntityIndex = mapper.subIndex(element, localKey.subEntity(), localKey.codim());
506 localSubEntityMeshSize_[i] = subEntityMeshSize[subEntityIndex];
509 localBasis_.bind(element, localSubEntityMeshSize_);
510 localInterpolation_.bind(element, localSubEntityMeshSize_);
511 type_ = element.type();
515 LocalSubEntityMeshSize localSubEntityMeshSize_;
516 LocalBasis localBasis_;
517 LocalCoefficients localCoefficients_;
518 LocalInterpolation localInterpolation_;
540template<
typename GV,
bool reduced>
541class CubicHermiteNode :
544 static const int gridDim = GV::dimension;
546 using CubeFiniteElement = Impl::CubicHermiteLocalFiniteElement<typename GV::ctype, double, gridDim, reduced>;
547 using SubEntityMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GV>;
548 using Base = LeafBasisNode;
552 using size_type = std::size_t;
553 using Element =
typename GV::template Codim<0>::Entity;
554 using FiniteElement = CubeFiniteElement;
557 CubicHermiteNode(
const SubEntityMapper& subEntityMapper,
const std::vector<double>& subEntityMeshSize) :
560 subEntityMapper_(&subEntityMapper),
561 subEntityMeshSize_(&subEntityMeshSize)
565 const Element& element()
const
574 const FiniteElement& finiteElement()
const
576 return finiteElement_;
580 void bind(
const Element& e)
583 finiteElement_.bind(*element_, *subEntityMapper_, *subEntityMeshSize_);
584 this->setSize(finiteElement_.size());
589 FiniteElement finiteElement_;
590 const Element* element_;
591 const SubEntityMapper* subEntityMapper_;
592 const std::vector<double>* subEntityMeshSize_;
604template<
typename GV,
bool reduced = false>
609 using SubEntityMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GV>;
611 static constexpr auto cubicHermiteMapperLayout(Dune::GeometryType type,
int gridDim) {
614 if ((type.isTriangle()) and (not reduced))
620 static constexpr auto subEntityLayout(Dune::GeometryType type,
int gridDim) {
621 return (cubicHermiteMapperLayout(type, gridDim) > 0);
627 using Node = CubicHermiteNode<GridView, reduced>;
630 Base(gv, cubicHermiteMapperLayout),
631 subEntityMapper_(gv, subEntityLayout)
633 static_assert(GridView::dimension<=2,
"CubicHermitePreBasis is only implemented for dim<=2");
634 subEntityMeshSize_ = Impl::computeAverageSubEntityMeshSize(subEntityMapper_);
640 subEntityMapper_.update(gv);
641 subEntityMeshSize_ = Impl::computeAverageSubEntityMeshSize(subEntityMapper_);
644 Node makeNode()
const
646 return Node{subEntityMapper_, subEntityMeshSize_};
650 std::vector<double> subEntityMeshSize_;
651 SubEntityMapper subEntityMapper_;
656namespace BasisFactory {
663template<
class Dummy=
void>
666 return [](
const auto& gridView) {
676template<
class Dummy=
void>
679 return [](
const auto& gridView) {
Pre-basis for a CubicHermite basis.
Definition: cubichermitebasis.hh:606
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
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:664
auto reducedCubicHermite()
Create a pre-basis factory that can create a CubicHermite pre-basis.
Definition: cubichermitebasis.hh:677
Definition: polynomial.hh:17