dune-functions  2.6-dev
pqknodalbasis.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_PQKNODALBASIS_HH
4 #define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_PQKNODALBASIS_HH
5 
6 #include <array>
7 #include <dune/common/exceptions.hh>
8 
9 #include <dune/localfunctions/lagrange/pqkfactory.hh>
10 
11 #include <dune/typetree/leafnode.hh>
12 
16 
17 
18 namespace Dune {
19 namespace Functions {
20 
21 
22 
23 // *****************************************************************************
24 // This is the reusable part of the basis. It contains
25 //
26 // PQkPreBasis
27 // PQkNodeIndexSet
28 // PQkNode
29 //
30 // The pre-basis allows to create the others and is the owner of possible shared
31 // state. These three components do _not_ depend on the global basis or index
32 // set and can be used without a global basis.
33 // *****************************************************************************
34 
35 template<typename GV, int k, typename TP>
36 class PQkNode;
37 
38 template<typename GV, int k, class MI, class TP>
40 
41 template<typename GV, int k, class MI>
43 
58 template<typename GV, int k, class MI>
59 class PQkPreBasis
60 {
61  static const int dim = GV::dimension;
62 
63 public:
64 
66  using GridView = GV;
67 
69  using size_type = std::size_t;
70 
71 private:
72 
73  template<typename, int, class, class>
74  friend class PQkNodeIndexSet;
75 
76  // Precompute the number of dofs per entity type
77  const static size_type dofsPerVertex =
78  k == 0 ? (dim == 0 ? 1 : 0) : 1;
79  const static size_type dofsPerEdge =
80  k == 0 ? (dim == 1 ? 1 : 0) : k-1;
81  const static size_type dofsPerTriangle =
82  k == 0 ? (dim == 2 ? 1 : 0) : (k-1)*(k-2)/2;
83  const static size_type dofsPerQuad =
84  k == 0 ? (dim == 2 ? 1 : 0) : (k-1)*(k-1);
85  const static size_type dofsPerTetrahedron =
86  k == 0 ? (dim == 3 ? 1 : 0) : (k-3)*(k-2)*(k-1)/6;
87  const static size_type dofsPerPrism =
88  k == 0 ? (dim == 3 ? 1 : 0) : (k-1)*(k-1)*(k-2)/2;
89  const static size_type dofsPerHexahedron =
90  k == 0 ? (dim == 3 ? 1 : 0) : (k-1)*(k-1)*(k-1);
91  const static size_type dofsPerPyramid =
92  k == 0 ? (dim == 3 ? 1 : 0) : (k-2)*(k-1)*(2*k-3)/6;
93 
94 public:
95 
97  template<class TP>
99 
101  template<class TP>
103 
105  using MultiIndex = MI;
106 
108  using SizePrefix = Dune::ReservedVector<size_type, 1>;
109 
111  PQkPreBasis(const GridView& gv) :
112  gridView_(gv)
113  {}
114 
117  {
118  vertexOffset_ = 0;
119  edgeOffset_ = vertexOffset_ + dofsPerVertex * ((size_type)gridView_.size(dim));
120  triangleOffset_ = edgeOffset_ + dofsPerEdge * ((size_type) gridView_.size(dim-1));
121 
122  quadrilateralOffset_ = triangleOffset_ + dofsPerTriangle * ((size_type)gridView_.size(Dune::GeometryTypes::triangle));
123 
124  if (dim==3) {
125  tetrahedronOffset_ = quadrilateralOffset_ + dofsPerQuad * ((size_type)gridView_.size(Dune::GeometryTypes::quadrilateral));
126 
127  prismOffset_ = tetrahedronOffset_ + dofsPerTetrahedron * ((size_type)gridView_.size(Dune::GeometryTypes::tetrahedron));
128 
129  hexahedronOffset_ = prismOffset_ + dofsPerPrism * ((size_type)gridView_.size(Dune::GeometryTypes::prism));
130 
131  pyramidOffset_ = hexahedronOffset_ + dofsPerHexahedron * ((size_type)gridView_.size(Dune::GeometryTypes::hexahedron));
132  }
133  }
134 
136  const GridView& gridView() const
137  {
138  return gridView_;
139  }
140 
142  void update (const GridView& gv)
143  {
144  gridView_ = gv;
145  }
146 
157  template<class TP>
158  Node<TP> node(const TP& tp) const
159  {
160  return Node<TP>{tp};
161  }
162 
172  template<class TP>
174  {
175  return IndexSet<TP>{*this};
176  }
177 
179  size_type size() const
180  {
181  switch (dim)
182  {
183  case 1:
184  return dofsPerVertex * ((size_type)gridView_.size(dim))
185  + dofsPerEdge*((size_type)gridView_.size(dim-1));
186  case 2:
187  {
188  return dofsPerVertex * ((size_type)gridView_.size(dim))
189  + dofsPerEdge * ((size_type)gridView_.size(dim-1))
190  + dofsPerTriangle * ((size_type)gridView_.size(Dune::GeometryTypes::triangle))
191  + dofsPerQuad * ((size_type)gridView_.size(Dune::GeometryTypes::quadrilateral));
192  }
193  case 3:
194  {
195  return dofsPerVertex * ((size_type)gridView_.size(dim))
196  + dofsPerEdge * ((size_type)gridView_.size(dim-1))
197  + dofsPerTriangle * ((size_type)gridView_.size(Dune::GeometryTypes::triangle))
198  + dofsPerQuad * ((size_type)gridView_.size(Dune::GeometryTypes::quadrilateral))
199  + dofsPerTetrahedron * ((size_type)gridView_.size(Dune::GeometryTypes::tetrahedron))
200  + dofsPerPyramid * ((size_type)gridView_.size(Dune::GeometryTypes::pyramid))
201  + dofsPerPrism * ((size_type)gridView_.size(Dune::GeometryTypes::prism))
202  + dofsPerHexahedron * ((size_type)gridView_.size(Dune::GeometryTypes::hexahedron));
203  }
204  }
205  DUNE_THROW(Dune::NotImplemented, "No size method for " << dim << "d grids available yet!");
206  }
207 
209  size_type size(const SizePrefix prefix) const
210  {
211  assert(prefix.size() == 0 || prefix.size() == 1);
212  return (prefix.size() == 0) ? size() : 0;
213  }
214 
217  {
218  return size();
219  }
220 
223  {
225  }
226 
227 protected:
229 
238 
239 };
240 
241 
242 
243 template<typename GV, int k, typename TP>
244 class PQkNode :
245  public LeafBasisNode<std::size_t, TP>
246 {
247  static const int dim = GV::dimension;
248  static const int maxSize = StaticPower<(k+1),GV::dimension>::power;
249 
250  using Base = LeafBasisNode<std::size_t,TP>;
251  using FiniteElementCache = typename Dune::PQkLocalFiniteElementCache<typename GV::ctype, double, dim, k>;
252 
253 public:
254 
255  using size_type = std::size_t;
256  using TreePath = TP;
257  using Element = typename GV::template Codim<0>::Entity;
258  using FiniteElement = typename FiniteElementCache::FiniteElementType;
259 
260  PQkNode(const TreePath& treePath) :
261  Base(treePath),
262  finiteElement_(nullptr),
263  element_(nullptr)
264  {}
265 
267  const Element& element() const
268  {
269  return *element_;
270  }
271 
277  {
278  return *finiteElement_;
279  }
280 
282  void bind(const Element& e)
283  {
284  element_ = &e;
285  finiteElement_ = &(cache_.get(element_->type()));
286  this->setSize(finiteElement_->size());
287  }
288 
289 protected:
290 
291  FiniteElementCache cache_;
294 };
295 
296 
297 
298 template<typename GV, int k, class MI, class TP>
299 class PQkNodeIndexSet
300 {
301  enum {dim = GV::dimension};
302 
303 public:
304 
305  using size_type = std::size_t;
306 
308  using MultiIndex = MI;
309 
311 
312  using Node = typename PreBasis::template Node<TP>;
313 
314  PQkNodeIndexSet(const PreBasis& preBasis) :
315  preBasis_(&preBasis),
316  node_(nullptr)
317  {}
318 
324  void bind(const Node& node)
325  {
326  node_ = &node;
327  }
328 
331  void unbind()
332  {
333  node_ = nullptr;
334  }
335 
338  size_type size() const
339  {
340  assert(node_ != nullptr);
341  return node_->finiteElement().size();
342  }
343 
345  template<typename It>
346  It indices(It it) const
347  {
348  assert(node_ != nullptr);
349  for (size_type i = 0, end = node_->finiteElement().size() ; i < end ; ++it, ++i)
350  {
351  Dune::LocalKey localKey = node_->finiteElement().localCoefficients().localKey(i);
352  const auto& gridIndexSet = preBasis_->gridView().indexSet();
353  const auto& element = node_->element();
354 
355  // The dimension of the entity that the current dof is related to
356  auto dofDim = dim - localKey.codim();
357 
358  if (dofDim==0) { // vertex dof
359  *it = {{ (size_type)(gridIndexSet.subIndex(element,localKey.subEntity(),dim)) }};
360  continue;
361  }
362 
363  if (dofDim==1)
364  { // edge dof
365  if (dim==1) // element dof -- any local numbering is fine
366  {
367  *it = {{ preBasis_->edgeOffset_
368  + preBasis_->dofsPerEdge * ((size_type)gridIndexSet.subIndex(element,0,0))
369  + localKey.index() }};
370  continue;
371  }
372  else
373  {
374  const auto refElement
375  = Dune::referenceElement<double,dim>(element.type());
376 
377  // we have to reverse the numbering if the local triangle edge is
378  // not aligned with the global edge
379  auto v0 = (size_type)gridIndexSet.subIndex(element,refElement.subEntity(localKey.subEntity(),localKey.codim(),0,dim),dim);
380  auto v1 = (size_type)gridIndexSet.subIndex(element,refElement.subEntity(localKey.subEntity(),localKey.codim(),1,dim),dim);
381  bool flip = (v0 > v1);
382  *it = {{ (flip)
383  ? preBasis_->edgeOffset_
384  + preBasis_->dofsPerEdge*((size_type)gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()))
385  + (preBasis_->dofsPerEdge-1)-localKey.index()
386  : preBasis_->edgeOffset_
387  + preBasis_->dofsPerEdge*((size_type)gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()))
388  + localKey.index() }};
389  continue;
390  }
391  }
392 
393  if (dofDim==2)
394  {
395  if (dim==2) // element dof -- any local numbering is fine
396  {
397  if (element.type().isTriangle())
398  {
399  const int interiorLagrangeNodesPerTriangle = (k-1)*(k-2)/2;
400  *it = {{ preBasis_->triangleOffset_ + interiorLagrangeNodesPerTriangle*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
401  continue;
402  }
403  else if (element.type().isQuadrilateral())
404  {
405  const int interiorLagrangeNodesPerQuadrilateral = (k-1)*(k-1);
406  *it = {{ preBasis_->quadrilateralOffset_ + interiorLagrangeNodesPerQuadrilateral*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
407  continue;
408  }
409  else
410  DUNE_THROW(Dune::NotImplemented, "2d elements have to be triangles or quadrilaterals");
411  } else
412  {
413  const auto refElement
414  = Dune::referenceElement<double,dim>(element.type());
415 
416  if (k>3)
417  DUNE_THROW(Dune::NotImplemented, "PQkNodalBasis for 3D grids is only implemented if k<=3");
418 
419  if (k==3 and !refElement.type(localKey.subEntity(), localKey.codim()).isTriangle())
420  DUNE_THROW(Dune::NotImplemented, "PQkNodalBasis for 3D grids with k==3 is only implemented if the grid is a simplex grid");
421 
422  *it = {{ preBasis_->triangleOffset_ + ((size_type)gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim())) }};
423  continue;
424  }
425  }
426 
427  if (dofDim==3)
428  {
429  if (dim==3) // element dof -- any local numbering is fine
430  {
431  if (element.type().isTetrahedron())
432  {
433  *it = {{ preBasis_->tetrahedronOffset_ + PreBasis::dofsPerTetrahedron*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
434  continue;
435  }
436  else if (element.type().isHexahedron())
437  {
438  *it = {{ preBasis_->hexahedronOffset_ + PreBasis::dofsPerHexahedron*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
439  continue;
440  }
441  else if (element.type().isPrism())
442  {
443  *it = {{ preBasis_->prismOffset_ + PreBasis::dofsPerPrism*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
444  continue;
445  }
446  else if (element.type().isPyramid())
447  {
448  *it = {{ preBasis_->pyramidOffset_ + PreBasis::dofsPerPyramid*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
449  continue;
450  }
451  else
452  DUNE_THROW(Dune::NotImplemented, "3d elements have to be tetrahedra, hexahedra, prisms, or pyramids");
453  } else
454  DUNE_THROW(Dune::NotImplemented, "Grids of dimension larger than 3 are no supported");
455  }
456  DUNE_THROW(Dune::NotImplemented, "Grid contains elements not supported for the PQkNodalBasis");
457  }
458  return it;
459  }
460 
461 protected:
463 
464  const Node* node_;
465 };
466 
467 
468 
469 namespace BasisBuilder {
470 
471 namespace Imp {
472 
473 template<std::size_t k>
474 class PQkPreBasisFactory
475 {
476 public:
477  static const std::size_t requiredMultiIndexSize = 1;
478 
479  template<class MultiIndex, class GridView>
480  auto makePreBasis(const GridView& gridView) const
481  {
483  }
484 
485 };
486 
487 } // end namespace BasisBuilder::Imp
488 
489 
490 
498 template<std::size_t k>
499 auto pq()
500 {
501  return Imp::PQkPreBasisFactory<k>();
502 }
503 
504 } // end namespace BasisBuilder
505 
506 
507 
508 // *****************************************************************************
509 // This is the actual global basis implementation based on the reusable parts.
510 // *****************************************************************************
511 
527 template<typename GV, int k>
529 
530 
531 
532 } // end namespace Functions
533 } // end namespace Dune
534 
535 
536 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_PQKNODALBASIS_HH
Dune::ReservedVector< size_type, 1 > SizePrefix
Type used for prefixes handed to the size() method.
Definition: pqknodalbasis.hh:108
const Element & element() const
Return current element, throw if unbound.
Definition: pqknodalbasis.hh:267
const FiniteElement * finiteElement_
Definition: pqknodalbasis.hh:292
TP TreePath
Definition: nodes.hh:126
IndexSet< TP > indexSet() const
Create tree node index set with given root tree path.
Definition: pqknodalbasis.hh:173
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition: pqknodalbasis.hh:276
Definition: pqknodalbasis.hh:39
std::size_t size_type
Definition: nodes.hh:127
GridView gridView_
Definition: pqknodalbasis.hh:228
size_type vertexOffset_
Definition: pqknodalbasis.hh:230
PQkNodeIndexSet(const PreBasis &preBasis)
Definition: pqknodalbasis.hh:314
size_type triangleOffset_
Definition: pqknodalbasis.hh:232
auto pq()
Create a pre-basis factory that can create a PQ_k pre-basis.
Definition: pqknodalbasis.hh:499
It indices(It it) const
Maps from subtree index set [0..size-1] to a globally unique multi index in global basis...
Definition: pqknodalbasis.hh:346
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: pqknodalbasis.hh:308
typename FiniteElementCache::FiniteElementType FiniteElement
Definition: pqknodalbasis.hh:258
Definition: pqknodalbasis.hh:36
size_type size() const
Size of subtree rooted in this node (element-local)
Definition: pqknodalbasis.hh:338
typename GV::template Codim< 0 >::Entity Element
Definition: pqknodalbasis.hh:257
FiniteElementCache cache_
Definition: pqknodalbasis.hh:291
const PreBasis * preBasis_
Definition: pqknodalbasis.hh:462
void bind(const Node &node)
Bind the view to a grid element.
Definition: pqknodalbasis.hh:324
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:42
Definition: nodes.hh:189
size_type hexahedronOffset_
Definition: pqknodalbasis.hh:237
void initializeIndices()
Initialize the global indices.
Definition: pqknodalbasis.hh:116
std::size_t size_type
Definition: pqknodalbasis.hh:305
size_type size() const
Same as size(prefix) with empty prefix.
Definition: pqknodalbasis.hh:179
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: pqknodalbasis.hh:136
size_type prismOffset_
Definition: pqknodalbasis.hh:236
auto power(ChildPreBasisFactory &&childPreBasisFactory, const IndexMergingStrategy &ims)
Create a pre-basis factory that can build a PowerPreBasis.
Definition: powerbasis.hh:493
size_type size(const SizePrefix prefix) const
Return number of possible values for next position in multi index.
Definition: pqknodalbasis.hh:209
size_type pyramidOffset_
Definition: pqknodalbasis.hh:235
PQMultiIndex MultiIndex
Type used for global numbering of the basis vectors.
Definition: pqknodalbasis.hh:105
const Node * node_
Definition: pqknodalbasis.hh:464
PQkNode(const TreePath &treePath)
Definition: pqknodalbasis.hh:260
size_type dimension() const
Get the total dimension of the space spanned by this basis.
Definition: pqknodalbasis.hh:216
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: pqknodalbasis.hh:142
std::size_t size_type
Type used for indices and size information.
Definition: pqknodalbasis.hh:69
size_type tetrahedronOffset_
Definition: pqknodalbasis.hh:234
size_type edgeOffset_
Definition: pqknodalbasis.hh:231
Definition: polynomial.hh:7
PQkPreBasis(const GridView &gv)
Constructor for a given grid view object.
Definition: pqknodalbasis.hh:111
A pre-basis for PQ-lagrange bases with given order.
Definition: pqknodalbasis.hh:42
typename PreBasis::template Node< TP > Node
Definition: pqknodalbasis.hh:312
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition: pqknodalbasis.hh:222
const Element * element_
Definition: pqknodalbasis.hh:293
Node< TP > node(const TP &tp) const
Create tree node with given root tree path.
Definition: pqknodalbasis.hh:158
GV GridView
The grid view that the FE basis is defined on.
Definition: pqknodalbasis.hh:66
void bind(const Element &e)
Bind to element.
Definition: pqknodalbasis.hh:282
size_type quadrilateralOffset_
Definition: pqknodalbasis.hh:233
void unbind()
Unbind the view.
Definition: pqknodalbasis.hh:331