DUNE PDELab (git)

lagrangedgbasis.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_LAGRANGEDGBASIS_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEDGBASIS_HH
9
10#include <array>
12#include <dune/common/math.hh>
13
14#include <dune/functions/functionspacebases/nodes.hh>
15#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
16#include <dune/functions/functionspacebases/lagrangebasis.hh>
17#include <dune/functions/functionspacebases/leafprebasismixin.hh>
18
19
20
21
22namespace Dune {
23namespace Functions {
24
25
26
27// *****************************************************************************
28// This is the reusable part of the basis. It contains
29//
30// LagrangeDGPreBasis
31// LagrangeDGNode
32//
33// The pre-basis allows to create the others and is the owner of possible shared
34// state. These components do _not_ depend on the global basis and local view
35// and can be used without a global basis.
36// *****************************************************************************
37
38template<typename GV, int k>
39using LagrangeDGNode = LagrangeNode<GV, k>;
40
41template<typename GV, int k>
42class LagrangeDGPreBasis :
43 public LeafPreBasisMixin< LagrangeDGPreBasis<GV,k> >
44{
45 static const int dim = GV::dimension;
46
47public:
48
50 using GridView = GV;
51 using size_type = std::size_t;
52
53
54 // Precompute the number of dofs per entity type
55 const static int dofsPerEdge = k+1;
56 const static int dofsPerTriangle = (k+1)*(k+2)/2;
57 const static int dofsPerQuad = (k+1)*(k+1);
58 const static int dofsPerTetrahedron = (k+1)*(k+2)*(k+3)/6;
59 const static int dofsPerPrism = (k+1)*(k+1)*(k+2)/2;
60 const static int dofsPerHexahedron = (k+1)*(k+1)*(k+1);
61 const static int dofsPerPyramid = (k+1)*(k+2)*(2*k+3)/6;
62
63
64 using Node = LagrangeDGNode<GV, k>;
65
67 LagrangeDGPreBasis(const GridView& gv) :
68 gridView_(gv)
69 {}
70
71
72 void initializeIndices()
73 {
74 switch (dim)
75 {
76 case 1:
77 {
78 break;
79 }
80 case 2:
81 {
82 quadrilateralOffset_ = dofsPerTriangle * gridView_.size(Dune::GeometryTypes::triangle);
83 break;
84 }
85 case 3:
86 {
87 prismOffset_ = dofsPerTetrahedron * gridView_.size(Dune::GeometryTypes::tetrahedron);
88
89 hexahedronOffset_ = prismOffset_ + dofsPerPrism * gridView_.size(Dune::GeometryTypes::prism);
90
91 pyramidOffset_ = hexahedronOffset_ + dofsPerHexahedron * gridView_.size(Dune::GeometryTypes::hexahedron);
92 break;
93 }
94 }
95 }
96
99 const GridView& gridView() const
100 {
101 return gridView_;
102 }
103
104 void update(const GridView& gv)
105 {
106 gridView_ = gv;
107 }
108
112 Node makeNode() const
113 {
114 return Node{};
115 }
116
117 size_type dimension() const
118 {
119 switch (dim)
120 {
121 case 1:
122 return dofsPerEdge*gridView_.size(0);
123 case 2:
124 {
125 return dofsPerTriangle*gridView_.size(Dune::GeometryTypes::triangle) + dofsPerQuad*gridView_.size(Dune::GeometryTypes::quadrilateral);
126 }
127 case 3:
128 {
129 return dofsPerTetrahedron*gridView_.size(Dune::GeometryTypes::tetrahedron)
130 + dofsPerPyramid*gridView_.size(Dune::GeometryTypes::pyramid)
131 + dofsPerPrism*gridView_.size(Dune::GeometryTypes::prism)
132 + dofsPerHexahedron*gridView_.size(Dune::GeometryTypes::hexahedron);
133 }
134 }
135 DUNE_THROW(Dune::NotImplemented, "No size method for " << dim << "d grids available yet!");
136 }
137
138 size_type maxNodeSize() const
139 {
140 return Dune::power(k+1, int(GV::dimension));
141 }
142
143 template<typename It>
144 It indices(const Node& node, It it) const
145 {
146 const auto& gridIndexSet = gridView().indexSet();
147 const auto& element = node.element();
148
149 size_type offset = 0;
150 if constexpr (dim==1)
151 offset = dofsPerEdge*gridIndexSet.subIndex(element,0,0);
152 else if constexpr (dim==2)
153 {
154 if (element.type().isTriangle())
155 offset = dofsPerTriangle*gridIndexSet.subIndex(element,0,0);
156 else if (element.type().isQuadrilateral())
157 offset = quadrilateralOffset_ + dofsPerQuad*gridIndexSet.subIndex(element,0,0);
158 else
159 DUNE_THROW(Dune::NotImplemented, "2d elements have to be triangles or quadrilaterals");
160 }
161 else if constexpr (dim==3)
162 {
163 if (element.type().isTetrahedron())
164 offset = dofsPerTetrahedron*gridIndexSet.subIndex(element,0,0);
165 else if (element.type().isPrism())
166 offset = prismOffset_ + dofsPerPrism*gridIndexSet.subIndex(element,0,0);
167 else if (element.type().isHexahedron())
168 offset = hexahedronOffset_ + dofsPerHexahedron*gridIndexSet.subIndex(element,0,0);
169 else if (element.type().isPyramid())
170 offset = pyramidOffset_ + dofsPerPyramid*gridIndexSet.subIndex(element,0,0);
171 else
172 DUNE_THROW(Dune::NotImplemented, "3d elements have to be tetrahedrons, prisms, hexahedrons or pyramids");
173 }
174 else
175 DUNE_THROW(Dune::NotImplemented, "No index method for " << dim << "d grids available yet!");
176 for (size_type i = 0, end = node.size() ; i < end ; ++i, ++it)
177 *it = {offset + i};
178 return it;
179 }
180
182 unsigned int order() const
183 {
184 return k;
185 }
186
187protected:
188 GridView gridView_;
189
190 size_t quadrilateralOffset_;
191 size_t pyramidOffset_;
192 size_t prismOffset_;
193 size_t hexahedronOffset_;
194};
195
196
197
198namespace BasisFactory {
199
207template<std::size_t k>
209{
210 return [](const auto& gridView) {
211 return LagrangeDGPreBasis<std::decay_t<decltype(gridView)>, k>(gridView);
212 };
213}
214
215} // end namespace BasisFactory
216
217
218
219// *****************************************************************************
220// This is the actual global basis implementation based on the reusable parts.
221// *****************************************************************************
222
230template<typename GV, int k>
232
233
234
235} // end namespace Functions
236} // end namespace Dune
237
238
239#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEDGBASIS_HH
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:50
Default exception for dummy implementations.
Definition: exceptions.hh:263
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
auto lagrangeDG()
Create a pre-basis factory that can create a LagrangeDG pre-basis.
Definition: lagrangedgbasis.hh:208
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
Some useful basic math stuff.
Dune namespace.
Definition: alignedallocator.hh:13
constexpr Base power(Base m, Exponent p)
Power method for integer exponents.
Definition: math.hh:75
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)