3#ifndef DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGESIMPLEX_HH
4#define DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGESIMPLEX_HH
14#include <dune/geometry/referenceelements.hh>
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>
21namespace Dune {
namespace Impl
33 template<
class D,
class R,
unsigned int dim,
unsigned int k>
34 class LagrangeSimplexLocalBasis
37 using Traits = LocalBasisTraits<D,dim,FieldVector<D,dim>,R,1,FieldVector<R,1>,FieldMatrix<R,1,dim> >;
43 static constexpr unsigned int size ()
50 std::vector<typename Traits::RangeType>& out)
const
65 for (
size_t i=0; i<dim; i++)
75 auto lagrangeNode = [](
unsigned int i) {
return ((D)i)/k; };
79 for (
unsigned int i=0; i<size(); i++)
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));
93 for (
unsigned int j=0; j<=k; j++)
94 for (
unsigned int i=0; i<=k-j; i++)
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));
110 DUNE_THROW(NotImplemented,
"LagrangeSimplexLocalBasis for k>=2 only implemented for dim==1 or dim==3");
117 for (i[2] = 0; i[2] <= k; ++i[2])
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])
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])
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];
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];
149 std::vector<typename Traits::JacobianType>& out)
const
156 std::fill(out[0][0].begin(), out[0][0].end(), 0);
163 std::fill(out[0][0].begin(), out[0][0].end(), -1);
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);
172 auto lagrangeNode = [](
unsigned int i) {
return ((D)i)/k; };
177 for (
unsigned int i=0; i<=k; i++)
182 for (
unsigned int a=0; a<i; a++)
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;
192 for (
unsigned int c=i+1; c<=k; c++)
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;
209 for (
unsigned int j=0; j<=k; j++)
210 for (
unsigned int i=0; i<=k-j; i++)
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++)
220 for (
unsigned int alpha=0; alpha<i; alpha++)
222 product *= D(1)/(lagrangeNode(i)-lagrangeNode(alpha));
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;
229 for (
unsigned int c=i+j+1; c<=k; c++)
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++)
236 product *= -D(1)/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
238 product *= (lagrangeNode(gamma)-x[0]-x[1])/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
239 out[n][0][0] += product;
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++)
250 for (
unsigned int beta=0; beta<j; beta++)
252 product *= D(1)/(lagrangeNode(j)-lagrangeNode(beta));
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;
259 for (
unsigned int c=i+j+1; c<=k; c++)
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++)
266 product *= -D(1)/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
268 product *= (lagrangeNode(gamma)-x[0]-x[1])/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
269 out[n][0][1] += product;
279 DUNE_THROW(NotImplemented,
"LagrangeSimplexLocalBasis only implemented for dim==3!");
287 for (i[2] = 0; i[2] <= k; ++i[2])
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])
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])
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];
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)
312 for (
unsigned int l = 0; l < i[3]; ++l)
319 for (
unsigned int j = 0; j < i[3]; ++j)
320 factor[3] *= kx3 - j;
321 for (
unsigned int m = 0; m < 3; ++m)
324 for (
unsigned int j = 0; j < i[m]; ++j)
327 for (
unsigned int p = 0; p < 3; ++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);
335 out[n][0][m] += prod;
350 void partial(
const std::array<unsigned int,dim>& order,
352 std::vector<typename Traits::RangeType>& out)
const
358 if (totalOrder == 0) {
359 evaluateFunction(in, out);
373 auto direction = std::find(order.begin(), order.end(), 1);
376 for (
unsigned int i=0; i<dim; i++)
377 out[i+1] = (i==(direction-order.begin()));
380 std::fill(out.begin(), out.end(), 0);
386 auto lagrangeNode = [](
unsigned int i) {
return ((D)i)/k; };
389 auto lagrangianFactor = [&lagrangeNode]
394 return (x[0]-lagrangeNode(no))/(lagrangeNode(i)-lagrangeNode(no));
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));
402 auto lagrangianFactorDerivative = [&lagrangeNode]
403 (
const int direction,
const int no,
const int i,
const int j,
const typename Traits::DomainType&)
408 return (direction == 0) ? T(1.0/(lagrangeNode(i)-lagrangeNode(no))) : T(0);
411 return (direction == 0) ? T(0) : T(1.0/(lagrangeNode(j)-lagrangeNode(no-i)));
413 return -1.0/(lagrangeNode(no+1)-lagrangeNode(i)-lagrangeNode(j));
418 int direction = std::find(order.begin(), order.end(), 1)-order.begin();
421 for (
unsigned int j=0; j<=k; j++)
423 for (
unsigned int i=0; i<=k-j; i++, n++)
426 for (
unsigned int no1=0; no1 < k; no1++)
428 R factor = lagrangianFactorDerivative(direction, no1, i, j, in);
429 for (
unsigned int no2=0; no2 < k; no2++)
431 factor *= lagrangianFactor(no2, i, j, in);
442 std::array<int,2> directions;
443 unsigned int counter = 0;
444 auto nonconstOrder = order;
445 for (
int i=0; i<2; i++)
447 while (nonconstOrder[i])
449 directions[counter++] = i;
456 for (
unsigned int j=0; j<=k; j++)
458 for (
unsigned int i=0; i<=k-j; i++, n++)
462 for (
unsigned int no1=0; no1 < k; no1++)
464 R factor1 = lagrangianFactorDerivative(directions[0], no1, i, j, in);
465 for (
unsigned int no2=0; no2 < k; no2++)
469 R factor2 = factor1*lagrangianFactorDerivative(directions[1], no2, i, j, in);
470 for (
unsigned int no3=0; no3 < k; no3++)
472 if (no3 == no1 || no3 == no2)
474 factor2 *= lagrangianFactor(no3, i, j, in);
488 DUNE_THROW(NotImplemented,
"Desired derivative order is not implemented");
492 static constexpr unsigned int order ()
503 template<
unsigned int dim,
unsigned int k>
504 class LagrangeSimplexLocalCoefficients
508 LagrangeSimplexLocalCoefficients ()
513 localKeys_[0] = LocalKey(0,0,0);
519 for (std::size_t i=0; i<size(); i++)
520 localKeys_[i] = LocalKey(i,dim,0);
527 localKeys_[0] = LocalKey(0,1,0);
528 for (
unsigned int i=1; i<k; i++)
529 localKeys_[i] = LocalKey(0,0,i-1);
530 localKeys_.back() = LocalKey(1,1,0);
538 for (
unsigned int j=0; j<=k; j++)
539 for (
unsigned int i=0; i<=k-j; i++)
543 localKeys_[n++] = LocalKey(0,2,0);
548 localKeys_[n++] = LocalKey(1,2,0);
553 localKeys_[n++] = LocalKey(2,2,0);
558 localKeys_[n++] = LocalKey(0,1,i-1);
563 localKeys_[n++] = LocalKey(1,1,j-1);
568 localKeys_[n++] = LocalKey(2,1,j-1);
571 localKeys_[n++] = LocalKey(0,0,c++);
578 std::array<unsigned int, dim+1> vertexMap;
579 for (
unsigned int i=0; i<=dim; i++)
581 generateLocalKeys(vertexMap);
584 DUNE_THROW(NotImplemented,
"LagrangeSimplexLocalCoefficients only implemented for k<=1 or dim<=3!");
593 LagrangeSimplexLocalCoefficients (
const std::array<unsigned int, dim+1> vertexMap)
596 if (dim!=2 && dim!=3)
597 DUNE_THROW(NotImplemented,
"LagrangeSimplexLocalCoefficients only implemented for dim==2 and dim==3!");
599 generateLocalKeys(vertexMap);
603 template<
class VertexMap>
604 LagrangeSimplexLocalCoefficients(
const VertexMap &vertexmap)
607 if (dim!=2 && dim!=3)
608 DUNE_THROW(NotImplemented,
"LagrangeSimplexLocalCoefficients only implemented for dim==2 and dim==3!");
610 std::array<unsigned int, dim+1> vertexmap_array;
611 std::copy(vertexmap, vertexmap + dim + 1, vertexmap_array.begin());
612 generateLocalKeys(vertexmap_array);
616 static constexpr std::size_t size ()
622 const LocalKey& localKey (std::size_t i)
const
624 return localKeys_[i];
628 std::vector<LocalKey> localKeys_;
630 void generateLocalKeys(
const std::array<unsigned int, dim+1> vertexMap)
634 localKeys_[0] = LocalKey(0,0,0);
643 for (
unsigned int j=0; j<=k; j++)
644 for (
unsigned int i=0; i<=k-j; i++)
648 localKeys_[n++] = LocalKey(0,2,0);
653 localKeys_[n++] = LocalKey(1,2,0);
658 localKeys_[n++] = LocalKey(2,2,0);
663 localKeys_[n++] = LocalKey(0,1,i-1);
668 localKeys_[n++] = LocalKey(1,1,j-1);
673 localKeys_[n++] = LocalKey(2,1,j-1);
676 localKeys_[n++] = LocalKey(0,0,c++);
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());
692 DUNE_THROW(NotImplemented,
"LagrangeSimplexLocalCoefficients only implemented for dim==3!");
694 unsigned int subindex[16];
695 unsigned int codim_count[4] = {0};
696 for (
unsigned int m = 1; m < 16; ++m)
698 unsigned int codim = !(m&1) + !(m&2) + !(m&4) + !(m&8);
699 subindex[m] = codim_count[codim]++;
702 int a1 = (3*k + 12)*k + 11;
704 unsigned int dof_count[16] = {0};
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])
710 i[0] = k - i[1] - i[2] - i[3];
712 unsigned int entity = 0;
713 unsigned int codim = 0;
714 for (
unsigned int m = 0; m < 4; ++m)
716 j[m] = i[vertexMap[m]];
717 entity += !!j[m] << m;
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]++);
731 template<
class LocalBasis>
732 class LagrangeSimplexLocalInterpolation
734 static const int kdiv = (LocalBasis::order() == 0 ? 1 : LocalBasis::order());
744 template<
typename F,
typename C>
745 void interpolate (
const F& ff, std::vector<C>& out)
const
747 constexpr auto dim = LocalBasis::Traits::dimDomain;
748 constexpr auto k = LocalBasis::order();
749 using D =
typename LocalBasis::Traits::DomainFieldType;
751 typename LocalBasis::Traits::DomainType x;
752 auto&& f = Impl::makeFunctionWithCallOperator<typename LocalBasis::Traits::DomainType>(ff);
754 out.resize(LocalBasis::size());
768 std::fill(x.begin(), x.end(), 0);
772 for (
int i=0; i<dim; i++)
774 for (
int j=0; j<dim; j++)
784 for (
unsigned int i=0; i<k+1; i++)
795 for (
unsigned int j=0; j<=k; j++)
796 for (
unsigned int i=0; i<=k-j; i++)
798 x = { ((D)i)/k, ((D)j)/k };
806 DUNE_THROW(NotImplemented,
"LagrangeSimplexLocalInterpolation only implemented for dim<=3!");
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++)
813 x[0] = ((D)i0)/((D)kdiv);
814 x[1] = ((D)i1)/((D)kdiv);
815 x[2] = ((D)i2)/((D)kdiv);
834 template<
class D,
class R,
int d,
int k>
841 Impl::LagrangeSimplexLocalCoefficients<d,k>,
842 Impl::LagrangeSimplexLocalInterpolation<Impl::LagrangeSimplexLocalBasis<D,R,d,k> > >;
851 template<
typename VertexMap>
853 : coefficients_(vertexmap)
867 return coefficients_;
874 return interpolation_;
878 static constexpr std::size_t
size ()
880 return Impl::LagrangeSimplexLocalBasis<D,R,d,k>::size();
891 Impl::LagrangeSimplexLocalBasis<D,R,d,k> basis_;
892 Impl::LagrangeSimplexLocalCoefficients<d,k> coefficients_;
893 Impl::LagrangeSimplexLocalInterpolation<Impl::LagrangeSimplexLocalBasis<D,R,d,k> > interpolation_;
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:123
Lagrange finite element for simplices with arbitrary compile-time dimension and polynomial order.
Definition: lagrangesimplex.hh:836
const Traits::LocalInterpolationType & localInterpolation() const
Returns object that evaluates degrees of freedom.
Definition: lagrangesimplex.hh:872
const Traits::LocalBasisType & localBasis() const
Returns the local basis, i.e., the set of shape functions.
Definition: lagrangesimplex.hh:858
LagrangeSimplexLocalFiniteElement()
Definition: lagrangesimplex.hh:845
static constexpr std::size_t size()
The number of shape functions.
Definition: lagrangesimplex.hh:878
LagrangeSimplexLocalFiniteElement(const VertexMap &vertexmap)
Definition: lagrangesimplex.hh:852
const Traits::LocalCoefficientsType & localCoefficients() const
Returns the assignment of the degrees of freedom to the element subentities.
Definition: lagrangesimplex.hh:865
static constexpr GeometryType type()
The reference element that the local finite element is defined on.
Definition: lagrangesimplex.hh:885
Definition of the DUNE_DEPRECATED macro for the case that config.h is not available.
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.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
constexpr GeometryType simplex(unsigned int dim)
Returns a GeometryType representing a simplex of dimension dim.
Definition: type.hh:461
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:289
Some useful basic math stuff.
Dune namespace.
Definition: alignedallocator.hh:11
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
static const ReferenceElement & simplex()
get simplex reference elements
Definition: referenceelements.hh:202
D DomainType
domain type
Definition: localbasis.hh:43
R RangeType
range type
Definition: localbasis.hh:55
traits helper struct
Definition: localfiniteelementtraits.hh:11
LB LocalBasisType
Definition: localfiniteelementtraits.hh:14
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:18
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:22