dune-localfunctions  2.8.0
nedelec1stkindsimplex.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_LOCALFUNCTIONS_NEDELEC_NEDELEC1STKINDSIMPLEX_HH
4 #define DUNE_LOCALFUNCTIONS_NEDELEC_NEDELEC1STKINDSIMPLEX_HH
5 
6 #include <numeric>
7 
8 #include <dune/common/fmatrix.hh>
9 #include <dune/common/fvector.hh>
10 
11 #include <dune/geometry/referenceelements.hh>
12 #include <dune/geometry/type.hh>
13 
16 #include <dune/localfunctions/common/localinterpolation.hh> // For deprecated makeFunctionWithCallOperator
18 
19 namespace Dune
20 {
21 namespace Impl
22 {
33  template<class D, class R, int dim, int k>
34  class Nedelec1stKindSimplexLocalBasis
35  {
36  // Number of edges of the reference simplex
37  constexpr static std::size_t numberOfEdges = dim*(dim+1)/2;
38 
39  public:
40  using Traits = LocalBasisTraits<D,dim,FieldVector<D,dim>,
41  R,dim,FieldVector<R,dim>,
42  FieldMatrix<R,dim,dim> >;
43 
50  Nedelec1stKindSimplexLocalBasis()
51  {
52  std::fill(edgeOrientation_.begin(), edgeOrientation_.end(), 1.0);
53  }
54 
57  Nedelec1stKindSimplexLocalBasis(std::bitset<numberOfEdges> edgeOrientation)
58  : Nedelec1stKindSimplexLocalBasis()
59  {
60  for (std::size_t i=0; i<edgeOrientation_.size(); i++)
61  edgeOrientation_[i] *= edgeOrientation[i] ? -1.0 : 1.0;
62  }
63 
65  static constexpr unsigned int size()
66  {
67  static_assert(dim==2 || dim==3, "Nedelec shape functions are implemented only for 2d and 3d simplices.");
68  if (dim==2)
69  return k * (k+2);
70  if (dim==3)
71  return k * (k+2) * (k+3) / 2;
72  }
73 
79  void evaluateFunction(const typename Traits::DomainType& in,
80  std::vector<typename Traits::RangeType>& out) const
81  {
82  static_assert(k==1, "Evaluating Nédélec shape functions is implemented only for first order.");
83  out.resize(size());
84 
85  if (dim==2)
86  {
87  // First-order Nédélec shape functions on a triangle are of the form
88  //
89  // (a1, a2) + b(-x2, x1)^T, a_1, a_2, b \in R
90  out[0] = {D(1) - in[1], in[0]};
91  out[1] = {in[1], -in[0]+D(1)};
92  out[2] = {-in[1], in[0]};
93  }
94 
95  if constexpr (dim==3)
96  {
97  // First-order Nédélec shape functions on a tetrahedron are of the form
98  //
99  // a + b \times x, a, b \in R^3
100  //
101  // The following coefficients create the six basis vectors
102  // that are dual to the edge degrees of freedom:
103  //
104  // a[0] = { 1, 0, 0} b[0] = { 0, -1, 1}
105  // a[1] = { 0, 1, 0} b[1] = { 1, 0, -1}
106  // a[2] = { 0, 0, 0} b[2] = { 0, 0, 1}
107  // a[3] = { 0, 0, 1} b[3] = {-1, 1, 0}
108  // a[4] = { 0, 0, 0} b[4] = { 0, -1, 0}
109  // a[5] = { 0, 0, 0} b[5] = { 1, 0, 0}
110  //
111  // The following implementation uses these values, and simply
112  // skips all the zeros.
113 
114  out[0] = { 1 - in[1] - in[2], in[0] , in[0] };
115  out[1] = { in[1] , 1 - in[0] - in[2], in[1]};
116  out[2] = { - in[1] , in[0] , 0 };
117  out[3] = { in[2], in[2], 1 - in[0] - in[1]};
118  out[4] = { -in[2], 0 , in[0] };
119  out[5] = { 0 , -in[2], in[1]};
120  }
121 
122  for (std::size_t i=0; i<out.size(); i++)
123  out[i] *= edgeOrientation_[i];
124  }
125 
131  void evaluateJacobian(const typename Traits::DomainType& in,
132  std::vector<typename Traits::JacobianType>& out) const
133  {
134  out.resize(size());
135  if (dim==2)
136  {
137  out[0][0] = { 0, -1};
138  out[0][1] = { 1, 0};
139 
140  out[1][0] = { 0, 1};
141  out[1][1] = {-1, 0};
142 
143  out[2][0] = { 0, -1};
144  out[2][1] = { 1, 0};
145  }
146  if (dim==3)
147  {
148  out[0][0] = { 0,-1,-1};
149  out[0][1] = { 1, 0, 0};
150  out[0][2] = { 1, 0, 0};
151 
152  out[1][0] = { 0, 1, 0};
153  out[1][1] = {-1, 0, -1};
154  out[1][2] = { 0, 1, 0};
155 
156  out[2][0] = { 0, -1, 0};
157  out[2][1] = { 1, 0, 0};
158  out[2][2] = { 0, 0, 0};
159 
160  out[3][0] = { 0, 0, 1};
161  out[3][1] = { 0, 0, 1};
162  out[3][2] = {-1, -1, 0};
163 
164  out[4][0] = { 0, 0, -1};
165  out[4][1] = { 0, 0, 0};
166  out[4][2] = { 1, 0, 0};
167 
168  out[5][0] = { 0, 0, 0};
169  out[5][1] = { 0, 0, -1};
170  out[5][2] = { 0, 1, 0};
171  }
172 
173  for (std::size_t i=0; i<out.size(); i++)
174  out[i] *= edgeOrientation_[i];
175 
176  }
177 
184  void partial(const std::array<unsigned int, dim>& order,
185  const typename Traits::DomainType& in,
186  std::vector<typename Traits::RangeType>& out) const
187  {
188  auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
189  if (totalOrder == 0) {
190  evaluateFunction(in, out);
191  } else if (totalOrder == 1) {
192  auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
193  out.resize(size());
194 
195  if (dim==2)
196  {
197  if (direction==0)
198  {
199  out[0] = {0, 1};
200  out[1] = {0, -1};
201  out[2] = {0, 1};
202  }
203  else
204  {
205  out[0] = {-1, 0};
206  out[1] = { 1, 0};
207  out[2] = {-1, 0};
208  }
209  }
210 
211  if (dim==3)
212  {
213  switch (direction)
214  {
215  case 0:
216  out[0] = { 0, 1, 1};
217  out[1] = { 0,-1, 0};
218  out[2] = { 0, 1, 0};
219  out[3] = { 0, 0,-1};
220  out[4] = { 0, 0, 1};
221  out[5] = { 0, 0, 0};
222  break;
223 
224  case 1:
225  out[0] = {-1, 0, 0};
226  out[1] = { 1, 0, 1};
227  out[2] = {-1, 0, 0};
228  out[3] = { 0, 0,-1};
229  out[4] = { 0, 0, 0};
230  out[5] = { 0, 0, 1};
231  break;
232 
233  case 2:
234  out[0] = {-1, 0, 0};
235  out[1] = { 0,-1, 0};
236  out[2] = { 0, 0, 0};
237  out[3] = { 1, 1, 0};
238  out[4] = {-1, 0, 0};
239  out[5] = { 0,-1, 0};
240  break;
241  }
242  }
243 
244  for (std::size_t i=0; i<out.size(); i++)
245  out[i] *= edgeOrientation_[i];
246 
247  } else {
248  out.resize(size());
249  for (std::size_t i = 0; i < size(); ++i)
250  for (std::size_t j = 0; j < dim; ++j)
251  out[i][j] = 0;
252  }
253 
254  }
255 
257  unsigned int order() const
258  {
259  return k;
260  }
261 
262  private:
263 
264  // Orientations of the simplex edges
265  std::array<R,numberOfEdges> edgeOrientation_;
266  };
267 
268 
273  template <int dim, int k>
274  class Nedelec1stKindSimplexLocalCoefficients
275  {
276  public:
278  Nedelec1stKindSimplexLocalCoefficients ()
279  : localKey_(size())
280  {
281  static_assert(k==1, "Only first-order Nédélec local coefficients are implemented.");
282  // Assign all degrees of freedom to edges
283  // TODO: This is correct only for first-order Nédélec elements
284  for (std::size_t i=0; i<size(); i++)
285  localKey_[i] = LocalKey(i,dim-1,0);
286  }
287 
289  std::size_t size() const
290  {
291  static_assert(dim==2 || dim==3, "Nédélec shape functions are implemented only for 2d and 3d simplices.");
292  return (dim==2) ? k * (k+2)
293  : k * (k+2) * (k+3) / 2;
294  }
295 
298  const LocalKey& localKey (std::size_t i) const
299  {
300  return localKey_[i];
301  }
302 
303  private:
304  std::vector<LocalKey> localKey_;
305  };
306 
311  template<class LB>
312  class Nedelec1stKindSimplexLocalInterpolation
313  {
314  static constexpr auto dim = LB::Traits::dimDomain;
315  static constexpr auto size = LB::size();
316 
317  // Number of edges of the reference simplex
318  constexpr static std::size_t numberOfEdges = dim*(dim+1)/2;
319 
320  public:
321 
323  Nedelec1stKindSimplexLocalInterpolation (std::bitset<numberOfEdges> s = 0)
324  {
325  auto refElement = Dune::referenceElement<double,dim>(GeometryTypes::simplex(dim));
326 
327  for (std::size_t i=0; i<numberOfEdges; i++)
328  m_[i] = refElement.position(i,dim-1);
329 
330  for (std::size_t i=0; i<numberOfEdges; i++)
331  {
332  auto vertexIterator = refElement.subEntities(i,dim-1,dim).begin();
333  auto v0 = *vertexIterator;
334  auto v1 = *(++vertexIterator);
335  // By default, edges point from the vertex with the smaller index
336  // to the vertex with the larger index.
337  if (v0>v1)
338  std::swap(v0,v1);
339  edge_[i] = refElement.position(v1,dim) - refElement.position(v0,dim);
340  edge_[i] *= (s[i]) ? -1.0 : 1.0;
341  }
342  }
343 
349  template<typename F, typename C>
350  void interpolate (const F& ff, std::vector<C>& out) const
351  {
352  out.resize(size);
353  auto&& f = Impl::makeFunctionWithCallOperator<typename LB::Traits::DomainType>(ff);
354 
355  for (std::size_t i=0; i<size; i++)
356  {
357  auto y = f(m_[i]);
358  out[i] = 0.0;
359  for (int j=0; j<dim; j++)
360  out[i] += y[j]*edge_[i][j];
361  }
362  }
363 
364  private:
365  // Edge midpoints of the reference simplex
366  std::array<typename LB::Traits::DomainType, numberOfEdges> m_;
367  // Edges of the reference simplex
368  std::array<typename LB::Traits::DomainType, numberOfEdges> edge_;
369  };
370 
371 }
372 
373 
399  template<class D, class R, int dim, int k>
401  {
402  public:
404  Impl::Nedelec1stKindSimplexLocalCoefficients<dim,k>,
405  Impl::Nedelec1stKindSimplexLocalInterpolation<Impl::Nedelec1stKindSimplexLocalBasis<D,R,dim,k> > >;
406 
407  static_assert(dim==2 || dim==3, "Nedelec elements are only implemented for 2d and 3d elements.");
408  static_assert(k==1, "Nedelec elements of the first kind are currently only implemented for order k==1.");
409 
413 
419  Nedelec1stKindSimplexLocalFiniteElement (std::bitset<dim*(dim+1)/2> s) :
420  basis_(s),
421  interpolation_(s)
422  {}
423 
424  const typename Traits::LocalBasisType& localBasis () const
425  {
426  return basis_;
427  }
428 
430  {
431  return coefficients_;
432  }
433 
435  {
436  return interpolation_;
437  }
438 
439  static constexpr unsigned int size ()
440  {
441  return Traits::LocalBasisType::size();
442  }
443 
444  static constexpr GeometryType type ()
445  {
446  return GeometryTypes::simplex(dim);
447  }
448 
449  private:
450  typename Traits::LocalBasisType basis_;
451  typename Traits::LocalCoefficientsType coefficients_;
452  typename Traits::LocalInterpolationType interpolation_;
453  };
454 
455 }
456 
457 #endif
Definition: bdfmcube.hh:16
D DomainType
domain type
Definition: common/localbasis.hh:43
traits helper struct
Definition: localfiniteelementtraits.hh:11
LB LocalBasisType
Definition: localfiniteelementtraits.hh:14
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:18
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:22
Nédélec elements of the first kind for simplex elements.
Definition: nedelec1stkindsimplex.hh:401
static constexpr unsigned int size()
Definition: nedelec1stkindsimplex.hh:439
Nedelec1stKindSimplexLocalFiniteElement()=default
Default constructor.
Nedelec1stKindSimplexLocalFiniteElement(std::bitset< dim *(dim+1)/2 > s)
Constructor with explicitly given edge orientations.
Definition: nedelec1stkindsimplex.hh:419
static constexpr GeometryType type()
Definition: nedelec1stkindsimplex.hh:444
const Traits::LocalBasisType & localBasis() const
Definition: nedelec1stkindsimplex.hh:424
const Traits::LocalInterpolationType & localInterpolation() const
Definition: nedelec1stkindsimplex.hh:434
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: nedelec1stkindsimplex.hh:429