dune-pdelab  2.7-git
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::LagrangeSimplexLocalFiniteElement<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::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 
93 
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 
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  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>
289  : public fem::PkLocalFiniteElementMapBase<GV,D,R,k,GV::dimension>
290  {
291 
292  public:
293 
295  static constexpr int dimension = GV::dimension;
296 
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::PDELab::SimpleLocalFiniteElementMap
simple implementation where all entities have the same finite element
Definition: finiteelementmap.hh:98
p
const P & p
Definition: constraints.hh:148
Dune::PDELab::PkLocalFiniteElementMap::PkLocalFiniteElementMap
PkLocalFiniteElementMap(const GV &gv)
Definition: pkfem.hh:297
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 2 >::hasDOFs
static constexpr bool hasDOFs(int codim)
Definition: pkfem.hh:129
Dune::PDELab::fem::PkLocalFiniteElementMapBase
Definition: pkfem.hh:21
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 3 >::maxLocalSize
static constexpr std::size_t maxLocalSize()
Definition: pkfem.hh:266
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 3 >::fixedSize
static constexpr bool fixedSize()
Definition: pkfem.hh:230
Dune::PDELab::LocalFiniteElementMapInterface
interface for a finite element map
Definition: finiteelementmap.hh:42
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 3 >::find
const Traits::FiniteElementType & find(const Entity &e) const
get local basis functions for entity
Definition: pkfem.hh:205
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 3 >::size
static constexpr std::size_t size(GeometryType gt)
Definition: pkfem.hh:253
Dune
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 1 >::fixedSize
static constexpr bool fixedSize()
Definition: pkfem.hh:38
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 2 >::fixedSize
static constexpr bool fixedSize()
Definition: pkfem.hh:124
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 2 >::find
const Traits::FiniteElementType & find(const Entity &e) const
get local basis functions for entity
Definition: pkfem.hh:108
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 1 >::size
static constexpr std::size_t size(GeometryType gt)
Definition: pkfem.hh:57
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 1 >::PkLocalFiniteElementMapBase
PkLocalFiniteElementMapBase(const GV &gv)
Definition: pkfem.hh:35
Dune::PDELab::PkLocalFiniteElementMap::dimension
static constexpr int dimension
The dimension of the finite elements returned by this map.
Definition: pkfem.hh:295
e
const Entity & e
Definition: localfunctionspace.hh:121
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 1 >::maxLocalSize
static constexpr std::size_t maxLocalSize()
Definition: pkfem.hh:66
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 2 >::Traits
LocalFiniteElementMapTraits< FE > Traits
export type of the signature
Definition: pkfem.hh:92
finiteelementmap.hh
Dune::PDELab::PkLocalFiniteElementMap
Definition: pkfem.hh:288
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 3 >::PkLocalFiniteElementMapBase
PkLocalFiniteElementMapBase(const GV &gv)
Definition: pkfem.hh:188
Dune::PDELab::InvalidGeometryType
FiniteElementMap exception raised when trying to obtain a finite element for an unsupported GeometryT...
Definition: finiteelementmap.hh:23
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 3 >::Traits
LocalFiniteElementMapTraits< FE > Traits
export type of the signature
Definition: pkfem.hh:186
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 2 >::maxLocalSize
static constexpr std::size_t maxLocalSize()
Definition: pkfem.hh:156
Dune::PDELab::LocalFiniteElementMapTraits
collect types exported by a finite element map
Definition: finiteelementmap.hh:38
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 2 >::PkLocalFiniteElementMapBase
PkLocalFiniteElementMapBase(const GV &gv)
Definition: pkfem.hh:94
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 2 >::size
static constexpr std::size_t size(GeometryType gt)
Definition: pkfem.hh:145
Dune::PDELab::FiniteElementMapTraits::FiniteElementType
T FiniteElementType
Type of finite element from local functions.
Definition: finiteelementmap.hh:30
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 1 >::hasDOFs
static constexpr bool hasDOFs(int codim)
Definition: pkfem.hh:43
Dune::PDELab::fem::PkLocalFiniteElementMapBase< GV, D, R, k, 3 >::hasDOFs
static constexpr bool hasDOFs(int codim)
Definition: pkfem.hh:235