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
23namespace 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
827namespace 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::localInterpolation
const Traits::LocalInterpolationType & localInterpolation() const
Returns object that evaluates degrees of freedom.
Definition: lagrangesimplex.hh:874
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::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::localCoefficients
const Traits::LocalCoefficientsType & localCoefficients() const
Returns the assignment of the degrees of freedom to the element subentities.
Definition: lagrangesimplex.hh:867
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
static constexpr 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.111.3 (Jul 15, 22:36, 2024)