DUNE PDELab (git)

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
14#include <dune/localfunctions/lagrange/equidistantpoints.hh>
15#include <dune/localfunctions/lagrange/lagrangelfecache.hh>
16
17#include <dune/functions/functionspacebases/nodes.hh>
18#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
19#include <dune/functions/functionspacebases/leafprebasismixin.hh>
20
21
22namespace Dune {
23namespace Functions {
24
25// *****************************************************************************
26// This is the reusable part of the LagrangeBasis. It contains
27//
28// LagrangePreBasis
29// LagrangeNode
30//
31// The pre-basis allows to create the others and is the owner of possible shared
32// state. These components do _not_ depend on the global basis and local view
33// and can be used without a global basis.
34// *****************************************************************************
35
36template<typename GV, int k, typename R=double>
37class LagrangeNode;
38
39template<typename GV, int k, typename R=double>
40class LagrangePreBasis;
41
42
43
58template<typename GV, int k, typename R>
60 public LeafPreBasisMixin< LagrangePreBasis<GV,k,R> >
61{
62 static const int dim = GV::dimension;
63 static const bool useDynamicOrder = (k<0);
64
65public:
66
68 using GridView = GV;
69
71 using size_type = std::size_t;
72
74 using Node = LagrangeNode<GV, k, R>;
75
78 : LagrangePreBasis(gv, std::numeric_limits<unsigned int>::max())
79 {}
80
82 LagrangePreBasis(const GridView& gv, unsigned int order) :
83 gridView_(gv), order_(order)
84 {
85 if (!useDynamicOrder && order!=std::numeric_limits<unsigned int>::max())
86 DUNE_THROW(RangeError, "Template argument k has to be -1 when supplying a run-time order!");
87
88 for (int i=0; i<=dim; i++)
89 {
90 dofsPerCube_[i] = computeDofsPerCube(i);
91 dofsPerSimplex_[i] = computeDofsPerSimplex(i);
92 }
93 dofsPerPrism_ = computeDofsPerPrism();
94 dofsPerPyramid_ = computeDofsPerPyramid();
95 }
96
99 {
100 vertexOffset_ = 0;
101 edgeOffset_ = vertexOffset_ + dofsPerCube(0) * ((size_type)gridView_.size(dim));
102
103 if (dim>=2)
104 {
105 triangleOffset_ = edgeOffset_ + dofsPerCube(1) * ((size_type) gridView_.size(dim-1));
106
107 quadrilateralOffset_ = triangleOffset_ + dofsPerSimplex(2) * ((size_type)gridView_.size(Dune::GeometryTypes::triangle));
108 }
109
110 if (dim==3) {
111 tetrahedronOffset_ = quadrilateralOffset_ + dofsPerCube(2) * ((size_type)gridView_.size(Dune::GeometryTypes::quadrilateral));
112
113 prismOffset_ = tetrahedronOffset_ + dofsPerSimplex(3) * ((size_type)gridView_.size(Dune::GeometryTypes::tetrahedron));
114
115 hexahedronOffset_ = prismOffset_ + dofsPerPrism() * ((size_type)gridView_.size(Dune::GeometryTypes::prism));
116
117 pyramidOffset_ = hexahedronOffset_ + dofsPerCube(3) * ((size_type)gridView_.size(Dune::GeometryTypes::hexahedron));
118 }
119 }
120
122 const GridView& gridView() const
123 {
124 return gridView_;
125 }
126
128 void update (const GridView& gv)
129 {
130 gridView_ = gv;
131 }
132
137 {
138 return Node{order_};
139 }
140
143 {
144 switch (dim)
145 {
146 case 1:
147 return dofsPerCube(0) * ((size_type)gridView_.size(dim))
148 + dofsPerCube(1) * ((size_type)gridView_.size(dim-1));
149 case 2:
150 {
151 return dofsPerCube(0) * ((size_type)gridView_.size(dim))
152 + dofsPerCube(1) * ((size_type)gridView_.size(dim-1))
155 }
156 case 3:
157 {
158 return dofsPerCube(0) * ((size_type)gridView_.size(dim))
159 + dofsPerCube(1) * ((size_type)gridView_.size(dim-1))
163 + dofsPerPyramid() * ((size_type)gridView_.size(Dune::GeometryTypes::pyramid))
164 + dofsPerPrism() * ((size_type)gridView_.size(Dune::GeometryTypes::prism))
165 + dofsPerCube(3) * ((size_type)gridView_.size(Dune::GeometryTypes::hexahedron));
166 }
167 }
168 DUNE_THROW(Dune::NotImplemented, "No size method for " << dim << "d grids available yet!");
169 }
170
173 {
174 // That cast to unsigned int is necessary because GV::dimension is an enum,
175 // which is not recognized by the power method as an integer type...
176 return power(order()+1, (unsigned int)GV::dimension);
177 }
178
179 template<typename It>
180 It indices(const Node& node, It it) const
181 {
182 for (size_type i = 0, end = node.finiteElement().size() ; i < end ; ++it, ++i)
183 {
184 Dune::LocalKey localKey = node.finiteElement().localCoefficients().localKey(i);
185 const auto& gridIndexSet = gridView().indexSet();
186 const auto& element = node.element();
187
188 // The dimension of the entity that the current dof is related to
189 auto dofDim = dim - localKey.codim();
190
191 // Test for a vertex dof
192 // The test for k==1 is redundant, but having it here allows the compiler to conclude
193 // at compile-time that the dofDim==0 case is the only one that will ever happen.
194 // This leads to measurable speed-up: see
195 // https://gitlab.dune-project.org/staging/dune-functions/issues/30
196 if (k==1 || dofDim==0) {
197 *it = {{ (size_type)(gridIndexSet.subIndex(element,localKey.subEntity(),dim)) }};
198 continue;
199 }
200
201 if (dofDim==1)
202 { // edge dof
203 if (dim==1) // element dof -- any local numbering is fine
204 {
205 *it = {{ edgeOffset_
206 + dofsPerCube(1) * ((size_type)gridIndexSet.subIndex(element,0,0))
207 + localKey.index() }};
208 continue;
209 }
210 else
211 {
212 const auto refElement
213 = Dune::referenceElement<double,dim>(element.type());
214
215 // We have to reverse the numbering if the local element edge is
216 // not aligned with the global edge.
217 auto v0 = (size_type)gridIndexSet.subIndex(element,refElement.subEntity(localKey.subEntity(),localKey.codim(),0,dim),dim);
218 auto v1 = (size_type)gridIndexSet.subIndex(element,refElement.subEntity(localKey.subEntity(),localKey.codim(),1,dim),dim);
219 bool flip = (v0 > v1);
220 *it = {{ (flip)
221 ? edgeOffset_
222 + dofsPerCube(1)*((size_type)gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()))
223 + (dofsPerCube(1)-1)-localKey.index()
224 : edgeOffset_
225 + dofsPerCube(1)*((size_type)gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()))
226 + localKey.index() }};
227 continue;
228 }
229 }
230
231 if (dofDim==2)
232 {
233 if (dim==2) // element dof -- any local numbering is fine
234 {
235 if (element.type().isTriangle())
236 {
237 *it = {{ triangleOffset_ + dofsPerSimplex(2)*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
238 continue;
239 }
240 else if (element.type().isQuadrilateral())
241 {
242 *it = {{ quadrilateralOffset_ + dofsPerCube(2)*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
243 continue;
244 }
245 else
246 DUNE_THROW(Dune::NotImplemented, "2d elements have to be triangles or quadrilaterals");
247 } else
248 {
249 const auto refElement
250 = Dune::referenceElement<double,dim>(element.type());
251
252 if (order()>3)
253 DUNE_THROW(Dune::NotImplemented, "LagrangeBasis for 3D grids is only implemented if k<=3");
254
255 if (order()==3 and !refElement.type(localKey.subEntity(), localKey.codim()).isTriangle())
256 DUNE_THROW(Dune::NotImplemented, "LagrangeBasis for 3D grids with k==3 is only implemented if the grid is a simplex grid");
257
258 *it = {{ triangleOffset_ + ((size_type)gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim())) }};
259 continue;
260 }
261 }
262
263 if (dofDim==3)
264 {
265 if (dim==3) // element dof -- any local numbering is fine
266 {
267 if (element.type().isTetrahedron())
268 {
269 *it = {{ tetrahedronOffset_ + dofsPerSimplex(3)*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
270 continue;
271 }
272 else if (element.type().isHexahedron())
273 {
274 *it = {{ hexahedronOffset_ + dofsPerCube(3)*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
275 continue;
276 }
277 else if (element.type().isPrism())
278 {
279 *it = {{ prismOffset_ + dofsPerPrism()*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
280 continue;
281 }
282 else if (element.type().isPyramid())
283 {
284 *it = {{ pyramidOffset_ + dofsPerPyramid()*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
285 continue;
286 }
287 else
288 DUNE_THROW(Dune::NotImplemented, "3d elements have to be tetrahedra, hexahedra, prisms, or pyramids");
289 } else
290 DUNE_THROW(Dune::NotImplemented, "Grids of dimension larger than 3 are no supported");
291 }
292 DUNE_THROW(Dune::NotImplemented, "Grid contains elements not supported for the LagrangeBasis");
293 }
294 return it;
295 }
296
298 unsigned int order() const
299 {
300 return (useDynamicOrder) ? order_ : k;
301 }
302
303protected:
304 GridView gridView_;
305
306 // Run-time order, only valid if k<0
307 unsigned int order_;
308
310 size_type dofsPerSimplex(std::size_t simplexDim) const
311 {
312 return useDynamicOrder ? dofsPerSimplex_[simplexDim] : computeDofsPerSimplex(simplexDim);
313 }
314
316 size_type dofsPerCube(std::size_t cubeDim) const
317 {
318 return useDynamicOrder ? dofsPerCube_[cubeDim] : computeDofsPerCube(cubeDim);
319 }
320
321 size_type dofsPerPrism() const
322 {
323 return useDynamicOrder ? dofsPerPrism_ : computeDofsPerPrism();
324 }
325
326 size_type dofsPerPyramid() const
327 {
328 return useDynamicOrder ? dofsPerPyramid_ : computeDofsPerPyramid();
329 }
330
332 size_type computeDofsPerSimplex(std::size_t simplexDim) const
333 {
334 return order() == 0 ? (dim == simplexDim ? 1 : 0) : Dune::binomial(std::size_t(order()-1),simplexDim);
335 }
336
338 size_type computeDofsPerCube(std::size_t cubeDim) const
339 {
340 return order() == 0 ? (dim == cubeDim ? 1 : 0) : Dune::power(order()-1, cubeDim);
341 }
342
343 size_type computeDofsPerPrism() const
344 {
345 return order() == 0 ? (dim == 3 ? 1 : 0) : (order()-1)*(order()-1)*(order()-2)/2;
346 }
347
348 size_type computeDofsPerPyramid() const
349 {
350 return order() == 0 ? (dim == 3 ? 1 : 0) : (order()-2)*(order()-1)*(2*order()-3)/6;
351 }
352
353 // When the order is given at run-time, the following numbers are pre-computed:
354 std::array<size_type,dim+1> dofsPerSimplex_;
355 std::array<size_type,dim+1> dofsPerCube_;
356 size_type dofsPerPrism_;
357 size_type dofsPerPyramid_;
358
359 size_type vertexOffset_;
360 size_type edgeOffset_;
361 size_type triangleOffset_;
362 size_type quadrilateralOffset_;
363 size_type tetrahedronOffset_;
364 size_type pyramidOffset_;
365 size_type prismOffset_;
366 size_type hexahedronOffset_;
367
368};
369
370
371
372template<typename GV, int k, typename R>
373class LagrangeNode :
374 public LeafBasisNode
375{
376 // Stores LocalFiniteElement implementations with run-time order as a function of GeometryType
377 template<typename Domain, typename Range, int dim>
378 class LagrangeRunTimeLFECache
379 {
380 public:
381 using FiniteElementType = LagrangeLocalFiniteElement<EquidistantPointSet,dim,Domain,Range>;
382
383 const FiniteElementType& get(GeometryType type)
384 {
385 auto i = data_.find(type);
386 if (i==data_.end())
387 i = data_.emplace(type,FiniteElementType(type,order_)).first;
388 return (*i).second;
389 }
390
391 std::map<GeometryType, FiniteElementType> data_;
392 unsigned int order_;
393 };
394
395 static constexpr int dim = GV::dimension;
396 static constexpr bool useDynamicOrder = (k<0);
397
398 using FiniteElementCache = std::conditional_t<(useDynamicOrder),
399 LagrangeRunTimeLFECache<typename GV::ctype, R, dim>,
400 LagrangeLocalFiniteElementCache<typename GV::ctype, R, dim, std::max(k,0)>
401 >;
402
403public:
404
405 using size_type = std::size_t;
406 using Element = typename GV::template Codim<0>::Entity;
407 using FiniteElement = typename FiniteElementCache::FiniteElementType;
408
410 LagrangeNode() :
411 finiteElement_(nullptr),
412 element_(nullptr)
413 {}
414
416 LagrangeNode(unsigned int order) :
417 order_(order),
418 finiteElement_(nullptr),
419 element_(nullptr)
420 {
421 // Only the cache for the run-time-order case (i.e., k<0), has the 'order_' member
422 if constexpr (useDynamicOrder)
423 cache_.order_ = order;
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 (useDynamicOrder) ? order_ : k;
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
526template<typename GV, int k=-1, typename R=double>
528
529
530
531
532
533} // end namespace Functions
534} // end namespace Dune
535
536
537#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEBASIS_HH
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:50
A pre-basis for a PQ-lagrange bases with given order.
Definition: lagrangebasis.hh:61
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:338
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:332
LagrangeNode< GV, k, R > Node
Template mapping root tree path to type of created tree node.
Definition: lagrangebasis.hh:74
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:310
std::size_t size_type
Type used for indices and size information.
Definition: lagrangebasis.hh:71
void initializeIndices()
Initialize the global indices.
Definition: lagrangebasis.hh:98
LagrangePreBasis(const GridView &gv, unsigned int order)
Constructor for a given grid view object and run-time order.
Definition: lagrangebasis.hh:82
size_type dimension() const
Get the total dimension of the space spanned by this basis.
Definition: lagrangebasis.hh:142
LagrangePreBasis(const GridView &gv)
Constructor for a given grid view object with compile-time order.
Definition: lagrangebasis.hh:77
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: lagrangebasis.hh:128
GV GridView
The grid view that the FE basis is defined on.
Definition: lagrangebasis.hh:68
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: lagrangebasis.hh:122
Node makeNode() const
Create tree node.
Definition: lagrangebasis.hh:136
unsigned int order() const
Polynomial order used in the local Lagrange finite-elements.
Definition: lagrangebasis.hh:298
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition: lagrangebasis.hh:172
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:316
A generic MixIn class for PreBasis.
Definition: leafprebasismixin.hh:36
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:263
Default exception class for range errors.
Definition: exceptions.hh:254
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
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:484
Convenience header that includes all implementations of Lagrange finite elements.
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:131
constexpr Base power(Base m, Exponent p)
Power method for integer exponents.
Definition: math.hh:75
LocalFiniteElementVariantCache< Impl::ImplementedLagrangeFiniteElements< D, R, dim, order > > LagrangeLocalFiniteElementCache
A cache that stores all available Pk/Qk like local finite elements for the given dimension and order.
Definition: lagrangelfecache.hh:118
constexpr auto get(std::integer_sequence< T, II... >, std::integral_constant< std::size_t, pos >={})
Return the entry at position pos of the given sequence.
Definition: integersequence.hh:22
STL namespace.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)