5#ifndef DUNE_HIERARCHICAL_SIMPLEX_P2_WITH_ELEMENT_BUBBLE_LOCALBASIS_HH
6#define DUNE_HIERARCHICAL_SIMPLEX_P2_WITH_ELEMENT_BUBBLE_LOCALBASIS_HH
22#include <dune/geometry/referenceelement.hh>
24#include <dune/localfunctions/common/localbasis.hh>
25#include <dune/localfunctions/common/localkey.hh>
44 template<
class D,
class R,
int dim>
47 template<
class,
int>
friend class HierarchicalSimplexP2WithElementBubbleLocalInterpolation;
59 static constexpr int numVertices = dim+1;
62 static constexpr int numEdges = (dim > 1 ? ((dim+1)*dim / 2) : 0);
66 static constexpr It evaluateVertexFunctions (
const DomainType& in, It outIt)
69 for (
int i = 0; i < dim; ++i)
72 for (
int i = 0; i < dim; ++i)
79 static constexpr It evaluateAllFunctions (
const DomainType& in, It outIt)
81 It vertexValues = outIt;
82 outIt = evaluateVertexFunctions(in, outIt);
84 if constexpr(dim > 1) {
86 for (
int i = 0; i < numEdges; ++i) {
87 const int v0 = refElem.subEntity(i,dim-1,0,dim);
88 const int v1 = refElem.subEntity(i,dim-1,1,dim);
89 *outIt++ = 4 * vertexValues[v0] * vertexValues[v1];
94 *outIt =
power(dim+1, dim+1);
95 for (
int i = 0; i < numVertices; ++i)
96 *outIt *= vertexValues[i];
105 static constexpr std::size_t
size () noexcept
107 return numVertices + numEdges + 1;
112 std::vector<RangeType>& out)
115 evaluateAllFunctions(in,out.
begin());
120 std::vector<JacobianType>& out)
126 for (
int i = 0; i < dim; ++i) {
128 for (
int j = 0; j < dim; ++j)
129 out[j+1][0][i] = (i == j);
134 std::array<RangeType,numVertices> shapeValues;
135 evaluateVertexFunctions(in, shapeValues.
begin());
138 if constexpr(dim > 1) {
140 for (
int i = 0; i < numEdges; ++i,++n) {
141 const int v0 = refElem.subEntity(i,dim-1,0,dim);
142 const int v1 = refElem.subEntity(i,dim-1,1,dim);
143 for (
int j = 0; j < dim; ++j)
144 out[n][0][j] = 4 * (out[v0][0][j] * shapeValues[v1] + shapeValues[v0] * out[v1][0][j]);
149 for (
int i = 0; i < dim; ++i) {
150 out[n][0][i] =
power(dim+1, dim+1) * (tmp - in[i]);
151 for (
int j = i+1; j < dim+i; ++j)
152 out[n][0][i] *= in[j % dim];
157 static constexpr void partial (
const std::array<unsigned int, dim>&
order,
159 std::vector<RangeType>& out)
161 unsigned int totalOrder = 0;
162 for (
int i = 0; i < dim; ++i)
163 totalOrder +=
order[i];
165 switch (totalOrder) {
172 for (
int i = 0; i < dim; ++i)
177 for (
int i = 0; i < dim; ++i) {
179 for (
int j = 0; j < dim; ++j)
180 out[j+1] = (dim == j);
185 std::array<RangeType,numVertices> shapeValues;
186 evaluateVertexFunctions(in, shapeValues.
begin());
189 if constexpr(dim > 1) {
191 for (
int i = 0; i < numEdges; ++i,++n) {
192 const int v0 = refElem.subEntity(i,dim-1,0,dim);
193 const int v1 = refElem.subEntity(i,dim-1,1,dim);
194 out[n] = 4 * (out[v0] * shapeValues[v1] + shapeValues[v0] * out[v1]);
199 out[n] =
power(dim+1, dim+1) * (tmp - in[d]);
200 for (
int j = d+1; j < dim+d; ++j)
201 out[n] *= in[j % dim];
204 throw std::runtime_error(
"Desired derivative order is not implemented");
209 static constexpr unsigned int order () noexcept
231 static constexpr int numVertices = dim+1;
234 static constexpr int numEdges = (dim > 1 ? ((dim+1)*dim / 2) : 0);
241 for (
int i = 0; i < numVertices; ++i)
243 if constexpr(dim > 1) {
244 for (
int i = 0; i < numEdges; ++i)
251 static constexpr std::size_t
size () noexcept
253 return numVertices + numEdges + 1;
263 std::array<LocalKey, numVertices+numEdges+1> li_;
269 template<
class LB,
int dim>
270 class HierarchicalSimplexP2WithElementBubbleLocalInterpolation
272 using LocalBasis = LB;
273 using DomainType =
typename LB::Traits::DomainType;
274 using RangeType =
typename LB::Traits::RangeType;
277 static constexpr int numVertices = dim+1;
280 static constexpr int numEdges = (dim > 1 ? ((dim+1)*dim / 2) : 0);
290 template<
class F,
class C,
291 class R = std::invoke_result_t<F, DomainType>,
292 std::enable_if_t<std::is_convertible_v<R, C>,
int> = 0>
293 static constexpr void interpolate (
const F& f, std::vector<C>& out)
297 out.resize(LB::size());
301 assert(numVertices == refElem.size(dim));
302 for (
int i = 0; i < numVertices; ++i)
303 out[n++] = f(refElem.position(i,dim));
305 std::array<RangeType,LB::size()> shapeValues;
308 if constexpr(dim > 1) {
309 assert(numEdges == refElem.size(dim-1));
310 for (
int i = 0; i < numEdges; ++i) {
311 R y = f(refElem.position(i,dim-1));
312 LB::evaluateVertexFunctions(refElem.position(i,dim-1), shapeValues.begin());
313 for (
int j = 0; j < numVertices; ++j)
314 y -= out[j]*shapeValues[j];
320 R y = f(refElem.position(0,0));
321 LB::evaluateAllFunctions(refElem.position(0,0), shapeValues.begin());
322 for (
int j = 0; j < numVertices+numEdges; ++j)
323 y -= out[j]*shapeValues[j];
Iterator begin()
begin iterator
Definition: densevector.hh:347
A dense n x m matrix.
Definition: fmatrix.hh:117
P1 basis in dim-d enriched by quadratic edge bubble functions and an element bubble function of order...
Definition: hierarchicalsimplexp2withelementbubble.hh:46
static constexpr void evaluateFunction(const DomainType &in, std::vector< RangeType > &out)
Evaluate all shape functions.
Definition: hierarchicalsimplexp2withelementbubble.hh:111
static constexpr void partial(const std::array< unsigned int, dim > &order, const DomainType &in, std::vector< RangeType > &out)
Evaluate partial derivatives of all shape functions.
Definition: hierarchicalsimplexp2withelementbubble.hh:157
static constexpr unsigned int order() noexcept
Polynomial order of the shape functions (4 in this case)
Definition: hierarchicalsimplexp2withelementbubble.hh:209
static constexpr void evaluateJacobian(const DomainType &in, std::vector< JacobianType > &out)
Evaluate Jacobian of all shape functions.
Definition: hierarchicalsimplexp2withelementbubble.hh:119
static constexpr std::size_t size() noexcept
Returns number of shape functions.
Definition: hierarchicalsimplexp2withelementbubble.hh:105
The local keys of the hierarchical basis functions with element bubble.
Definition: hierarchicalsimplexp2withelementbubble.hh:229
HierarchicalSimplexP2WithElementBubbleLocalCoefficients() noexcept
Default constructor, initializes the local keys.
Definition: hierarchicalsimplexp2withelementbubble.hh:238
static constexpr std::size_t size() noexcept
Returns number of coefficients.
Definition: hierarchicalsimplexp2withelementbubble.hh:251
const LocalKey & localKey(std::size_t i) const noexcept
Returns the i'th local key.
Definition: hierarchicalsimplexp2withelementbubble.hh:257
Describe position of one degree of freedom.
Definition: localkey.hh:24
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
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
Some useful basic math stuff.
Dune namespace.
Definition: alignedallocator.hh:13
constexpr Base power(Base m, Exponent p)
Power method for integer exponents.
Definition: math.hh:75
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:35