DUNE PDELab (unstable)

lagrangebasis.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3
4// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
5// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
6
7#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEBASIS_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEBASIS_HH
9
10#include <type_traits>
12
13#include <dune/localfunctions/lagrange/cache.hh>
14
15#include <dune/functions/functionspacebases/nodes.hh>
16#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
17#include <dune/functions/functionspacebases/leafprebasismixin.hh>
18
20
21namespace Dune {
22namespace Functions {
23
24// *****************************************************************************
25// This is the reusable part of the LagrangeBasis. It contains
26//
27// LagrangePreBasis
28// LagrangeNode
29//
30// The pre-basis allows to create the others and is the owner of possible shared
31// state. These components do _not_ depend on the global basis and local view
32// and can be used without a global basis.
33// *****************************************************************************
34
35template<typename GV, int k, typename R=double>
36class LagrangeNode;
37
38template<typename GV, int k, typename R=double>
39class LagrangePreBasis;
40
41
42
62template<typename GV, int k, typename R>
64 public LeafPreBasisMixin< LagrangePreBasis<GV,k,R> >
65{
66 static const int dim = GV::dimension;
67 static const bool useDynamicOrder = (k<0);
68
69public:
70
72 using GridView = GV;
73
75 using size_type = std::size_t;
76
78 using Node = LagrangeNode<GV, k, R>;
79
82 : LagrangePreBasis(gv, std::numeric_limits<unsigned int>::max())
83 {}
84
86 LagrangePreBasis(const GridView& gv, unsigned int order) :
87 gridView_(gv), order_(order)
88 {
89 if (!useDynamicOrder && order!=std::numeric_limits<unsigned int>::max())
90 DUNE_THROW(RangeError, "Template argument k has to be -1 when supplying a run-time order!");
91
92 for (int i=0; i<=dim; i++)
93 {
94 dofsPerCube_[i] = computeDofsPerCube(i);
95 dofsPerSimplex_[i] = computeDofsPerSimplex(i);
96 }
97 dofsPerPrism_ = computeDofsPerPrism();
98 dofsPerPyramid_ = computeDofsPerPyramid();
99 }
100
103 {
104 vertexOffset_ = 0;
105 edgeOffset_ = vertexOffset_ + dofsPerCube(0) * ((size_type)gridView_.size(dim));
106
107 if (dim>=2)
108 {
109 triangleOffset_ = edgeOffset_ + dofsPerCube(1) * ((size_type) gridView_.size(dim-1));
110
111 quadrilateralOffset_ = triangleOffset_ + dofsPerSimplex(2) * ((size_type)gridView_.size(Dune::GeometryTypes::triangle));
112 }
113
114 if (dim==3) {
115 tetrahedronOffset_ = quadrilateralOffset_ + dofsPerCube(2) * ((size_type)gridView_.size(Dune::GeometryTypes::quadrilateral));
116
117 prismOffset_ = tetrahedronOffset_ + dofsPerSimplex(3) * ((size_type)gridView_.size(Dune::GeometryTypes::tetrahedron));
118
119 hexahedronOffset_ = prismOffset_ + dofsPerPrism() * ((size_type)gridView_.size(Dune::GeometryTypes::prism));
120
121 pyramidOffset_ = hexahedronOffset_ + dofsPerCube(3) * ((size_type)gridView_.size(Dune::GeometryTypes::hexahedron));
122 }
123 }
124
126 const GridView& gridView() const
127 {
128 return gridView_;
129 }
130
132 void update (const GridView& gv)
133 {
134 gridView_ = gv;
135 }
136
141 {
142 return Node{order_};
143 }
144
147 {
148 switch (dim)
149 {
150 case 1:
151 return dofsPerCube(0) * ((size_type)gridView_.size(dim))
152 + dofsPerCube(1) * ((size_type)gridView_.size(dim-1));
153 case 2:
154 {
155 return dofsPerCube(0) * ((size_type)gridView_.size(dim))
156 + dofsPerCube(1) * ((size_type)gridView_.size(dim-1))
159 }
160 case 3:
161 {
162 return dofsPerCube(0) * ((size_type)gridView_.size(dim))
163 + dofsPerCube(1) * ((size_type)gridView_.size(dim-1))
167 + dofsPerPyramid() * ((size_type)gridView_.size(Dune::GeometryTypes::pyramid))
168 + dofsPerPrism() * ((size_type)gridView_.size(Dune::GeometryTypes::prism))
169 + dofsPerCube(3) * ((size_type)gridView_.size(Dune::GeometryTypes::hexahedron));
170 }
171 }
172 DUNE_THROW(Dune::NotImplemented, "No size method for " << dim << "d grids available yet!");
173 }
174
177 {
178 // That cast to unsigned int is necessary because GV::dimension is an enum,
179 // which is not recognized by the power method as an integer type...
180 return power(order()+1, (unsigned int)GV::dimension);
181 }
182
183 template<typename It>
184 It indices(const Node& node, It it) const
185 {
186 for (size_type i = 0, end = node.finiteElement().size() ; i < end ; ++it, ++i)
187 {
188 Dune::LocalKey localKey = node.finiteElement().localCoefficients().localKey(i);
189 const auto& gridIndexSet = gridView().indexSet();
190 const auto& element = node.element();
191
192 // The dimension of the entity that the current dof is related to
193 auto dofDim = dim - localKey.codim();
194
195 // Test for a vertex dof
196 // The test for k==1 is redundant, but having it here allows the compiler to conclude
197 // at compile-time that the dofDim==0 case is the only one that will ever happen.
198 // This leads to measurable speed-up: see
199 // https://gitlab.dune-project.org/staging/dune-functions/issues/30
200 if (k==1 || dofDim==0) {
201 *it = {{ (size_type)(gridIndexSet.subIndex(element,localKey.subEntity(),dim)) }};
202 continue;
203 }
204
205 if (dofDim==1)
206 { // edge dof
207 if (dim==1) // element dof -- any local numbering is fine
208 {
209 *it = {{ edgeOffset_
210 + dofsPerCube(1) * ((size_type)gridIndexSet.subIndex(element,0,0))
211 + localKey.index() }};
212 continue;
213 }
214 else
215 {
216 const auto refElement
217 = Dune::referenceElement<double,dim>(element.type());
218
219 // We have to reverse the numbering if the local element edge is
220 // not aligned with the global edge.
221 auto v0 = (size_type)gridIndexSet.subIndex(element,refElement.subEntity(localKey.subEntity(),localKey.codim(),0,dim),dim);
222 auto v1 = (size_type)gridIndexSet.subIndex(element,refElement.subEntity(localKey.subEntity(),localKey.codim(),1,dim),dim);
223 bool flip = (v0 > v1);
224 *it = {{ (flip)
225 ? edgeOffset_
226 + dofsPerCube(1)*((size_type)gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()))
227 + (dofsPerCube(1)-1)-localKey.index()
228 : edgeOffset_
229 + dofsPerCube(1)*((size_type)gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()))
230 + localKey.index() }};
231 continue;
232 }
233 }
234
235 if (dofDim==2)
236 {
237 if (dim==2) // element dof -- any local numbering is fine
238 {
239 if (element.type().isTriangle())
240 {
241 *it = {{ triangleOffset_ + dofsPerSimplex(2)*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
242 continue;
243 }
244 else if (element.type().isQuadrilateral())
245 {
246 *it = {{ quadrilateralOffset_ + dofsPerCube(2)*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
247 continue;
248 }
249 else
250 DUNE_THROW(Dune::NotImplemented, "2d elements have to be triangles or quadrilaterals");
251 } else
252 {
253 const auto refElement
254 = Dune::referenceElement<double,dim>(element.type());
255
256 if (order()>3)
257 DUNE_THROW(Dune::NotImplemented, "LagrangeBasis for 3D grids is only implemented if k<=3");
258
259 if (order()==3 and !refElement.type(localKey.subEntity(), localKey.codim()).isTriangle())
260 DUNE_THROW(Dune::NotImplemented, "LagrangeBasis for 3D grids with k==3 is only implemented if the grid is a simplex grid");
261
262 *it = {{ triangleOffset_ + ((size_type)gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim())) }};
263 continue;
264 }
265 }
266
267 if (dofDim==3)
268 {
269 if (dim==3) // element dof -- any local numbering is fine
270 {
271 if (element.type().isTetrahedron())
272 {
273 *it = {{ tetrahedronOffset_ + dofsPerSimplex(3)*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
274 continue;
275 }
276 else if (element.type().isHexahedron())
277 {
278 *it = {{ hexahedronOffset_ + dofsPerCube(3)*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
279 continue;
280 }
281 else if (element.type().isPrism())
282 {
283 *it = {{ prismOffset_ + dofsPerPrism()*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
284 continue;
285 }
286 else if (element.type().isPyramid())
287 {
288 *it = {{ pyramidOffset_ + dofsPerPyramid()*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
289 continue;
290 }
291 else
292 DUNE_THROW(Dune::NotImplemented, "3d elements have to be tetrahedra, hexahedra, prisms, or pyramids");
293 } else
294 DUNE_THROW(Dune::NotImplemented, "Grids of dimension larger than 3 are no supported");
295 }
296 DUNE_THROW(Dune::NotImplemented, "Grid contains elements not supported for the LagrangeBasis");
297 }
298 return it;
299 }
300
302 unsigned int order() const
303 {
304 return (useDynamicOrder) ? order_ : k;
305 }
306
307protected:
308 GridView gridView_;
309
310 // Run-time order, only valid if k<0
311 unsigned int order_;
312
314 size_type dofsPerSimplex(std::size_t simplexDim) const
315 {
316 return useDynamicOrder ? dofsPerSimplex_[simplexDim] : computeDofsPerSimplex(simplexDim);
317 }
318
320 size_type dofsPerCube(std::size_t cubeDim) const
321 {
322 return useDynamicOrder ? dofsPerCube_[cubeDim] : computeDofsPerCube(cubeDim);
323 }
324
325 size_type dofsPerPrism() const
326 {
327 return useDynamicOrder ? dofsPerPrism_ : computeDofsPerPrism();
328 }
329
330 size_type dofsPerPyramid() const
331 {
332 return useDynamicOrder ? dofsPerPyramid_ : computeDofsPerPyramid();
333 }
334
336 size_type computeDofsPerSimplex(std::size_t simplexDim) const
337 {
338 return order() == 0 ? (dim == simplexDim ? 1 : 0) : Dune::binomial(std::size_t(order()-1),simplexDim);
339 }
340
342 size_type computeDofsPerCube(std::size_t cubeDim) const
343 {
344 return order() == 0 ? (dim == cubeDim ? 1 : 0) : Dune::power(order()-1, cubeDim);
345 }
346
347 size_type computeDofsPerPrism() const
348 {
349 return order() == 0 ? (dim == 3 ? 1 : 0) : (order()-1)*(order()-1)*(order()-2)/2;
350 }
351
352 size_type computeDofsPerPyramid() const
353 {
354 return order() == 0 ? (dim == 3 ? 1 : 0) : (order()-2)*(order()-1)*(2*order()-3)/6;
355 }
356
357 // When the order is given at run-time, the following numbers are pre-computed:
358 std::array<size_type,dim+1> dofsPerSimplex_;
359 std::array<size_type,dim+1> dofsPerCube_;
360 size_type dofsPerPrism_;
361 size_type dofsPerPyramid_;
362
363 size_type vertexOffset_;
364 size_type edgeOffset_;
365 size_type triangleOffset_;
366 size_type quadrilateralOffset_;
367 size_type tetrahedronOffset_;
368 size_type pyramidOffset_;
369 size_type prismOffset_;
370 size_type hexahedronOffset_;
371
372};
373
374
375
376template<typename GV, int k, typename R>
377class LagrangeNode :
378 public LeafBasisNode
379{
380 static constexpr int dim = GV::dimension;
381 static constexpr bool useDynamicOrder = (k<0);
382
383 // Compute the GeometryType id in case the grid has only a single GeometryType
384 static constexpr GeometryType::Id geometryTypeId()
385 {
388 else
389 return GeometryType::Id(~0u);
390 }
391
392 // Select the static LFECache if k >= 0, else the dynamic LFECache
393 using FiniteElementCache = std::conditional_t<(useDynamicOrder),
394 DynamicLagrangeLocalFiniteElementCache<typename GV::ctype,R,dim>,
395 StaticLagrangeLocalFiniteElementCache<geometryTypeId(),typename GV::ctype,R,dim,std::max(k,0)>
396 >;
397
398 // Construct the FiniteElementCache depending on whether the order is dynamic or static
399 static auto makeFiniteElementCache(unsigned int order)
400 {
401 if constexpr (useDynamicOrder)
402 return FiniteElementCache{order};
403 else
404 return FiniteElementCache{};
405 }
406
407public:
408
409 using size_type = std::size_t;
410 using Element = typename GV::template Codim<0>::Entity;
411 using FiniteElement = typename FiniteElementCache::FiniteElementType;
412
414 LagrangeNode() :
415 LagrangeNode(k)
416 {}
417
419 LagrangeNode(unsigned int order) :
420 order_(order),
421 cache_(makeFiniteElementCache(order)),
422 finiteElement_(nullptr),
423 element_(nullptr)
424 {}
425
427 const Element& element() const
428 {
429 return *element_;
430 }
431
436 const FiniteElement& finiteElement() const
437 {
438 return *finiteElement_;
439 }
440
442 void bind(const Element& e)
443 {
444 element_ = &e;
445 finiteElement_ = &(cache_.get(element_->type()));
446 this->setSize(finiteElement_->size());
447 }
448
449protected:
450
451 unsigned int order() const
452 {
453 return order_;
454 }
455
456 // Run-time order, only valid if k<0
457 unsigned int order_;
458
459 FiniteElementCache cache_;
460 const FiniteElement* finiteElement_;
461 const Element* element_;
462};
463
464
465
466namespace BasisFactory {
467
476template<std::size_t k, typename R=double>
478{
479 return [](const auto& gridView) {
480 return LagrangePreBasis<std::decay_t<decltype(gridView)>, k, R>(gridView);
481 };
482}
483
491template<typename R=double>
492auto lagrange(int order)
493{
494 return [=](const auto& gridView) {
495 return LagrangePreBasis<std::decay_t<decltype(gridView)>, -1, R>(gridView, order);
496 };
497}
498
499} // end namespace BasisFactory
500
501
502
531template<typename GV, int k=-1, typename R=double>
533
534
535
536
537
538} // end namespace Functions
539} // end namespace Dune
540
541
542#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEBASIS_HH
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:51
A pre-basis for a PQ-lagrange bases with given order.
Definition: lagrangebasis.hh:65
size_type computeDofsPerCube(std::size_t cubeDim) const
Number of degrees of freedom assigned to a cube (without the ones assigned to its faces!...
Definition: lagrangebasis.hh:342
size_type computeDofsPerSimplex(std::size_t simplexDim) const
Number of degrees of freedom assigned to a simplex (without the ones assigned to its faces!...
Definition: lagrangebasis.hh:336
LagrangeNode< GV, k, R > Node
Template mapping root tree path to type of created tree node.
Definition: lagrangebasis.hh:78
size_type dofsPerSimplex(std::size_t simplexDim) const
Number of degrees of freedom assigned to a simplex (without the ones assigned to its faces!...
Definition: lagrangebasis.hh:314
std::size_t size_type
Type used for indices and size information.
Definition: lagrangebasis.hh:75
void initializeIndices()
Initialize the global indices.
Definition: lagrangebasis.hh:102
LagrangePreBasis(const GridView &gv, unsigned int order)
Constructor for a given grid view object and run-time order.
Definition: lagrangebasis.hh:86
size_type dimension() const
Get the total dimension of the space spanned by this basis.
Definition: lagrangebasis.hh:146
LagrangePreBasis(const GridView &gv)
Constructor for a given grid view object with compile-time order.
Definition: lagrangebasis.hh:81
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: lagrangebasis.hh:132
GV GridView
The grid view that the FE basis is defined on.
Definition: lagrangebasis.hh:72
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: lagrangebasis.hh:126
Node makeNode() const
Create tree node.
Definition: lagrangebasis.hh:140
unsigned int order() const
Polynomial order used in the local Lagrange finite-elements.
Definition: lagrangebasis.hh:302
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition: lagrangebasis.hh:176
size_type dofsPerCube(std::size_t cubeDim) const
Number of degrees of freedom assigned to a cube (without the ones assigned to its faces!...
Definition: lagrangebasis.hh:320
A generic MixIn class for PreBasis.
Definition: leafprebasismixin.hh:36
IdType Id
An integral id representing a GeometryType.
Definition: type.hh:181
Grid view abstract base class.
Definition: gridview.hh:66
Describe position of one degree of freedom.
Definition: localkey.hh:24
constexpr unsigned int index() const noexcept
Return offset within subentity.
Definition: localkey.hh:70
constexpr unsigned int codim() const noexcept
Return codim of associated entity.
Definition: localkey.hh:63
constexpr unsigned int subEntity() const noexcept
Return number of associated subentity.
Definition: localkey.hh:56
Default exception for dummy implementations.
Definition: exceptions.hh:357
Default exception class for range errors.
Definition: exceptions.hh:348
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,...)
Definition: exceptions.hh:314
auto lagrange()
Create a pre-basis factory that can create a Lagrange pre-basis.
Definition: lagrangebasis.hh:477
constexpr GeometryType prism
GeometryType representing a 3D prism.
Definition: type.hh:528
constexpr GeometryType triangle
GeometryType representing a triangle.
Definition: type.hh:504
constexpr GeometryType quadrilateral
GeometryType representing a quadrilateral (a square).
Definition: type.hh:510
constexpr GeometryType hexahedron
GeometryType representing a hexahedron.
Definition: type.hh:534
constexpr GeometryType pyramid
GeometryType representing a 3D pyramid.
Definition: type.hh:522
constexpr GeometryType tetrahedron
GeometryType representing a tetrahedron.
Definition: type.hh:516
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:485
Dune namespace.
Definition: alignedallocator.hh:13
static constexpr T binomial(const T &n, const T &k) noexcept
calculate the binomial coefficient n over k as a constexpr
Definition: math.hh:113
constexpr Base power(Base m, Exponent p)
Power method for integer exponents.
Definition: math.hh:75
STL namespace.
Specialize with 'true' for if the codimension 0 entity of the grid has only one possible geometry typ...
Definition: capabilities.hh:27
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden & Uni Heidelberg  |  generated with Hugo v0.111.3 (Sep 28, 22:46, 2025)