2 #ifndef DUNE_PDELAB_FINITEELEMENTMAP_PKFEM_HH
3 #define DUNE_PDELAB_FINITEELEMENTMAP_PKFEM_HH
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>
20 template<
typename GV,
typename D,
typename R,
unsigned int k,
unsigned int d>
28 template<
typename GV,
typename D,
typename R,
unsigned int k>
52 assert(k >= 0 and k <= 1);
57 static constexpr std::size_t
size(GeometryType gt)
59 if (gt == GeometryTypes::vertex)
61 if (gt == GeometryTypes::line)
62 return k > 0 ? k - 1 : 1;
78 template<
typename GV,
typename D,
typename R,
unsigned int k>
81 Dune::PkLocalFiniteElement<D,R,2,k>
83 PkLocalFiniteElementMapBase<GV,D,R,k,2>
87 typedef Dune::PkLocalFiniteElement<D,R,2,k> FE;
98 unsigned int p[3] = {0,1,2};
99 for (
int i = 0; i < 6; ++i)
102 std::next_permutation(
p,
p+3);
107 template<
typename Entity>
111 if (!
e.type().isSimplex())
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);
119 unsigned int n0_compressed = (n0 > n1) + (n0 > n2);
121 return _variant[2 * n0_compressed + (n1 > n2)];
138 return k > 2 || k == 0;
140 assert(
false &&
"Invalid codim specified!");
145 static constexpr std::size_t
size(GeometryType gt)
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);
158 return (k+1)*(k+2)/2;
162 std::array<FE,6> _variant;
172 template<
typename GV,
typename D,
typename R,
unsigned int k>
175 Dune::PkLocalFiniteElement<D,R,3,k>
177 PkLocalFiniteElementMapBase<GV,D,R,k,3>
181 typedef Dune::PkLocalFiniteElement<D,R,3,k> FE;
191 std::fill(_perm_index.begin(),_perm_index.end(),0);
195 unsigned int vertexmap[4];
196 for(vertexmap[0] = 0; vertexmap[0] < 4; ++vertexmap[0])
198 for(vertexmap[1] = 0; vertexmap[1] < 4; ++vertexmap[1])
200 if (vertexmap[0] == vertexmap[1])
202 for(vertexmap[2] = 0; vertexmap[2] < 4; ++vertexmap[2])
204 if (vertexmap[0] == vertexmap[2] ||
205 vertexmap[1] == vertexmap[2])
207 vertexmap[3] = 6 - vertexmap[0] - vertexmap[1] - vertexmap[2];
208 _variant[n] = FE(vertexmap);
209 _perm_index[compressPerm(vertexmap)] = n++;
216 template<
typename Entity>
220 if (!
e.type().isSimplex())
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);
230 for (
unsigned int i = 0; i < 4; ++i)
233 for (
unsigned int j = 0; j < 4; ++j)
234 if ((min_index < 0 || vertexmap[j] < vertexmap[min_index]) && vertexmap[j] >= i)
236 assert(min_index >= 0);
237 vertexmap[min_index] = i;
239 return _variant[_perm_index[compressPerm(vertexmap)]];
258 return k == 0 || k > 3;
260 assert(
false &&
"Invalid codim specified!");
265 static constexpr std::size_t
size(GeometryType gt)
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;
280 return (k+1)*(k+2)*(k+3)/6;
285 unsigned int compressPerm(
const unsigned int vertexmap[4])
const
287 return vertexmap[0] + (vertexmap[1]<<2) + (vertexmap[2]<<4) + (vertexmap[3]<<6);
290 std::array<FE,24> _variant;
291 std::array<unsigned int,256> _perm_index;
299 template<
typename GV,
typename D,
typename R,
unsigned int k>
310 : fem::PkLocalFiniteElementMapBase<GV,D,R,k,GV::
dimension>(gv)
319 #endif // DUNE_PDELAB_FINITEELEMENTMAP_PKFEM_HH