DUNE-FEM (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
22namespace 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
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:114
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_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::Fem::interpolate
static void interpolate(const GridFunction &u, DiscreteFunction &v)
perform native interpolation of a discrete function space
Definition: interpolate.hh:54
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
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: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.111.3 (Nov 12, 23:30, 2024)