5#ifndef DUNE_LOCALFUNCTIONS_NEDELEC_NEDELEC1STKINDCUBE_HH
6#define DUNE_LOCALFUNCTIONS_NEDELEC_NEDELEC1STKINDCUBE_HH
14#include <dune/geometry/referenceelements.hh>
17#include <dune/localfunctions/common/localbasis.hh>
18#include <dune/localfunctions/common/localfiniteelementtraits.hh>
19#include <dune/localfunctions/common/localkey.hh>
37 template<
class D,
class R,
int dim,
int k>
38 class Nedelec1stKindCubeLocalBasis
41 constexpr static std::size_t numberOfEdges =
power(2,dim-1)*dim;
44 using Traits = LocalBasisTraits<D,dim,FieldVector<D,dim>,
45 R,dim,FieldVector<R,dim>,
46 FieldMatrix<R,dim,dim> >;
54 Nedelec1stKindCubeLocalBasis()
56 std::fill(edgeOrientation_.begin(), edgeOrientation_.end(), 1.0);
61 Nedelec1stKindCubeLocalBasis(std::bitset<numberOfEdges> edgeOrientation)
62 : Nedelec1stKindCubeLocalBasis()
64 for (std::size_t i=0; i<edgeOrientation_.size(); i++)
65 edgeOrientation_[i] *= edgeOrientation[i] ? -1.0 : 1.0;
69 static constexpr unsigned int size()
71 static_assert(dim==2 || dim==3,
"Nedelec shape functions are implemented only for 2d and 3d cubes.");
75 return 3*k * (k+1) * (k+1);
84 std::vector<typename Traits::RangeType>& out)
const
86 static_assert(k==1,
"Evaluating Nédélec shape functions is implemented only for first order.");
103 out[0] = { 0, D(1) - in[0]};
104 out[1] = { 0, in[0]};
105 out[2] = { D(1) - in[1], 0};
106 out[3] = { in[1], 0};
109 if constexpr (dim==3)
128 for (std::size_t i=0; i<out.size(); i++)
131 out[0][2] = { 1 - in[0] - in[1] + in[0]*in[1]};
132 out[1][2] = { in[0] - in[0]*in[1]};
133 out[2][2] = { in[1] - in[0]*in[1]};
134 out[3][2] = { in[0]*in[1]};
136 out[4][1] = { 1 - in[0] - in[2] + in[0]*in[2]};
137 out[5][1] = { in[0] - in[0]*in[2]};
138 out[8][1] = { in[2] - in[0]*in[2]};
139 out[9][1] = { in[0]*in[2]};
141 out[6][0] = { 1 - in[1] - in[2] + in[1]*in[2]};
142 out[7][0] = { in[1] - in[1]*in[2]};
143 out[10][0] = { in[2] - in[1]*in[2]};
144 out[11][0] = { in[1]*in[2]};
147 for (std::size_t i=0; i<out.size(); i++)
148 out[i] *= edgeOrientation_[i];
157 std::vector<typename Traits::JacobianType>& out)
const
162 for (std::size_t i=0; i<out.size(); i++)
163 for (std::size_t j=0; j<dim; j++)
166 out[0][1] = { -1, 0};
169 out[2][0] = { 0, -1};
174 for (std::size_t i=0; i<out.size(); i++)
175 for(std::size_t j=0;j<dim; j++)
179 out[0][2] = {-1 +in[1], -1 + in[0], 0};
180 out[1][2] = { 1 -in[1], - in[0], 0};
181 out[2][2] = { -in[1], 1 - in[0], 0};
182 out[3][2] = { in[1], in[0], 0};
184 out[4][1] = {-1 +in[2], 0, -1 + in[0]};
185 out[5][1] = { 1 -in[2], 0, - in[0]};
186 out[8][1] = { -in[2], 0, 1 - in[0]};
187 out[9][1] = { in[2], 0, in[0]};
189 out[6][0] = { 0, -1 + in[2], -1 + in[1]};
190 out[7][0] = { 0, 1 - in[2], - in[1]};
191 out[10][0] = { 0, - in[2], 1 - in[1]};
192 out[11][0] = { 0, in[2], in[1]};
196 for (std::size_t i=0; i<out.size(); i++)
197 out[i] *= edgeOrientation_[i];
207 void partial(
const std::array<unsigned int, dim>& order,
209 std::vector<typename Traits::RangeType>& out)
const
212 if (totalOrder == 0) {
213 evaluateFunction(in, out);
214 }
else if (totalOrder == 1) {
215 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
241 out[0] = { 0, 0, -1 +in[1]};
242 out[1] = { 0, 0, 1 -in[1]};
243 out[2] = { 0, 0, -in[1]};
244 out[3] = { 0, 0, in[1]};
246 out[4] = { 0, -1 +in[2], 0};
247 out[5] = { 0, 1 -in[2], 0};
248 out[8] = { 0, -in[2], 0};
249 out[9] = { 0, in[2], 0};
258 out[0] = { 0, 0, -1 + in[0]};
259 out[1] = { 0, 0, - in[0]};
260 out[2] = { 0, 0, 1 - in[0]};
261 out[3] = { 0, 0, in[0]};
268 out[6] = { -1 + in[2], 0, 0};
269 out[7] = { 1 - in[2], 0, 0};
270 out[10] = { - in[2], 0, 0};
271 out[11] = { in[2], 0, 0};
280 out[4] = { 0, -1 + in[0], 0};
281 out[5] = { 0, - in[0], 0};
282 out[8] = { 0, 1 - in[0], 0};
283 out[9] = { 0, in[0], 0};
285 out[6] = { -1 + in[1], 0, 0};
286 out[7] = { - in[1], 0, 0};
287 out[10] = { 1 - in[1], 0, 0};
288 out[11] = { in[1], 0, 0};
293 for (std::size_t i=0; i<out.size(); i++)
294 out[i] *= edgeOrientation_[i];
296 }
else if (totalOrder == 2) {
300 for (std::size_t i = 0; i <
size(); ++i)
301 for (std::size_t j = 0; j < dim; ++j)
306 for(
size_t i=0; i<out.size(); i++)
310 if( order[0] == 1 and order[1]==1)
313 out[1] = { 0, 0, -1};
314 out[2] = { 0, 0, -1};
319 if( order[0] == 1 and order[2]==1)
322 out[5] = { 0, -1, 0};
323 out[8] = { 0, -1, 0};
328 if( order[1] == 1 and order[2]==1)
331 out[7] = { -1, 0, 0};
332 out[10] = { -1, 0, 0};
333 out[11] = { 1, 0, 0};
336 for (std::size_t i=0; i<out.size(); i++)
337 out[i] *= edgeOrientation_[i];
343 for (std::size_t i = 0; i <
size(); ++i)
344 for (std::size_t j = 0; j < dim; ++j)
351 unsigned int order()
const
362 std::array<R,numberOfEdges> edgeOrientation_;
370 template <
int dim,
int k>
371 class Nedelec1stKindCubeLocalCoefficients
375 Nedelec1stKindCubeLocalCoefficients ()
378 static_assert(k==1,
"Only first-order Nédélec local coefficients are implemented.");
381 for (std::size_t i=0; i<
size(); i++)
382 localKey_[i] = LocalKey(i,dim-1,0);
386 std::size_t
size()
const
388 static_assert(dim==2 || dim==3,
"Nédélec shape functions are implemented only for 2d and 3d cubes.");
389 return (dim==2) ? 2*k * (k+1)
390 : 3*k * (k+1) * (k+1);
395 const LocalKey& localKey (std::size_t i)
const
401 std::vector<LocalKey> localKey_;
409 class Nedelec1stKindCubeLocalInterpolation
411 static constexpr auto dim = LB::Traits::dimDomain;
412 static constexpr auto size = LB::size();
415 constexpr static std::size_t numberOfEdges =
power(2,dim-1)*dim;
420 Nedelec1stKindCubeLocalInterpolation (std::bitset<numberOfEdges> s = 0)
424 for (std::size_t i=0; i<numberOfEdges; i++)
425 m_[i] = refElement.position(i,dim-1);
427 for (std::size_t i=0; i<numberOfEdges; i++)
429 auto vertexIterator = refElement.subEntities(i,dim-1,dim).begin();
430 auto v0 = *vertexIterator;
431 auto v1 = *(++vertexIterator);
437 edge_[i] = refElement.position(v1,dim) - refElement.position(v0,dim);
438 edge_[i] *= (s[i]) ? -1.0 : 1.0;
447 template<
typename F,
typename C>
448 void interpolate (
const F& f, std::vector<C>& out)
const
452 for (std::size_t i=0; i<
size; i++)
456 for (
int j=0; j<dim; j++)
457 out[i] += y[j]*edge_[i][j];
463 std::array<typename LB::Traits::DomainType, numberOfEdges> m_;
465 std::array<typename LB::Traits::DomainType, numberOfEdges> edge_;
494 template<
class D,
class R,
int dim,
int k>
499 Impl::Nedelec1stKindCubeLocalCoefficients<dim,k>,
500 Impl::Nedelec1stKindCubeLocalInterpolation<Impl::Nedelec1stKindCubeLocalBasis<D,R,dim,k> > >;
502 static_assert(dim==2 || dim==3,
"Nedelec elements are only implemented for 2d and 3d elements.");
503 static_assert(k==1,
"Nedelec elements of the first kind are currently only implemented for order k==1.");
526 return coefficients_;
531 return interpolation_;
534 static constexpr unsigned int size ()
536 return Traits::LocalBasisType::size();
Nédélec elements of the first kind for cube elements.
Definition: nedelec1stkindcube.hh:496
Nedelec1stKindCubeLocalFiniteElement()=default
Default constructor.
Nedelec1stKindCubeLocalFiniteElement(std::bitset< power(2, dim-1) *dim > s)
Constructor with explicitly given edge orientations.
Definition: nedelec1stkindcube.hh:514
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:132
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
static void interpolate(const GridFunction &u, DiscreteFunction &v)
perform native interpolation of a discrete function space
Definition: interpolate.hh:54
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:462
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
Some useful basic math stuff.
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
constexpr Base power(Base m, Exponent p)
Power method for integer exponents.
Definition: math.hh:75
D DomainType
domain type
Definition: localbasis.hh:43
traits helper struct
Definition: localfiniteelementtraits.hh:13
LB LocalBasisType
Definition: localfiniteelementtraits.hh:16
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:20
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:24
A unique label for each type of element that can occur in a grid.