dune-functions  2.8.0
lagrangebasis.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 #ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEBASIS_HH
4 #define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEBASIS_HH
5 
6 #include <type_traits>
7 #include <dune/common/exceptions.hh>
8 
9 #include <dune/localfunctions/lagrange.hh>
10 #include <dune/localfunctions/lagrange/equidistantpoints.hh>
11 #include <dune/localfunctions/lagrange/pqkfactory.hh>
12 
16 
17 
18 namespace Dune {
19 namespace 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 
32 template<typename GV, int k, typename R=double>
33 class LagrangeNode;
34 
35 template<typename GV, int k, class MI, typename R=double>
36 class LagrangePreBasis;
37 
38 
39 
55 template<typename GV, int k, class MI, typename R>
57 {
58  static const int dim = GV::dimension;
59  static const bool useDynamicOrder = (k<0);
60 
61 public:
62 
64  using GridView = GV;
65 
67  using size_type = std::size_t;
68 
71 
73  using IndexSet = Impl::DefaultNodeIndexSet<LagrangePreBasis>;
74 
76  using MultiIndex = MI;
77 
79  using SizePrefix = Dune::ReservedVector<size_type, 1>;
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  {
97  }
100  }
101 
104  {
105  vertexOffset_ = 0;
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 
141  Node makeNode() const
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 
161  size_type size() const
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))
172  + dofsPerSimplex(2) * ((size_type)gridView_.size(Dune::GeometryTypes::triangle))
173  + dofsPerCube(2) * ((size_type)gridView_.size(Dune::GeometryTypes::quadrilateral));
174  }
175  case 3:
176  {
177  return dofsPerCube(0) * ((size_type)gridView_.size(dim))
178  + dofsPerCube(1) * ((size_type)gridView_.size(dim-1))
179  + dofsPerSimplex(2) * ((size_type)gridView_.size(Dune::GeometryTypes::triangle))
180  + dofsPerCube(2) * ((size_type)gridView_.size(Dune::GeometryTypes::quadrilateral))
181  + dofsPerSimplex(3) * ((size_type)gridView_.size(Dune::GeometryTypes::tetrahedron))
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 
329 protected:
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 
353  {
354  return useDynamicOrder ? dofsPerPrism_ : computeDofsPerPrism();
355  }
356 
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 
375  {
376  return order() == 0 ? (dim == 3 ? 1 : 0) : (order()-1)*(order()-1)*(order()-2)/2;
377  }
378 
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_;
389 
398 
399 };
400 
401 
402 
403 template<typename GV, int k, typename R>
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 
434 public:
435 
436  using size_type = std::size_t;
437  using Element = typename GV::template Codim<0>::Entity;
438  using FiniteElement = typename FiniteElementCache::FiniteElementType;
439 
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 
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 
480 protected:
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_;
493 };
494 
495 
496 
497 namespace BasisFactory {
498 
499 namespace Imp {
500 
501 template<int k, typename R=double>
502 class LagrangePreBasisFactory
503 {
504  static const bool useDynamicOrder = (k<0);
505 public:
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)
523  }
524 
525 private:
526  unsigned int order_;
527 };
528 
529 } // end namespace BasisFactory::Imp
530 
531 
532 
541 template<std::size_t k, typename R=double>
542 auto lagrange()
543 {
544  return Imp::LagrangePreBasisFactory<k,R>();
545 }
546 
554 template<typename R=double>
555 auto lagrange(int order)
556 {
557  return Imp::LagrangePreBasisFactory<-1,R>(order);
558 }
559 
560 } // end namespace BasisFactory
561 
562 
563 
587 template<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
auto lagrange()
Create a pre-basis factory that can create a Lagrange pre-basis.
Definition: lagrangebasis.hh:542
auto power(ChildPreBasisFactory &&childPreBasisFactory, const IndexMergingStrategy &ims)
Create a pre-basis factory that can build a PowerPreBasis.
Definition: powerbasis.hh:425
Definition: polynomial.hh:10
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:47
Definition: lagrangebasis.hh:406
LagrangeNode(unsigned int order)
Constructor with a run-time order.
Definition: lagrangebasis.hh:447
const Element & element() const
Return current element, throw if unbound.
Definition: lagrangebasis.hh:458
unsigned int order() const
Definition: lagrangebasis.hh:482
const Element * element_
Definition: lagrangebasis.hh:492
FiniteElementCache cache_
Definition: lagrangebasis.hh:490
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition: lagrangebasis.hh:467
typename FiniteElementCache::FiniteElementType FiniteElement
Definition: lagrangebasis.hh:438
void bind(const Element &e)
Bind to element.
Definition: lagrangebasis.hh:473
typename GV::template Codim< 0 >::Entity Element
Definition: lagrangebasis.hh:437
const FiniteElement * finiteElement_
Definition: lagrangebasis.hh:491
unsigned int order_
Definition: lagrangebasis.hh:488
std::size_t size_type
Definition: lagrangebasis.hh:436
LagrangeNode()
Constructor without order (uses the compile-time value)
Definition: lagrangebasis.hh:441
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
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: lagrangebasis.hh:127
size_type dimension() const
Get the total dimension of the space spanned by this basis.
Definition: lagrangebasis.hh:198
const unsigned int order_
Definition: lagrangebasis.hh:338
size_type vertexOffset_
Definition: lagrangebasis.hh:390
size_type hexahedronOffset_
Definition: lagrangebasis.hh:397
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 quadrilateralOffset_
Definition: lagrangebasis.hh:393
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
size_type triangleOffset_
Definition: lagrangebasis.hh:392
It indices(const Node &node, It it) const
Definition: lagrangebasis.hh:212
LagrangePreBasis(const GridView &gv)
Constructor for a given grid view object with compile-time order.
Definition: lagrangebasis.hh:82
size_type dofsPerPrism_
Definition: lagrangebasis.hh:387
GV GridView
The grid view that the FE basis is defined on.
Definition: lagrangebasis.hh:64
std::size_t size_type
Type used for indices and size information.
Definition: lagrangebasis.hh:67
size_type dofsPerPrism() const
Definition: lagrangebasis.hh:352
IndexSet makeIndexSet() const
Create tree node index set.
Definition: lagrangebasis.hh:155
GridView gridView_
Definition: lagrangebasis.hh:330
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: lagrangebasis.hh:133
size_type pyramidOffset_
Definition: lagrangebasis.hh:395
size_type dofsPerPyramid_
Definition: lagrangebasis.hh:388
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: lagrangebasis.hh:76
size_type edgeOffset_
Definition: lagrangebasis.hh:391
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
size_type tetrahedronOffset_
Definition: lagrangebasis.hh:394
size_type computeDofsPerPyramid() const
Definition: lagrangebasis.hh:379
unsigned int order() const
Definition: lagrangebasis.hh:332
size_type dofsPerPyramid() const
Definition: lagrangebasis.hh:357
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 prismOffset_
Definition: lagrangebasis.hh:396
size_type computeDofsPerPrism() const
Definition: lagrangebasis.hh:374
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
std::array< size_type, dim+1 > dofsPerSimplex_
Definition: lagrangebasis.hh:385
Dune::ReservedVector< size_type, 1 > SizePrefix
Type used for prefixes handed to the size() method.
Definition: lagrangebasis.hh:79
std::array< size_type, dim+1 > dofsPerCube_
Definition: lagrangebasis.hh:386
void setSize(const size_type size)
Definition: nodes.hh:162
Definition: nodes.hh:184