Dune Core Modules (2.7.1)

lagrangesimplex.hh
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_LAGRANGE_LAGRANGESIMPLEX_HH
4 #define DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGESIMPLEX_HH
5 
6 #include <array>
7 #include <numeric>
8 
9 #include <dune/common/deprecated.hh>
10 #include <dune/common/fmatrix.hh>
11 #include <dune/common/fvector.hh>
12 #include <dune/common/math.hh>
13 
14 #include <dune/geometry/referenceelements.hh>
15 
16 #include <dune/localfunctions/common/localbasis.hh>
17 #include <dune/localfunctions/common/localfiniteelementtraits.hh>
18 #include <dune/localfunctions/common/localinterpolation.hh>
19 #include <dune/localfunctions/common/localkey.hh>
20 
21 namespace Dune { namespace Impl
22 {
33  template<class D, class R, unsigned int dim, unsigned int k>
34  class LagrangeSimplexLocalBasis
35  {
36  public:
37  using Traits = LocalBasisTraits<D,dim,FieldVector<D,dim>,R,1,FieldVector<R,1>,FieldMatrix<R,1,dim> >;
38 
43  static constexpr unsigned int size ()
44  {
45  return binomial(k+dim,dim);
46  }
47 
49  void evaluateFunction(const typename Traits::DomainType& x,
50  std::vector<typename Traits::RangeType>& out) const
51  {
52  out.resize(size());
53 
54  // Specialization for zero-order case
55  if (k==0)
56  {
57  out[0] = 1;
58  return;
59  }
60 
61  // Specialization for first-order case
62  if (k==1)
63  {
64  out[0] = 1.0;
65  for (size_t i=0; i<dim; i++)
66  {
67  out[0] -= x[i];
68  out[i+1] = x[i];
69  }
70  return;
71  }
72 
73  assert(k>=2);
74 
75  auto lagrangeNode = [](unsigned int i) { return ((D)i)/k; };
76 
77  if (dim==1)
78  {
79  for (unsigned int i=0; i<size(); i++)
80  {
81  out[i] = 1.0;
82  for (unsigned int alpha=0; alpha<i; alpha++)
83  out[i] *= (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
84  for (unsigned int gamma=i+1; gamma<=k; gamma++)
85  out[i] *= (x[0]-lagrangeNode(gamma))/(lagrangeNode(i)-lagrangeNode(gamma));
86  }
87  return;
88  }
89 
90  if (dim==2)
91  {
92  int n=0;
93  for (unsigned int j=0; j<=k; j++)
94  for (unsigned int i=0; i<=k-j; i++)
95  {
96  out[n] = 1.0;
97  for (unsigned int alpha=0; alpha<i; alpha++)
98  out[n] *= (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
99  for (unsigned int beta=0; beta<j; beta++)
100  out[n] *= (x[1]-lagrangeNode(beta))/(lagrangeNode(j)-lagrangeNode(beta));
101  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
102  out[n] *= (lagrangeNode(gamma)-x[0]-x[1])/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
103  n++;
104  }
105 
106  return;
107  }
108 
109  if (dim!=3)
110  DUNE_THROW(NotImplemented, "LagrangeSimplexLocalBasis for k>=2 only implemented for dim==1 or dim==3");
111 
112  typename Traits::DomainType kx = x;
113  kx *= k;
114  unsigned int n = 0;
115  unsigned int i[4];
116  R factor[4];
117  for (i[2] = 0; i[2] <= k; ++i[2])
118  {
119  factor[2] = 1.0;
120  for (unsigned int j = 0; j < i[2]; ++j)
121  factor[2] *= (kx[2]-j) / (i[2]-j);
122  for (i[1] = 0; i[1] <= k - i[2]; ++i[1])
123  {
124  factor[1] = 1.0;
125  for (unsigned int j = 0; j < i[1]; ++j)
126  factor[1] *= (kx[1]-j) / (i[1]-j);
127  for (i[0] = 0; i[0] <= k - i[1] - i[2]; ++i[0])
128  {
129  factor[0] = 1.0;
130  for (unsigned int j = 0; j < i[0]; ++j)
131  factor[0] *= (kx[0]-j) / (i[0]-j);
132  i[3] = k - i[0] - i[1] - i[2];
133  D kx3 = k - kx[0] - kx[1] - kx[2];
134  factor[3] = 1.0;
135  for (unsigned int j = 0; j < i[3]; ++j)
136  factor[3] *= (kx3-j) / (i[3]-j);
137  out[n++] = factor[0] * factor[1] * factor[2] * factor[3];
138  }
139  }
140  }
141  }
142 
148  void evaluateJacobian(const typename Traits::DomainType& x,
149  std::vector<typename Traits::JacobianType>& out) const
150  {
151  out.resize(size());
152 
153  // Specialization for k==0
154  if (k==0)
155  {
156  std::fill(out[0][0].begin(), out[0][0].end(), 0);
157  return;
158  }
159 
160  // Specialization for k==1
161  if (k==1)
162  {
163  std::fill(out[0][0].begin(), out[0][0].end(), -1);
164 
165  for (unsigned int i=0; i<dim; i++)
166  for (unsigned int j=0; j<dim; j++)
167  out[i+1][0][j] = (i==j);
168 
169  return;
170  }
171 
172  auto lagrangeNode = [](unsigned int i) { return ((D)i)/k; };
173 
174  // Specialization for dim==1
175  if (dim==1)
176  {
177  for (unsigned int i=0; i<=k; i++)
178  {
179  // x_0 derivative
180  out[i][0][0] = 0.0;
181  R factor=1.0;
182  for (unsigned int a=0; a<i; a++)
183  {
184  R product=factor;
185  for (unsigned int alpha=0; alpha<i; alpha++)
186  product *= (alpha==a) ? 1.0/(lagrangeNode(i)-lagrangeNode(alpha))
187  : (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
188  for (unsigned int gamma=i+1; gamma<=k; gamma++)
189  product *= (lagrangeNode(gamma)-x[0])/(lagrangeNode(gamma)-lagrangeNode(i));
190  out[i][0][0] += product;
191  }
192  for (unsigned int c=i+1; c<=k; c++)
193  {
194  R product=factor;
195  for (unsigned int alpha=0; alpha<i; alpha++)
196  product *= (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
197  for (unsigned int gamma=i+1; gamma<=k; gamma++)
198  product *= (gamma==c) ? -1.0/(lagrangeNode(gamma)-lagrangeNode(i))
199  : (lagrangeNode(gamma)-x[0])/(lagrangeNode(gamma)-lagrangeNode(i));
200  out[i][0][0] += product;
201  }
202  }
203  return;
204  }
205 
206  if (dim==2)
207  {
208  int n=0;
209  for (unsigned int j=0; j<=k; j++)
210  for (unsigned int i=0; i<=k-j; i++)
211  {
212  // x_0 derivative
213  out[n][0][0] = 0.0;
214  R factor=1.0;
215  for (unsigned int beta=0; beta<j; beta++)
216  factor *= (x[1]-lagrangeNode(beta))/(lagrangeNode(j)-lagrangeNode(beta));
217  for (unsigned int a=0; a<i; a++)
218  {
219  R product=factor;
220  for (unsigned int alpha=0; alpha<i; alpha++)
221  if (alpha==a)
222  product *= D(1)/(lagrangeNode(i)-lagrangeNode(alpha));
223  else
224  product *= (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
225  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
226  product *= (lagrangeNode(gamma)-x[0]-x[1])/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
227  out[n][0][0] += product;
228  }
229  for (unsigned int c=i+j+1; c<=k; c++)
230  {
231  R product=factor;
232  for (unsigned int alpha=0; alpha<i; alpha++)
233  product *= (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
234  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
235  if (gamma==c)
236  product *= -D(1)/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
237  else
238  product *= (lagrangeNode(gamma)-x[0]-x[1])/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
239  out[n][0][0] += product;
240  }
241 
242  // x_1 derivative
243  out[n][0][1] = 0.0;
244  factor = 1.0;
245  for (unsigned int alpha=0; alpha<i; alpha++)
246  factor *= (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
247  for (unsigned int b=0; b<j; b++)
248  {
249  R product=factor;
250  for (unsigned int beta=0; beta<j; beta++)
251  if (beta==b)
252  product *= D(1)/(lagrangeNode(j)-lagrangeNode(beta));
253  else
254  product *= (x[1]-lagrangeNode(beta))/(lagrangeNode(j)-lagrangeNode(beta));
255  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
256  product *= (lagrangeNode(gamma)-x[0]-x[1])/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
257  out[n][0][1] += product;
258  }
259  for (unsigned int c=i+j+1; c<=k; c++)
260  {
261  R product=factor;
262  for (unsigned int beta=0; beta<j; beta++)
263  product *= (x[1]-lagrangeNode(beta))/(lagrangeNode(j)-lagrangeNode(beta));
264  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
265  if (gamma==c)
266  product *= -D(1)/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
267  else
268  product *= (lagrangeNode(gamma)-x[0]-x[1])/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
269  out[n][0][1] += product;
270  }
271 
272  n++;
273  }
274 
275  return;
276  }
277 
278  if (dim!=3)
279  DUNE_THROW(NotImplemented, "LagrangeSimplexLocalBasis only implemented for dim==3!");
280 
281  // Specialization for arbitrary order and dim==3
282  typename Traits::DomainType kx = x;
283  kx *= k;
284  unsigned int n = 0;
285  unsigned int i[4];
286  R factor[4];
287  for (i[2] = 0; i[2] <= k; ++i[2])
288  {
289  factor[2] = 1.0;
290  for (unsigned int j = 0; j < i[2]; ++j)
291  factor[2] *= (kx[2]-j) / (i[2]-j);
292  for (i[1] = 0; i[1] <= k - i[2]; ++i[1])
293  {
294  factor[1] = 1.0;
295  for (unsigned int j = 0; j < i[1]; ++j)
296  factor[1] *= (kx[1]-j) / (i[1]-j);
297  for (i[0] = 0; i[0] <= k - i[1] - i[2]; ++i[0])
298  {
299  factor[0] = 1.0;
300  for (unsigned int j = 0; j < i[0]; ++j)
301  factor[0] *= (kx[0]-j) / (i[0]-j);
302  i[3] = k - i[0] - i[1] - i[2];
303  D kx3 = k - kx[0] - kx[1] - kx[2];
304  R sum3 = 0.0;
305  factor[3] = 1.0;
306  for (unsigned int j = 0; j < i[3]; ++j)
307  factor[3] /= i[3] - j;
308  R prod_all = factor[0] * factor[1] * factor[2] * factor[3];
309  for (unsigned int j = 0; j < i[3]; ++j)
310  {
311  R prod = prod_all;
312  for (unsigned int l = 0; l < i[3]; ++l)
313  if (j == l)
314  prod *= -R(k);
315  else
316  prod *= kx3 - l;
317  sum3 += prod;
318  }
319  for (unsigned int j = 0; j < i[3]; ++j)
320  factor[3] *= kx3 - j;
321  for (unsigned int m = 0; m < 3; ++m)
322  {
323  out[n][0][m] = sum3;
324  for (unsigned int j = 0; j < i[m]; ++j)
325  {
326  R prod = factor[3];
327  for (unsigned int p = 0; p < 3; ++p)
328  {
329  if (m == p)
330  for (unsigned int l = 0; l < i[p]; ++l)
331  prod *= (j==l) ? R(k) / (i[p]-l) : (kx[p]-l) / (i[p]-l);
332  else
333  prod *= factor[p];
334  }
335  out[n][0][m] += prod;
336  }
337  }
338  n++;
339  }
340  }
341  }
342  }
343 
350  void partial(const std::array<unsigned int,dim>& order,
351  const typename Traits::DomainType& in,
352  std::vector<typename Traits::RangeType>& out) const
353  {
354  auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
355 
356  out.resize(size());
357 
358  if (totalOrder == 0) {
359  evaluateFunction(in, out);
360  return;
361  }
362 
363  if (k==0)
364  {
365  out[0] = 0;
366  return;
367  }
368 
369  if (k==1)
370  {
371  if (totalOrder==1)
372  {
373  auto direction = std::find(order.begin(), order.end(), 1);
374 
375  out[0] = -1;
376  for (unsigned int i=0; i<dim; i++)
377  out[i+1] = (i==(direction-order.begin()));
378  }
379  else // all higher order derivatives are zero
380  std::fill(out.begin(), out.end(), 0);
381  return;
382  }
383 
384  if (dim==2)
385  {
386  auto lagrangeNode = [](unsigned int i) { return ((D)i)/k; };
387 
388  // Helper method: Return a single Lagrangian factor of l_ij evaluated at x
389  auto lagrangianFactor = [&lagrangeNode]
390  (const int no, const int i, const int j, const typename Traits::DomainType& x)
391  -> typename Traits::RangeType
392  {
393  if ( no < i)
394  return (x[0]-lagrangeNode(no))/(lagrangeNode(i)-lagrangeNode(no));
395  if (no < i+j)
396  return (x[1]-lagrangeNode(no-i))/(lagrangeNode(j)-lagrangeNode(no-i));
397  return (lagrangeNode(no+1)-x[0]-x[1])/(lagrangeNode(no+1)-lagrangeNode(i)-lagrangeNode(j));
398  };
399 
400  // Helper method: Return the derivative of a single Lagrangian factor of l_ij evaluated at x
401  // direction: Derive in x-direction if this is 0, otherwise derive in y direction
402  auto lagrangianFactorDerivative = [&lagrangeNode]
403  (const int direction, const int no, const int i, const int j, const typename Traits::DomainType& x)
404  -> typename Traits::RangeType
405  {
406  using T = typename Traits::RangeType;
407  if ( no < i)
408  return (direction == 0) ? T(1.0/(lagrangeNode(i)-lagrangeNode(no))) : T(0);
409 
410  if (no < i+j)
411  return (direction == 0) ? T(0) : T(1.0/(lagrangeNode(j)-lagrangeNode(no-i)));
412 
413  return -1.0/(lagrangeNode(no+1)-lagrangeNode(i)-lagrangeNode(j));
414  };
415 
416  if (totalOrder==1)
417  {
418  int direction = std::find(order.begin(), order.end(), 1)-order.begin();
419 
420  int n=0;
421  for (unsigned int j=0; j<=k; j++)
422  {
423  for (unsigned int i=0; i<=k-j; i++, n++)
424  {
425  out[n] = 0.0;
426  for (unsigned int no1=0; no1 < k; no1++)
427  {
428  R factor = lagrangianFactorDerivative(direction, no1, i, j, in);
429  for (unsigned int no2=0; no2 < k; no2++)
430  if (no1 != no2)
431  factor *= lagrangianFactor(no2, i, j, in);
432 
433  out[n] += factor;
434  }
435  }
436  }
437  return;
438  }
439 
440  if (totalOrder==2)
441  {
442  std::array<int,2> directions;
443  unsigned int counter = 0;
444  auto nonconstOrder = order; // need a copy that I can modify
445  for (int i=0; i<2; i++)
446  {
447  while (nonconstOrder[i])
448  {
449  directions[counter++] = i;
450  nonconstOrder[i]--;
451  }
452  }
453 
454  //f = prod_{i} f_i -> dxa dxb f = sum_{i} {dxa f_i sum_{k \neq i} dxb f_k prod_{l \neq k,i} f_l
455  int n=0;
456  for (unsigned int j=0; j<=k; j++)
457  {
458  for (unsigned int i=0; i<=k-j; i++, n++)
459  {
460  R res = 0.0;
461 
462  for (unsigned int no1=0; no1 < k; no1++)
463  {
464  R factor1 = lagrangianFactorDerivative(directions[0], no1, i, j, in);
465  for (unsigned int no2=0; no2 < k; no2++)
466  {
467  if (no1 == no2)
468  continue;
469  R factor2 = factor1*lagrangianFactorDerivative(directions[1], no2, i, j, in);
470  for (unsigned int no3=0; no3 < k; no3++)
471  {
472  if (no3 == no1 || no3 == no2)
473  continue;
474  factor2 *= lagrangianFactor(no3, i, j, in);
475  }
476  res += factor2;
477  }
478  }
479  out[n] = res;
480  }
481  }
482 
483  return;
484  } // totalOrder==2
485 
486  } // dim==2
487 
488  DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
489  }
490 
492  static constexpr unsigned int order ()
493  {
494  return k;
495  }
496  };
497 
503  template<unsigned int dim, unsigned int k>
504  class LagrangeSimplexLocalCoefficients
505  {
506  public:
508  LagrangeSimplexLocalCoefficients ()
509  : localKeys_(size())
510  {
511  if (k==0)
512  {
513  localKeys_[0] = LocalKey(0,0,0);
514  return;
515  }
516 
517  if (k==1)
518  {
519  for (std::size_t i=0; i<size(); i++)
520  localKeys_[i] = LocalKey(i,dim,0);
521  return;
522  }
523 
524  if (dim==1)
525  {
526  // Order is at least 2 here
527  localKeys_[0] = LocalKey(0,1,0); // vertex dof
528  for (unsigned int i=1; i<k; i++)
529  localKeys_[i] = LocalKey(0,0,i-1); // element dofs
530  localKeys_.back() = LocalKey(1,1,0); // vertex dof
531  return;
532  }
533 
534  if (dim==2)
535  {
536  int n=0;
537  int c=0;
538  for (unsigned int j=0; j<=k; j++)
539  for (unsigned int i=0; i<=k-j; i++)
540  {
541  if (i==0 && j==0)
542  {
543  localKeys_[n++] = LocalKey(0,2,0);
544  continue;
545  }
546  if (i==k && j==0)
547  {
548  localKeys_[n++] = LocalKey(1,2,0);
549  continue;
550  }
551  if (i==0 && j==k)
552  {
553  localKeys_[n++] = LocalKey(2,2,0);
554  continue;
555  }
556  if (j==0)
557  {
558  localKeys_[n++] = LocalKey(0,1,i-1);
559  continue;
560  }
561  if (i==0)
562  {
563  localKeys_[n++] = LocalKey(1,1,j-1);
564  continue;
565  }
566  if (i+j==k)
567  {
568  localKeys_[n++] = LocalKey(2,1,j-1);
569  continue;
570  }
571  localKeys_[n++] = LocalKey(0,0,c++);
572  }
573  return;
574  }
575 
576  if (dim==3)
577  {
578  std::array<unsigned int, dim+1> vertexMap;
579  for (unsigned int i=0; i<=dim; i++)
580  vertexMap[i] = i;
581  generateLocalKeys(vertexMap);
582  return;
583  }
584  DUNE_THROW(NotImplemented, "LagrangeSimplexLocalCoefficients only implemented for k<=1 or dim<=3!");
585  }
586 
593  LagrangeSimplexLocalCoefficients (const std::array<unsigned int, dim+1> vertexMap)
594  : localKeys_(size())
595  {
596  if (dim!=2 && dim!=3)
597  DUNE_THROW(NotImplemented, "LagrangeSimplexLocalCoefficients only implemented for dim==2 and dim==3!");
598 
599  generateLocalKeys(vertexMap);
600  }
601 
602 
603  template<class VertexMap>
604  LagrangeSimplexLocalCoefficients(const VertexMap &vertexmap)
605  : localKeys_(size())
606  {
607  if (dim!=2 && dim!=3)
608  DUNE_THROW(NotImplemented, "LagrangeSimplexLocalCoefficients only implemented for dim==2 and dim==3!");
609 
610  std::array<unsigned int, dim+1> vertexmap_array;
611  std::copy(vertexmap, vertexmap + dim + 1, vertexmap_array.begin());
612  generateLocalKeys(vertexmap_array);
613  }
614 
616  static constexpr std::size_t size ()
617  {
618  return binomial(k+dim,dim);
619  }
620 
622  const LocalKey& localKey (std::size_t i) const
623  {
624  return localKeys_[i];
625  }
626 
627  private:
628  std::vector<LocalKey> localKeys_;
629 
630  void generateLocalKeys(const std::array<unsigned int, dim+1> vertexMap)
631  {
632  if (k==0)
633  {
634  localKeys_[0] = LocalKey(0,0,0);
635  return;
636  }
637 
638  if (dim==2)
639  {
640  // Create default assignment
641  int n=0;
642  int c=0;
643  for (unsigned int j=0; j<=k; j++)
644  for (unsigned int i=0; i<=k-j; i++)
645  {
646  if (i==0 && j==0)
647  {
648  localKeys_[n++] = LocalKey(0,2,0);
649  continue;
650  }
651  if (i==k && j==0)
652  {
653  localKeys_[n++] = LocalKey(1,2,0);
654  continue;
655  }
656  if (i==0 && j==k)
657  {
658  localKeys_[n++] = LocalKey(2,2,0);
659  continue;
660  }
661  if (j==0)
662  {
663  localKeys_[n++] = LocalKey(0,1,i-1);
664  continue;
665  }
666  if (i==0)
667  {
668  localKeys_[n++] = LocalKey(1,1,j-1);
669  continue;
670  }
671  if (i+j==k)
672  {
673  localKeys_[n++] = LocalKey(2,1,j-1);
674  continue;
675  }
676  localKeys_[n++] = LocalKey(0,0,c++);
677  }
678 
679  // Flip edge orientations, if requested
680  bool flip[3];
681  flip[0] = vertexMap[0] > vertexMap[1];
682  flip[1] = vertexMap[0] > vertexMap[2];
683  flip[2] = vertexMap[1] > vertexMap[2];
684  for (std::size_t i=0; i<size(); i++)
685  if (localKeys_[i].codim()==1 && flip[localKeys_[i].subEntity()])
686  localKeys_[i].index(k-2-localKeys_[i].index());
687 
688  return;
689  }
690 
691  if (dim!=3)
692  DUNE_THROW(NotImplemented, "LagrangeSimplexLocalCoefficients only implemented for dim==3!");
693 
694  unsigned int subindex[16];
695  unsigned int codim_count[4] = {0};
696  for (unsigned int m = 1; m < 16; ++m)
697  {
698  unsigned int codim = !(m&1) + !(m&2) + !(m&4) + !(m&8);
699  subindex[m] = codim_count[codim]++;
700  }
701 
702  int a1 = (3*k + 12)*k + 11;
703  int a2 = -3*k - 6;
704  unsigned int dof_count[16] = {0};
705  unsigned int i[4];
706  for (i[3] = 0; i[3] <= k; ++i[3])
707  for (i[2] = 0; i[2] <= k - i[3]; ++i[2])
708  for (i[1] = 0; i[1] <= k - i[2] - i[3]; ++i[1])
709  {
710  i[0] = k - i[1] - i[2] - i[3];
711  unsigned int j[4];
712  unsigned int entity = 0;
713  unsigned int codim = 0;
714  for (unsigned int m = 0; m < 4; ++m)
715  {
716  j[m] = i[vertexMap[m]];
717  entity += !!j[m] << m;
718  codim += !j[m];
719  }
720  int local_index = j[3]*(a1 + (a2 + j[3])*j[3])/6
721  + j[2]*(2*(k - j[3]) + 3 - j[2])/2 + j[1];
722  localKeys_[local_index] = LocalKey(subindex[entity], codim, dof_count[entity]++);
723  }
724  }
725  };
726 
731  template<class LocalBasis>
732  class LagrangeSimplexLocalInterpolation
733  {
734  static const int kdiv = (LocalBasis::order() == 0 ? 1 : LocalBasis::order());
735  public:
736 
744  template<typename F, typename C>
745  void interpolate (const F& ff, std::vector<C>& out) const
746  {
747  constexpr auto dim = LocalBasis::Traits::dimDomain;
748  constexpr auto k = LocalBasis::order();
749  using D = typename LocalBasis::Traits::DomainFieldType;
750 
751  typename LocalBasis::Traits::DomainType x;
752  auto&& f = Impl::makeFunctionWithCallOperator<typename LocalBasis::Traits::DomainType>(ff);
753 
754  out.resize(LocalBasis::size());
755 
756  // Specialization for zero-order case
757  if (k==0)
758  {
759  auto center = ReferenceElements<D,dim>::simplex().position(0,0);
760  out[0] = f(center);
761  return;
762  }
763 
764  // Specialization for first-order case
765  if (k==1)
766  {
767  // vertex 0
768  std::fill(x.begin(), x.end(), 0);
769  out[0] = f(x);
770 
771  // remaining vertices
772  for (int i=0; i<dim; i++)
773  {
774  for (int j=0; j<dim; j++)
775  x[j] = (i==j);
776 
777  out[i+1] = f(x);
778  }
779  return;
780  }
781 
782  if (dim==1)
783  {
784  for (unsigned int i=0; i<k+1; i++)
785  {
786  x[0] = ((D)i)/k;
787  out[i] = f(x);
788  }
789  return;
790  }
791 
792  if (dim==2)
793  {
794  int n=0;
795  for (unsigned int j=0; j<=k; j++)
796  for (unsigned int i=0; i<=k-j; i++)
797  {
798  x = { ((D)i)/k, ((D)j)/k };
799  out[n] = f(x);
800  n++;
801  }
802  return;
803  }
804 
805  if (dim!=3)
806  DUNE_THROW(NotImplemented, "LagrangeSimplexLocalInterpolation only implemented for dim<=3!");
807 
808  int n=0;
809  for (int i2 = 0; i2 <= (int)k; i2++)
810  for (int i1 = 0; i1 <= (int)k-i2; i1++)
811  for (int i0 = 0; i0 <= (int)k-i1-i2; i0++)
812  {
813  x[0] = ((D)i0)/((D)kdiv);
814  x[1] = ((D)i1)/((D)kdiv);
815  x[2] = ((D)i2)/((D)kdiv);
816  out[n] = f(x);
817  n++;
818  }
819  }
820 
821  };
822 
823 } } // namespace Dune::Impl
824 
825 namespace Dune
826 {
834  template<class D, class R, int d, int k>
835  class LagrangeSimplexLocalFiniteElement
836  {
837  public:
840  using Traits = LocalFiniteElementTraits<Impl::LagrangeSimplexLocalBasis<D,R,d,k>,
841  Impl::LagrangeSimplexLocalCoefficients<d,k>,
842  Impl::LagrangeSimplexLocalInterpolation<Impl::LagrangeSimplexLocalBasis<D,R,d,k> > >;
843 
845  LagrangeSimplexLocalFiniteElement() {}
846 
851  template<typename VertexMap>
852  LagrangeSimplexLocalFiniteElement(const VertexMap& vertexmap)
853  : coefficients_(vertexmap)
854  {}
855 
858  const typename Traits::LocalBasisType& localBasis () const
859  {
860  return basis_;
861  }
862 
865  const typename Traits::LocalCoefficientsType& localCoefficients () const
866  {
867  return coefficients_;
868  }
869 
872  const typename Traits::LocalInterpolationType& localInterpolation () const
873  {
874  return interpolation_;
875  }
876 
878  static constexpr std::size_t size ()
879  {
880  return Impl::LagrangeSimplexLocalBasis<D,R,d,k>::size();
881  }
882 
885  static constexpr GeometryType type ()
886  {
887  return GeometryTypes::simplex(d);
888  }
889 
890  private:
891  Impl::LagrangeSimplexLocalBasis<D,R,d,k> basis_;
892  Impl::LagrangeSimplexLocalCoefficients<d,k> coefficients_;
893  Impl::LagrangeSimplexLocalInterpolation<Impl::LagrangeSimplexLocalBasis<D,R,d,k> > interpolation_;
894  };
895 
896 } // namespace Dune
897 
898 #endif // DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGESIMPLEX_HH
Dune::GeometryType
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:279
Dune::LagrangeSimplexLocalFiniteElement
Lagrange finite element for simplices with arbitrary compile-time dimension and polynomial order.
Definition: lagrangesimplex.hh:836
Dune::LagrangeSimplexLocalFiniteElement::localBasis
const Traits::LocalBasisType & localBasis() const
Returns the local basis, i.e., the set of shape functions.
Definition: lagrangesimplex.hh:858
Dune::LagrangeSimplexLocalFiniteElement::localInterpolation
const Traits::LocalInterpolationType & localInterpolation() const
Returns object that evaluates degrees of freedom.
Definition: lagrangesimplex.hh:872
Dune::LagrangeSimplexLocalFiniteElement::localCoefficients
const Traits::LocalCoefficientsType & localCoefficients() const
Returns the assignment of the degrees of freedom to the element subentities.
Definition: lagrangesimplex.hh:865
Dune::LagrangeSimplexLocalFiniteElement::LagrangeSimplexLocalFiniteElement
LagrangeSimplexLocalFiniteElement()
Definition: lagrangesimplex.hh:845
Dune::LagrangeSimplexLocalFiniteElement::size
static constexpr std::size_t size()
The number of shape functions.
Definition: lagrangesimplex.hh:878
Dune::LagrangeSimplexLocalFiniteElement::LagrangeSimplexLocalFiniteElement
LagrangeSimplexLocalFiniteElement(const VertexMap &vertexmap)
Definition: lagrangesimplex.hh:852
Dune::LagrangeSimplexLocalFiniteElement::type
static constexpr GeometryType type()
The reference element that the local finite element is defined on.
Definition: lagrangesimplex.hh:885
deprecated.hh
Definition of the DUNE_DEPRECATED macro for the case that config.h is not available.
fmatrix.hh
Implements a matrix constructed from a given type representing a field and compile-time given number ...
fvector.hh
Implements a vector constructed from a given type representing a field and a compile-time given size.
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune::GeometryTypes::simplex
constexpr GeometryType simplex(unsigned int dim)
Returns a GeometryType representing a simplex of dimension dim.
Definition: type.hh:766
Dune::Hybrid::accumulate
T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:290
math.hh
Some useful basic math stuff.
Dune
Dune namespace.
Definition: alignedallocator.hh:14
Dune::binomial
constexpr static T binomial(const T &n, const T &k) noexcept
calculate the binomial coefficient n over k as a constexpr
Definition: math.hh:128
Dune::Geo::ReferenceElements::simplex
static const ReferenceElement & simplex()
get simplex reference elements
Definition: referenceelements.hh:202
Dune::LocalBasisTraits::DomainType
D DomainType
domain type
Definition: localbasis.hh:43
Dune::LocalBasisTraits::RangeType
R RangeType
range type
Definition: localbasis.hh:55
Dune::LocalFiniteElementTraits
traits helper struct
Definition: localfiniteelementtraits.hh:11
Dune::LocalFiniteElementTraits::LocalBasisType
LB LocalBasisType
Definition: localfiniteelementtraits.hh:14
Dune::LocalFiniteElementTraits::LocalCoefficientsType
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:18
Dune::LocalFiniteElementTraits::LocalInterpolationType
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:22
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.80.0 (May 16, 22:29, 2024)