From fc08ee91f8707ceed68750ca230a2e1fd7405f4d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= <graeser@dune-project.org>
Date: Wed, 19 Mar 2025 15:20:10 +0100
Subject: [PATCH 1/4] Unify LagrangeSimplexFiniteElement across dimensions
The `LagrangeSimplexFiniteElement` contains 5 different code paths
to implement the basis and its derivatives:
* A special implementation for order 0 in any dimension.
* A special implementation for order 1 in any dimension.
* An implementation for any order in dimension 1.
* An implementation for any order in dimension 2.
* An implementation for any order in dimension 3.
There are good reasons for the first two special cases.
Unfortunately the ones for dim=1,2,3 are all implemented
differently for no good reason (except history) makes the
code hard to read and understand.
This patch unifies the latter three implementations. While
there are still three branches, there is a very simple common
structure making the now making the code much simpler to understand.
The unified implementation basically follows the former
`evaluateFunction()` implementation of the 3d case:
1. Rescale the simplex by `k` such that Lagrange points have integer
coordinates.
2. Transform to barycentric coordinates in the scaled simplex.
3. Loop over Lagrange points in lexicographic order.
4. Compute basis functions as products of univariate
Lagrange-polynomials wrt. barycentric coordinates.
This approach is now also used for `evaluateJacobian()`
and first order derivatives in `partial()` extending the
support of the latter to dim=1 and dim=2.
For now the special implementation of second order
partial derivatives in 2d is kept, but could be updated
accordingly.
---
.../lagrange/lagrangesimplex.hh | 421 +++++++++---------
1 file changed, 200 insertions(+), 221 deletions(-)
diff --git a/dune/localfunctions/lagrange/lagrangesimplex.hh b/dune/localfunctions/lagrange/lagrangesimplex.hh
index b59a922b..e85f7d55 100644
--- a/dune/localfunctions/lagrange/lagrangesimplex.hh
+++ b/dune/localfunctions/lagrange/lagrangesimplex.hh
@@ -33,6 +33,65 @@ namespace Dune { namespace Impl
template<class D, class R, unsigned int dim, unsigned int k>
class LagrangeSimplexLocalBasis
{
+
+ // Compute the rescaled barycentric coordinates of x.
+ // We rescale the simplex by k and then compute the
+ // barycentric coordinates with respect to the points
+ // p_i = e_i (for i=0,...,dim-1) and p_dim=0.
+ // Notice that then the Lagrange points have the barycentric
+ // coordinates (i_0,...,i_d) where i_j are all non-negative
+ // integers satisfying the constraint sum i_j = k.
+ constexpr auto barycentric(const auto& x) const
+ {
+ auto b = std::array<R,dim+1>{};
+ b[dim] = k;
+ for(auto i : Dune::range(dim))
+ {
+ b[i] = k*x[i];
+ b[dim] -= b[i];
+ }
+ return b;
+ }
+
+ // Evaluate univariate Lagrange factor Li(t) = (t-0)/(i-0) * ... * (t-(i-1))/(i-(i-1))
+ static constexpr R lagrangeFactor(unsigned int i, const R& t)
+ {
+ auto y = R(1);
+ for (unsigned int j = 0; j < i; ++j)
+ y *= (t-j) / (i-j);
+ return y;
+ }
+
+ // Evaluate univariate Lagrange factor Li(t) = (t-0)/(i-0) * ... * (t-(i-1))/(i-(i-1))
+ // and its derivatives up to given order. Currently only maxDerivativeOrder 0 and 1
+ // are supported.
+ static constexpr auto lagrangeFactor(unsigned int i, const R& t, unsigned int maxDerivativeOrder)
+ {
+ auto y = std::array<R, 2>{R(1), R(0)};
+ if (i==0)
+ return y;
+
+ // Split of the zero-th factor, since we can reuse it for the derivative.
+ auto tmp = R(1);
+ for (unsigned int j = 1; j < i; ++j)
+ tmp *= (t-j) / (i-j);
+ y[0] = (t/i) * tmp;
+
+ if (maxDerivativeOrder>0)
+ {
+ y[1] = R(1)/i * tmp;
+ for (unsigned int l = 1; l < i; ++l)
+ {
+ auto yl = R(1)/(i-l);
+ for (unsigned int j = 0; j < i; ++j)
+ if (j!=l)
+ yl *= (t-j) / (i-j);
+ y[1] += yl;
+ }
+ }
+ return y;
+ }
+
public:
using Traits = LocalBasisTraits<D,dim,FieldVector<D,dim>,R,1,FieldVector<R,1>,FieldMatrix<R,1,dim> >;
@@ -70,74 +129,65 @@ namespace Dune { namespace Impl
return;
}
- assert(k>=2);
+ // Compute rescaled barycentric coordinates of x
+ auto z = barycentric(x);
- auto lagrangeNode = [](unsigned int i) { return ((D)i)/k; };
+ // To improve code readability we introduce a short cut
+ constexpr auto d = dim;
if (dim==1)
{
- for (unsigned int i=0; i<size(); i++)
+ unsigned int n = 0;
+ for (auto i0 : Dune::range(k + 1))
{
- out[i] = 1.0;
- for (unsigned int alpha=0; alpha<i; alpha++)
- out[i] *= (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
- for (unsigned int gamma=i+1; gamma<=k; gamma++)
- out[i] *= (x[0]-lagrangeNode(gamma))/(lagrangeNode(i)-lagrangeNode(gamma));
+ auto id = k - i0;
+ auto Lz0 = lagrangeFactor(i0, z[0]);
+ auto Lzd = lagrangeFactor(id, z[d]);
+ out[n] = Lz0 * Lzd;
+ ++n;
}
return;
}
-
if (dim==2)
{
- int n=0;
- for (unsigned int j=0; j<=k; j++)
- for (unsigned int i=0; i<=k-j; i++)
+ unsigned int n=0;
+ for (auto i1 : Dune::range(k + 1))
+ {
+ auto Lz1 = lagrangeFactor(i1, z[1]);
+ for (auto i0 : Dune::range(k - i1 + 1))
{
- out[n] = 1.0;
- for (unsigned int alpha=0; alpha<i; alpha++)
- out[n] *= (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
- for (unsigned int beta=0; beta<j; beta++)
- out[n] *= (x[1]-lagrangeNode(beta))/(lagrangeNode(j)-lagrangeNode(beta));
- for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
- out[n] *= (lagrangeNode(gamma)-x[0]-x[1])/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
- n++;
+ auto id = k - i0 - i1;
+ auto Lz0 = lagrangeFactor(i0, z[0]);
+ auto Lzd = lagrangeFactor(id, z[d]);
+ out[n] = Lz0 * Lz1 * Lzd;
+ ++n;
}
-
+ }
return;
}
-
- if (dim!=3)
- DUNE_THROW(NotImplemented, "LagrangeSimplexLocalBasis for k>=2 only implemented for dim==1 or dim==3");
-
- typename Traits::DomainType kx = x;
- kx *= k;
- unsigned int n = 0;
- unsigned int i[4];
- R factor[4];
- for (i[2] = 0; i[2] <= k; ++i[2])
+ if (dim==3)
{
- factor[2] = 1.0;
- for (unsigned int j = 0; j < i[2]; ++j)
- factor[2] *= (kx[2]-j) / (i[2]-j);
- for (i[1] = 0; i[1] <= k - i[2]; ++i[1])
+ unsigned int n = 0;
+ for (auto i2 : Dune::range(k + 1))
{
- factor[1] = 1.0;
- for (unsigned int j = 0; j < i[1]; ++j)
- factor[1] *= (kx[1]-j) / (i[1]-j);
- for (i[0] = 0; i[0] <= k - i[1] - i[2]; ++i[0])
+ auto Lz2 = lagrangeFactor(i2, z[2]);
+ for (auto i1 : Dune::range(k - i2 + 1))
{
- factor[0] = 1.0;
- for (unsigned int j = 0; j < i[0]; ++j)
- factor[0] *= (kx[0]-j) / (i[0]-j);
- i[3] = k - i[0] - i[1] - i[2];
- D kx3 = k - kx[0] - kx[1] - kx[2];
- factor[3] = 1.0;
- for (unsigned int j = 0; j < i[3]; ++j)
- factor[3] *= (kx3-j) / (i[3]-j);
- out[n++] = factor[0] * factor[1] * factor[2] * factor[3];
+ auto Lz1 = lagrangeFactor(i1, z[1]);
+ for (auto i0 : Dune::range(k - i2 - i1 + 1))
+ {
+ auto id = k - i0 - i1 -i2;
+ auto Lz0 = lagrangeFactor(i0, z[0]);
+ auto Lzd = lagrangeFactor(id, z[d]);
+ out[n] = Lz0 * Lz1 * Lz2 * Lzd;
+ ++n;
+ }
}
}
+ return;
}
+
+ DUNE_THROW(NotImplemented, "LagrangeSimplexLocalBasis for k>=2 only implemented for dim<=3");
}
/** \brief Evaluate Jacobian of all shape functions
@@ -169,176 +219,68 @@ namespace Dune { namespace Impl
return;
}
- auto lagrangeNode = [](unsigned int i) { return ((D)i)/k; };
+ // Compute rescaled barycentric coordinates of x
+ auto z = barycentric(x);
+
+ // To improve code readability we introduce a short cut
+ constexpr auto d = dim;
- // Specialization for dim==1
if (dim==1)
{
- for (unsigned int i=0; i<=k; i++)
+ unsigned int n = 0;
+ for (auto i0 : Dune::range(k + 1))
{
- // x_0 derivative
- out[i][0][0] = 0.0;
- R factor=1.0;
- for (unsigned int a=0; a<i; a++)
- {
- R product=factor;
- for (unsigned int alpha=0; alpha<i; alpha++)
- product *= (alpha==a) ? 1.0/(lagrangeNode(i)-lagrangeNode(alpha))
- : (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
- for (unsigned int gamma=i+1; gamma<=k; gamma++)
- product *= (lagrangeNode(gamma)-x[0])/(lagrangeNode(gamma)-lagrangeNode(i));
- out[i][0][0] += product;
- }
- for (unsigned int c=i+1; c<=k; c++)
- {
- R product=factor;
- for (unsigned int alpha=0; alpha<i; alpha++)
- product *= (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
- for (unsigned int gamma=i+1; gamma<=k; gamma++)
- product *= (gamma==c) ? -1.0/(lagrangeNode(gamma)-lagrangeNode(i))
- : (lagrangeNode(gamma)-x[0])/(lagrangeNode(gamma)-lagrangeNode(i));
- out[i][0][0] += product;
- }
+ auto id = k - i0;
+ auto Lz0 = lagrangeFactor(i0, z[0], 1);
+ auto Lzd = lagrangeFactor(id, z[d], 1);
+ out[n][0][0] = (Lz0[1] * Lzd[0] - Lz0[0] * Lzd[1])*k;
+ ++n;
}
return;
}
-
if (dim==2)
{
- int n=0;
- for (unsigned int j=0; j<=k; j++)
- for (unsigned int i=0; i<=k-j; i++)
+ unsigned int n=0;
+ for (auto i1 : Dune::range(k + 1))
+ {
+ auto Lz1 = lagrangeFactor(i1, z[1], 1);
+ for (auto i0 : Dune::range(k - i1 + 1))
{
- // x_0 derivative
- out[n][0][0] = 0.0;
- R factor=1.0;
- for (unsigned int beta=0; beta<j; beta++)
- factor *= (x[1]-lagrangeNode(beta))/(lagrangeNode(j)-lagrangeNode(beta));
- for (unsigned int a=0; a<i; a++)
- {
- R product=factor;
- for (unsigned int alpha=0; alpha<i; alpha++)
- if (alpha==a)
- product *= D(1)/(lagrangeNode(i)-lagrangeNode(alpha));
- else
- product *= (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
- for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
- product *= (lagrangeNode(gamma)-x[0]-x[1])/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
- out[n][0][0] += product;
- }
- for (unsigned int c=i+j+1; c<=k; c++)
- {
- R product=factor;
- for (unsigned int alpha=0; alpha<i; alpha++)
- product *= (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
- for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
- if (gamma==c)
- product *= -D(1)/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
- else
- product *= (lagrangeNode(gamma)-x[0]-x[1])/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
- out[n][0][0] += product;
- }
-
- // x_1 derivative
- out[n][0][1] = 0.0;
- factor = 1.0;
- for (unsigned int alpha=0; alpha<i; alpha++)
- factor *= (x[0]-lagrangeNode(alpha))/(lagrangeNode(i)-lagrangeNode(alpha));
- for (unsigned int b=0; b<j; b++)
- {
- R product=factor;
- for (unsigned int beta=0; beta<j; beta++)
- if (beta==b)
- product *= D(1)/(lagrangeNode(j)-lagrangeNode(beta));
- else
- product *= (x[1]-lagrangeNode(beta))/(lagrangeNode(j)-lagrangeNode(beta));
- for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
- product *= (lagrangeNode(gamma)-x[0]-x[1])/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
- out[n][0][1] += product;
- }
- for (unsigned int c=i+j+1; c<=k; c++)
- {
- R product=factor;
- for (unsigned int beta=0; beta<j; beta++)
- product *= (x[1]-lagrangeNode(beta))/(lagrangeNode(j)-lagrangeNode(beta));
- for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
- if (gamma==c)
- product *= -D(1)/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
- else
- product *= (lagrangeNode(gamma)-x[0]-x[1])/(lagrangeNode(gamma)-lagrangeNode(i)-lagrangeNode(j));
- out[n][0][1] += product;
- }
-
- n++;
+ auto id = k - i0 - i1;
+ auto Lz0 = lagrangeFactor(i0, z[0], 1);
+ auto Lzd = lagrangeFactor(id, z[d], 1);
+ out[n][0][0] = (Lz0[1] * Lz1[0] * Lzd[0] - Lz0[0] * Lz1[0] * Lzd[1])*k;
+ out[n][0][1] = (Lz0[0] * Lz1[1] * Lzd[0] - Lz0[0] * Lz1[0] * Lzd[1])*k;
+ ++n;
}
-
+ }
return;
}
-
- if (dim!=3)
- DUNE_THROW(NotImplemented, "LagrangeSimplexLocalBasis only implemented for dim==3!");
-
- // Specialization for arbitrary order and dim==3
- typename Traits::DomainType kx = x;
- kx *= k;
- unsigned int n = 0;
- unsigned int i[4];
- R factor[4];
- for (i[2] = 0; i[2] <= k; ++i[2])
+ if (dim==3)
{
- factor[2] = 1.0;
- for (unsigned int j = 0; j < i[2]; ++j)
- factor[2] *= (kx[2]-j) / (i[2]-j);
- for (i[1] = 0; i[1] <= k - i[2]; ++i[1])
+ unsigned int n = 0;
+ for (auto i2 : Dune::range(k + 1))
{
- factor[1] = 1.0;
- for (unsigned int j = 0; j < i[1]; ++j)
- factor[1] *= (kx[1]-j) / (i[1]-j);
- for (i[0] = 0; i[0] <= k - i[1] - i[2]; ++i[0])
+ auto Lz2 = lagrangeFactor(i2, z[2], 1);
+ for (auto i1 : Dune::range(k - i2 + 1))
{
- factor[0] = 1.0;
- for (unsigned int j = 0; j < i[0]; ++j)
- factor[0] *= (kx[0]-j) / (i[0]-j);
- i[3] = k - i[0] - i[1] - i[2];
- D kx3 = k - kx[0] - kx[1] - kx[2];
- R sum3 = 0.0;
- factor[3] = 1.0;
- for (unsigned int j = 0; j < i[3]; ++j)
- factor[3] /= i[3] - j;
- R prod_all = factor[0] * factor[1] * factor[2] * factor[3];
- for (unsigned int j = 0; j < i[3]; ++j)
+ auto Lz1 = lagrangeFactor(i1, z[1], 1);
+ for (auto i0 : Dune::range(k - i2 - i1 + 1))
{
- R prod = prod_all;
- for (unsigned int l = 0; l < i[3]; ++l)
- if (j == l)
- prod *= -R(k);
- else
- prod *= kx3 - l;
- sum3 += prod;
+ auto id = k - i0 - i1 -i2;
+ auto Lz0 = lagrangeFactor(i0, z[0], 1);
+ auto Lzd = lagrangeFactor(id, z[d], 1);
+ out[n][0][0] = (Lz0[1] * Lz1[0] * Lz2[0] * Lzd[0] - Lz0[0] * Lz1[0] * Lz2[0] * Lzd[1])*k;
+ out[n][0][1] = (Lz0[0] * Lz1[1] * Lz2[0] * Lzd[0] - Lz0[0] * Lz1[0] * Lz2[0] * Lzd[1])*k;
+ out[n][0][2] = (Lz0[0] * Lz1[0] * Lz2[1] * Lzd[0] - Lz0[0] * Lz1[0] * Lz2[0] * Lzd[1])*k;
+ ++n;
}
- for (unsigned int j = 0; j < i[3]; ++j)
- factor[3] *= kx3 - j;
- for (unsigned int m = 0; m < 3; ++m)
- {
- out[n][0][m] = sum3;
- for (unsigned int j = 0; j < i[m]; ++j)
- {
- R prod = factor[3];
- for (unsigned int p = 0; p < 3; ++p)
- {
- if (m == p)
- for (unsigned int l = 0; l < i[p]; ++l)
- prod *= (j==l) ? R(k) / (i[p]-l) : R(kx[p]-l) / (i[p]-l);
- else
- prod *= factor[p];
- }
- out[n][0][m] += prod;
- }
- }
- n++;
}
}
+ return;
}
+
+ DUNE_THROW(NotImplemented, "LagrangeSimplexLocalBasis for k>=2 only implemented for dim<=3");
}
/** \brief Evaluate partial derivatives of any order of all shape functions
@@ -381,6 +323,67 @@ namespace Dune { namespace Impl
return;
}
+ if (totalOrder==1)
+ {
+ // Compute rescaled barycentric coordinates of x
+ auto z = barycentric(in);
+
+ // To improve code readability we introduce a short cut
+ constexpr auto d = dim;
+
+ if (dim==1)
+ {
+ unsigned int n = 0;
+ for (auto i0 : Dune::range(k + 1))
+ {
+ auto id = k - i0;
+ auto Lz0 = lagrangeFactor(i0, z[0], 1);
+ auto Lzd = lagrangeFactor(id, z[d], 1);
+ out[n] = (Lz0[1] * Lzd[0] - Lz0[0] * Lzd[1])*k;
+ ++n;
+ }
+ return;
+ }
+ if (dim==2)
+ {
+ unsigned int n=0;
+ for (auto i1 : Dune::range(k + 1))
+ {
+ auto Lz1 = lagrangeFactor(i1, z[1], order[1]);
+ for (auto i0 : Dune::range(k - i1 + 1))
+ {
+ auto id = k - i0 - i1;
+ auto Lz0 = lagrangeFactor(i0, z[0], order[0]);
+ auto Lzd = lagrangeFactor(id, z[d], 1);
+ out[n] = (Lz0[order[0]] * Lz1[order[1]] * Lzd[0] - Lz0[0] * Lz1[0] * Lzd[1])*k;
+ ++n;
+ }
+ }
+ return;
+ }
+ if (dim==3)
+ {
+ unsigned int n = 0;
+ for (auto i2 : Dune::range(k + 1))
+ {
+ auto Lz2 = lagrangeFactor(i2, z[2], order[2]);
+ for (auto i1 : Dune::range(k - i2 + 1))
+ {
+ auto Lz1 = lagrangeFactor(i1, z[1], order[1]);
+ for (auto i0 : Dune::range(k - i2 - i1 + 1))
+ {
+ auto id = k - i0 - i1 -i2;
+ auto Lz0 = lagrangeFactor(i0, z[0], order[0]);
+ auto Lzd = lagrangeFactor(id, z[d], 1);
+ out[n] = (Lz0[order[0]] * Lz1[order[1]] * Lz2[order[2]] * Lzd[0] - Lz0[0] * Lz1[0] * Lz2[0] * Lzd[1])*k;
+ ++n;
+ }
+ }
+ }
+ return;
+ }
+ }
+
if (dim==2)
{
auto lagrangeNode = [](unsigned int i) { return ((D)i)/k; };
@@ -413,30 +416,6 @@ namespace Dune { namespace Impl
return -1.0/(lagrangeNode(no+1)-lagrangeNode(i)-lagrangeNode(j));
};
- if (totalOrder==1)
- {
- int direction = std::find(order.begin(), order.end(), 1)-order.begin();
-
- int n=0;
- for (unsigned int j=0; j<=k; j++)
- {
- for (unsigned int i=0; i<=k-j; i++, n++)
- {
- out[n] = 0.0;
- for (unsigned int no1=0; no1 < k; no1++)
- {
- R factor = lagrangianFactorDerivative(direction, no1, i, j, in);
- for (unsigned int no2=0; no2 < k; no2++)
- if (no1 != no2)
- factor *= lagrangianFactor(no2, i, j, in);
-
- out[n] += factor;
- }
- }
- }
- return;
- }
-
if (totalOrder==2)
{
std::array<int,2> directions;
--
GitLab
From 853580c47ed15b035264f8fa4fc91a3abed2fb73 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= <graeser@dune-project.org>
Date: Fri, 21 Mar 2025 13:14:54 +0100
Subject: [PATCH 2/4] Precompute barycentric factors in LagrangeSimplexBasis
Instead of evaluating the barycentric factors during the loop
over the Lagrange nodes, it is simpler and more efficient to
do this once in advance. This also adds the prcisice formula
for the basis functions to the doxygen documentation.
---
.../lagrange/lagrangesimplex.hh | 227 +++++++++---------
1 file changed, 109 insertions(+), 118 deletions(-)
diff --git a/dune/localfunctions/lagrange/lagrangesimplex.hh b/dune/localfunctions/lagrange/lagrangesimplex.hh
index e85f7d55..d3d96811 100644
--- a/dune/localfunctions/lagrange/lagrangesimplex.hh
+++ b/dune/localfunctions/lagrange/lagrangesimplex.hh
@@ -23,7 +23,7 @@ namespace Dune { namespace Impl
/** \brief Lagrange shape functions of arbitrary order on the reference simplex
Lagrange shape functions of arbitrary order have the property that
- \f$\hat\phi^i(x_j) = \delta_{i,j}\f$ for certain points \f$x_j\f$.
+ \f$\hat\phi^i(x^{(j)}) = \delta_{i,j}\f$ for certain points \f$x^{(j)}\f$.
\tparam D Type to represent the field in the domain
\tparam R Type to represent the field in the range
@@ -53,43 +53,32 @@ namespace Dune { namespace Impl
return b;
}
- // Evaluate univariate Lagrange factor Li(t) = (t-0)/(i-0) * ... * (t-(i-1))/(i-(i-1))
- static constexpr R lagrangeFactor(unsigned int i, const R& t)
+ // Evaluate the univariate Lagrange polynomials L_i(t) for i=0,...,k where
+ //
+ // L_i(t) = (t-0)/(i-0) * ... * (t-(i-1))/(i-(i-1))
+ // = (t-0)*...*(t-(i-1))/(i!);
+ static constexpr void evaluateLagrangePolynomials(const R& t, auto& L)
{
- auto y = R(1);
- for (unsigned int j = 0; j < i; ++j)
- y *= (t-j) / (i-j);
- return y;
+ L[0] = 1;
+ for (auto i : Dune::range(k))
+ L[i+1] = L[i] * (t - i) / (i+1);
}
- // Evaluate univariate Lagrange factor Li(t) = (t-0)/(i-0) * ... * (t-(i-1))/(i-(i-1))
- // and its derivatives up to given order. Currently only maxDerivativeOrder 0 and 1
- // are supported.
- static constexpr auto lagrangeFactor(unsigned int i, const R& t, unsigned int maxDerivativeOrder)
+ // Evaluate the univariate Lagrange polynomial derivatives L_i(t) for i=0,...,k
+ // up to given maxDerivativeOrder. Currently only maxDerivativeOrder 0 and 1 are supported.
+ static constexpr void evaluateLagrangePolynomialDerivative(const R& t, auto& LL, unsigned int maxDerivativeOrder)
{
- auto y = std::array<R, 2>{R(1), R(0)};
- if (i==0)
- return y;
-
- // Split of the zero-th factor, since we can reuse it for the derivative.
- auto tmp = R(1);
- for (unsigned int j = 1; j < i; ++j)
- tmp *= (t-j) / (i-j);
- y[0] = (t/i) * tmp;
-
+ auto& L = LL[0];
+ L[0] = 1;
+ for (auto i : Dune::range(k))
+ L[i+1] = L[i] * (t - i) / (i+1);
+ auto& DL = LL[1];
if (maxDerivativeOrder>0)
{
- y[1] = R(1)/i * tmp;
- for (unsigned int l = 1; l < i; ++l)
- {
- auto yl = R(1)/(i-l);
- for (unsigned int j = 0; j < i; ++j)
- if (j!=l)
- yl *= (t-j) / (i-j);
- y[1] += yl;
- }
+ DL[0] = 0;
+ for (auto i : Dune::range(k))
+ DL[i+1] = (DL[i] * (t - i) + L[i]) / (i+1);
}
- return y;
}
public:
@@ -132,58 +121,35 @@ namespace Dune { namespace Impl
// Compute rescaled barycentric coordinates of x
auto z = barycentric(x);
- // To improve code readability we introduce a short cut
- constexpr auto d = dim;
+ auto L = std::array<std::array<R,k+1>, dim+1>();
+ for (auto j : Dune::range(dim+1))
+ evaluateLagrangePolynomials(z[j], L[j]);
if (dim==1)
{
unsigned int n = 0;
for (auto i0 : Dune::range(k + 1))
- {
- auto id = k - i0;
- auto Lz0 = lagrangeFactor(i0, z[0]);
- auto Lzd = lagrangeFactor(id, z[d]);
- out[n] = Lz0 * Lzd;
- ++n;
- }
+ for (auto i1 : std::array{k - i0})
+ out[n++] = L[0][i0] * L[1][i1];
return;
}
if (dim==2)
{
unsigned int n=0;
for (auto i1 : Dune::range(k + 1))
- {
- auto Lz1 = lagrangeFactor(i1, z[1]);
for (auto i0 : Dune::range(k - i1 + 1))
- {
- auto id = k - i0 - i1;
- auto Lz0 = lagrangeFactor(i0, z[0]);
- auto Lzd = lagrangeFactor(id, z[d]);
- out[n] = Lz0 * Lz1 * Lzd;
- ++n;
- }
- }
+ for (auto i2 : std::array{k - i1 - i0})
+ out[n++] = L[0][i0] * L[1][i1] * L[2][i2];
return;
}
if (dim==3)
{
unsigned int n = 0;
for (auto i2 : Dune::range(k + 1))
- {
- auto Lz2 = lagrangeFactor(i2, z[2]);
for (auto i1 : Dune::range(k - i2 + 1))
- {
- auto Lz1 = lagrangeFactor(i1, z[1]);
for (auto i0 : Dune::range(k - i2 - i1 + 1))
- {
- auto id = k - i0 - i1 -i2;
- auto Lz0 = lagrangeFactor(i0, z[0]);
- auto Lzd = lagrangeFactor(id, z[d]);
- out[n] = Lz0 * Lz1 * Lz2 * Lzd;
- ++n;
- }
- }
- }
+ for (auto i3 : std::array{k - i2 - i1 - i0})
+ out[n++] = L[0][i0] * L[1][i1] * L[2][i2] * L[3][i3];
return;
}
@@ -222,19 +188,21 @@ namespace Dune { namespace Impl
// Compute rescaled barycentric coordinates of x
auto z = barycentric(x);
- // To improve code readability we introduce a short cut
- constexpr auto d = dim;
+ // L[j][m][i] is the m-th derivative of the i-th Lagrange polynomial at z[j]
+ auto L = std::array<std::array<std::array<R,k+1>, 2>, dim+1>();
+ for (auto j : Dune::range(dim+1))
+ evaluateLagrangePolynomialDerivative(z[j], L[j], 1);
if (dim==1)
{
unsigned int n = 0;
for (auto i0 : Dune::range(k + 1))
{
- auto id = k - i0;
- auto Lz0 = lagrangeFactor(i0, z[0], 1);
- auto Lzd = lagrangeFactor(id, z[d], 1);
- out[n][0][0] = (Lz0[1] * Lzd[0] - Lz0[0] * Lzd[1])*k;
- ++n;
+ for (auto i1 : std::array{k-i0})
+ {
+ out[n][0][0] = (L[0][1][i0] * L[1][0][i1] - L[0][0][i0] * L[1][1][i1])*k;
+ ++n;
+ }
}
return;
}
@@ -243,15 +211,14 @@ namespace Dune { namespace Impl
unsigned int n=0;
for (auto i1 : Dune::range(k + 1))
{
- auto Lz1 = lagrangeFactor(i1, z[1], 1);
for (auto i0 : Dune::range(k - i1 + 1))
{
- auto id = k - i0 - i1;
- auto Lz0 = lagrangeFactor(i0, z[0], 1);
- auto Lzd = lagrangeFactor(id, z[d], 1);
- out[n][0][0] = (Lz0[1] * Lz1[0] * Lzd[0] - Lz0[0] * Lz1[0] * Lzd[1])*k;
- out[n][0][1] = (Lz0[0] * Lz1[1] * Lzd[0] - Lz0[0] * Lz1[0] * Lzd[1])*k;
- ++n;
+ for (auto i2 : std::array{k - i1 - i0})
+ {
+ out[n][0][0] = (L[0][1][i0] * L[1][0][i1] * L[2][0][i2] - L[0][0][i0] * L[1][0][i1] * L[2][1][i2])*k;
+ out[n][0][1] = (L[0][0][i0] * L[1][1][i1] * L[2][0][i2] - L[0][0][i0] * L[1][0][i1] * L[2][1][i2])*k;
+ ++n;
+ }
}
}
return;
@@ -261,22 +228,21 @@ namespace Dune { namespace Impl
unsigned int n = 0;
for (auto i2 : Dune::range(k + 1))
{
- auto Lz2 = lagrangeFactor(i2, z[2], 1);
for (auto i1 : Dune::range(k - i2 + 1))
{
- auto Lz1 = lagrangeFactor(i1, z[1], 1);
for (auto i0 : Dune::range(k - i2 - i1 + 1))
{
- auto id = k - i0 - i1 -i2;
- auto Lz0 = lagrangeFactor(i0, z[0], 1);
- auto Lzd = lagrangeFactor(id, z[d], 1);
- out[n][0][0] = (Lz0[1] * Lz1[0] * Lz2[0] * Lzd[0] - Lz0[0] * Lz1[0] * Lz2[0] * Lzd[1])*k;
- out[n][0][1] = (Lz0[0] * Lz1[1] * Lz2[0] * Lzd[0] - Lz0[0] * Lz1[0] * Lz2[0] * Lzd[1])*k;
- out[n][0][2] = (Lz0[0] * Lz1[0] * Lz2[1] * Lzd[0] - Lz0[0] * Lz1[0] * Lz2[0] * Lzd[1])*k;
- ++n;
+ for (auto i3 : std::array{k - i2 - i1 - i0})
+ {
+ out[n][0][0] = (L[0][1][i0] * L[1][0][i1] * L[2][0][i2] * L[3][0][i3] - L[0][0][i0] * L[1][0][i1] * L[2][0][i2] * L[3][1][i3])*k;
+ out[n][0][1] = (L[0][0][i0] * L[1][1][i1] * L[2][0][i2] * L[3][0][i3] - L[0][0][i0] * L[1][0][i1] * L[2][0][i2] * L[3][1][i3])*k;
+ out[n][0][2] = (L[0][0][i0] * L[1][0][i1] * L[2][1][i2] * L[3][0][i3] - L[0][0][i0] * L[1][0][i1] * L[2][0][i2] * L[3][1][i3])*k;
+ ++n;
+ }
}
}
}
+
return;
}
@@ -328,58 +294,37 @@ namespace Dune { namespace Impl
// Compute rescaled barycentric coordinates of x
auto z = barycentric(in);
- // To improve code readability we introduce a short cut
- constexpr auto d = dim;
+ // L[j][m][i] is the m-th derivative of the i-th Lagrange polynomial at z[j]
+ auto L = std::array<std::array<std::array<R,k+1>, 2>, dim+1>();
+ for (auto j : Dune::range(dim))
+ evaluateLagrangePolynomialDerivative(z[j], L[j], order[j]);
+ evaluateLagrangePolynomialDerivative(z[dim], L[dim], 1);
if (dim==1)
{
unsigned int n = 0;
for (auto i0 : Dune::range(k + 1))
- {
- auto id = k - i0;
- auto Lz0 = lagrangeFactor(i0, z[0], 1);
- auto Lzd = lagrangeFactor(id, z[d], 1);
- out[n] = (Lz0[1] * Lzd[0] - Lz0[0] * Lzd[1])*k;
- ++n;
- }
+ for (auto i1 : std::array{k - i0})
+ out[n++] = (L[0][1][i0] * L[1][0][i1] - L[0][0][i0] * L[1][1][i1])*k;
return;
}
if (dim==2)
{
unsigned int n=0;
for (auto i1 : Dune::range(k + 1))
- {
- auto Lz1 = lagrangeFactor(i1, z[1], order[1]);
for (auto i0 : Dune::range(k - i1 + 1))
- {
- auto id = k - i0 - i1;
- auto Lz0 = lagrangeFactor(i0, z[0], order[0]);
- auto Lzd = lagrangeFactor(id, z[d], 1);
- out[n] = (Lz0[order[0]] * Lz1[order[1]] * Lzd[0] - Lz0[0] * Lz1[0] * Lzd[1])*k;
- ++n;
- }
- }
+ for (auto i2 : std::array{k - i1 - i0})
+ out[n++] = (L[0][order[0]][i0] * L[1][order[1]][i1] * L[2][0][i2] - L[0][0][i0] * L[1][0][i1] * L[2][1][i2])*k;
return;
}
if (dim==3)
{
unsigned int n = 0;
for (auto i2 : Dune::range(k + 1))
- {
- auto Lz2 = lagrangeFactor(i2, z[2], order[2]);
for (auto i1 : Dune::range(k - i2 + 1))
- {
- auto Lz1 = lagrangeFactor(i1, z[1], order[1]);
for (auto i0 : Dune::range(k - i2 - i1 + 1))
- {
- auto id = k - i0 - i1 -i2;
- auto Lz0 = lagrangeFactor(i0, z[0], order[0]);
- auto Lzd = lagrangeFactor(id, z[d], 1);
- out[n] = (Lz0[order[0]] * Lz1[order[1]] * Lz2[order[2]] * Lzd[0] - Lz0[0] * Lz1[0] * Lz2[0] * Lzd[1])*k;
- ++n;
- }
- }
- }
+ for (auto i3 : std::array{k - i2 - i1 - i0})
+ out[n++] = (L[0][order[0]][i0] * L[1][order[1]][i1] * L[2][order[2]][i2] * L[3][0][i3] - L[0][0][i0] * L[1][0][i1] * L[2][0][i2] * L[3][1][i3])*k;
return;
}
}
@@ -808,6 +753,52 @@ namespace Dune
* \tparam R type used for function values
* \tparam d dimension of the reference element
* \tparam k polynomial order
+ *
+ * The Lagrange basis functions \f$\phi_i\f$ of order \f$k>1\f$ on the unit simplex
+ * \f$G = \{ x \in [0,1]^{d} \,|\, \sum_{j=1}^d x_j \leq 1\}\f$ are implemented as
+ * \f$\phi_i(x) = \hat{\phi}_i(kx)\f$ where \f$\hat{\phi}_i\f$ is the \f$i\f$-th
+ * basis function on the scaled simplex
+ * \f$\hat{G} = kG = \{ x \in [0,k]^{d} \,|\, \sum_{j=1}^d x_i \leq k\}\f$.
+ * On \f$\hat{G}\f$ the uniform Lagrange points of order \f$k\f$ are
+ * exactly the points \f$(i_1,\dots,i_d) \in \mathbb{N}_0^d \cap \hat{G}\f$
+ * with non-negative integer coordinates in \f$\hat{G}\f$.
+ * Using the lexicographic enumeration
+ * of those points we can identify each \f$i=(i_1,...,i_d)\f$ with the
+ * flat index \f$i\f$ of the Lagrange node and associated basis function.
+
+ * Since we can map any point \f$\hat{x} \in \hat{G}\f$ to its barycentric coordinates
+ * \f$(\hat{x}_1, \dots, \hat{x}_d, k-\sum_{j=1}^d \hat{x}_j)\f$ we define for any such point
+ * its auxiliary \f$(d+1)\f$-th coordinate \f$\hat{x}_{d+1} = k-\sum_{j=1}^d \hat{x}_j\f$
+ * and use this in particular for Lagrange points \f$i=(i_1,...,i_d)\f$.
+ * Then the Lagrange basis function on \f$\hat{G}\f$ associated to the
+ * Lagrange point \f$i=(i_1,\dots,i_d)\f$ is given by
+ * \f[
+ * \hat{\phi}_i(\hat{x}) = \prod_{j=1}^{d+1} L_{i_j}(\hat{x}_j)
+ * \f]
+ * where we used the barycentric coordinates of \f$\hat{x}\f$ and \f$i\f$ and the
+ * univariate Lagrange polynomials
+ * \f[
+ * L_n(t) = \prod_{m=0}^{n-1} \frac{t-m}{n-m} = \frac{1}{n!}\prod_{m=0}^{n-1}(t-m).
+ * \f]
+ * This factorization can be interpreted geometrically. To this end note that
+ * any component \f$1,\dots,d+1\f$ of the barycentric coordinates is associated
+ * to a vertex of the simplex and thus to the opposite facet. Furthermore the Lagrange
+ * nodes are all located on hyperplanes parallel to those facets. In the factorized
+ * form the term \f$L_{i_j}(\hat{x}_j)\f$ guarantees that \f$\hat{\phi}_i\f$ is zero
+ * on any of these hyperplane that is parallel to the facet associated with \f$j\f$
+ * and lying between the Lagrange node \f$i\f$ and this facet
+ * (excluding \f$i\f$ and including the facet).
+ *
+ * Using the factorized form, evaluating all basis functions \f$\hat{\phi}_i\f$ at a point \f$\hat{x}\f$
+ * requires to first evaluate \f$\alpha_{m,j}=L_m(\hat{x}_j)\f$ for all \f$j\in \{1,\dots,d+1\}\f$
+ * and \f$m \in \{0,\dots,k\}\f$ and then computing the products
+ * \f$\hat{\phi}_i(\hat{x})=\prod_{j=1}^{d+1} \alpha_{i_j,j}\f$ for all admissible multi indices
+ * (or, equivalently, Lagrange nodes in barycentric coordinates)
+ * \f$(i_1,\dots,i_{d+1})\f$ with \f$\sum_{j=1}^{d+1} i_j = k\f$.
+ * The evaluation of the univariate Lagrange polynomials can be done efficiently using the recursion
+ * \f[
+ * L_0(t) = 1, \qquad L_{n+1}(t) = L_n(t)\frac{t-n}{n+1} \qquad n\geq 0.
+ * \f]
*/
template<class D, class R, int d, int k>
class LagrangeSimplexLocalFiniteElement
--
GitLab
From 6886e465e4dba1ad6370be4f61cfdfd7bd7b4197 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= <graeser@dune-project.org>
Date: Fri, 21 Mar 2025 13:28:45 +0100
Subject: [PATCH 3/4] Implement arbitrary order partial derivatives in
LagrangeSimplexLFE
We can implement arbitrary order derivatives by a simple recursion.
Notice that the effort grows exponentially with the total order because
the number of terms from the product rule is 2^totalOrder. Furthermore
the stack-storage grows linearly with the dynamic total order. Hence
we need to use a dynamic-to-static switch using an explicit enumeration
of the supported static derivative orders. Since derivatives of order >k
are zero, the number of these cases is known at compile time.
---
.../lagrange/lagrangesimplex.hh | 201 ++++++++----------
1 file changed, 93 insertions(+), 108 deletions(-)
diff --git a/dune/localfunctions/lagrange/lagrangesimplex.hh b/dune/localfunctions/lagrange/lagrangesimplex.hh
index d3d96811..42b91ac7 100644
--- a/dune/localfunctions/lagrange/lagrangesimplex.hh
+++ b/dune/localfunctions/lagrange/lagrangesimplex.hh
@@ -7,10 +7,14 @@
#include <array>
#include <numeric>
+#include <algorithm>
+#include <dune/common/exceptions.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
+#include <dune/common/hybridutilities.hh>
#include <dune/common/math.hh>
+#include <dune/common/rangeutilities.hh>
#include <dune/geometry/referenceelements.hh>
@@ -65,22 +69,76 @@ namespace Dune { namespace Impl
}
// Evaluate the univariate Lagrange polynomial derivatives L_i(t) for i=0,...,k
- // up to given maxDerivativeOrder. Currently only maxDerivativeOrder 0 and 1 are supported.
+ // up to given maxDerivativeOrder.
static constexpr void evaluateLagrangePolynomialDerivative(const R& t, auto& LL, unsigned int maxDerivativeOrder)
{
auto& L = LL[0];
L[0] = 1;
for (auto i : Dune::range(k))
L[i+1] = L[i] * (t - i) / (i+1);
- auto& DL = LL[1];
- if (maxDerivativeOrder>0)
+ for(auto j : Dune::range(maxDerivativeOrder))
{
- DL[0] = 0;
+ auto& F = LL[j];
+ auto& DF = LL[j+1];
+ DF[0] = 0;
for (auto i : Dune::range(k))
- DL[i+1] = (DL[i] * (t - i) + L[i]) / (i+1);
+ DF[i+1] = (DF[i] * (t - i) + (j+1)*F[i]) / (i+1);
}
}
+ using BarycentricMultiIndex = std::array<unsigned int,dim+1>;
+
+
+ // This computed the required partial derivative given by the multi-index
+ // beta of a product of a function given as a product of dim+1 derivatives
+ // of univariate Lagrange polynomials of the dim+1 barycentric coordinates.
+ // The polynomials in the product are specified as follows:
+ //
+ // The table L contains all required derivatives of all univariate
+ // polynomials evaluated at all barycentric coordinates. The two
+ // multi-indices i and alpha encode that the polynomial for the
+ // j-th barycentric coordinate is the alpha-j-th derivative of
+ // the i_j-the Lagrange polynomial.
+ //
+ // Hence this method computes D_beta f(x) where f(x) is the product
+ // \f$f(x) = \prod_{j=0}^{d} L_{i_j}^{(alpha_j)}(x_j) \f$.
+ static constexpr R barycentricDerivative(
+ BarycentricMultiIndex beta,
+ const auto&L,
+ const BarycentricMultiIndex& i,
+ const BarycentricMultiIndex& alpha = {})
+ {
+ // If there are unprocessed derivatives left we search the first unprocessed
+ // partial derivative direction j and compute it using the product and chain rule.
+ // The remaining partial derivatives are forwarded to the recursion.
+ // Notice that the partial derivative of the last barycentric component
+ // wrt to the j-th component is -1 which is responsible for the second
+ // term in the sum. Furthermore we get the factor k due to the scaling of
+ // the simplex.
+ for(auto j : Dune::range(dim))
+ {
+ if (beta[j] > 0)
+ {
+ auto leftDerivatives = alpha;
+ auto rightDerivatives = alpha;
+ leftDerivatives[j]++;
+ rightDerivatives.back()++;
+ beta[j]--;
+ return (barycentricDerivative(beta, L, i, leftDerivatives) - barycentricDerivative(beta, L, i, rightDerivatives))*k;
+ }
+ }
+ // If there are no unprocessed derivatives we can simply evaluate
+ // the product of the derivatives of the Lagrange polynomials with
+ // given indices i_j and derivative orders alpha_j
+ // Evaluate the product of the univariate Lagrange polynomials
+ // with given indices and orders.
+ auto y = R(1);
+ for(auto j : Dune::range(dim+1))
+ y *= L[j][alpha[j]][i[j]];
+ return y;
+ }
+
+
public:
using Traits = LocalBasisTraits<D,dim,FieldVector<D,dim>,R,1,FieldVector<R,1>,FieldMatrix<R,1,dim> >;
@@ -259,157 +317,84 @@ namespace Dune { namespace Impl
const typename Traits::DomainType& in,
std::vector<typename Traits::RangeType>& out) const
{
- auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
+ auto totalOrder = std::accumulate(order.begin(), order.end(), 0u);
out.resize(size());
- if (totalOrder == 0) {
+ // Derivative order zero corresponds to the function evaluation.
+ if (totalOrder == 0)
+ {
evaluateFunction(in, out);
return;
}
- if (k==0)
+ // Derivatives of order >k are all zero.
+ if (totalOrder > k)
{
- out[0] = 0;
+ for(auto& out_i : out)
+ out_i = 0;
return;
}
+ // It remains to cover the cases 0 < totalOrder<= k.
+
if (k==1)
{
if (totalOrder==1)
{
auto direction = std::find(order.begin(), order.end(), 1);
-
out[0] = -1;
for (unsigned int i=0; i<dim; i++)
out[i+1] = (i==(direction-order.begin()));
}
- else // all higher order derivatives are zero
- std::fill(out.begin(), out.end(), 0);
return;
}
- if (totalOrder==1)
- {
+ // Since the required stack storage depends on the dynamic total order,
+ // we need to do a dynamic to static dispatch by enumerating all supported
+ // static orders.
+ auto supportedStaticOrders = Dune::range(Dune::index_constant<1>{}, Dune::index_constant<k+1>{});
+ return Dune::Hybrid::switchCases(supportedStaticOrders, totalOrder, [&](auto staticTotalOrder) {
+
// Compute rescaled barycentric coordinates of x
auto z = barycentric(in);
// L[j][m][i] is the m-th derivative of the i-th Lagrange polynomial at z[j]
- auto L = std::array<std::array<std::array<R,k+1>, 2>, dim+1>();
+ auto L = std::array<std::array<std::array<R, k+1>, staticTotalOrder+1>, dim+1>();
for (auto j : Dune::range(dim))
evaluateLagrangePolynomialDerivative(z[j], L[j], order[j]);
- evaluateLagrangePolynomialDerivative(z[dim], L[dim], 1);
+ evaluateLagrangePolynomialDerivative(z[dim], L[dim], totalOrder);
- if (dim==1)
+ auto barycentricOrder = BarycentricMultiIndex{};
+ for (auto j : Dune::range(dim))
+ barycentricOrder[j] = order[j];
+ barycentricOrder[dim] = 0;
+
+ if constexpr (dim==1)
{
unsigned int n = 0;
for (auto i0 : Dune::range(k + 1))
for (auto i1 : std::array{k - i0})
- out[n++] = (L[0][1][i0] * L[1][0][i1] - L[0][0][i0] * L[1][1][i1])*k;
- return;
+ out[n++] = barycentricDerivative(barycentricOrder, L, BarycentricMultiIndex{i0, i1});
}
- if (dim==2)
+ if constexpr (dim==2)
{
unsigned int n=0;
for (auto i1 : Dune::range(k + 1))
for (auto i0 : Dune::range(k - i1 + 1))
for (auto i2 : std::array{k - i1 - i0})
- out[n++] = (L[0][order[0]][i0] * L[1][order[1]][i1] * L[2][0][i2] - L[0][0][i0] * L[1][0][i1] * L[2][1][i2])*k;
- return;
+ out[n++] = barycentricDerivative(barycentricOrder, L, BarycentricMultiIndex{i0, i1, i2});
}
- if (dim==3)
+ if constexpr (dim==3)
{
unsigned int n = 0;
for (auto i2 : Dune::range(k + 1))
for (auto i1 : Dune::range(k - i2 + 1))
for (auto i0 : Dune::range(k - i2 - i1 + 1))
for (auto i3 : std::array{k - i2 - i1 - i0})
- out[n++] = (L[0][order[0]][i0] * L[1][order[1]][i1] * L[2][order[2]][i2] * L[3][0][i3] - L[0][0][i0] * L[1][0][i1] * L[2][0][i2] * L[3][1][i3])*k;
- return;
+ out[n++] = barycentricDerivative(barycentricOrder, L, BarycentricMultiIndex{i0, i1, i2, i3});
}
- }
-
- if (dim==2)
- {
- auto lagrangeNode = [](unsigned int i) { return ((D)i)/k; };
-
- // Helper method: Return a single Lagrangian factor of l_ij evaluated at x
- auto lagrangianFactor = [&lagrangeNode]
- (const int no, const int i, const int j, const typename Traits::DomainType& x)
- -> typename Traits::RangeType
- {
- if ( no < i)
- return (x[0]-lagrangeNode(no))/(lagrangeNode(i)-lagrangeNode(no));
- if (no < i+j)
- return (x[1]-lagrangeNode(no-i))/(lagrangeNode(j)-lagrangeNode(no-i));
- return (lagrangeNode(no+1)-x[0]-x[1])/(lagrangeNode(no+1)-lagrangeNode(i)-lagrangeNode(j));
- };
-
- // Helper method: Return the derivative of a single Lagrangian factor of l_ij evaluated at x
- // direction: Derive in x-direction if this is 0, otherwise derive in y direction
- auto lagrangianFactorDerivative = [&lagrangeNode]
- (const int direction, const int no, const int i, const int j, const typename Traits::DomainType&)
- -> typename Traits::RangeType
- {
- using T = typename Traits::RangeType;
- if ( no < i)
- return (direction == 0) ? T(1.0/(lagrangeNode(i)-lagrangeNode(no))) : T(0);
-
- if (no < i+j)
- return (direction == 0) ? T(0) : T(1.0/(lagrangeNode(j)-lagrangeNode(no-i)));
-
- return -1.0/(lagrangeNode(no+1)-lagrangeNode(i)-lagrangeNode(j));
- };
-
- if (totalOrder==2)
- {
- std::array<int,2> directions;
- unsigned int counter = 0;
- auto nonconstOrder = order; // need a copy that I can modify
- for (int i=0; i<2; i++)
- {
- while (nonconstOrder[i])
- {
- directions[counter++] = i;
- nonconstOrder[i]--;
- }
- }
-
- //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
- int n=0;
- for (unsigned int j=0; j<=k; j++)
- {
- for (unsigned int i=0; i<=k-j; i++, n++)
- {
- R res = 0.0;
-
- for (unsigned int no1=0; no1 < k; no1++)
- {
- R factor1 = lagrangianFactorDerivative(directions[0], no1, i, j, in);
- for (unsigned int no2=0; no2 < k; no2++)
- {
- if (no1 == no2)
- continue;
- R factor2 = factor1*lagrangianFactorDerivative(directions[1], no2, i, j, in);
- for (unsigned int no3=0; no3 < k; no3++)
- {
- if (no3 == no1 || no3 == no2)
- continue;
- factor2 *= lagrangianFactor(no3, i, j, in);
- }
- res += factor2;
- }
- }
- out[n] = res;
- }
- }
-
- return;
- } // totalOrder==2
-
- } // dim==2
-
- DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
+ });
}
//! \brief Polynomial order of the shape functions
--
GitLab
From d6a607f5456780f2e9179739e55c24c01eecd2d3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= <graeser@dune-project.org>
Date: Thu, 20 Mar 2025 11:19:09 +0100
Subject: [PATCH 4/4] [test] Test 2st order partial() for
LagrangeSimplexFiniteElement
The implementation supports arbitrary order, but the test is only
implemented for order up to 2.
---
dune/localfunctions/test/lagrangeshapefunctiontest.cc | 4 +---
1 file changed, 1 insertion(+), 3 deletions(-)
diff --git a/dune/localfunctions/test/lagrangeshapefunctiontest.cc b/dune/localfunctions/test/lagrangeshapefunctiontest.cc
index c5310c97..555c7c9d 100644
--- a/dune/localfunctions/test/lagrangeshapefunctiontest.cc
+++ b/dune/localfunctions/test/lagrangeshapefunctiontest.cc
@@ -178,9 +178,7 @@ DUNE_NO_DEPRECATED_END
{
Dune::Hybrid::forEach(std::make_index_sequence<5>{},[&](auto order)
{
- // In general we implement partial derivatives up to order 1
- // For some special cases, we also have 2nd order partial derivatives
- auto diffOrder = (dim==2 or order<2) ? 2 : 0;
+ auto diffOrder = 2;
auto lfe = LagrangeSimplexLocalFiniteElement<double,double,dim,order>();
testSuite.subTest(testVirtualLFE(lfe, DisableNone, diffOrder));
});
--
GitLab