DUNE PDELab (2.7)

pkfem.hh
1// -*- tab-width: 4; indent-tabs-mode: nil -*-
2#ifndef DUNE_PDELAB_FINITEELEMENTMAP_PKFEM_HH
3#define DUNE_PDELAB_FINITEELEMENTMAP_PKFEM_HH
4
5#include <algorithm>
6#include <array>
7
8#include <dune/common/deprecated.hh>
9#include <dune/common/exceptions.hh>
10#include <dune/geometry/type.hh>
11#include <dune/localfunctions/lagrange/pk.hh>
12
13#include "finiteelementmap.hh"
14
15namespace Dune {
16 namespace PDELab {
17
18 namespace fem {
19
20 template<typename GV, typename D, typename R, unsigned int k, unsigned int d>
21 class PkLocalFiniteElementMapBase;
22
23
24 // ********************************************************************************
25 // 1D version
26 // ********************************************************************************
27
28 template<typename GV, typename D, typename R, unsigned int k>
29 class PkLocalFiniteElementMapBase<GV,D,R,k,1>
30 : public SimpleLocalFiniteElementMap<Dune::LagrangeSimplexLocalFiniteElement<D,R,1,k>,1>
31 {
32
33 public:
34
35 PkLocalFiniteElementMapBase(const GV& gv)
36 {}
37
38 static constexpr bool fixedSize()
39 {
40 return true;
41 }
42
43 static constexpr bool hasDOFs(int codim)
44 {
45 switch (codim)
46 {
47 case 1: // vertex
48 return k > 0;
49 case 0: // line
50 return k != 1;
51 default:
52 assert(k >= 0 and k <= 1);
53 }
54 return false;
55 }
56
57 static constexpr std::size_t size(GeometryType gt)
58 {
59 if (gt == GeometryTypes::vertex)
60 return k > 0 ? 1 : 0;
61 if (gt == GeometryTypes::line)
62 return k > 0 ? k - 1 : 1;
63 return 0;
64 }
65
66 static constexpr std::size_t maxLocalSize()
67 {
68 return k + 1;
69 }
70
71 };
72
73
74 // ********************************************************************************
75 // 2D version
76 // ********************************************************************************
77
78 template<typename GV, typename D, typename R, unsigned int k>
79 class PkLocalFiniteElementMapBase<GV,D,R,k,2>
80 : public LocalFiniteElementMapInterface<LocalFiniteElementMapTraits<
81 Dune::LagrangeSimplexLocalFiniteElement<D,R,2,k>
82 >,
83 PkLocalFiniteElementMapBase<GV,D,R,k,2>
84 >
85 {
86
87 typedef Dune::LagrangeSimplexLocalFiniteElement<D,R,2,k> FE;
88
89 public:
90
92 typedef LocalFiniteElementMapTraits<FE> Traits;
93
94 PkLocalFiniteElementMapBase(const GV& gv)
95 : _gv(gv)
96 {
97 // generate permutations
98 std::array<unsigned int, 3> p = {0,1,2};
99 for (int i = 0; i < 6; ++i)
100 {
101 _variant[i] = FE(p);
102 std::next_permutation(p.begin(), p.end());
103 }
104 }
105
107 template<typename Entity>
108 const typename Traits::FiniteElementType& find (const Entity& e) const
109 {
110
111 if (!e.type().isSimplex())
112 DUNE_THROW(InvalidGeometryType,"PkLocalFiniteElementMap only works for simplices!");
113
114 const typename GV::IndexSet& is = _gv.indexSet();
115 unsigned int n0 = is.subIndex(e,0,2);
116 unsigned int n1 = is.subIndex(e,1,2);
117 unsigned int n2 = is.subIndex(e,2,2);
118 // compress first number to [0,2]
119 unsigned int n0_compressed = (n0 > n1) + (n0 > n2);
120 // translate to permutation index
121 return _variant[2 * n0_compressed + (n1 > n2)];
122 }
123
124 static constexpr bool fixedSize()
125 {
126 return true;
127 }
128
129 static constexpr bool hasDOFs(int codim)
130 {
131 switch (codim)
132 {
133 case 2: // vertex
134 return k > 0;
135 case 1: // line
136 return k > 1;
137 case 0: // triangle
138 return k > 2 || k == 0;
139 default:
140 assert(false && "Invalid codim specified!");
141 }
142 return false;
143 }
144
145 static constexpr std::size_t size(GeometryType gt)
146 {
147 if (gt == GeometryTypes::vertex)
148 return k > 0 ? 1 : 0;
149 if (gt == GeometryTypes::line)
150 return k > 1 ? k - 1 : 0;
151 if (gt == GeometryTypes::triangle)
152 return k > 2 ? (k-2)*(k-1)/2 : (k == 0);
153 return 0;
154 }
155
156 static constexpr std::size_t maxLocalSize()
157 {
158 return (k+1)*(k+2)/2;
159 }
160
161 private:
162 std::array<FE,6> _variant;
163 GV _gv;
164
165 };
166
167
168 // ********************************************************************************
169 // 3D version
170 // ********************************************************************************
171
172 template<typename GV, typename D, typename R, unsigned int k>
173 class PkLocalFiniteElementMapBase<GV,D,R,k,3>
174 : public LocalFiniteElementMapInterface<LocalFiniteElementMapTraits<
175 Dune::LagrangeSimplexLocalFiniteElement<D,R,3,k>
176 >,
177 PkLocalFiniteElementMapBase<GV,D,R,k,3>
178 >
179 {
180
181 typedef Dune::LagrangeSimplexLocalFiniteElement<D,R,3,k> FE;
182
183 public:
184
186 typedef LocalFiniteElementMapTraits<FE> Traits;
187
188 PkLocalFiniteElementMapBase(const GV& gv)
189 : _gv(gv)
190 {
191 std::fill(_perm_index.begin(),_perm_index.end(),0);
192
193 // create all variants by iterating over all permutations
194 unsigned int n = 0;
195 std::array<unsigned int, 4> vertexmap = {0, 1, 2, 3};
196 do {
197 _variant[n] = FE(vertexmap);
198 _perm_index[compressPerm(vertexmap)] = n++;
199 }
200 while(std::next_permutation(vertexmap.begin(), vertexmap.end()));
201 }
202
204 template<typename Entity>
205 const typename Traits::FiniteElementType& find (const Entity& e) const
206 {
207
208 if (!e.type().isSimplex())
209 DUNE_THROW(InvalidGeometryType,"PkLocalFiniteElementMap only works for simplices!");
210
211 // get the vertex indices
212 const typename GV::IndexSet& is = _gv.indexSet();
213 std::array<unsigned int, 4> vertexmap;
214 for (unsigned int i = 0; i < 4; ++i)
215 vertexmap[i] = is.subIndex(e,i,3);
216
217 // reduce the indices to the interval 0..3
218 for (unsigned int i = 0; i < 4; ++i)
219 {
220 int min_index = -1;
221 for (unsigned int j = 0; j < 4; ++j)
222 if ((min_index < 0 || vertexmap[j] < vertexmap[min_index]) && vertexmap[j] >= i)
223 min_index = j;
224 assert(min_index >= 0);
225 vertexmap[min_index] = i;
226 }
227 return _variant[_perm_index[compressPerm(vertexmap)]];
228 }
229
230 static constexpr bool fixedSize()
231 {
232 return true;
233 }
234
235 static constexpr bool hasDOFs(int codim)
236 {
237 switch (codim)
238 {
239 case 3: // vertex
240 return k > 0;
241 case 2: // line
242 return k > 1;
243 case 1: // triangle
244 return k > 2;
245 case 0: // tetrahedron
246 return k == 0 || k > 3;
247 default:
248 assert(false && "Invalid codim specified!");
249 }
250 return false;
251 }
252
253 static constexpr std::size_t size(GeometryType gt)
254 {
255 if (gt == GeometryTypes::vertex)
256 return k > 0 ? 1 : 0;
257 if (gt == GeometryTypes::line)
258 return k > 1 ? k - 1 : 0;
259 if (gt == GeometryTypes::triangle)
260 return k > 2 ? (k-2)*(k-1)/2 : 0;
261 if (gt == GeometryTypes::tetrahedron)
262 return k == 0 ? 1 : (k-3)*(k-2)*(k-1)/6;
263 return 0;
264 }
265
266 static constexpr std::size_t maxLocalSize()
267 {
268 return (k+1)*(k+2)*(k+3)/6;
269 }
270
271 private:
272
273 unsigned int compressPerm(const std::array<unsigned int, 4> vertexmap) const
274 {
275 return vertexmap[0] + (vertexmap[1]<<2) + (vertexmap[2]<<4) + (vertexmap[3]<<6);
276 }
277
278 std::array<FE,24> _variant;
279 std::array<unsigned int,256> _perm_index;
280 GV _gv;
281
282 };
283
284 } // namespace fem
285
286
287 template<typename GV, typename D, typename R, unsigned int k>
288 class PkLocalFiniteElementMap
289 : public fem::PkLocalFiniteElementMapBase<GV,D,R,k,GV::dimension>
290 {
291
292 public:
293
295 static constexpr int dimension = GV::dimension;
296
297 PkLocalFiniteElementMap(const GV& gv)
298 : fem::PkLocalFiniteElementMapBase<GV,D,R,k,GV::dimension>(gv)
299 {}
300
301 };
302
303
304 }
305}
306
307#endif // DUNE_PDELAB_FINITEELEMENTMAP_PKFEM_HH
Dune::LagrangeSimplexLocalFiniteElement
Lagrange finite element for simplices with arbitrary compile-time dimension and polynomial order.
Definition: lagrangesimplex.hh:836
deprecated.hh
Definition of the DUNE_DEPRECATED macro for the case that config.h is not available.
exceptions.hh
A few common exception classes.
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune::FloatCmp::gt
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:156
Dune::GeometryTypes::line
constexpr GeometryType line
GeometryType representing a line.
Definition: type.hh:804
Dune::GeometryTypes::triangle
constexpr GeometryType triangle
GeometryType representing a triangle.
Definition: type.hh:810
Dune::GeometryTypes::tetrahedron
constexpr GeometryType tetrahedron
GeometryType representing a tetrahedron.
Definition: type.hh:822
Dune::GeometryTypes::vertex
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:798
Dune
Dune namespace.
Definition: alignedallocator.hh:14
type.hh
A unique label for each type of element that can occur in a grid.
fixedSize
std::size_t fixedSize
The number of data items per index if it is fixed, 0 otherwise.
Definition: variablesizecommunicator.hh:245
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 15, 22:36, 2024)