dune-functions  2.5.0
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 // PQkNodeFactory
27 // PQkNodeIndexSet
28 // PQkNode
29 //
30 // The factory 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 
44 
45 
60 template<typename GV, int k, class MI>
61 class PQkNodeFactory
62 {
63  static const int dim = GV::dimension;
64 
65 public:
66 
68  using GridView = GV;
69 
71  using size_type = std::size_t;
72 
73 private:
74 
75  template<typename, int, class, class>
76  friend class PQkNodeIndexSet;
77 
78  // Precompute the number of dofs per entity type
79  const static size_type dofsPerVertex =
80  k == 0 ? (dim == 0 ? 1 : 0) : 1;
81  const static size_type dofsPerEdge =
82  k == 0 ? (dim == 1 ? 1 : 0) : k-1;
83  const static size_type dofsPerTriangle =
84  k == 0 ? (dim == 2 ? 1 : 0) : (k-1)*(k-2)/2;
85  const static size_type dofsPerQuad =
86  k == 0 ? (dim == 2 ? 1 : 0) : (k-1)*(k-1);
87  const static size_type dofsPerTetrahedron =
88  k == 0 ? (dim == 3 ? 1 : 0) : (k-3)*(k-2)*(k-1)/6;
89  const static size_type dofsPerPrism =
90  k == 0 ? (dim == 3 ? 1 : 0) : (k-1)*(k-1)*(k-2)/2;
91  const static size_type dofsPerHexahedron =
92  k == 0 ? (dim == 3 ? 1 : 0) : (k-1)*(k-1)*(k-1);
93  const static size_type dofsPerPyramid =
94  k == 0 ? (dim == 3 ? 1 : 0) : (k-2)*(k-1)*(2*k-3)/6;
95 
96 public:
97 
99  template<class TP>
101 
103  template<class TP>
105 
107  using MultiIndex = MI;
108 
110  using SizePrefix = Dune::ReservedVector<size_type, 1>;
111 
114  gridView_(gv)
115  {}
116 
119  {
120  vertexOffset_ = 0;
121  edgeOffset_ = vertexOffset_ + dofsPerVertex * ((size_type)gridView_.size(dim));
122  triangleOffset_ = edgeOffset_ + dofsPerEdge * ((size_type) gridView_.size(dim-1));
123 
124  GeometryType triangle;
125  triangle.makeTriangle();
126  quadrilateralOffset_ = triangleOffset_ + dofsPerTriangle * ((size_type)gridView_.size(triangle));
127 
128  Dune::GeometryType quadrilateral;
129  quadrilateral.makeQuadrilateral();
130  if (dim==3) {
131  tetrahedronOffset_ = quadrilateralOffset_ + dofsPerQuad * ((size_type)gridView_.size(quadrilateral));
132 
133  GeometryType tetrahedron;
134  tetrahedron.makeSimplex(3);
135  prismOffset_ = tetrahedronOffset_ + dofsPerTetrahedron * ((size_type)gridView_.size(tetrahedron));
136 
137  GeometryType prism;
138  prism.makePrism();
139  hexahedronOffset_ = prismOffset_ + dofsPerPrism * ((size_type)gridView_.size(prism));
140 
141  GeometryType hexahedron;
142  hexahedron.makeCube(3);
143  pyramidOffset_ = hexahedronOffset_ + dofsPerHexahedron * ((size_type)gridView_.size(hexahedron));
144  }
145  }
146 
148  const GridView& gridView() const
149  {
150  return gridView_;
151  }
152 
154  void update (const GridView& gv)
155  {
156  gridView_ = gv;
157  }
158 
169  template<class TP>
170  Node<TP> node(const TP& tp) const
171  {
172  return Node<TP>{tp};
173  }
174 
184  template<class TP>
186  {
187  return IndexSet<TP>{*this};
188  }
189 
191  size_type size() const
192  {
193  switch (dim)
194  {
195  case 1:
196  return dofsPerVertex * ((size_type)gridView_.size(dim))
197  + dofsPerEdge*((size_type)gridView_.size(dim-1));
198  case 2:
199  {
200  GeometryType triangle, quad;
201  triangle.makeTriangle();
202  quad.makeQuadrilateral();
203 
204  return dofsPerVertex * ((size_type)gridView_.size(dim))
205  + dofsPerEdge * ((size_type)gridView_.size(dim-1))
206  + dofsPerTriangle * ((size_type)gridView_.size(triangle))
207  + dofsPerQuad * ((size_type)gridView_.size(quad));
208  }
209  case 3:
210  {
211  GeometryType triangle, quad, tetrahedron, pyramid, prism, hexahedron;
212  triangle.makeTriangle();
213  quad.makeQuadrilateral();
214  tetrahedron.makeTetrahedron();
215  pyramid.makePyramid();
216  prism.makePrism();
217  hexahedron.makeCube(3);
218  return dofsPerVertex * ((size_type)gridView_.size(dim))
219  + dofsPerEdge * ((size_type)gridView_.size(dim-1))
220  + dofsPerTriangle * ((size_type)gridView_.size(triangle))
221  + dofsPerQuad * ((size_type)gridView_.size(quad))
222  + dofsPerTetrahedron * ((size_type)gridView_.size(tetrahedron))
223  + dofsPerPyramid * ((size_type)gridView_.size(pyramid))
224  + dofsPerPrism * ((size_type)gridView_.size(prism))
225  + dofsPerHexahedron * ((size_type)gridView_.size(hexahedron));
226  }
227  }
228  DUNE_THROW(Dune::NotImplemented, "No size method for " << dim << "d grids available yet!");
229  }
230 
232  size_type size(const SizePrefix prefix) const
233  {
234  assert(prefix.size() == 0 || prefix.size() == 1);
235  return (prefix.size() == 0) ? size() : 0;
236  }
237 
240  {
241  return size();
242  }
243 
246  {
248  }
249 
250 protected:
252 
261 
262 };
263 
264 
265 
266 template<typename GV, int k, typename TP>
267 class PQkNode :
268  public LeafBasisNode<std::size_t, TP>
269 {
270  static const int dim = GV::dimension;
271  static const int maxSize = StaticPower<(k+1),GV::dimension>::power;
272 
273  using Base = LeafBasisNode<std::size_t,TP>;
274  using FiniteElementCache = typename Dune::PQkLocalFiniteElementCache<typename GV::ctype, double, dim, k>;
275 
276 public:
277 
278  using size_type = std::size_t;
279  using TreePath = TP;
280  using Element = typename GV::template Codim<0>::Entity;
281  using FiniteElement = typename FiniteElementCache::FiniteElementType;
282 
283  PQkNode(const TreePath& treePath) :
284  Base(treePath),
285  finiteElement_(nullptr),
286  element_(nullptr)
287  {}
288 
290  const Element& element() const
291  {
292  return *element_;
293  }
294 
300  {
301  return *finiteElement_;
302  }
303 
305  void bind(const Element& e)
306  {
307  element_ = &e;
308  finiteElement_ = &(cache_.get(element_->type()));
309  this->setSize(finiteElement_->size());
310  }
311 
312 protected:
313 
314  FiniteElementCache cache_;
317 };
318 
319 
320 
321 template<typename GV, int k, class MI, class TP>
322 class PQkNodeIndexSet
323 {
324  enum {dim = GV::dimension};
325 
326 public:
327 
328  using size_type = std::size_t;
329 
331  using MultiIndex = MI;
332 
334 
335  using Node = typename NodeFactory::template Node<TP>;
336 
337  PQkNodeIndexSet(const NodeFactory& nodeFactory) :
338  nodeFactory_(&nodeFactory),
339  node_(nullptr)
340  {}
341 
347  void bind(const Node& node)
348  {
349  node_ = &node;
350  }
351 
354  void unbind()
355  {
356  node_ = nullptr;
357  }
358 
361  size_type size() const
362  {
363  assert(node_ != nullptr);
364  return node_->finiteElement().size();
365  }
366 
369  {
370  assert(node_ != nullptr);
371  Dune::LocalKey localKey = node_->finiteElement().localCoefficients().localKey(i);
372  const auto& gridIndexSet = nodeFactory_->gridView().indexSet();
373  const auto& element = node_->element();
374 
375  // The dimension of the entity that the current dof is related to
376  auto dofDim = dim - localKey.codim();
377 
378  if (dofDim==0) { // vertex dof
379  return {{ (size_type)(gridIndexSet.subIndex(element,localKey.subEntity(),dim)) }};
380  }
381 
382  if (dofDim==1)
383  { // edge dof
384  if (dim==1) // element dof -- any local numbering is fine
385  return {{ nodeFactory_->edgeOffset_
386  + nodeFactory_->dofsPerEdge * ((size_type)gridIndexSet.subIndex(element,0,0))
387  + localKey.index() }};
388  else
389  {
390  const Dune::ReferenceElement<double,dim>& refElement
391  = Dune::ReferenceElements<double,dim>::general(element.type());
392 
393  // we have to reverse the numbering if the local triangle edge is
394  // not aligned with the global edge
395  auto v0 = (size_type)gridIndexSet.subIndex(element,refElement.subEntity(localKey.subEntity(),localKey.codim(),0,dim),dim);
396  auto v1 = (size_type)gridIndexSet.subIndex(element,refElement.subEntity(localKey.subEntity(),localKey.codim(),1,dim),dim);
397  bool flip = (v0 > v1);
398  return {{ (flip)
399  ? nodeFactory_->edgeOffset_
400  + nodeFactory_->dofsPerEdge*((size_type)gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()))
401  + (nodeFactory_->dofsPerEdge-1)-localKey.index()
402  : nodeFactory_->edgeOffset_
403  + nodeFactory_->dofsPerEdge*((size_type)gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()))
404  + localKey.index() }};
405  }
406  }
407 
408  if (dofDim==2)
409  {
410  if (dim==2) // element dof -- any local numbering is fine
411  {
412  if (element.type().isTriangle())
413  {
414  const int interiorLagrangeNodesPerTriangle = (k-1)*(k-2)/2;
415  return {{ nodeFactory_->triangleOffset_ + interiorLagrangeNodesPerTriangle*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
416  }
417  else if (element.type().isQuadrilateral())
418  {
419  const int interiorLagrangeNodesPerQuadrilateral = (k-1)*(k-1);
420  return {{ nodeFactory_->quadrilateralOffset_ + interiorLagrangeNodesPerQuadrilateral*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
421  }
422  else
423  DUNE_THROW(Dune::NotImplemented, "2d elements have to be triangles or quadrilaterals");
424  } else
425  {
426  const Dune::ReferenceElement<double,dim>& refElement
427  = Dune::ReferenceElements<double,dim>::general(element.type());
428 
429  if (k>3)
430  DUNE_THROW(Dune::NotImplemented, "PQkNodalBasis for 3D grids is only implemented if k<=3");
431 
432  if (k==3 and !refElement.type(localKey.subEntity(), localKey.codim()).isTriangle())
433  DUNE_THROW(Dune::NotImplemented, "PQkNodalBasis for 3D grids with k==3 is only implemented if the grid is a simplex grid");
434 
435  return {{ nodeFactory_->triangleOffset_ + ((size_type)gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim())) }};
436  }
437  }
438 
439  if (dofDim==3)
440  {
441  if (dim==3) // element dof -- any local numbering is fine
442  {
443  if (element.type().isTetrahedron())
444  return {{ nodeFactory_->tetrahedronOffset_ + NodeFactory::dofsPerTetrahedron*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
445  else if (element.type().isHexahedron())
446  return {{ nodeFactory_->hexahedronOffset_ + NodeFactory::dofsPerHexahedron*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
447  else if (element.type().isPrism())
448  return {{ nodeFactory_->prismOffset_ + NodeFactory::dofsPerPrism*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
449  else if (element.type().isPyramid())
450  return {{ nodeFactory_->pyramidOffset_ + NodeFactory::dofsPerPyramid*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
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 
459 protected:
461 
462  const Node* node_;
463 };
464 
465 
466 
467 namespace BasisBuilder {
468 
469 namespace Imp {
470 
471 template<std::size_t k>
473 {
474  static const std::size_t requiredMultiIndexSize=1;
475 
476  template<class MultiIndex, class GridView>
477  auto build(const GridView& gridView)
479  {
480  return {gridView};
481  }
482 };
483 
484 } // end namespace BasisBuilder::Imp
485 
486 
487 
495 template<std::size_t k>
497 {
498  return{};
499 }
500 
501 } // end namespace BasisBuilder
502 
503 
504 
505 // *****************************************************************************
506 // This is the actual global basis implementation based on the reusable parts.
507 // *****************************************************************************
508 
524 template<typename GV, int k>
526 
527 
528 
529 } // end namespace Functions
530 } // end namespace Dune
531 
532 
533 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_PQKNODALBASIS_HH
Definition: pqknodalbasis.hh:39
size_type prismOffset_
Definition: pqknodalbasis.hh:259
Imp::PQkNodeFactoryBuilder< k > pq()
Create a factory builder that can build a PQkNodeFactory.
Definition: pqknodalbasis.hh:496
std::size_t size_type
Definition: nodes.hh:127
size_type tetrahedronOffset_
Definition: pqknodalbasis.hh:257
const Element * element_
Definition: pqknodalbasis.hh:316
Global basis for given node factory.
Definition: defaultglobalbasis.hh:42
Definition: polynomial.hh:7
Definition: nodes.hh:189
size_type edgeOffset_
Definition: pqknodalbasis.hh:254
PQkNodeFactory(const GridView &gv)
Constructor for a given grid view object.
Definition: pqknodalbasis.hh:113
void unbind()
Unbind the view.
Definition: pqknodalbasis.hh:354
typename GV::template Codim< 0 >::Entity Element
Definition: pqknodalbasis.hh:280
size_type pyramidOffset_
Definition: pqknodalbasis.hh:258
size_type size() const
Size of subtree rooted in this node (element-local)
Definition: pqknodalbasis.hh:361
const Element & element() const
Return current element, throw if unbound.
Definition: pqknodalbasis.hh:290
void bind(const Node &node)
Bind the view to a grid element.
Definition: pqknodalbasis.hh:347
FiniteElementCache cache_
Definition: pqknodalbasis.hh:314
size_type quadrilateralOffset_
Definition: pqknodalbasis.hh:256
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition: pqknodalbasis.hh:299
Definition: pqknodalbasis.hh:36
const NodeFactory * nodeFactory_
Definition: pqknodalbasis.hh:460
size_type size() const
Same as size(prefix) with empty prefix.
Definition: pqknodalbasis.hh:191
typename NodeFactory::template Node< TP > Node
Definition: pqknodalbasis.hh:335
typename FiniteElementCache::FiniteElementType FiniteElement
Definition: pqknodalbasis.hh:281
size_type dimension() const
Get the total dimension of the space spanned by this basis.
Definition: pqknodalbasis.hh:239
std::size_t size_type
Type used for indices and size information.
Definition: pqknodalbasis.hh:71
Dune::ReservedVector< size_type, 1 > SizePrefix
Type used for prefixes handed to the size() method.
Definition: pqknodalbasis.hh:110
TP TreePath
Definition: nodes.hh:126
size_type hexahedronOffset_
Definition: pqknodalbasis.hh:260
const Node * node_
Definition: pqknodalbasis.hh:462
MultiIndex index(size_type i) const
Maps from subtree index set [0..size-1] to a globally unique multi index in global basis...
Definition: pqknodalbasis.hh:368
PQkNodeIndexSet(const NodeFactory &nodeFactory)
Definition: pqknodalbasis.hh:337
A factory for PQ-lagrange bases with given order.
Definition: pqknodalbasis.hh:42
IndexSet< TP > indexSet() const
Create tree node index set with given root tree path.
Definition: pqknodalbasis.hh:185
size_type vertexOffset_
Definition: pqknodalbasis.hh:253
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: pqknodalbasis.hh:148
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition: pqknodalbasis.hh:245
Node< TP > node(const TP &tp) const
Create tree node with given root tree path.
Definition: pqknodalbasis.hh:170
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: pqknodalbasis.hh:154
size_type size(const SizePrefix prefix) const
Return number of possible values for next position in multi index.
Definition: pqknodalbasis.hh:232
Imp::PowerNodeFactoryBuilder< k, IndexMergingStrategy, SubFactoryTag > power(SubFactoryTag &&tag, const IndexMergingStrategy &ims)
Create a factory builder that can build a PowerNodeFactory.
Definition: powerbasis.hh:424
const FiniteElement * finiteElement_
Definition: pqknodalbasis.hh:315
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: pqknodalbasis.hh:331
PQkNode(const TreePath &treePath)
Definition: pqknodalbasis.hh:283
auto build(const GridView &gridView) -> PQkNodeFactory< GridView, k, MultiIndex >
Definition: pqknodalbasis.hh:477
PQMultiIndex MultiIndex
Type used for global numbering of the basis vectors.
Definition: pqknodalbasis.hh:107
size_type triangleOffset_
Definition: pqknodalbasis.hh:255
GridView gridView_
Definition: pqknodalbasis.hh:251
GV GridView
The grid view that the FE basis is defined on.
Definition: pqknodalbasis.hh:68
void bind(const Element &e)
Bind to element.
Definition: pqknodalbasis.hh:305
std::size_t size_type
Definition: pqknodalbasis.hh:328
void initializeIndices()
Initialize the global indices.
Definition: pqknodalbasis.hh:118