Dune Core Modules (2.8.0)

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#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEBASIS_HH
4#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEBASIS_HH
5
6#include <type_traits>
8
10#include <dune/localfunctions/lagrange/equidistantpoints.hh>
11#include <dune/localfunctions/lagrange/pqkfactory.hh>
12
13#include <dune/functions/functionspacebases/nodes.hh>
14#include <dune/functions/functionspacebases/flatmultiindex.hh>
15#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
16
17
18namespace Dune {
19namespace Functions {
20
21// *****************************************************************************
22// This is the reusable part of the LagrangeBasis. It contains
23//
24// LagrangePreBasis
25// LagrangeNode
26//
27// The pre-basis allows to create the others and is the owner of possible shared
28// state. These components do _not_ depend on the global basis and local view
29// and can be used without a global basis.
30// *****************************************************************************
31
32template<typename GV, int k, typename R=double>
33class LagrangeNode;
34
35template<typename GV, int k, class MI, typename R=double>
36class LagrangePreBasis;
37
38
39
55template<typename GV, int k, class MI, typename R>
57{
58 static const int dim = GV::dimension;
59 static const bool useDynamicOrder = (k<0);
60
61public:
62
64 using GridView = GV;
65
67 using size_type = std::size_t;
68
70 using Node = LagrangeNode<GV, k, R>;
71
73 using IndexSet = Impl::DefaultNodeIndexSet<LagrangePreBasis>;
74
76 using MultiIndex = MI;
77
80
83 : LagrangePreBasis(gv, std::numeric_limits<unsigned int>::max())
84 {}
85
87 LagrangePreBasis(const GridView& gv, unsigned int order) :
88 gridView_(gv), order_(order)
89 {
90 if (!useDynamicOrder && order!=std::numeric_limits<unsigned int>::max())
91 DUNE_THROW(RangeError, "Template argument k has to be -1 when supplying a run-time order!");
92
93 for (int i=0; i<=dim; i++)
94 {
95 dofsPerCube_[i] = computeDofsPerCube(i);
96 dofsPerSimplex_[i] = computeDofsPerSimplex(i);
97 }
98 dofsPerPrism_ = computeDofsPerPrism();
99 dofsPerPyramid_ = computeDofsPerPyramid();
100 }
101
104 {
105 vertexOffset_ = 0;
106 edgeOffset_ = vertexOffset_ + dofsPerCube(0) * ((size_type)gridView_.size(dim));
107
108 if (dim>=2)
109 {
110 triangleOffset_ = edgeOffset_ + dofsPerCube(1) * ((size_type) gridView_.size(dim-1));
111
112 quadrilateralOffset_ = triangleOffset_ + dofsPerSimplex(2) * ((size_type)gridView_.size(Dune::GeometryTypes::triangle));
113 }
114
115 if (dim==3) {
116 tetrahedronOffset_ = quadrilateralOffset_ + dofsPerCube(2) * ((size_type)gridView_.size(Dune::GeometryTypes::quadrilateral));
117
118 prismOffset_ = tetrahedronOffset_ + dofsPerSimplex(3) * ((size_type)gridView_.size(Dune::GeometryTypes::tetrahedron));
119
120 hexahedronOffset_ = prismOffset_ + dofsPerPrism() * ((size_type)gridView_.size(Dune::GeometryTypes::prism));
121
122 pyramidOffset_ = hexahedronOffset_ + dofsPerCube(3) * ((size_type)gridView_.size(Dune::GeometryTypes::hexahedron));
123 }
124 }
125
127 const GridView& gridView() const
128 {
129 return gridView_;
130 }
131
133 void update (const GridView& gv)
134 {
135 gridView_ = gv;
136 }
137
142 {
143 return Node{order_};
144 }
145
153 [[deprecated("Warning: The IndexSet typedef and the makeIndexSet method are deprecated. "\
154 "As a replacement use the indices() method of the PreBasis directly.")]]
156 {
157 return IndexSet{*this};
158 }
159
162 {
163 switch (dim)
164 {
165 case 1:
166 return dofsPerCube(0) * ((size_type)gridView_.size(dim))
167 + dofsPerCube(1) * ((size_type)gridView_.size(dim-1));
168 case 2:
169 {
170 return dofsPerCube(0) * ((size_type)gridView_.size(dim))
171 + dofsPerCube(1) * ((size_type)gridView_.size(dim-1))
174 }
175 case 3:
176 {
177 return dofsPerCube(0) * ((size_type)gridView_.size(dim))
178 + dofsPerCube(1) * ((size_type)gridView_.size(dim-1))
182 + dofsPerPyramid() * ((size_type)gridView_.size(Dune::GeometryTypes::pyramid))
183 + dofsPerPrism() * ((size_type)gridView_.size(Dune::GeometryTypes::prism))
184 + dofsPerCube(3) * ((size_type)gridView_.size(Dune::GeometryTypes::hexahedron));
185 }
186 }
187 DUNE_THROW(Dune::NotImplemented, "No size method for " << dim << "d grids available yet!");
188 }
189
191 size_type size(const SizePrefix prefix) const
192 {
193 assert(prefix.size() == 0 || prefix.size() == 1);
194 return (prefix.size() == 0) ? size() : 0;
195 }
196
199 {
200 return size();
201 }
202
205 {
206 // That cast to unsigned int is necessary because GV::dimension is an enum,
207 // which is not recognized by the power method as an integer type...
208 return power(order()+1, (unsigned int)GV::dimension);
209 }
210
211 template<typename It>
212 It indices(const Node& node, It it) const
213 {
214 for (size_type i = 0, end = node.finiteElement().size() ; i < end ; ++it, ++i)
215 {
216 Dune::LocalKey localKey = node.finiteElement().localCoefficients().localKey(i);
217 const auto& gridIndexSet = gridView().indexSet();
218 const auto& element = node.element();
219
220 // The dimension of the entity that the current dof is related to
221 auto dofDim = dim - localKey.codim();
222
223 // Test for a vertex dof
224 // The test for k==1 is redundant, but having it here allows the compiler to conclude
225 // at compile-time that the dofDim==0 case is the only one that will ever happen.
226 // This leads to measurable speed-up: see
227 // https://gitlab.dune-project.org/staging/dune-functions/issues/30
228 if (k==1 || dofDim==0) {
229 *it = {{ (size_type)(gridIndexSet.subIndex(element,localKey.subEntity(),dim)) }};
230 continue;
231 }
232
233 if (dofDim==1)
234 { // edge dof
235 if (dim==1) // element dof -- any local numbering is fine
236 {
237 *it = {{ edgeOffset_
238 + dofsPerCube(1) * ((size_type)gridIndexSet.subIndex(element,0,0))
239 + localKey.index() }};
240 continue;
241 }
242 else
243 {
244 const auto refElement
245 = Dune::referenceElement<double,dim>(element.type());
246
247 // we have to reverse the numbering if the local triangle edge is
248 // not aligned with the global edge
249 auto v0 = (size_type)gridIndexSet.subIndex(element,refElement.subEntity(localKey.subEntity(),localKey.codim(),0,dim),dim);
250 auto v1 = (size_type)gridIndexSet.subIndex(element,refElement.subEntity(localKey.subEntity(),localKey.codim(),1,dim),dim);
251 bool flip = (v0 > v1);
252 *it = {{ (flip)
253 ? edgeOffset_
254 + dofsPerCube(1)*((size_type)gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()))
255 + (dofsPerCube(1)-1)-localKey.index()
256 : edgeOffset_
257 + dofsPerCube(1)*((size_type)gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()))
258 + localKey.index() }};
259 continue;
260 }
261 }
262
263 if (dofDim==2)
264 {
265 if (dim==2) // element dof -- any local numbering is fine
266 {
267 if (element.type().isTriangle())
268 {
269 *it = {{ triangleOffset_ + dofsPerSimplex(2)*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
270 continue;
271 }
272 else if (element.type().isQuadrilateral())
273 {
274 *it = {{ quadrilateralOffset_ + dofsPerCube(2)*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
275 continue;
276 }
277 else
278 DUNE_THROW(Dune::NotImplemented, "2d elements have to be triangles or quadrilaterals");
279 } else
280 {
281 const auto refElement
282 = Dune::referenceElement<double,dim>(element.type());
283
284 if (order()>3)
285 DUNE_THROW(Dune::NotImplemented, "LagrangeNodalBasis for 3D grids is only implemented if k<=3");
286
287 if (order()==3 and !refElement.type(localKey.subEntity(), localKey.codim()).isTriangle())
288 DUNE_THROW(Dune::NotImplemented, "LagrangeNodalBasis for 3D grids with k==3 is only implemented if the grid is a simplex grid");
289
290 *it = {{ triangleOffset_ + ((size_type)gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim())) }};
291 continue;
292 }
293 }
294
295 if (dofDim==3)
296 {
297 if (dim==3) // element dof -- any local numbering is fine
298 {
299 if (element.type().isTetrahedron())
300 {
301 *it = {{ tetrahedronOffset_ + dofsPerSimplex(3)*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
302 continue;
303 }
304 else if (element.type().isHexahedron())
305 {
306 *it = {{ hexahedronOffset_ + dofsPerCube(3)*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
307 continue;
308 }
309 else if (element.type().isPrism())
310 {
311 *it = {{ prismOffset_ + dofsPerPrism()*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
312 continue;
313 }
314 else if (element.type().isPyramid())
315 {
316 *it = {{ pyramidOffset_ + dofsPerPyramid()*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
317 continue;
318 }
319 else
320 DUNE_THROW(Dune::NotImplemented, "3d elements have to be tetrahedra, hexahedra, prisms, or pyramids");
321 } else
322 DUNE_THROW(Dune::NotImplemented, "Grids of dimension larger than 3 are no supported");
323 }
324 DUNE_THROW(Dune::NotImplemented, "Grid contains elements not supported for the LagrangeNodalBasis");
325 }
326 return it;
327 }
328
329protected:
330 GridView gridView_;
331
332 unsigned int order() const
333 {
334 return (useDynamicOrder) ? order_ : k;
335 }
336
337 // Run-time order, only valid if k<0
338 const unsigned int order_;
339
341 size_type dofsPerSimplex(std::size_t simplexDim) const
342 {
343 return useDynamicOrder ? dofsPerSimplex_[simplexDim] : computeDofsPerSimplex(simplexDim);
344 }
345
347 size_type dofsPerCube(std::size_t cubeDim) const
348 {
349 return useDynamicOrder ? dofsPerCube_[cubeDim] : computeDofsPerCube(cubeDim);
350 }
351
352 size_type dofsPerPrism() const
353 {
354 return useDynamicOrder ? dofsPerPrism_ : computeDofsPerPrism();
355 }
356
357 size_type dofsPerPyramid() const
358 {
359 return useDynamicOrder ? dofsPerPyramid_ : computeDofsPerPyramid();
360 }
361
363 size_type computeDofsPerSimplex(std::size_t simplexDim) const
364 {
365 return order() == 0 ? (dim == simplexDim ? 1 : 0) : Dune::binomial(std::size_t(order()-1),simplexDim);
366 }
367
369 size_type computeDofsPerCube(std::size_t cubeDim) const
370 {
371 return order() == 0 ? (dim == cubeDim ? 1 : 0) : Dune::power(order()-1, cubeDim);
372 }
373
374 size_type computeDofsPerPrism() const
375 {
376 return order() == 0 ? (dim == 3 ? 1 : 0) : (order()-1)*(order()-1)*(order()-2)/2;
377 }
378
379 size_type computeDofsPerPyramid() const
380 {
381 return order() == 0 ? (dim == 3 ? 1 : 0) : (order()-2)*(order()-1)*(2*order()-3)/6;
382 }
383
384 // When the order is given at run-time, the following numbers are pre-computed:
385 std::array<size_type,dim+1> dofsPerSimplex_;
386 std::array<size_type,dim+1> dofsPerCube_;
387 size_type dofsPerPrism_;
388 size_type dofsPerPyramid_;
389
390 size_type vertexOffset_;
391 size_type edgeOffset_;
392 size_type triangleOffset_;
393 size_type quadrilateralOffset_;
394 size_type tetrahedronOffset_;
395 size_type pyramidOffset_;
396 size_type prismOffset_;
397 size_type hexahedronOffset_;
398
399};
400
401
402
403template<typename GV, int k, typename R>
404class LagrangeNode :
405 public LeafBasisNode
406{
407 // Stores LocalFiniteElement implementations with run-time order as a function of GeometryType
408 template<typename Domain, typename Range, int dim>
409 class LagrangeRunTimeLFECache
410 {
411 public:
412 using FiniteElementType = LagrangeLocalFiniteElement<EquidistantPointSet,dim,Domain,Range>;
413
414 const FiniteElementType& get(GeometryType type)
415 {
416 auto i = data_.find(type);
417 if (i==data_.end())
418 i = data_.emplace(type,FiniteElementType(type,order_)).first;
419 return (*i).second;
420 }
421
422 std::map<GeometryType, FiniteElementType> data_;
423 unsigned int order_;
424 };
425
426 static constexpr int dim = GV::dimension;
427 static constexpr bool useDynamicOrder = (k<0);
428
429 using FiniteElementCache = typename std::conditional<(useDynamicOrder),
430 LagrangeRunTimeLFECache<typename GV::ctype, R, dim>,
431 PQkLocalFiniteElementCache<typename GV::ctype, R, dim, k>
432 >::type;
433
434public:
435
436 using size_type = std::size_t;
437 using Element = typename GV::template Codim<0>::Entity;
438 using FiniteElement = typename FiniteElementCache::FiniteElementType;
439
441 LagrangeNode() :
442 finiteElement_(nullptr),
443 element_(nullptr)
444 {}
445
447 LagrangeNode(unsigned int order) :
448 order_(order),
449 finiteElement_(nullptr),
450 element_(nullptr)
451 {
452 // Only the cache for the run-time-order case (i.e., k<0), has the 'order_' member
453 if constexpr (useDynamicOrder)
454 cache_.order_ = order;
455 }
456
458 const Element& element() const
459 {
460 return *element_;
461 }
462
467 const FiniteElement& finiteElement() const
468 {
469 return *finiteElement_;
470 }
471
473 void bind(const Element& e)
474 {
475 element_ = &e;
476 finiteElement_ = &(cache_.get(element_->type()));
477 this->setSize(finiteElement_->size());
478 }
479
480protected:
481
482 unsigned int order() const
483 {
484 return (useDynamicOrder) ? order_ : k;
485 }
486
487 // Run-time order, only valid if k<0
488 unsigned int order_;
489
490 FiniteElementCache cache_;
491 const FiniteElement* finiteElement_;
492 const Element* element_;
493};
494
495
496
497namespace BasisFactory {
498
499namespace Imp {
500
501template<int k, typename R=double>
502class LagrangePreBasisFactory
503{
504 static const bool useDynamicOrder = (k<0);
505public:
506 static const std::size_t requiredMultiIndexSize = 1;
507
508 // \brief Constructor for factory with compile-time order
509 LagrangePreBasisFactory() : order_(0)
510 {}
511
512 // \brief Constructor for factory with run-time order (template argument k is disregarded)
513 LagrangePreBasisFactory(unsigned int order)
514 : order_(order)
515 {}
516
517 template<class MultiIndex, class GridView>
518 auto makePreBasis(const GridView& gridView) const
519 {
520 return (useDynamicOrder)
521 ? LagrangePreBasis<GridView, k, MultiIndex, R>(gridView, order_)
522 : LagrangePreBasis<GridView, k, MultiIndex, R>(gridView);
523 }
524
525private:
526 unsigned int order_;
527};
528
529} // end namespace BasisFactory::Imp
530
531
532
541template<std::size_t k, typename R=double>
543{
544 return Imp::LagrangePreBasisFactory<k,R>();
545}
546
554template<typename R=double>
555auto lagrange(int order)
556{
557 return Imp::LagrangePreBasisFactory<-1,R>(order);
558}
559
560} // end namespace BasisFactory
561
562
563
587template<typename GV, int k=-1, typename R=double>
589
590
591
592
593
594} // end namespace Functions
595} // end namespace Dune
596
597
598#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEBASIS_HH
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:47
A pre-basis for a PQ-lagrange bases with given order.
Definition: lagrangebasis.hh:57
Node makeNode() const
Create tree node.
Definition: lagrangebasis.hh:141
size_type dimension() const
Get the total dimension of the space spanned by this basis.
Definition: lagrangebasis.hh:198
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:347
Impl::DefaultNodeIndexSet< LagrangePreBasis > IndexSet
Type of created tree node index set.
Definition: lagrangebasis.hh:73
size_type size() const
Same as size(prefix) with empty prefix.
Definition: lagrangebasis.hh:161
LagrangePreBasis(const GridView &gv, unsigned int order)
Constructor for a given grid view object and run-time order.
Definition: lagrangebasis.hh:87
LagrangePreBasis(const GridView &gv)
Constructor for a given grid view object with compile-time order.
Definition: lagrangebasis.hh:82
GV GridView
The grid view that the FE basis is defined on.
Definition: lagrangebasis.hh:64
LagrangeNode< GV, k, R > Node
Template mapping root tree path to type of created tree node.
Definition: lagrangebasis.hh:70
std::size_t size_type
Type used for indices and size information.
Definition: lagrangebasis.hh:67
IndexSet makeIndexSet() const
Create tree node index set.
Definition: lagrangebasis.hh:155
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: lagrangebasis.hh:133
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: lagrangebasis.hh:76
void initializeIndices()
Initialize the global indices.
Definition: lagrangebasis.hh:103
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:341
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:369
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: lagrangebasis.hh:127
size_type size(const SizePrefix prefix) const
Return number of possible values for next position in multi index.
Definition: lagrangebasis.hh:191
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition: lagrangebasis.hh:204
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:363
Describe position of one degree of freedom.
Definition: localkey.hh:21
unsigned int index() const
Return offset within subentity.
Definition: localkey.hh:66
unsigned int codim() const
Return codim of associated entity.
Definition: localkey.hh:60
unsigned int subEntity() const
Return number of associated subentity.
Definition: localkey.hh:54
Default exception for dummy implementations.
Definition: exceptions.hh:261
Default exception class for range errors.
Definition: exceptions.hh:252
A Vector class with statically reserved memory.
Definition: reservedvector.hh:43
size_type size() const
Returns number of elements in the vector.
Definition: reservedvector.hh:184
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
auto lagrange(int order)
Create a pre-basis factory that can create a Lagrange pre-basis with a run-time order.
Definition: lagrangebasis.hh:555
constexpr GeometryType prism
GeometryType representing a 3D prism.
Definition: type.hh:540
constexpr GeometryType triangle
GeometryType representing a triangle.
Definition: type.hh:516
constexpr GeometryType quadrilateral
GeometryType representing a quadrilateral (a square).
Definition: type.hh:522
constexpr GeometryType hexahedron
GeometryType representing a hexahedron.
Definition: type.hh:546
constexpr GeometryType pyramid
GeometryType representing a 3D pyramid.
Definition: type.hh:534
constexpr GeometryType tetrahedron
GeometryType representing a tetrahedron.
Definition: type.hh:528
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Convenience header that includes all implementations of Lagrange finite elements.
Dune namespace.
Definition: alignedallocator.hh:11
static constexpr T binomial(const T &n, const T &k) noexcept
calculate the binomial coefficient n over k as a constexpr
Definition: math.hh:128
constexpr Mantissa power(Mantissa m, Exponent p)
Power method for integer exponents.
Definition: math.hh:73
STL namespace.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 22, 23:30, 2024)