DUNE PDELab (git)

hierarchicalsimplexp2withelementbubble.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// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5#ifndef DUNE_HIERARCHICAL_SIMPLEX_P2_WITH_ELEMENT_BUBBLE_LOCALBASIS_HH
6#define DUNE_HIERARCHICAL_SIMPLEX_P2_WITH_ELEMENT_BUBBLE_LOCALBASIS_HH
7
12#include <array>
13#include <cassert>
14#include <numeric>
15#include <stdexcept>
16#include <vector>
17
20#include <dune/common/math.hh>
21
22#include <dune/geometry/referenceelement.hh>
23
24#include <dune/localfunctions/common/localbasis.hh>
25#include <dune/localfunctions/common/localkey.hh>
26
27namespace Dune
28{
44 template<class D, class R, int dim>
46 {
47 template<class,int> friend class HierarchicalSimplexP2WithElementBubbleLocalInterpolation;
48
51
54
57
58 // Number of vertices
59 static constexpr int numVertices = dim+1;
60
61 // Number of edges (or zero for dim==1)
62 static constexpr int numEdges = (dim > 1 ? ((dim+1)*dim / 2) : 0);
63
64 // helper function to evaluate the vertex basis functions
65 template <class It>
66 static constexpr It evaluateVertexFunctions (const DomainType& in, It outIt)
67 {
68 *outIt = 1;
69 for (int i = 0; i < dim; ++i)
70 *outIt -= in[i];
71 ++outIt;
72 for (int i = 0; i < dim; ++i)
73 *outIt++ = in[i];
74 return outIt;
75 }
76
77 // helper function to evaluate the basis functions
78 template <class It>
79 static constexpr It evaluateAllFunctions (const DomainType& in, It outIt)
80 {
81 It vertexValues = outIt;
82 outIt = evaluateVertexFunctions(in, outIt);
83
84 if constexpr(dim > 1) {
85 auto refElem = referenceElement<D,dim>(GeometryTypes::simplex(dim));
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];
90 }
91 }
92
93 // element bubble function
94 *outIt = power(dim+1, dim+1);
95 for (int i = 0; i < numVertices; ++i)
96 *outIt *= vertexValues[i];
97 return outIt;
98 }
99
100 public:
103
105 static constexpr std::size_t size () noexcept
106 {
107 return numVertices + numEdges + 1;
108 }
109
111 static constexpr void evaluateFunction (const DomainType& in,
112 std::vector<RangeType>& out)
113 {
114 out.resize(size());
115 evaluateAllFunctions(in,out.begin());
116 }
117
119 static constexpr void evaluateJacobian (const DomainType& in,
120 std::vector<JacobianType>& out)
121 {
122 out.resize(size());
123
124 // vertex basis functions
125 RangeType tmp = 1;
126 for (int i = 0; i < dim; ++i) {
127 out[0][0][i] = -1;
128 for (int j = 0; j < dim; ++j)
129 out[j+1][0][i] = (i == j);
130 tmp -= in[i];
131 }
132
133 int n = numVertices;
134 std::array<RangeType,numVertices> shapeValues;
135 evaluateVertexFunctions(in, shapeValues.begin());
136
137 // edge basis functions
138 if constexpr(dim > 1) {
139 auto refElem = referenceElement<D,dim>(GeometryTypes::simplex(dim));
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]);
145 }
146 }
147
148 // element bubble function
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];
153 }
154 }
155
157 static constexpr void partial (const std::array<unsigned int, dim>& order,
158 const DomainType& in,
159 std::vector<RangeType>& out)
160 {
161 unsigned int totalOrder = 0;
162 for (int i = 0; i < dim; ++i)
163 totalOrder += order[i];
164
165 switch (totalOrder) {
166 case 0:
167 evaluateFunction(in,out);
168 break;
169 case 1: {
170 out.resize(size());
171 int d = 0; // the direction of differentiation
172 for (int i = 0; i < dim; ++i)
173 d += i * order[i];
174
175 // vertex basis functions
176 RangeType tmp = 1;
177 for (int i = 0; i < dim; ++i) {
178 out[0] = -1;
179 for (int j = 0; j < dim; ++j)
180 out[j+1] = (dim == j);
181 tmp -= in[i];
182 }
183
184 int n = numVertices;
185 std::array<RangeType,numVertices> shapeValues;
186 evaluateVertexFunctions(in, shapeValues.begin());
187
188 // edge basis functions
189 if constexpr(dim > 1) {
190 auto refElem = referenceElement<D,dim>(GeometryTypes::simplex(dim));
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]);
195 }
196 }
197
198 // element bubble function
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];
202 } break;
203 default:
204 throw std::runtime_error("Desired derivative order is not implemented");
205 }
206 }
207
209 static constexpr unsigned int order () noexcept
210 {
211 return dim+1;
212 }
213 };
214
215
227 template <int dim>
229 {
230 // Number of vertices
231 static constexpr int numVertices = dim+1;
232
233 // Number of edges (or zero for dim==1)
234 static constexpr int numEdges = (dim > 1 ? ((dim+1)*dim / 2) : 0);
235
236 public:
239 {
240 int n = 0;
241 for (int i = 0; i < numVertices; ++i)
242 li_[n++] = LocalKey(i,dim,0); // Vertex
243 if constexpr(dim > 1) {
244 for (int i = 0; i < numEdges; ++i)
245 li_[n++] = LocalKey(i,dim-1,0); // Edges
246 }
247 li_[n++] = LocalKey(0,0,0); // Element
248 }
249
251 static constexpr std::size_t size () noexcept
252 {
253 return numVertices + numEdges + 1;
254 }
255
257 const LocalKey& localKey (std::size_t i) const noexcept
258 {
259 return li_[i];
260 }
261
262 private:
263 std::array<LocalKey, numVertices+numEdges+1> li_;
264 };
265
269 template<class LB, int dim>
270 class HierarchicalSimplexP2WithElementBubbleLocalInterpolation
271 {
272 using LocalBasis = LB;
273 using DomainType = typename LB::Traits::DomainType;
274 using RangeType = typename LB::Traits::RangeType;
275
276 // Number of vertices
277 static constexpr int numVertices = dim+1;
278
279 // Number of edges (or zero for dim==1)
280 static constexpr int numEdges = (dim > 1 ? ((dim+1)*dim / 2) : 0);
281
282 public:
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)
294 {
295 auto refElem = referenceElement<typename LB::Traits::DomainFieldType,dim>(GeometryTypes::simplex(dim));
296
297 out.resize(LB::size());
298 int n = 0;
299
300 // vertices
301 assert(numVertices == refElem.size(dim));
302 for (int i = 0; i < numVertices; ++i)
303 out[n++] = f(refElem.position(i,dim));
304
305 std::array<RangeType,LB::size()> shapeValues;
306
307 // edge bubbles
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];
315 out[n++] = y;
316 }
317 }
318
319 // element bubble
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];
324 out[n++] = y;
325 }
326 };
327
328} // end namespace Dune
329
330#endif // DUNE_HIERARCHICAL_SIMPLEX_P2_WITH_ELEMENT_BUBBLE_LOCALBASIS_HH
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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)