Dune Core Modules (unstable)

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 // SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4 // SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5 #ifndef DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGESIMPLEX_HH
6 #define DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGESIMPLEX_HH
7 
8 #include <array>
9 #include <numeric>
10 
11 #include <dune/common/deprecated.hh>
12 #include <dune/common/fmatrix.hh>
13 #include <dune/common/fvector.hh>
14 #include <dune/common/math.hh>
15 
16 #include <dune/geometry/referenceelements.hh>
17 
18 #include <dune/localfunctions/common/localbasis.hh>
19 #include <dune/localfunctions/common/localfiniteelementtraits.hh>
20 #include <dune/localfunctions/common/localkey.hh>
21 
22 namespace Dune { namespace Impl
23 {
34  template<class D, class R, unsigned int dim, unsigned int k>
35  class LagrangeSimplexLocalBasis
36  {
37  public:
38  using Traits = LocalBasisTraits<D,dim,FieldVector<D,dim>,R,1,FieldVector<R,1>,FieldMatrix<R,1,dim> >;
39 
44  static constexpr unsigned int size ()
45  {
46  return binomial(k+dim,dim);
47  }
48 
50  void evaluateFunction(const typename Traits::DomainType& x,
51  std::vector<typename Traits::RangeType>& out) const
52  {
53  out.resize(size());
54 
55  // Specialization for zero-order case
56  if (k==0)
57  {
58  out[0] = 1;
59  return;
60  }
61 
62  // Specialization for first-order case
63  if (k==1)
64  {
65  out[0] = 1.0;
66  for (size_t i=0; i<dim; i++)
67  {
68  out[0] -= x[i];
69  out[i+1] = x[i];
70  }
71  return;
72  }
73 
74  assert(k>=2);
75 
76  auto lagrangeNode = [](unsigned int i) { return ((D)i)/k; };
77 
78  if (dim==1)
79  {
80  for (unsigned int i=0; i<size(); i++)
81  {
82  out[i] = 1.0;
83  for (unsigned int alpha=0; alpha<i; alpha++)
84  out[i] *= (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
85  for (unsigned int gamma=i+1; gamma<=k; gamma++)
86  out[i] *= (x[0]-lagrangeNode(gamma))/(lagrangeNode(i)-lagrangeNode(gamma));
87  }
88  return;
89  }
90 
91  if (dim==2)
92  {
93  int n=0;
94  for (unsigned int j=0; j<=k; j++)
95  for (unsigned int i=0; i<=k-j; i++)
96  {
97  out[n] = 1.0;
98  for (unsigned int alpha=0; alpha<i; alpha++)
99  out[n] *= (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
100  for (unsigned int beta=0; beta<j; beta++)
101  out[n] *= (x[1]-lagrangeNode(beta))/(lagrangeNode(j)-lagrangeNode(beta));
102  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
103  out[n] *= (lagrangeNode(gamma)-x[0]-x[1])/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
104  n++;
105  }
106 
107  return;
108  }
109 
110  if (dim!=3)
111  DUNE_THROW(NotImplemented, "LagrangeSimplexLocalBasis for k>=2 only implemented for dim==1 or dim==3");
112 
113  typename Traits::DomainType kx = x;
114  kx *= k;
115  unsigned int n = 0;
116  unsigned int i[4];
117  R factor[4];
118  for (i[2] = 0; i[2] <= k; ++i[2])
119  {
120  factor[2] = 1.0;
121  for (unsigned int j = 0; j < i[2]; ++j)
122  factor[2] *= (kx[2]-j) / (i[2]-j);
123  for (i[1] = 0; i[1] <= k - i[2]; ++i[1])
124  {
125  factor[1] = 1.0;
126  for (unsigned int j = 0; j < i[1]; ++j)
127  factor[1] *= (kx[1]-j) / (i[1]-j);
128  for (i[0] = 0; i[0] <= k - i[1] - i[2]; ++i[0])
129  {
130  factor[0] = 1.0;
131  for (unsigned int j = 0; j < i[0]; ++j)
132  factor[0] *= (kx[0]-j) / (i[0]-j);
133  i[3] = k - i[0] - i[1] - i[2];
134  D kx3 = k - kx[0] - kx[1] - kx[2];
135  factor[3] = 1.0;
136  for (unsigned int j = 0; j < i[3]; ++j)
137  factor[3] *= (kx3-j) / (i[3]-j);
138  out[n++] = factor[0] * factor[1] * factor[2] * factor[3];
139  }
140  }
141  }
142  }
143 
149  void evaluateJacobian(const typename Traits::DomainType& x,
150  std::vector<typename Traits::JacobianType>& out) const
151  {
152  out.resize(size());
153 
154  // Specialization for k==0
155  if (k==0)
156  {
157  std::fill(out[0][0].begin(), out[0][0].end(), 0);
158  return;
159  }
160 
161  // Specialization for k==1
162  if (k==1)
163  {
164  std::fill(out[0][0].begin(), out[0][0].end(), -1);
165 
166  for (unsigned int i=0; i<dim; i++)
167  for (unsigned int j=0; j<dim; j++)
168  out[i+1][0][j] = (i==j);
169 
170  return;
171  }
172 
173  auto lagrangeNode = [](unsigned int i) { return ((D)i)/k; };
174 
175  // Specialization for dim==1
176  if (dim==1)
177  {
178  for (unsigned int i=0; i<=k; i++)
179  {
180  // x_0 derivative
181  out[i][0][0] = 0.0;
182  R factor=1.0;
183  for (unsigned int a=0; a<i; a++)
184  {
185  R product=factor;
186  for (unsigned int alpha=0; alpha<i; alpha++)
187  product *= (alpha==a) ? 1.0/(lagrangeNode(i)-lagrangeNode(alpha))
188  : (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
189  for (unsigned int gamma=i+1; gamma<=k; gamma++)
190  product *= (lagrangeNode(gamma)-x[0])/(lagrangeNode(gamma)-lagrangeNode(i));
191  out[i][0][0] += product;
192  }
193  for (unsigned int c=i+1; c<=k; c++)
194  {
195  R product=factor;
196  for (unsigned int alpha=0; alpha<i; alpha++)
197  product *= (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
198  for (unsigned int gamma=i+1; gamma<=k; gamma++)
199  product *= (gamma==c) ? -1.0/(lagrangeNode(gamma)-lagrangeNode(i))
200  : (lagrangeNode(gamma)-x[0])/(lagrangeNode(gamma)-lagrangeNode(i));
201  out[i][0][0] += product;
202  }
203  }
204  return;
205  }
206 
207  if (dim==2)
208  {
209  int n=0;
210  for (unsigned int j=0; j<=k; j++)
211  for (unsigned int i=0; i<=k-j; i++)
212  {
213  // x_0 derivative
214  out[n][0][0] = 0.0;
215  R factor=1.0;
216  for (unsigned int beta=0; beta<j; beta++)
217  factor *= (x[1]-lagrangeNode(beta))/(lagrangeNode(j)-lagrangeNode(beta));
218  for (unsigned int a=0; a<i; a++)
219  {
220  R product=factor;
221  for (unsigned int alpha=0; alpha<i; alpha++)
222  if (alpha==a)
223  product *= D(1)/(lagrangeNode(i)-lagrangeNode(alpha));
224  else
225  product *= (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
226  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
227  product *= (lagrangeNode(gamma)-x[0]-x[1])/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
228  out[n][0][0] += product;
229  }
230  for (unsigned int c=i+j+1; c<=k; c++)
231  {
232  R product=factor;
233  for (unsigned int alpha=0; alpha<i; alpha++)
234  product *= (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
235  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
236  if (gamma==c)
237  product *= -D(1)/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
238  else
239  product *= (lagrangeNode(gamma)-x[0]-x[1])/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
240  out[n][0][0] += product;
241  }
242 
243  // x_1 derivative
244  out[n][0][1] = 0.0;
245  factor = 1.0;
246  for (unsigned int alpha=0; alpha<i; alpha++)
247  factor *= (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
248  for (unsigned int b=0; b<j; b++)
249  {
250  R product=factor;
251  for (unsigned int beta=0; beta<j; beta++)
252  if (beta==b)
253  product *= D(1)/(lagrangeNode(j)-lagrangeNode(beta));
254  else
255  product *= (x[1]-lagrangeNode(beta))/(lagrangeNode(j)-lagrangeNode(beta));
256  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
257  product *= (lagrangeNode(gamma)-x[0]-x[1])/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
258  out[n][0][1] += product;
259  }
260  for (unsigned int c=i+j+1; c<=k; c++)
261  {
262  R product=factor;
263  for (unsigned int beta=0; beta<j; beta++)
264  product *= (x[1]-lagrangeNode(beta))/(lagrangeNode(j)-lagrangeNode(beta));
265  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
266  if (gamma==c)
267  product *= -D(1)/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
268  else
269  product *= (lagrangeNode(gamma)-x[0]-x[1])/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
270  out[n][0][1] += product;
271  }
272 
273  n++;
274  }
275 
276  return;
277  }
278 
279  if (dim!=3)
280  DUNE_THROW(NotImplemented, "LagrangeSimplexLocalBasis only implemented for dim==3!");
281 
282  // Specialization for arbitrary order and dim==3
283  typename Traits::DomainType kx = x;
284  kx *= k;
285  unsigned int n = 0;
286  unsigned int i[4];
287  R factor[4];
288  for (i[2] = 0; i[2] <= k; ++i[2])
289  {
290  factor[2] = 1.0;
291  for (unsigned int j = 0; j < i[2]; ++j)
292  factor[2] *= (kx[2]-j) / (i[2]-j);
293  for (i[1] = 0; i[1] <= k - i[2]; ++i[1])
294  {
295  factor[1] = 1.0;
296  for (unsigned int j = 0; j < i[1]; ++j)
297  factor[1] *= (kx[1]-j) / (i[1]-j);
298  for (i[0] = 0; i[0] <= k - i[1] - i[2]; ++i[0])
299  {
300  factor[0] = 1.0;
301  for (unsigned int j = 0; j < i[0]; ++j)
302  factor[0] *= (kx[0]-j) / (i[0]-j);
303  i[3] = k - i[0] - i[1] - i[2];
304  D kx3 = k - kx[0] - kx[1] - kx[2];
305  R sum3 = 0.0;
306  factor[3] = 1.0;
307  for (unsigned int j = 0; j < i[3]; ++j)
308  factor[3] /= i[3] - j;
309  R prod_all = factor[0] * factor[1] * factor[2] * factor[3];
310  for (unsigned int j = 0; j < i[3]; ++j)
311  {
312  R prod = prod_all;
313  for (unsigned int l = 0; l < i[3]; ++l)
314  if (j == l)
315  prod *= -R(k);
316  else
317  prod *= kx3 - l;
318  sum3 += prod;
319  }
320  for (unsigned int j = 0; j < i[3]; ++j)
321  factor[3] *= kx3 - j;
322  for (unsigned int m = 0; m < 3; ++m)
323  {
324  out[n][0][m] = sum3;
325  for (unsigned int j = 0; j < i[m]; ++j)
326  {
327  R prod = factor[3];
328  for (unsigned int p = 0; p < 3; ++p)
329  {
330  if (m == p)
331  for (unsigned int l = 0; l < i[p]; ++l)
332  prod *= (j==l) ? R(k) / (i[p]-l) : R(kx[p]-l) / (i[p]-l);
333  else
334  prod *= factor[p];
335  }
336  out[n][0][m] += prod;
337  }
338  }
339  n++;
340  }
341  }
342  }
343  }
344 
351  void partial(const std::array<unsigned int,dim>& order,
352  const typename Traits::DomainType& in,
353  std::vector<typename Traits::RangeType>& out) const
354  {
355  auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
356 
357  out.resize(size());
358 
359  if (totalOrder == 0) {
360  evaluateFunction(in, out);
361  return;
362  }
363 
364  if (k==0)
365  {
366  out[0] = 0;
367  return;
368  }
369 
370  if (k==1)
371  {
372  if (totalOrder==1)
373  {
374  auto direction = std::find(order.begin(), order.end(), 1);
375 
376  out[0] = -1;
377  for (unsigned int i=0; i<dim; i++)
378  out[i+1] = (i==(direction-order.begin()));
379  }
380  else // all higher order derivatives are zero
381  std::fill(out.begin(), out.end(), 0);
382  return;
383  }
384 
385  if (dim==2)
386  {
387  auto lagrangeNode = [](unsigned int i) { return ((D)i)/k; };
388 
389  // Helper method: Return a single Lagrangian factor of l_ij evaluated at x
390  auto lagrangianFactor = [&lagrangeNode]
391  (const int no, const int i, const int j, const typename Traits::DomainType& x)
392  -> typename Traits::RangeType
393  {
394  if ( no < i)
395  return (x[0]-lagrangeNode(no))/(lagrangeNode(i)-lagrangeNode(no));
396  if (no < i+j)
397  return (x[1]-lagrangeNode(no-i))/(lagrangeNode(j)-lagrangeNode(no-i));
398  return (lagrangeNode(no+1)-x[0]-x[1])/(lagrangeNode(no+1)-lagrangeNode(i)-lagrangeNode(j));
399  };
400 
401  // Helper method: Return the derivative of a single Lagrangian factor of l_ij evaluated at x
402  // direction: Derive in x-direction if this is 0, otherwise derive in y direction
403  auto lagrangianFactorDerivative = [&lagrangeNode]
404  (const int direction, const int no, const int i, const int j, const typename Traits::DomainType&)
405  -> typename Traits::RangeType
406  {
407  using T = typename Traits::RangeType;
408  if ( no < i)
409  return (direction == 0) ? T(1.0/(lagrangeNode(i)-lagrangeNode(no))) : T(0);
410 
411  if (no < i+j)
412  return (direction == 0) ? T(0) : T(1.0/(lagrangeNode(j)-lagrangeNode(no-i)));
413 
414  return -1.0/(lagrangeNode(no+1)-lagrangeNode(i)-lagrangeNode(j));
415  };
416 
417  if (totalOrder==1)
418  {
419  int direction = std::find(order.begin(), order.end(), 1)-order.begin();
420 
421  int n=0;
422  for (unsigned int j=0; j<=k; j++)
423  {
424  for (unsigned int i=0; i<=k-j; i++, n++)
425  {
426  out[n] = 0.0;
427  for (unsigned int no1=0; no1 < k; no1++)
428  {
429  R factor = lagrangianFactorDerivative(direction, no1, i, j, in);
430  for (unsigned int no2=0; no2 < k; no2++)
431  if (no1 != no2)
432  factor *= lagrangianFactor(no2, i, j, in);
433 
434  out[n] += factor;
435  }
436  }
437  }
438  return;
439  }
440 
441  if (totalOrder==2)
442  {
443  std::array<int,2> directions;
444  unsigned int counter = 0;
445  auto nonconstOrder = order; // need a copy that I can modify
446  for (int i=0; i<2; i++)
447  {
448  while (nonconstOrder[i])
449  {
450  directions[counter++] = i;
451  nonconstOrder[i]--;
452  }
453  }
454 
455  //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
456  int n=0;
457  for (unsigned int j=0; j<=k; j++)
458  {
459  for (unsigned int i=0; i<=k-j; i++, n++)
460  {
461  R res = 0.0;
462 
463  for (unsigned int no1=0; no1 < k; no1++)
464  {
465  R factor1 = lagrangianFactorDerivative(directions[0], no1, i, j, in);
466  for (unsigned int no2=0; no2 < k; no2++)
467  {
468  if (no1 == no2)
469  continue;
470  R factor2 = factor1*lagrangianFactorDerivative(directions[1], no2, i, j, in);
471  for (unsigned int no3=0; no3 < k; no3++)
472  {
473  if (no3 == no1 || no3 == no2)
474  continue;
475  factor2 *= lagrangianFactor(no3, i, j, in);
476  }
477  res += factor2;
478  }
479  }
480  out[n] = res;
481  }
482  }
483 
484  return;
485  } // totalOrder==2
486 
487  } // dim==2
488 
489  DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
490  }
491 
493  static constexpr unsigned int order ()
494  {
495  return k;
496  }
497  };
498 
504  template<unsigned int dim, unsigned int k>
505  class LagrangeSimplexLocalCoefficients
506  {
507  public:
509  LagrangeSimplexLocalCoefficients ()
510  : localKeys_(size())
511  {
512  if (k==0)
513  {
514  localKeys_[0] = LocalKey(0,0,0);
515  return;
516  }
517 
518  if (k==1)
519  {
520  for (std::size_t i=0; i<size(); i++)
521  localKeys_[i] = LocalKey(i,dim,0);
522  return;
523  }
524 
525  if (dim==1)
526  {
527  // Order is at least 2 here
528  localKeys_[0] = LocalKey(0,1,0); // vertex dof
529  for (unsigned int i=1; i<k; i++)
530  localKeys_[i] = LocalKey(0,0,i-1); // element dofs
531  localKeys_.back() = LocalKey(1,1,0); // vertex dof
532  return;
533  }
534 
535  if (dim==2)
536  {
537  int n=0;
538  int c=0;
539  for (unsigned int j=0; j<=k; j++)
540  for (unsigned int i=0; i<=k-j; i++)
541  {
542  if (i==0 && j==0)
543  {
544  localKeys_[n++] = LocalKey(0,2,0);
545  continue;
546  }
547  if (i==k && j==0)
548  {
549  localKeys_[n++] = LocalKey(1,2,0);
550  continue;
551  }
552  if (i==0 && j==k)
553  {
554  localKeys_[n++] = LocalKey(2,2,0);
555  continue;
556  }
557  if (j==0)
558  {
559  localKeys_[n++] = LocalKey(0,1,i-1);
560  continue;
561  }
562  if (i==0)
563  {
564  localKeys_[n++] = LocalKey(1,1,j-1);
565  continue;
566  }
567  if (i+j==k)
568  {
569  localKeys_[n++] = LocalKey(2,1,j-1);
570  continue;
571  }
572  localKeys_[n++] = LocalKey(0,0,c++);
573  }
574  return;
575  }
576 
577  if (dim==3)
578  {
579  std::array<unsigned int, dim+1> vertexMap;
580  for (unsigned int i=0; i<=dim; i++)
581  vertexMap[i] = i;
582  generateLocalKeys(vertexMap);
583  return;
584  }
585  DUNE_THROW(NotImplemented, "LagrangeSimplexLocalCoefficients only implemented for k<=1 or dim<=3!");
586  }
587 
594  LagrangeSimplexLocalCoefficients (const std::array<unsigned int, dim+1> vertexMap)
595  : localKeys_(size())
596  {
597  if (dim!=2 && dim!=3)
598  DUNE_THROW(NotImplemented, "LagrangeSimplexLocalCoefficients only implemented for dim==2 and dim==3!");
599 
600  generateLocalKeys(vertexMap);
601  }
602 
603 
604  template<class VertexMap>
605  LagrangeSimplexLocalCoefficients(const VertexMap &vertexmap)
606  : localKeys_(size())
607  {
608  if (dim!=2 && dim!=3)
609  DUNE_THROW(NotImplemented, "LagrangeSimplexLocalCoefficients only implemented for dim==2 and dim==3!");
610 
611  std::array<unsigned int, dim+1> vertexmap_array;
612  std::copy(vertexmap, vertexmap + dim + 1, vertexmap_array.begin());
613  generateLocalKeys(vertexmap_array);
614  }
615 
617  static constexpr std::size_t size ()
618  {
619  return binomial(k+dim,dim);
620  }
621 
623  const LocalKey& localKey (std::size_t i) const
624  {
625  return localKeys_[i];
626  }
627 
628  private:
629  std::vector<LocalKey> localKeys_;
630 
631  void generateLocalKeys(const std::array<unsigned int, dim+1> vertexMap)
632  {
633  if (k==0)
634  {
635  localKeys_[0] = LocalKey(0,0,0);
636  return;
637  }
638 
639  if (dim==2)
640  {
641  // Create default assignment
642  int n=0;
643  int c=0;
644  for (unsigned int j=0; j<=k; j++)
645  for (unsigned int i=0; i<=k-j; i++)
646  {
647  if (i==0 && j==0)
648  {
649  localKeys_[n++] = LocalKey(0,2,0);
650  continue;
651  }
652  if (i==k && j==0)
653  {
654  localKeys_[n++] = LocalKey(1,2,0);
655  continue;
656  }
657  if (i==0 && j==k)
658  {
659  localKeys_[n++] = LocalKey(2,2,0);
660  continue;
661  }
662  if (j==0)
663  {
664  localKeys_[n++] = LocalKey(0,1,i-1);
665  continue;
666  }
667  if (i==0)
668  {
669  localKeys_[n++] = LocalKey(1,1,j-1);
670  continue;
671  }
672  if (i+j==k)
673  {
674  localKeys_[n++] = LocalKey(2,1,j-1);
675  continue;
676  }
677  localKeys_[n++] = LocalKey(0,0,c++);
678  }
679 
680  // Flip edge orientations, if requested
681  bool flip[3];
682  flip[0] = vertexMap[0] > vertexMap[1];
683  flip[1] = vertexMap[0] > vertexMap[2];
684  flip[2] = vertexMap[1] > vertexMap[2];
685  for (std::size_t i=0; i<size(); i++)
686  if (localKeys_[i].codim()==1 && flip[localKeys_[i].subEntity()])
687  localKeys_[i].index(k-2-localKeys_[i].index());
688 
689  return;
690  }
691 
692  if (dim!=3)
693  DUNE_THROW(NotImplemented, "LagrangeSimplexLocalCoefficients only implemented for dim==3!");
694 
695  unsigned int subindex[16];
696  unsigned int codim_count[4] = {0};
697  for (unsigned int m = 1; m < 16; ++m)
698  {
699  unsigned int codim = !(m&1) + !(m&2) + !(m&4) + !(m&8);
700  subindex[m] = codim_count[codim]++;
701  }
702 
703  int a1 = (3*k + 12)*k + 11;
704  int a2 = -3*k - 6;
705  unsigned int dof_count[16] = {0};
706  unsigned int i[4];
707  for (i[3] = 0; i[3] <= k; ++i[3])
708  for (i[2] = 0; i[2] <= k - i[3]; ++i[2])
709  for (i[1] = 0; i[1] <= k - i[2] - i[3]; ++i[1])
710  {
711  i[0] = k - i[1] - i[2] - i[3];
712  unsigned int j[4];
713  unsigned int entity = 0;
714  unsigned int codim = 0;
715  for (unsigned int m = 0; m < 4; ++m)
716  {
717  j[m] = i[vertexMap[m]];
718  entity += !!j[m] << m;
719  codim += !j[m];
720  }
721  int local_index = j[3]*(a1 + (a2 + j[3])*j[3])/6
722  + j[2]*(2*(k - j[3]) + 3 - j[2])/2 + j[1];
723  localKeys_[local_index] = LocalKey(subindex[entity], codim, dof_count[entity]++);
724  }
725  }
726  };
727 
732  template<class LocalBasis>
733  class LagrangeSimplexLocalInterpolation
734  {
735  static const int kdiv = (LocalBasis::order() == 0 ? 1 : LocalBasis::order());
736  public:
737 
745  template<typename F, typename C>
746  void interpolate (const F& f, std::vector<C>& out) const
747  {
748  constexpr auto dim = LocalBasis::Traits::dimDomain;
749  constexpr auto k = LocalBasis::order();
750  using D = typename LocalBasis::Traits::DomainFieldType;
751 
752  typename LocalBasis::Traits::DomainType x;
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:114
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_NO_DEPRECATED_* macros.
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:218
Dune::GeometryTypes::simplex
constexpr GeometryType simplex(unsigned int dim)
Returns a GeometryType representing a simplex of dimension dim.
Definition: type.hh:453
Dune::Hybrid::accumulate
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
math.hh
Some useful basic math stuff.
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Dune::size
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
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:131
Dune::Geo::ReferenceElements::simplex
static const ReferenceElement & simplex()
get simplex reference elements
Definition: referenceelements.hh:162
Dune::LocalBasisTraits::DomainType
D DomainType
domain type
Definition: localbasis.hh:43
Dune::LocalBasisTraits::RangeType
R RangeType
range type
Definition: localbasis.hh:52
Dune::LocalFiniteElementTraits
traits helper struct
Definition: localfiniteelementtraits.hh:13
Dune::LocalFiniteElementTraits::LocalBasisType
LB LocalBasisType
Definition: localfiniteelementtraits.hh:16
Dune::LocalFiniteElementTraits::LocalCoefficientsType
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:20
Dune::LocalFiniteElementTraits::LocalInterpolationType
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:24
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.80.0 (Apr 27, 22:29, 2024)