dune-localfunctions  2.8.0
nedelec1stkindcube.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_NEDELEC1STKINDCUBE_HH
4 #define DUNE_LOCALFUNCTIONS_NEDELEC_NEDELEC1STKINDCUBE_HH
5 
6 #include <numeric>
7 
8 #include <dune/common/fmatrix.hh>
9 #include <dune/common/fvector.hh>
10 #include <dune/common/math.hh>
11 
12 #include <dune/geometry/referenceelements.hh>
13 #include <dune/geometry/type.hh>
14 
17 #include <dune/localfunctions/common/localinterpolation.hh> // For deprecated makeFunctionWithCallOperator
19 
20 namespace Dune
21 {
22 namespace Impl
23 {
34  template<class D, class R, int dim, int k>
35  class Nedelec1stKindCubeLocalBasis
36  {
37  // Number of edges of the reference cube
38  constexpr static std::size_t numberOfEdges = power(2,dim-1)*dim;
39 
40  public:
41  using Traits = LocalBasisTraits<D,dim,FieldVector<D,dim>,
42  R,dim,FieldVector<R,dim>,
43  FieldMatrix<R,dim,dim> >;
44 
51  Nedelec1stKindCubeLocalBasis()
52  {
53  std::fill(edgeOrientation_.begin(), edgeOrientation_.end(), 1.0);
54  }
55 
58  Nedelec1stKindCubeLocalBasis(std::bitset<numberOfEdges> edgeOrientation)
59  : Nedelec1stKindCubeLocalBasis()
60  {
61  for (std::size_t i=0; i<edgeOrientation_.size(); i++)
62  edgeOrientation_[i] *= edgeOrientation[i] ? -1.0 : 1.0;
63  }
64 
66  static constexpr unsigned int size()
67  {
68  static_assert(dim==2 || dim==3, "Nedelec shape functions are implemented only for 2d and 3d cubes.");
69  if (dim==2)
70  return 2*k * (k+1);
71  if (dim==3)
72  return 3*k * (k+1) * (k+1);
73  }
74 
80  void evaluateFunction(const typename Traits::DomainType& in,
81  std::vector<typename Traits::RangeType>& out) const
82  {
83  static_assert(k==1, "Evaluating Nédélec shape functions is implemented only for first order.");
84  out.resize(size());
85 
86  if (dim==2)
87  {
88  // First-order Nédélec shape functions on a square are of the form
89  //
90  // (a, b)^T + (c y, d x)^T, a, b, c, d \in R
91  //
92  // The following coefficients create the four basis vectors
93  // that are dual to the edge degrees of freedom:
94  //
95  // a[0] = 0 b[0] = 1 c[0] = 0 d[0] = -1
96  // a[1] = 0 b[1] = 0 c[1] = 0 d[1] = 1
97  // a[2] = 1 b[2] = 0 c[2] = 0 d[2] = -1
98  // a[3] = 0 b[3] = 0 c[3] = 0 d[3] = 1
99 
100  out[0] = { 0, D(1) - in[0]};
101  out[1] = { 0, in[0]};
102  out[2] = { D(1) - in[1], 0};
103  out[3] = { in[1], 0};
104  }
105 
106  if constexpr (dim==3)
107  {
108  // First-order Nédélec shape functions on a cube are of the form
109  //
110  // (e1 yz)
111  // a + b x + c y + d z + (e2 xz) , a, b, c, d \in R^3 and b[0]=c[1]=d[2]=0
112  // (e3 xy)
113  //
114  // The following coefficients create the twelve basis vectors
115  // that are dual to the edge degrees of freedom:
116  //
117  // a[0] = { 0, 0, 1} b[0] = { 0, 0, -1} c[0] = { 0, 0, -1} d[0] = { 0, 0, 0} e[0] = { 0, 0, 1}
118  // a[1] = { 0, 0, 0} b[1] = { 0, 0, 1} c[1] = { 0, 0, 0} d[1] = { 0, 0, 0} e[1] = { 0, 0, -1}
119  // a[2] = { 0, 0, 0} b[2] = { 0, 0, 0} c[2] = { 0, 0, 1} d[2] = { 0, 0, 0} e[2] = { 0, 0, -1}
120  // a[3] = { 0, 0, 0} b[3] = { 0, 0, 0} c[3] = { 0, 0, 0} d[3] = { 0, 0, 0} e[3] = { 0, 0, 1}
121  //
122  // The following implementation uses these values, and simply
123  // skips all the zeros.
124 
125  for (std::size_t i=0; i<out.size(); i++)
126  out[i] = {0,0,0};
127 
128  out[0][2] = { 1 - in[0] - in[1] + in[0]*in[1]};
129  out[1][2] = { in[0] - in[0]*in[1]};
130  out[2][2] = { in[1] - in[0]*in[1]};
131  out[3][2] = { in[0]*in[1]};
132 
133  out[4][1] = { 1 - in[0] - in[2] + in[0]*in[2]};
134  out[5][1] = { in[0] - in[0]*in[2]};
135  out[8][1] = { in[2] - in[0]*in[2]};
136  out[9][1] = { in[0]*in[2]};
137 
138  out[6][0] = { 1 - in[1] - in[2] + in[1]*in[2]};
139  out[7][0] = { in[1] - in[1]*in[2]};
140  out[10][0] = { in[2] - in[1]*in[2]};
141  out[11][0] = { in[1]*in[2]};
142  }
143 
144  for (std::size_t i=0; i<out.size(); i++)
145  out[i] *= edgeOrientation_[i];
146  }
147 
153  void evaluateJacobian(const typename Traits::DomainType& in,
154  std::vector<typename Traits::JacobianType>& out) const
155  {
156  out.resize(size());
157  if (dim==2)
158  {
159  for (std::size_t i=0; i<out.size(); i++)
160  for (std::size_t j=0; j<dim; j++)
161  out[i][j] = { 0, 0};
162 
163  out[0][1] = { -1, 0};
164  out[1][1] = { 1, 0};
165 
166  out[2][0] = { 0, -1};
167  out[3][0] = { 0, 1};
168  }
169  if (dim==3)
170  {
171  for (std::size_t i=0; i<out.size(); i++)
172  for(std::size_t j=0;j<dim; j++)
173  out[i][j] = {0,0,0};
174 
175 
176  out[0][2] = {-1 +in[1], -1 + in[0], 0};
177  out[1][2] = { 1 -in[1], - in[0], 0};
178  out[2][2] = { -in[1], 1 - in[0], 0};
179  out[3][2] = { in[1], in[0], 0};
180 
181  out[4][1] = {-1 +in[2], 0, -1 + in[0]};
182  out[5][1] = { 1 -in[2], 0, - in[0]};
183  out[8][1] = { -in[2], 0, 1 - in[0]};
184  out[9][1] = { in[2], 0, in[0]};
185 
186  out[6][0] = { 0, -1 + in[2], -1 + in[1]};
187  out[7][0] = { 0, 1 - in[2], - in[1]};
188  out[10][0] = { 0, - in[2], 1 - in[1]};
189  out[11][0] = { 0, in[2], in[1]};
190 
191  }
192 
193  for (std::size_t i=0; i<out.size(); i++)
194  out[i] *= edgeOrientation_[i];
195 
196  }
197 
204  void partial(const std::array<unsigned int, dim>& order,
205  const typename Traits::DomainType& in,
206  std::vector<typename Traits::RangeType>& out) const
207  {
208  auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
209  if (totalOrder == 0) {
210  evaluateFunction(in, out);
211  } else if (totalOrder == 1) {
212  auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
213  out.resize(size());
214 
215  if (dim==2)
216  {
217  if (direction==0)
218  {
219  out[0] = { 0, -1};
220  out[1] = { 0, 1};
221  out[2] = { 0, 0};
222  out[3] = { 0, 0};
223  }
224  else
225  {
226  out[0] = { 0, 0};
227  out[1] = { 0, 0};
228  out[2] = { -1, 0};
229  out[3] = { 1, 0};
230  }
231  }
232 
233  if (dim==3)
234  {
235  switch (direction)
236  {
237  case 0:
238  out[0] = { 0, 0, -1 +in[1]};
239  out[1] = { 0, 0, 1 -in[1]};
240  out[2] = { 0, 0, -in[1]};
241  out[3] = { 0, 0, in[1]};
242 
243  out[4] = { 0, -1 +in[2], 0};
244  out[5] = { 0, 1 -in[2], 0};
245  out[8] = { 0, -in[2], 0};
246  out[9] = { 0, in[2], 0};
247 
248  out[6] = {0,0,0};
249  out[7] = {0,0,0};
250  out[10] = {0,0,0};
251  out[11] = {0,0,0};
252  break;
253 
254  case 1:
255  out[0] = { 0, 0, -1 + in[0]};
256  out[1] = { 0, 0, - in[0]};
257  out[2] = { 0, 0, 1 - in[0]};
258  out[3] = { 0, 0, in[0]};
259 
260  out[4] = {0,0,0};
261  out[5] = {0,0,0};
262  out[8] = {0,0,0};
263  out[9] = {0,0,0};
264 
265  out[6] = { -1 + in[2], 0, 0};
266  out[7] = { 1 - in[2], 0, 0};
267  out[10] = { - in[2], 0, 0};
268  out[11] = { in[2], 0, 0};
269  break;
270 
271  case 2:
272  out[0] = {0,0,0};
273  out[1] = {0,0,0};
274  out[2] = {0,0,0};
275  out[3] = {0,0,0};
276 
277  out[4] = { 0, -1 + in[0], 0};
278  out[5] = { 0, - in[0], 0};
279  out[8] = { 0, 1 - in[0], 0};
280  out[9] = { 0, in[0], 0};
281 
282  out[6] = { -1 + in[1], 0, 0};
283  out[7] = { - in[1], 0, 0};
284  out[10] = { 1 - in[1], 0, 0};
285  out[11] = { in[1], 0, 0};
286  break;
287  }
288  }
289 
290  for (std::size_t i=0; i<out.size(); i++)
291  out[i] *= edgeOrientation_[i];
292 
293  } else if (totalOrder == 2) {
294  out.resize(size());
295 
296  if (dim==2)
297  for (std::size_t i = 0; i < size(); ++i)
298  for (std::size_t j = 0; j < dim; ++j)
299  out[i][j] = 0;
300 
301  if (dim==3)
302  {
303  for(size_t i=0; i<out.size(); i++)
304  out[i] = { 0, 0, 0};
305 
306  //case (1,1,0):
307  if( order[0] == 1 and order[1]==1)
308  {
309  out[0] = { 0, 0, 1};
310  out[1] = { 0, 0, -1};
311  out[2] = { 0, 0, -1};
312  out[3] = { 0, 0, 1};
313  }
314 
315  //case (1,0,1):
316  if( order[0] == 1 and order[2]==1)
317  {
318  out[4] = { 0, 1, 0};
319  out[5] = { 0, -1, 0};
320  out[8] = { 0, -1, 0};
321  out[9] = { 0, 1, 0};
322  }
323 
324  //case (0,1,1):
325  if( order[1] == 1 and order[2]==1)
326  {
327  out[6] = { 1, 0, 0};
328  out[7] = { -1, 0, 0};
329  out[10] = { -1, 0, 0};
330  out[11] = { 1, 0, 0};
331  }
332 
333  for (std::size_t i=0; i<out.size(); i++)
334  out[i] *= edgeOrientation_[i];
335  }
336 
337 
338  }else {
339  out.resize(size());
340  for (std::size_t i = 0; i < size(); ++i)
341  for (std::size_t j = 0; j < dim; ++j)
342  out[i][j] = 0;
343  }
344 
345  }
346 
348  unsigned int order() const
349  {
350  if (dim==2)
351  return 2*k-1;
352  if (dim==3)
353  return 3*k-1;
354  }
355 
356  private:
357 
358  // Orientations of the cube edges
359  std::array<R,numberOfEdges> edgeOrientation_;
360  };
361 
362 
367  template <int dim, int k>
368  class Nedelec1stKindCubeLocalCoefficients
369  {
370  public:
372  Nedelec1stKindCubeLocalCoefficients ()
373  : localKey_(size())
374  {
375  static_assert(k==1, "Only first-order Nédélec local coefficients are implemented.");
376  // Assign all degrees of freedom to edges
377  // TODO: This is correct only for first-order Nédélec elements
378  for (std::size_t i=0; i<size(); i++)
379  localKey_[i] = LocalKey(i,dim-1,0);
380  }
381 
383  std::size_t size() const
384  {
385  static_assert(dim==2 || dim==3, "Nédélec shape functions are implemented only for 2d and 3d cubes.");
386  return (dim==2) ? 2*k * (k+1)
387  : 3*k * (k+1) * (k+1);
388  }
389 
392  const LocalKey& localKey (std::size_t i) const
393  {
394  return localKey_[i];
395  }
396 
397  private:
398  std::vector<LocalKey> localKey_;
399  };
400 
405  template<class LB>
406  class Nedelec1stKindCubeLocalInterpolation
407  {
408  static constexpr auto dim = LB::Traits::dimDomain;
409  static constexpr auto size = LB::size();
410 
411  // Number of edges of the reference cube
412  constexpr static std::size_t numberOfEdges = power(2,dim-1)*dim;
413 
414  public:
415 
417  Nedelec1stKindCubeLocalInterpolation (std::bitset<numberOfEdges> s = 0)
418  {
419  auto refElement = Dune::referenceElement<double,dim>(GeometryTypes::cube(dim));
420 
421  for (std::size_t i=0; i<numberOfEdges; i++)
422  m_[i] = refElement.position(i,dim-1);
423 
424  for (std::size_t i=0; i<numberOfEdges; i++)
425  {
426  auto vertexIterator = refElement.subEntities(i,dim-1,dim).begin();
427  auto v0 = *vertexIterator;
428  auto v1 = *(++vertexIterator);
429 
430  // By default, edges point from the vertex with the smaller index
431  // to the vertex with the larger index.
432  if (v0>v1)
433  std::swap(v0,v1);
434  edge_[i] = refElement.position(v1,dim) - refElement.position(v0,dim);
435  edge_[i] *= (s[i]) ? -1.0 : 1.0;
436  }
437  }
438 
444  template<typename F, typename C>
445  void interpolate (const F& ff, std::vector<C>& out) const
446  {
447  out.resize(size);
448  auto&& f = Impl::makeFunctionWithCallOperator<typename LB::Traits::DomainType>(ff);
449 
450  for (std::size_t i=0; i<size; i++)
451  {
452  auto y = f(m_[i]);
453  out[i] = 0.0;
454  for (int j=0; j<dim; j++)
455  out[i] += y[j]*edge_[i][j];
456  }
457  }
458 
459  private:
460  // Edge midpoints of the reference cube
461  std::array<typename LB::Traits::DomainType, numberOfEdges> m_;
462  // Edges of the reference cube
463  std::array<typename LB::Traits::DomainType, numberOfEdges> edge_;
464  };
465 
466 }
467 
468 
492  template<class D, class R, int dim, int k>
494  {
495  public:
497  Impl::Nedelec1stKindCubeLocalCoefficients<dim,k>,
498  Impl::Nedelec1stKindCubeLocalInterpolation<Impl::Nedelec1stKindCubeLocalBasis<D,R,dim,k> > >;
499 
500  static_assert(dim==2 || dim==3, "Nedelec elements are only implemented for 2d and 3d elements.");
501  static_assert(k==1, "Nedelec elements of the first kind are currently only implemented for order k==1.");
502 
506 
512  Nedelec1stKindCubeLocalFiniteElement (std::bitset<power(2,dim-1)*dim> s) :
513  basis_(s),
514  interpolation_(s)
515  {}
516 
517  const typename Traits::LocalBasisType& localBasis () const
518  {
519  return basis_;
520  }
521 
523  {
524  return coefficients_;
525  }
526 
528  {
529  return interpolation_;
530  }
531 
532  static constexpr unsigned int size ()
533  {
534  return Traits::LocalBasisType::size();
535  }
536 
537  static constexpr GeometryType type ()
538  {
539  return GeometryTypes::cube(dim);
540  }
541 
542  private:
543  typename Traits::LocalBasisType basis_;
544  typename Traits::LocalCoefficientsType coefficients_;
545  typename Traits::LocalInterpolationType interpolation_;
546  };
547 
548 }
549 
550 #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 cube elements.
Definition: nedelec1stkindcube.hh:494
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: nedelec1stkindcube.hh:522
static constexpr unsigned int size()
Definition: nedelec1stkindcube.hh:532
Nedelec1stKindCubeLocalFiniteElement()=default
Default constructor.
static constexpr GeometryType type()
Definition: nedelec1stkindcube.hh:537
const Traits::LocalBasisType & localBasis() const
Definition: nedelec1stkindcube.hh:517
const Traits::LocalInterpolationType & localInterpolation() const
Definition: nedelec1stkindcube.hh:527
Nedelec1stKindCubeLocalFiniteElement(std::bitset< power(2, dim-1) *dim > s)
Constructor with explicitly given edge orientations.
Definition: nedelec1stkindcube.hh:512