DUNE PDELab (2.7)

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
21namespace 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
825namespace 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:280
Dune::LagrangeSimplexLocalFiniteElement
Lagrange finite element for simplices with arbitrary compile-time dimension and polynomial order.
Definition: lagrangesimplex.hh:836
Dune::LagrangeSimplexLocalFiniteElement::localInterpolation
const Traits::LocalInterpolationType & localInterpolation() const
Returns object that evaluates degrees of freedom.
Definition: lagrangesimplex.hh:872
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::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::localCoefficients
const Traits::LocalCoefficientsType & localCoefficients() const
Returns the assignment of the degrees of freedom to the element subentities.
Definition: lagrangesimplex.hh:865
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:767
Dune::PDELab::interpolate
void interpolate(const F &f, const GFS &gfs, XG &xg)
interpolation from a given grid function
Definition: interpolate.hh:174
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
static constexpr 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.111.3 (Jul 15, 22:36, 2024)