Dune Core Modules (2.9.0)

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