2#ifndef DUNE_PDELAB_FINITEELEMENTMAP_PKFEM_HH
3#define DUNE_PDELAB_FINITEELEMENTMAP_PKFEM_HH
11#include<dune/localfunctions/lagrange/lagrangesimplex.hh>
13#include "finiteelementmap.hh"
20 template<
typename GV,
typename D,
typename R,
unsigned int k,
unsigned int d>
21 class PkLocalFiniteElementMapBase;
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>
35 PkLocalFiniteElementMapBase(
const GV& gv)
43 static constexpr bool hasDOFs(
int codim)
52 assert(k >= 0 and k <= 1);
57 static constexpr std::size_t
size(GeometryType
gt)
62 return k > 0 ? k - 1 : 1;
66 static constexpr std::size_t maxLocalSize()
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>
83 PkLocalFiniteElementMapBase<GV,D,R,k,2>
92 typedef LocalFiniteElementMapTraits<FE> Traits;
94 PkLocalFiniteElementMapBase(
const GV& gv)
98 std::array<unsigned int, 3> p = {0,1,2};
99 for (
int i = 0; i < 6; ++i)
102 std::next_permutation(p.begin(), p.end());
107 template<
typename Entity>
108 const typename Traits::FiniteElementType& find (
const Entity& e)
const
111 if (!e.type().isSimplex())
112 DUNE_THROW(InvalidGeometryType,
"PkLocalFiniteElementMap only works for simplices!");
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)];
129 static constexpr bool hasDOFs(
int codim)
138 return k > 2 || k == 0;
140 assert(
false &&
"Invalid codim specified!");
145 static constexpr std::size_t
size(GeometryType
gt)
148 return k > 0 ? 1 : 0;
150 return k > 1 ? k - 1 : 0;
152 return k > 2 ? (k-2)*(k-1)/2 : (k == 0);
156 static constexpr std::size_t maxLocalSize()
158 return (k+1)*(k+2)/2;
162 std::array<FE,6> _variant;
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>
177 PkLocalFiniteElementMapBase<GV,D,R,k,3>
186 typedef LocalFiniteElementMapTraits<FE> Traits;
188 PkLocalFiniteElementMapBase(
const GV& gv)
191 std::fill(_perm_index.begin(),_perm_index.end(),0);
195 std::array<unsigned int, 4> vertexmap = {0, 1, 2, 3};
197 _variant[n] = FE(vertexmap);
198 _perm_index[compressPerm(vertexmap)] = n++;
200 while(std::next_permutation(vertexmap.begin(), vertexmap.end()));
204 template<
typename Entity>
205 const typename Traits::FiniteElementType& find (
const Entity& e)
const
208 if (!e.type().isSimplex())
209 DUNE_THROW(InvalidGeometryType,
"PkLocalFiniteElementMap only works for simplices!");
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);
218 for (
unsigned int i = 0; i < 4; ++i)
221 for (
unsigned int j = 0; j < 4; ++j)
222 if ((min_index < 0 || vertexmap[j] < vertexmap[min_index]) && vertexmap[j] >= i)
224 assert(min_index >= 0);
225 vertexmap[min_index] = i;
227 return _variant[_perm_index[compressPerm(vertexmap)]];
235 static constexpr bool hasDOFs(
int codim)
246 return k == 0 || k > 3;
248 assert(
false &&
"Invalid codim specified!");
253 static constexpr std::size_t
size(GeometryType
gt)
256 return k > 0 ? 1 : 0;
258 return k > 1 ? k - 1 : 0;
260 return k > 2 ? (k-2)*(k-1)/2 : 0;
262 return k == 0 ? 1 : (k-3)*(k-2)*(k-1)/6;
266 static constexpr std::size_t maxLocalSize()
268 return (k+1)*(k+2)*(k+3)/6;
273 unsigned int compressPerm(
const std::array<unsigned int, 4> vertexmap)
const
275 return vertexmap[0] + (vertexmap[1]<<2) + (vertexmap[2]<<4) + (vertexmap[3]<<6);
278 std::array<FE,24> _variant;
279 std::array<unsigned int,256> _perm_index;
287 template<
typename GV,
typename D,
typename R,
unsigned int k>
288 class PkLocalFiniteElementMap
289 :
public fem::PkLocalFiniteElementMapBase<GV,D,R,k,GV::dimension>
295 static constexpr int dimension = GV::dimension;
297 PkLocalFiniteElementMap(
const GV& gv)
298 : fem::PkLocalFiniteElementMapBase<GV,D,R,k,GV::dimension>(gv)
Lagrange finite element for simplices with arbitrary compile-time dimension and polynomial order.
Definition: lagrangesimplex.hh:836
Definition of the DUNE_NO_DEPRECATED_* macros.
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:158
constexpr GeometryType line
GeometryType representing a line.
Definition: type.hh:498
constexpr GeometryType triangle
GeometryType representing a triangle.
Definition: type.hh:504
constexpr GeometryType tetrahedron
GeometryType representing a tetrahedron.
Definition: type.hh:516
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:492
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
A unique label for each type of element that can occur in a grid.
std::size_t fixedSize
The number of data items per index if it is fixed, 0 otherwise.
Definition: variablesizecommunicator.hh:264