DUNE-FEM (unstable)

lagrangecube.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_LAGRANGECUBE_HH
6#define DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGECUBE_HH
7
8#include <array>
9#include <numeric>
10
13#include <dune/common/math.hh>
14
15#include <dune/geometry/referenceelements.hh>
16
17#include <dune/localfunctions/common/localbasis.hh>
18#include <dune/localfunctions/common/localfiniteelementtraits.hh>
19#include <dune/localfunctions/common/localkey.hh>
20
21namespace Dune { namespace Impl
22{
23 // Forward declaration
24 template<class LocalBasis>
25 class LagrangeCubeLocalInterpolation;
26
37 template<class D, class R, unsigned int dim, unsigned int k>
38 class LagrangeCubeLocalBasis
39 {
40 friend class LagrangeCubeLocalInterpolation<LagrangeCubeLocalBasis<D,R,dim,k> >;
41
42 // i-th Lagrange polynomial of degree k in one dimension
43 static R p(unsigned int i, D x)
44 {
45 R result(1.0);
46 for (unsigned int j=0; j<=k; j++)
47 if (j!=i) result *= (k*x-j)/((int)i-(int)j);
48 return result;
49 }
50
51 // derivative of ith Lagrange polynomial of degree k in one dimension
52 static R dp(unsigned int i, D x)
53 {
54 R result(0.0);
55
56 for (unsigned int j=0; j<=k; j++)
57 {
58 if (j!=i)
59 {
60 R prod( (k*1.0)/((int)i-(int)j) );
61 for (unsigned int l=0; l<=k; l++)
62 if (l!=i && l!=j)
63 prod *= (k*x-l)/((int)i-(int)l);
64 result += prod;
65 }
66 }
67 return result;
68 }
69
70 // Second derivative of j-th Lagrange polynomial of degree k in one dimension
71 // Formula and notation taken from https://en.wikipedia.org/wiki/Lagrange_polynomial#Derivatives
72 static R ddp(unsigned int j, D x)
73 {
74 R result(0.0);
75
76 for (unsigned int i=0; i<=k; i++)
77 {
78 if (i==j)
79 continue;
80
81 R sum(0);
82
83 for (unsigned int m=0; m<=k; m++)
84 {
85 if (m==i || m==j)
86 continue;
87
88 R prod( (k*1.0)/((int)j-(int)m) );
89 for (unsigned int l=0; l<=k; l++)
90 if (l!=i && l!=j && l!=m)
91 prod *= (k*x-l)/((int)j-(int)l);
92 sum += prod;
93 }
94
95 result += sum * ( (k*1.0)/((int)j-(int)i) );
96 }
97
98 return result;
99 }
100
101 // Return i as a d-digit number in the (k+1)-nary system
102 static std::array<unsigned int,dim> multiindex (unsigned int i)
103 {
104 std::array<unsigned int,dim> alpha;
105 for (unsigned int j=0; j<dim; j++)
106 {
107 alpha[j] = i % (k+1);
108 i = i/(k+1);
109 }
110 return alpha;
111 }
112
113 public:
114 using Traits = LocalBasisTraits<D,dim,FieldVector<D,dim>,R,1,FieldVector<R,1>,FieldMatrix<R,1,dim> >;
115
118 static constexpr unsigned int size ()
119 {
120 return power(k+1, dim);
121 }
122
124 void evaluateFunction(const typename Traits::DomainType& x,
125 std::vector<typename Traits::RangeType>& out) const
126 {
127 out.resize(size());
128
129 // Specialization for zero-order case
130 if (k==0)
131 {
132 out[0] = 1;
133 return;
134 }
135
136 if (k==1)
137 {
138 for (size_t i=0; i<size(); i++)
139 {
140 out[i] = 1;
141
142 for (unsigned int j=0; j<dim; j++)
143 // if j-th bit of i is set multiply with x[j], else with 1-x[j]
144 out[i] *= (i & (1<<j)) ? x[j] : 1-x[j];
145 }
146 return;
147 }
148
149 // General case
150 for (size_t i=0; i<size(); i++)
151 {
152 // convert index i to multiindex
153 std::array<unsigned int,dim> alpha(multiindex(i));
154
155 // initialize product
156 out[i] = 1.0;
157
158 // dimension by dimension
159 for (unsigned int j=0; j<dim; j++)
160 out[i] *= p(alpha[j],x[j]);
161 }
162 }
163
169 void evaluateJacobian(const typename Traits::DomainType& x,
170 std::vector<typename Traits::JacobianType>& out) const
171 {
172 out.resize(size());
173
174 // Specialization for k==0
175 if (k==0)
176 {
177 std::fill(out[0][0].begin(), out[0][0].end(), 0);
178 return;
179 }
180
181 // Specialization for k==1
182 if (k==1)
183 {
184 // Loop over all shape functions
185 for (size_t i=0; i<size(); i++)
186 {
187 // Loop over all coordinate directions
188 for (unsigned int j=0; j<dim; j++)
189 {
190 // Initialize: the overall expression is a product
191 // if j-th bit of i is set to 1, else -11
192 out[i][0][j] = (i & (1<<j)) ? 1 : -1;
193
194 for (unsigned int l=0; l<dim; l++)
195 {
196 if (j!=l)
197 // if l-th bit of i is set multiply with x[l], else with 1-x[l]
198 out[i][0][j] *= (i & (1<<l)) ? x[l] : 1-x[l];
199 }
200 }
201 }
202 return;
203 }
204
205 // The general case
206
207 // Loop over all shape functions
208 for (size_t i=0; i<size(); i++)
209 {
210 // convert index i to multiindex
211 std::array<unsigned int,dim> alpha(multiindex(i));
212
213 // Loop over all coordinate directions
214 for (unsigned int j=0; j<dim; j++)
215 {
216 // Initialize: the overall expression is a product
217 // if j-th bit of i is set to -1, else 1
218 out[i][0][j] = dp(alpha[j],x[j]);
219
220 // rest of the product
221 for (unsigned int l=0; l<dim; l++)
222 if (l!=j)
223 out[i][0][j] *= p(alpha[l],x[l]);
224 }
225 }
226 }
227
234 void partial(const std::array<unsigned int,dim>& order,
235 const typename Traits::DomainType& in,
236 std::vector<typename Traits::RangeType>& out) const
237 {
238 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
239
240 out.resize(size());
241
242 if (k==0)
243 {
244 out[0] = (totalOrder==0);
245 return;
246 }
247
248 if (k==1)
249 {
250 if (totalOrder == 0)
251 {
252 evaluateFunction(in, out);
253 }
254 else if (totalOrder == 1)
255 {
256 out.resize(size());
257
258 auto direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
259 if (direction >= dim)
260 DUNE_THROW(RangeError, "Direction of partial derivative not found!");
261
262 // Loop over all shape functions
263 for (std::size_t i = 0; i < size(); ++i)
264 {
265 // Initialize: the overall expression is a product
266 // if j-th bit of i is set to 1, otherwise to -1
267 out[i] = (i & (1<<direction)) ? 1 : -1;
268
269 for (unsigned int j = 0; j < dim; ++j)
270 {
271 if (direction != j)
272 // if j-th bit of i is set multiply with in[j], else with 1-in[j]
273 out[i] *= (i & (1<<j)) ? in[j] : 1-in[j];
274 }
275 }
276 }
277 else if (totalOrder == 2)
278 {
279
280 for (size_t i=0; i<size(); i++)
281 {
282 // convert index i to multiindex
283 std::array<unsigned int,dim> alpha(multiindex(i));
284
285 // Initialize: the overall expression is a product
286 out[i][0] = 1.0;
287
288 // rest of the product
289 for (std::size_t l=0; l<dim; l++)
290 {
291 switch (order[l])
292 {
293 case 0:
294 out[i][0] *= p(alpha[l],in[l]);
295 break;
296 case 1:
297 //std::cout << "dp: " << dp(alpha[l],in[l]) << std::endl;
298 out[i][0] *= dp(alpha[l],in[l]);
299 break;
300 case 2:
301 out[i][0] *= 0;
302 break;
303 default:
304 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
305 }
306 }
307 }
308 }
309 else
310 DUNE_THROW(NotImplemented, "Partial derivative of order " << totalOrder << " is not implemented!");
311
312 return;
313 }
314
315 // The case k>1
316
317 // Loop over all shape functions
318 for (size_t i=0; i<size(); i++)
319 {
320 // convert index i to multiindex
321 std::array<unsigned int,dim> alpha(multiindex(i));
322
323 // Initialize: the overall expression is a product
324 out[i][0] = 1.0;
325
326 // rest of the product
327 for (std::size_t l=0; l<dim; l++)
328 {
329 switch (order[l])
330 {
331 case 0:
332 out[i][0] *= p(alpha[l],in[l]);
333 break;
334 case 1:
335 out[i][0] *= dp(alpha[l],in[l]);
336 break;
337 case 2:
338 out[i][0] *= ddp(alpha[l],in[l]);
339 break;
340 default:
341 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
342 }
343 }
344 }
345 }
346
348 static constexpr unsigned int order ()
349 {
350 return k;
351 }
352 };
353
359 template<unsigned int dim, unsigned int k>
360 class LagrangeCubeLocalCoefficients
361 {
362 // Return i as a d-digit number in the (k+1)-nary system
363 static std::array<unsigned int,dim> multiindex (unsigned int i)
364 {
365 std::array<unsigned int,dim> alpha;
366 for (unsigned int j=0; j<dim; j++)
367 {
368 alpha[j] = i % (k+1);
369 i = i/(k+1);
370 }
371 return alpha;
372 }
373
375 void setup1d(std::vector<unsigned int>& subEntity)
376 {
377 assert(k>0);
378
379 unsigned lastIndex=0;
380
381 /* edge and vertex numbering
382 0----0----1
383 */
384
385 // edge (0)
386 subEntity[lastIndex++] = 0; // corner 0
387 for (unsigned i = 0; i < k - 1; ++i)
388 subEntity[lastIndex++] = 0; // inner dofs of element (0)
389
390 subEntity[lastIndex++] = 1; // corner 1
391
392 assert(power(k+1,dim)==lastIndex);
393 }
394
395 void setup2d(std::vector<unsigned int>& subEntity)
396 {
397 assert(k>0);
398
399 unsigned lastIndex=0;
400
401 // LocalKey: entity number, entity codim, dof indices within each entity
402 /* edge and vertex numbering
403 2----3----3
404 | |
405 | |
406 0 1
407 | |
408 | |
409 0----2----1
410 */
411
412 // lower edge (2)
413 subEntity[lastIndex++] = 0; // corner 0
414 for (unsigned i = 0; i < k - 1; ++i)
415 subEntity[lastIndex++] = 2; // inner dofs of lower edge (2)
416
417 subEntity[lastIndex++] = 1; // corner 1
418
419 // iterate from bottom to top over inner edge dofs
420 for (unsigned e = 0; e < k - 1; ++e) {
421 subEntity[lastIndex++] = 0; // left edge (0)
422 for (unsigned i = 0; i < k - 1; ++i)
423 subEntity[lastIndex++] = 0; // face dofs
424 subEntity[lastIndex++] = 1; // right edge (1)
425 }
426
427 // upper edge (3)
428 subEntity[lastIndex++] = 2; // corner 2
429 for (unsigned i = 0; i < k - 1; ++i)
430 subEntity[lastIndex++] = 3; // inner dofs of upper edge (3)
431
432 subEntity[lastIndex++] = 3; // corner 3
433
434 assert(power(k+1,dim)==lastIndex);
435 }
436
437 void setup3d(std::vector<unsigned int>& subEntity)
438 {
439 assert(k>0);
440
441 unsigned lastIndex=0;
442#ifndef NDEBUG
443 const unsigned numIndices = power(k+1,dim);
444 const unsigned numFaceIndices = power(k+1,dim-1);
445#endif
446 const unsigned numInnerEdgeDofs = k-1;
447
448 // LocalKey: entity number, entity codim, dof indices within each entity
449 /* edge and vertex numbering
450
451 6---(11)--7 6---------7
452 /| /| /| (5) /|
453 (8)| (9)| / | top / |
454 / (2) / (3) / |(3)bac/k |
455 4---(10)--5 | 4---------5 |
456 | | | | left|(0)| |(1)|right
457 | 2--(7)|---3 | 2-----|---3
458 (0) / (1) / |(2)front | /
459 |(4) |(5) | / (4) | /
460 |/ |/ |/ bottom |/
461 0---(6)---1 0---------1
462 */
463
464 // bottom face (4)
465 lastIndex=0;
466 // lower edge (6)
467 subEntity[lastIndex++] = 0; // corner 0
468 for (unsigned i = 0; i < numInnerEdgeDofs; ++i)
469 subEntity[lastIndex++] = 6; // inner dofs of lower edge (6)
470
471 subEntity[lastIndex++] = 1; // corner 1
472
473 // iterate from bottom to top over inner edge dofs
474 for (unsigned e = 0; e < numInnerEdgeDofs; ++e) {
475 subEntity[lastIndex++] = 4; // left edge (4)
476 for (unsigned i = 0; i < numInnerEdgeDofs; ++i)
477 subEntity[lastIndex++] = 4; // inner face dofs
478 subEntity[lastIndex++] = 5; // right edge (5)
479 }
480
481 // upper edge (7)
482 subEntity[lastIndex++] = 2; // corner 2
483 for (unsigned i = 0; i < k - 1; ++i)
484 subEntity[lastIndex++] = 7; // inner dofs of upper edge (7)
485 subEntity[lastIndex++] = 3; // corner 3
486
487 assert(numFaceIndices==lastIndex); // added 1 face so far
489
491 for(unsigned f = 0; f < numInnerEdgeDofs; ++f) {
492
493 // lower edge (connecting edges 0 and 1)
494 subEntity[lastIndex++] = 0; // dof on edge 0
495 for (unsigned i = 0; i < numInnerEdgeDofs; ++i)
496 subEntity[lastIndex++] = 2; // dof in front face
497 subEntity[lastIndex++] = 1; // dof on edge 1
498
499 // iterate from bottom to top over inner edge dofs
500 for (unsigned e = 0; e < numInnerEdgeDofs; ++e) {
501 subEntity[lastIndex++] = 0; // on left face (0)
502 for (unsigned i = 0; i < numInnerEdgeDofs; ++i)
503 subEntity[lastIndex++] = 0; // volume dofs
504 subEntity[lastIndex++] = 1; // right face (1)
505 }
506
507 // upper edge (connecting edges 0 and 1)
508 subEntity[lastIndex++] = 2; // dof on edge 2
509 for (unsigned i = 0; i < numInnerEdgeDofs; ++i)
510 subEntity[lastIndex++] = 3; // dof on rear face (3)
511 subEntity[lastIndex++] = 3; // dof on edge 3
512
513 assert(lastIndex==(f+1+1)*numFaceIndices);
514 }
515
517 // lower edge (10)
518 subEntity[lastIndex++] = 4; // corner 4
519 for (unsigned i = 0; i < k - 1; ++i)
520 subEntity[lastIndex++] = 10; // inner dofs on lower edge (10)
521 subEntity[lastIndex++] = 5; // corner 5
522
523 // iterate from bottom to top over inner edge dofs
524 for (unsigned e = 0; e < k - 1; ++e) {
525 subEntity[lastIndex++] = 8; // left edge (8)
526 for (unsigned i = 0; i < k - 1; ++i)
527 subEntity[lastIndex++] = 5; // face dofs
528 subEntity[lastIndex++] = 9; // right edge (9)
529 }
530
531 // upper edge (11)
532 subEntity[lastIndex++] = 6; // corner 6
533 for (unsigned i = 0; i < k - 1; ++i)
534 subEntity[lastIndex++] = 11; // inner dofs of upper edge (11)
535 subEntity[lastIndex++] = 7; // corner 7
536
537 assert(numIndices==lastIndex);
538 }
539
540 public:
542 LagrangeCubeLocalCoefficients ()
543 : localKeys_(size())
544 {
545 if (k==0)
546 {
547 localKeys_[0] = LocalKey(0,0,0);
548 return;
549 }
550
551 if (k==1)
552 {
553 for (std::size_t i=0; i<size(); i++)
554 localKeys_[i] = LocalKey(i,dim,0);
555 return;
556 }
557
558 // Now: the general case
559
560 // Set up array of codimension-per-dof-number
561 std::vector<unsigned int> codim(size());
562
563 for (std::size_t i=0; i<codim.size(); i++)
564 {
565 codim[i] = 0;
566
567 // Codimension gets increased by 1 for each coordinate direction
568 // where dof is on boundary
569 std::array<unsigned int,dim> mIdx = multiindex(i);
570 for (unsigned int j=0; j<dim; j++)
571 if (mIdx[j]==0 or mIdx[j]==k)
572 codim[i]++;
573 }
574
575 // Set up index vector (the index of the dof in the set of dofs of a given subentity)
576 // Algorithm: the 'index' has the same ordering as the dof number 'i'.
577 // To make it consecutive we interpret 'i' in the (k+1)-adic system, omit all digits
578 // that correspond to axes where the dof is on the element boundary, and transform the
579 // rest to the (k-1)-adic system.
580 std::vector<unsigned int> index(size());
581
582 for (std::size_t i=0; i<size(); i++)
583 {
584 index[i] = 0;
585
586 std::array<unsigned int,dim> mIdx = multiindex(i);
587
588 for (int j=dim-1; j>=0; j--)
589 if (mIdx[j]>0 && mIdx[j]<k)
590 index[i] = (k-1)*index[i] + (mIdx[j]-1);
591 }
592
593 // Set up entity and dof numbers for each (supported) dimension separately
594 std::vector<unsigned int> subEntity(size());
595
596 if (dim==1) {
597
598 setup1d(subEntity);
599
600 } else if (dim==2) {
601
602 setup2d(subEntity);
603
604 } else if (dim==3) {
605
606 setup3d(subEntity);
607
608 } else
609 DUNE_THROW(Dune::NotImplemented, "LagrangeCubeLocalCoefficients for order " << k << " and dim == " << dim);
610
611 for (size_t i=0; i<size(); i++)
612 localKeys_[i] = LocalKey(subEntity[i], codim[i], index[i]);
613 }
614
616 static constexpr std::size_t size ()
617 {
618 return power(k+1,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
635 template<class LocalBasis>
636 class LagrangeCubeLocalInterpolation
637 {
638 public:
639
647 template<typename F, typename C>
648 void interpolate (const F& f, std::vector<C>& out) const
649 {
650 constexpr auto dim = LocalBasis::Traits::dimDomain;
651 constexpr auto k = LocalBasis::order();
652 using D = typename LocalBasis::Traits::DomainFieldType;
653
654 typename LocalBasis::Traits::DomainType x;
655
656 out.resize(LocalBasis::size());
657
658 // Specialization for zero-order case
659 if (k==0)
660 {
661 auto center = ReferenceElements<D,dim>::cube().position(0,0);
662 out[0] = f(center);
663 return;
664 }
665
666 // Specialization for first-order case
667 if (k==1)
668 {
669 for (unsigned int i=0; i<LocalBasis::size(); i++)
670 {
671 // Generate coordinate of the i-th corner of the reference cube
672 for (int j=0; j<dim; j++)
673 x[j] = (i & (1<<j)) ? 1.0 : 0.0;
674
675 out[i] = f(x);
676 }
677 return;
678 }
679
680 // The general case
681 for (unsigned int i=0; i<LocalBasis::size(); i++)
682 {
683 // convert index i to multiindex
684 std::array<unsigned int,dim> alpha(LocalBasis::multiindex(i));
685
686 // Generate coordinate of the i-th Lagrange point
687 for (unsigned int j=0; j<dim; j++)
688 x[j] = (1.0*alpha[j])/k;
689
690 out[i] = f(x);
691 }
692 }
693
694 };
695
696} } // namespace Dune::Impl
697
698namespace Dune
699{
707 template<class D, class R, int dim, int k>
709 {
710 public:
714 Impl::LagrangeCubeLocalCoefficients<dim,k>,
715 Impl::LagrangeCubeLocalInterpolation<Impl::LagrangeCubeLocalBasis<D,R,dim,k> > >;
716
719 const typename Traits::LocalBasisType& localBasis () const
720 {
721 return basis_;
722 }
723
727 {
728 return coefficients_;
729 }
730
734 {
735 return interpolation_;
736 }
737
739 static constexpr std::size_t size ()
740 {
741 return power(k+1,dim);
742 }
743
746 static constexpr GeometryType type ()
747 {
748 return GeometryTypes::cube(dim);
749 }
750
751 private:
752 Impl::LagrangeCubeLocalBasis<D,R,dim,k> basis_;
753 Impl::LagrangeCubeLocalCoefficients<dim,k> coefficients_;
754 Impl::LagrangeCubeLocalInterpolation<Impl::LagrangeCubeLocalBasis<D,R,dim,k> > interpolation_;
755 };
756
757} // namespace Dune
758
759#endif // DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGECUBE_HH
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
Lagrange finite element for cubes with arbitrary compile-time dimension and polynomial order.
Definition: lagrangecube.hh:709
const Traits::LocalBasisType & localBasis() const
Returns the local basis, i.e., the set of shape functions.
Definition: lagrangecube.hh:719
const Traits::LocalInterpolationType & localInterpolation() const
Returns object that evaluates degrees of freedom.
Definition: lagrangecube.hh:733
static constexpr GeometryType type()
The reference element that the local finite element is defined on.
Definition: lagrangecube.hh:746
static constexpr std::size_t size()
The number of shape functions.
Definition: lagrangecube.hh:739
const Traits::LocalCoefficientsType & localCoefficients() const
Returns the assignment of the degrees of freedom to the element subentities.
Definition: lagrangecube.hh:726
Default exception for dummy implementations.
Definition: exceptions.hh:263
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
static void interpolate(const GridFunction &u, DiscreteFunction &v)
perform native interpolation of a discrete function space
Definition: interpolate.hh:54
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:462
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
Some useful basic math stuff.
Dune namespace.
Definition: alignedallocator.hh:13
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
constexpr Base power(Base m, Exponent p)
Power method for integer exponents.
Definition: math.hh:75
static const ReferenceElement & cube()
get hypercube reference elements
Definition: referenceelements.hh:168
D DomainType
domain type
Definition: localbasis.hh:43
traits helper struct
Definition: localfiniteelementtraits.hh:13
LB LocalBasisType
Definition: localfiniteelementtraits.hh:16
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:20
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:24
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)