5#ifndef DUNE_LOCALFUNCTIONS_NEDELEC_NEDELEC1STKINDCUBE_HH
6#define DUNE_LOCALFUNCTIONS_NEDELEC_NEDELEC1STKINDCUBE_HH
10#include <dune/common/fmatrix.hh>
11#include <dune/common/fvector.hh>
12#include <dune/common/math.hh>
14#include <dune/geometry/referenceelements.hh>
15#include <dune/geometry/type.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
211 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
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)
422 auto refElement = Dune::referenceElement<double,dim>(GeometryTypes::cube(dim));
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();
539 static constexpr GeometryType
type ()
541 return GeometryTypes::cube(dim);
Definition bdfmcube.hh:18
D DomainType
domain type
Definition common/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
Nédélec elements of the first kind for cube elements.
Definition nedelec1stkindcube.hh:496
const Traits::LocalInterpolationType & localInterpolation() const
Definition nedelec1stkindcube.hh:529
static constexpr unsigned int size()
Definition nedelec1stkindcube.hh:534
Nedelec1stKindCubeLocalFiniteElement()=default
Default constructor.
static constexpr GeometryType type()
Definition nedelec1stkindcube.hh:539
const Traits::LocalCoefficientsType & localCoefficients() const
Definition nedelec1stkindcube.hh:524
const Traits::LocalBasisType & localBasis() const
Definition nedelec1stkindcube.hh:519
Nedelec1stKindCubeLocalFiniteElement(std::bitset< power(2, dim-1) *dim > s)
Constructor with explicitly given edge orientations.
Definition nedelec1stkindcube.hh:514