dune-pdelab  2.0.0
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 
7 #include <dune/common/deprecated.hh>
8 #include <dune/common/array.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> >
31  {
32 
33  public:
34 
36  {}
37 
38  bool fixedSize() const
39  {
40  return true;
41  }
42 
43  std::size_t size(GeometryType gt) const
44  {
45  if (gt.isVertex())
46  return k > 0 ? 1 : 0;
47  if (gt.isLine())
48  return k > 0 ? k - 1 : 1;
49  return 0;
50  }
51 
52  std::size_t maxLocalSize() const
53  {
54  return k + 1;
55  }
56 
57  };
58 
59 
60  // ********************************************************************************
61  // 2D version
62  // ********************************************************************************
63 
64  template<typename GV, typename D, typename R, unsigned int k>
65  class PkLocalFiniteElementMapBase<GV,D,R,k,2>
66  : public LocalFiniteElementMapInterface<LocalFiniteElementMapTraits<
67  Dune::PkLocalFiniteElement<D,R,2,k>
68  >,
69  PkLocalFiniteElementMapBase<GV,D,R,k,2>
70  >
71  {
72 
73  typedef Dune::PkLocalFiniteElement<D,R,2,k> FE;
74 
75  public:
76 
79 
81  : _gv(gv)
82  {
83  // generate permutations
84  unsigned int p[3] = {0,1,2};
85  for (int i = 0; i < 6; ++i)
86  {
87  _variant[i] = FE(p);
88  std::next_permutation(p,p+3);
89  }
90  }
91 
93  template<typename Entity>
94  const typename Traits::FiniteElementType& find (const Entity& e) const
95  {
96 
97  if (!e.type().isSimplex())
98  DUNE_THROW(InvalidGeometryType,"PkLocalFiniteElementMap only works for simplices!");
99 
100  const typename GV::IndexSet& is = _gv.indexSet();
101  unsigned int n0 = is.subIndex(e,0,2);
102  unsigned int n1 = is.subIndex(e,1,2);
103  unsigned int n2 = is.subIndex(e,2,2);
104  // compress first number to [0,2]
105  unsigned int n0_compressed = (n0 > n1) + (n0 > n2);
106  // translate to permutation index
107  return _variant[2 * n0_compressed + (n1 > n2)];
108  }
109 
110  bool fixedSize() const
111  {
112  return true;
113  }
114 
115  std::size_t size(GeometryType gt) const
116  {
117  if (gt.isVertex())
118  return k > 0 ? 1 : 0;
119  if (gt.isLine())
120  return k > 1 ? k - 1 : 0;
121  if (gt.isTriangle())
122  return k > 2 ? (k-2)*(k-1)/2 : (k == 0);
123  return 0;
124  }
125 
126  std::size_t maxLocalSize() const
127  {
128  return (k+1)*(k+2)/2;
129  }
130 
131  private:
132  array<FE,6> _variant;
133  GV _gv;
134 
135  };
136 
137 
138  // ********************************************************************************
139  // 3D version
140  // ********************************************************************************
141 
142  template<typename GV, typename D, typename R, unsigned int k>
143  class PkLocalFiniteElementMapBase<GV,D,R,k,3>
144  : public LocalFiniteElementMapInterface<LocalFiniteElementMapTraits<
145  Dune::PkLocalFiniteElement<D,R,3,k>
146  >,
147  PkLocalFiniteElementMapBase<GV,D,R,k,3>
148  >
149  {
150 
151  typedef Dune::PkLocalFiniteElement<D,R,3,k> FE;
152 
153  public:
154 
157 
159  : _gv(gv)
160  {
161  std::fill(_perm_index.begin(),_perm_index.end(),0);
162 
163  // create all variants by iterating over all permutations
164  unsigned int n = 0;
165  unsigned int vertexmap[4];
166  for(vertexmap[0] = 0; vertexmap[0] < 4; ++vertexmap[0])
167  {
168  for(vertexmap[1] = 0; vertexmap[1] < 4; ++vertexmap[1])
169  {
170  if (vertexmap[0] == vertexmap[1])
171  continue;
172  for(vertexmap[2] = 0; vertexmap[2] < 4; ++vertexmap[2])
173  {
174  if (vertexmap[0] == vertexmap[2] ||
175  vertexmap[1] == vertexmap[2])
176  continue;
177  vertexmap[3] = 6 - vertexmap[0] - vertexmap[1] - vertexmap[2];
178  _variant[n] = FE(vertexmap);
179  _perm_index[compressPerm(vertexmap)] = n++;
180  }
181  }
182  }
183  }
184 
186  template<typename Entity>
187  const typename Traits::FiniteElementType& find (const Entity& e) const
188  {
189 
190  if (!e.type().isSimplex())
191  DUNE_THROW(InvalidGeometryType,"PkLocalFiniteElementMap only works for simplices!");
192 
193  // get the vertex indices
194  const typename GV::IndexSet& is = _gv.indexSet();
195  unsigned int vertexmap[4];
196  for (unsigned int i = 0; i < 4; ++i)
197  vertexmap[i] = is.subIndex(e,i,3);
198 
199  // reduce the indices to the interval 0..3
200  for (unsigned int i = 0; i < 4; ++i)
201  {
202  int min_index = -1;
203  for (unsigned int j = 0; j < 4; ++j)
204  if ((min_index < 0 || vertexmap[j] < vertexmap[min_index]) && vertexmap[j] >= i)
205  min_index = j;
206  assert(min_index >= 0);
207  vertexmap[min_index] = i;
208  }
209  return _variant[_perm_index[compressPerm(vertexmap)]];
210  }
211 
212  bool fixedSize() const
213  {
214  return true;
215  }
216 
217  std::size_t size(GeometryType gt) const
218  {
219  if (gt.isVertex())
220  return k > 0 ? 1 : 0;
221  if (gt.isLine())
222  return k > 1 ? k - 1 : 0;
223  if (gt.isTriangle())
224  return k > 2 ? (k-2)*(k-1)/2 : 0;
225  if (gt.isTetrahedron())
226  return k == 0 ? 1 : (k-3)*(k-2)*(k-1)/6;
227  return 0;
228  }
229 
230  std::size_t maxLocalSize() const
231  {
232  return (k+1)*(k+2)*(k+3)/6;
233  }
234 
235  private:
236 
237  unsigned int compressPerm(const unsigned int vertexmap[4]) const
238  {
239  return vertexmap[0] + (vertexmap[1]<<2) + (vertexmap[2]<<4) + (vertexmap[3]<<6);
240  }
241 
242  array<FE,24> _variant;
243  array<unsigned int,256> _perm_index;
244  GV _gv;
245 
246  };
247 
248 
249 #ifndef DOXYGEN
250 
251  template<unsigned int d>
252  struct DUNE_DEPRECATED_MSG("The fifth template parameter for PkLocalFiniteElementMap (the dimension) is deprecated") warn_deprecated_pk_interface
253  {
254  warn_deprecated_pk_interface() DUNE_DEPRECATED_MSG("The fifth template parameter for PkLocalFiniteElementMap (the dimension) is deprecated")
255  {}
256  };
257 
258  template<>
259  struct warn_deprecated_pk_interface<0>
260  {
261  warn_deprecated_pk_interface()
262  {}
263  };
264 
265 #endif // DOXYGEN
266 
267 
268  } // namespace fem
269 
270 
271  template<typename GV, typename D, typename R, unsigned int k, unsigned int d = 0>
273  : public fem::PkLocalFiniteElementMapBase<GV,D,R,k,GV::dimension>
274 #ifndef DOXYGEN
275  , public fem::warn_deprecated_pk_interface<d>
276 #endif
277  {
278 
279  public:
280 
282  : fem::PkLocalFiniteElementMapBase<GV,D,R,k,GV::dimension>(gv)
283 #ifndef DOXYGEN
284  , fem::warn_deprecated_pk_interface<d>()
285 #endif
286  {}
287 
288  };
289 
290 
291  }
292 }
293 
294 #endif // DUNE_PDELAB_FINITEELEMENTMAP_PKFEM_HH
std::size_t maxLocalSize() const
Definition: pkfem.hh:126
FiniteElementMap exception raised when trying to obtain a finite element for an unsupported GeometryT...
Definition: finiteelementmap.hh:23
PkLocalFiniteElementMapBase(const GV &gv)
Definition: pkfem.hh:35
const Traits::FiniteElementType & find(const Entity &e) const
get local basis functions for entity
Definition: pkfem.hh:94
PkLocalFiniteElementMapBase(const GV &gv)
Definition: pkfem.hh:80
std::size_t size(GeometryType gt) const
Definition: pkfem.hh:43
std::size_t size(GeometryType gt) const
Definition: pkfem.hh:217
interface for a finite element map
Definition: finiteelementmap.hh:42
LocalFiniteElementMapTraits< FE > Traits
export type of the signature
Definition: pkfem.hh:156
const Traits::FiniteElementType & find(const Entity &e) const
get local basis functions for entity
Definition: pkfem.hh:187
LocalFiniteElementMapTraits< FE > Traits
export type of the signature
Definition: pkfem.hh:78
T FiniteElementType
Type of finite element from local functions.
Definition: finiteelementmap.hh:30
std::size_t size(GeometryType gt) const
Definition: pkfem.hh:115
std::size_t maxLocalSize() const
Definition: pkfem.hh:52
simple implementation where all entities have the same finite element
Definition: finiteelementmap.hh:107
PkLocalFiniteElementMap(const GV &gv)
Definition: pkfem.hh:281
R p(int i, D x)
Definition: qkdg.hh:62
collect types exported by a finite element map
Definition: finiteelementmap.hh:38
const E & e
Definition: interpolate.hh:172
std::size_t maxLocalSize() const
Definition: pkfem.hh:230
PkLocalFiniteElementMapBase(const GV &gv)
Definition: pkfem.hh:158