dune-pdelab  2.5-dev
pkfem.hh
Go to the documentation of this file.
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 
15 namespace Dune {
16  namespace PDELab {
17 
18  namespace fem {
19 
20  template<typename GV, typename D, typename R, unsigned int k, unsigned int d>
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::PkLocalFiniteElement<D,R,1,k>,1>
31  {
32 
33  public:
34 
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::PkLocalFiniteElement<D,R,2,k>
82  >,
83  PkLocalFiniteElementMapBase<GV,D,R,k,2>
84  >
85  {
86 
87  typedef Dune::PkLocalFiniteElement<D,R,2,k> FE;
88 
89  public:
90 
93 
95  : _gv(gv)
96  {
97  // generate permutations
98  unsigned int p[3] = {0,1,2};
99  for (int i = 0; i < 6; ++i)
100  {
101  _variant[i] = FE(p);
102  std::next_permutation(p,p+3);
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::PkLocalFiniteElement<D,R,3,k>
176  >,
177  PkLocalFiniteElementMapBase<GV,D,R,k,3>
178  >
179  {
180 
181  typedef Dune::PkLocalFiniteElement<D,R,3,k> FE;
182 
183  public:
184 
187 
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  unsigned int vertexmap[4];
196  for(vertexmap[0] = 0; vertexmap[0] < 4; ++vertexmap[0])
197  {
198  for(vertexmap[1] = 0; vertexmap[1] < 4; ++vertexmap[1])
199  {
200  if (vertexmap[0] == vertexmap[1])
201  continue;
202  for(vertexmap[2] = 0; vertexmap[2] < 4; ++vertexmap[2])
203  {
204  if (vertexmap[0] == vertexmap[2] ||
205  vertexmap[1] == vertexmap[2])
206  continue;
207  vertexmap[3] = 6 - vertexmap[0] - vertexmap[1] - vertexmap[2];
208  _variant[n] = FE(vertexmap);
209  _perm_index[compressPerm(vertexmap)] = n++;
210  }
211  }
212  }
213  }
214 
216  template<typename Entity>
217  const typename Traits::FiniteElementType& find (const Entity& e) const
218  {
219 
220  if (!e.type().isSimplex())
221  DUNE_THROW(InvalidGeometryType,"PkLocalFiniteElementMap only works for simplices!");
222 
223  // get the vertex indices
224  const typename GV::IndexSet& is = _gv.indexSet();
225  unsigned int vertexmap[4];
226  for (unsigned int i = 0; i < 4; ++i)
227  vertexmap[i] = is.subIndex(e,i,3);
228 
229  // reduce the indices to the interval 0..3
230  for (unsigned int i = 0; i < 4; ++i)
231  {
232  int min_index = -1;
233  for (unsigned int j = 0; j < 4; ++j)
234  if ((min_index < 0 || vertexmap[j] < vertexmap[min_index]) && vertexmap[j] >= i)
235  min_index = j;
236  assert(min_index >= 0);
237  vertexmap[min_index] = i;
238  }
239  return _variant[_perm_index[compressPerm(vertexmap)]];
240  }
241 
242  static constexpr bool fixedSize()
243  {
244  return true;
245  }
246 
247  static constexpr bool hasDOFs(int codim)
248  {
249  switch (codim)
250  {
251  case 3: // vertex
252  return k > 0;
253  case 2: // line
254  return k > 1;
255  case 1: // triangle
256  return k > 2;
257  case 0: // tetrahedron
258  return k == 0 || k > 3;
259  default:
260  assert(false && "Invalid codim specified!");
261  }
262  return false;
263  }
264 
265  static constexpr std::size_t size(GeometryType gt)
266  {
267  if (gt == GeometryTypes::vertex)
268  return k > 0 ? 1 : 0;
269  if (gt == GeometryTypes::line)
270  return k > 1 ? k - 1 : 0;
271  if (gt == GeometryTypes::triangle)
272  return k > 2 ? (k-2)*(k-1)/2 : 0;
273  if (gt == GeometryTypes::tetrahedron)
274  return k == 0 ? 1 : (k-3)*(k-2)*(k-1)/6;
275  return 0;
276  }
277 
278  static constexpr std::size_t maxLocalSize()
279  {
280  return (k+1)*(k+2)*(k+3)/6;
281  }
282 
283  private:
284 
285  unsigned int compressPerm(const unsigned int vertexmap[4]) const
286  {
287  return vertexmap[0] + (vertexmap[1]<<2) + (vertexmap[2]<<4) + (vertexmap[3]<<6);
288  }
289 
290  std::array<FE,24> _variant;
291  std::array<unsigned int,256> _perm_index;
292  GV _gv;
293 
294  };
295 
296  } // namespace fem
297 
298 
299  template<typename GV, typename D, typename R, unsigned int k>
301  : public fem::PkLocalFiniteElementMapBase<GV,D,R,k,GV::dimension>
302  {
303 
304  public:
305 
307  static constexpr int dimension = GV::dimension;
308 
310  : fem::PkLocalFiniteElementMapBase<GV,D,R,k,GV::dimension>(gv)
311  {}
312 
313  };
314 
315 
316  }
317 }
318 
319 #endif // DUNE_PDELAB_FINITEELEMENTMAP_PKFEM_HH
static constexpr std::size_t maxLocalSize()
Definition: pkfem.hh:66
const Entity & e
Definition: localfunctionspace.hh:120
LocalFiniteElementMapTraits< FE > Traits
export type of the signature
Definition: pkfem.hh:186
FiniteElementMap exception raised when trying to obtain a finite element for an unsupported GeometryT...
Definition: finiteelementmap.hh:23
const Traits::FiniteElementType & find(const Entity &e) const
get local basis functions for entity
Definition: pkfem.hh:217
static constexpr bool fixedSize()
Definition: pkfem.hh:38
static constexpr bool hasDOFs(int codim)
Definition: pkfem.hh:247
collect types exported by a finite element map
Definition: finiteelementmap.hh:38
static constexpr std::size_t maxLocalSize()
Definition: pkfem.hh:278
PkLocalFiniteElementMap(const GV &gv)
Definition: pkfem.hh:309
const Traits::FiniteElementType & find(const Entity &e) const
get local basis functions for entity
Definition: pkfem.hh:108
T FiniteElementType
Type of finite element from local functions.
Definition: finiteelementmap.hh:30
static constexpr std::size_t size(GeometryType gt)
Definition: pkfem.hh:145
PkLocalFiniteElementMapBase(const GV &gv)
Definition: pkfem.hh:94
static constexpr bool hasDOFs(int codim)
Definition: pkfem.hh:129
static constexpr std::size_t size(GeometryType gt)
Definition: pkfem.hh:57
PkLocalFiniteElementMapBase(const GV &gv)
Definition: pkfem.hh:188
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
static constexpr std::size_t size(GeometryType gt)
Definition: pkfem.hh:265
static constexpr std::size_t maxLocalSize()
Definition: pkfem.hh:156
LocalFiniteElementMapTraits< FE > Traits
export type of the signature
Definition: pkfem.hh:92
static constexpr bool fixedSize()
Definition: pkfem.hh:242
PkLocalFiniteElementMapBase(const GV &gv)
Definition: pkfem.hh:35
static constexpr bool fixedSize()
Definition: pkfem.hh:124
const P & p
Definition: constraints.hh:147
simple implementation where all entities have the same finite element
Definition: finiteelementmap.hh:111
static constexpr bool hasDOFs(int codim)
Definition: pkfem.hh:43
interface for a finite element map
Definition: finiteelementmap.hh:42