Skip to content

Local functions

Local functions are functions that are bound to single grid elements. Therefore, they are constructed from some local basis, a coefficient vector, and the geometry of the grid element. Since the implementation is quite involved, localfefunctions do not reside at Ikarus but in the separate module dune-localfefunctions.

Usually, local functions are need to be evaluated in the local coordinate system ξTrefRn :

f:ξnRm

where Tref is the reference element, e.g., for a hypercube, Tref=[0,1]d.

Interface

Local functions provide the following interface:

LocalFunction(const Dune::CachedLocalBasis<DuneBasis>& p_basis, const CoeffContainer& coeffs_,
              const std::shared_ptr<const Geometry>& geo,
              Dune::template index_constant<ID> = Dune::template index_constant<std::size_t(0)>{}); 

FunctionReturnType evaluate(const DomainType& local);
FunctionReturnType evaluate(const unsigned int& integrationPointIndex);

auto evaluateDerivative(const DomainType& local,...);
auto evaluateDerivative(const unsigned int& integrationPointIndex,...);
auto viewOverIntegrationPoints(); 

template<std::size_t ID=0>
constexpr int order(Dune::index_constant<ID> ); 

template<std::size_t ID=0>
auto basis(Dune::index_constant<ID> ); 

template<std::size_t ID=0>
auto coefficientsRef(Dune::index_constant<ID>); 

template <typename IntegrationRule, typename... Ints>
void bind(IntegrationRule&& p_rule, Derivatives<Ints...>&& ints); 

auto clone (); 

template<typename ScalarType, std::size_t ID=0>
auto rebindClone (ScalarType, Dune::index_constant<ID>); 

The "..." in the evaluateDerivative function call refers to several possible variadic templates. The implementation looks like the following:

using namespace Dune::DerivativeDirections;
localFunction.bind(rule, bindDerivatives(0,1));
for(auto& [gpIndex, gp] : localFunction.viewOverIntegrationPoints()){
  localFunction.evaluateDerivative(gpIndex, wrt(spatialAll)); 
  localFunction.evaluateDerivative(gpIndex, wrt(spatialAll), on(gridElement)); 
}
using namespace Dune::DerivativeDirections;
for(auto& gp : rule){
  localFunction.evaluateDerivative(gp.position(), wrt(spatialAll)); 
  localFunction.evaluateDerivative(gp.position(), wrt(spatialAll), on(gridElement)); 
}

where the first call implements

gradξf:ξRm×d.

The second one takes into account the fact that the local function is defined in some physical space X with the coordinate x. Therefore, it transforms the Jacobian from the reference element gradξ to the Jacobian on the grid element gradx. This behavior is activated, if on(gridElement) is passed; otherwise, if the derivatives on the reference element is needed, pass on (referenceElement). Here, gridElement and referenceElement are global constants in the namespace Dune::DerivativeDirections.

Thus, if on(gridElement) is passed, the local function usually implements

gradx=gradξJ1

where J is the Jacobian of the mapping from the reference element Tref to the element living in physical space T. For details, see page 22 of the Dune book1.

Instead of passing spatialAll to wrt(..), there are other helper functions such as:

localFunction.evaluateDerivative(gpIndex, wrt(spatial(0))); 
localFunction.evaluateDerivative(gpIndex, wrt(spatial(1))); 

which can also be combined with on(...).

Derivatives w.r.t. coefficients

localFunction.evaluateDerivative(gpIndex, wrt(coeff(j)));

which evaluates the first derivative for a vector space valued function,e.g., for f(ξ)=I=1nNI(ξ)xI, we arrive at a matrix A such that

[A]ij=Aij=fi(ξ)xj

and the second derivative

localFunction.evaluateDerivative(gpIndex, wrt(coeff(j,k)), along(q));
[B]jk=Bjk=qiAijk=2(qifi(ξ))xjxk

where q is an arbitrary vector of the same size as f, i.e., it is the direction of the derivative in this case. A and B are simply then the returned matrices and do not have any special meaning. If a vector is not passed while evaluating the second derivative, the result would be a third-order tensor for a vector-valued function f. As a result, a direction derivative in the direction given by along(q) is computed to return a matrix B in this case. This helps for both readability and performance. See the example later for more details.

Derivatives w.r.t. coefficients and spatial derivatives

Spatial derivatives and derivatives w.r.t. the coefficients can be combined. Therefore, it is legal to call

auto B = localFunction.evaluateDerivative(gpIndex, wrt(coeff(j,k),spatialAll), along(Q));
auto b1 = localFunction.evaluateDerivative(gpIndex, wrt(coeff(j,k),spatial(0)), along(q));
auto b2 = localFunction.evaluateDerivative(gpIndex, wrt(coeff(j,k),spatial(1)), along(q));

Warning

The order of spatial and coefficient derivatives does not matter. The returned value is always rearranged so that the first derivative is the spatial one.

The first line is then equivalent to

[B]jk=Bjk=QilAiljk=2([gradξf(ξ)]ilQil)xjxk.

For the second and third line, we have

b0,jk=2([gradξ0f(ξ)]iqi)xjxk,b1,jk=2([gradξ1f(ξ)]iqi)xjxk.

These objects are also returned when the second and third lines above are used.

All of these function calls, once again, can be combined with on(gridElement) as shown below:

localFunction.evaluateDerivative(gpIndex, wrt(coeff(j,k),spatialAll), along(Q), on(gridElement));

which computes

2([gradxf(ξ)]ilQil)xjxk.

Warning

Currently, only first-order spatial derivatives and second-order derivatives w.r.t. the coefficients are supported.

Example: Dirichlet energy

This example shows how the energy, gradient, and Hessian of a dirichlet energy can be calculated. E(u)=12Ω||gradxu(x)||2dx

If the energy is to be minimized w.r.t. the coefficients of the nodes, the energy, gradient, and Hessian w.r.t. the coefficients are to be calculated. Of course, this depends on the optimization algorithms, but for now, the general case where all three are needed is considered.

auto dirichletEnergy() {
  double energy = 0;
  // bind localBasis to some integration rule
  // create uNodalCoeffs
  Ikarus::StandardLocalFunction uFunction(localBasis, uNodalCoeffs, sharedGeometry);
  for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
    //.. calculate the inverse Jacobian of the geometry
    const auto gradu = uFunction.evaluateDerivative(gpIndex, wrt(spatialAll), on(gridElement));
    energy+= 0.5 * (gradu.transpose() * gradu).trace() * gp.weight() * sharedGeometry->integrationElement(gp.position());
  }
}
auto gradientDirichletEnergy(Eigen::VectorXd& g) {
  //.bind localBasis to some integration rule
  // create uNodalCoeffs
  constexpr int size =  // spatial size of u
      Ikarus::StandardLocalFunction uFunction(localBasis, uNodalCoeffs, sharedGeometry);
  for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
    //.. calculate the inverse Jacobian of the geometry
    const auto gradu = uFunction.evaluateDerivative(gpIndex, wrt(spatialAll), on(gridElement));
    for (auto i : fe.size()) { //loop over coeffs, i.e.nodes of the finite element
      const auto graduDCoeffs
          = uFunction.evaluateDerivative(gpIndex, wrt(spatialAll, coeff(i)), on(gridElement));
      Eigen::Vector<double, size> tmp;
      tmp.setZero();
      for (int k = 0; k < gridDimension; ++k)
        tmp += graduDCoeffs[k] * gradu.col(k);  
      g.segment<size>(i * size) += tmp * gp.weight() * sharedGeometry->integrationElement(gp.position());
    }
  }
}
auto hessianDirichletEnergy(Matrix& h) {
  //... bind localBasis to some integration rule
  // and create uNodalCoeffs
  constexpr int size =  // spatial size of u
      Ikarus::StandardLocalFunction uFunction(localBasis, uNodalCoeffs, sharedGeometry);
  for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
    //.. calculate the inverse Jacobian of the geometry
    for (auto i : loop over coeffs, i.e.nodes of the finite element) {
      const auto graduDCoeffsI
          = uFunction.evaluateDerivative(gpIndex, wrt(spatialAll, coeff(i)), on(gridElement));
      for (auto j : fe.size()) { //loop over coeffs, i.e.nodes of the finite element
        const auto graduDCoeffsJ
          = uFunction.evaluateDerivative(gpIndex, wrt(spatialAll, coeffs), on(gridElement), coeffIndices(j));
        Eigen::Matrix<double, size, size> tmp;
        tmp.setZero();
        for (int k = 0; k < gridDimension; ++k)
          tmp += graduDCoeffsI[k] * graduDCoeffsJ[k];
        h.block<size, size>(i * size, j * size) += tmp * gp.weight() * sharedGeometry->integrationElement(gp.position());
      }
    }
  }
}

Implementations

In the following, the local functions that are currently available are summarized. The ansatz functions are denoted as Ni(ξ) in the table below.

Name Interpolation formula Note Header
Standard x=i=1nNi(ξ)xi - standardLocalFunction.hh
Projection-Based2 x=P(i=1nNi(ξ)xi) This is one version of geometric finite elements. These are finite elements suited
for interpolation on manifolds. Here P:RmM is an operator that projects
the usual linear interpolation onto a manifold.
projectionBasedLocalFunction.hh

How to implement your own local functions

To implement your own local function, the file ikarus/localFunctions/impl/localFunctionTemplate.hh is made available.

The file can be copied, and then you can rename the class to your preferred name and implement the functions mentioned below. If a particular function is not required, it has to be deleted explicitly. Then, if someone calls that function, it returns a zero matrix or vector of appropriate size.

These functions are all templated with DomainTypeOrIntegrationPointIndex, which is an integration point index or position. Additionally, On<TransformArgs> specifies whether the function should be evaluated on the reference element or the grid element (see above).

FunctionReturnType evaluateFunctionImpl(const DomainTypeOrIntegrationPointIndex& ipIndexOrPosition,
                                        const On<TransformArgs>&) const; 

Jacobian evaluateDerivativeWRTSpaceAllImpl(const DomainTypeOrIntegrationPointIndex& ipIndexOrPosition,
                                           const On<TransformArgs>& transArgs) const; 

JacobianColType evaluateDerivativeWRTSpaceSingleImpl(const DomainTypeOrIntegrationPointIndex& ipIndexOrPosition,
                                                     int spaceIndex, const On<TransformArgs>& transArgs) const; 


CoeffDerivMatrix evaluateDerivativeWRTCoeffsImpl(const DomainTypeOrIntegrationPointIndex& ipIndexOrPosition,
                                                 int coeffsIndex, const On<TransformArgs>& transArgs) const; 

CoeffDerivMatrix evaluateSecondDerivativeWRTCoeffsImpl(const DomainTypeOrIntegrationPointIndex& ipIndexOrPosition,
                                                       const std::array<size_t, 2>& coeffsIndex,
                                                       const Along<AlongArgs...>& alongArgs,
                                                       const On<TransformArgs>& transArgs) const; 

std::array<CoeffDerivEukRieMatrix, gridDim> evaluateDerivativeWRTCoeffsANDSpatialImpl(
    const DomainTypeOrIntegrationPointIndex& ipIndexOrPosition, int coeffsIndex,
    const On<TransformArgs>& transArgs) const; 


CoeffDerivEukRieMatrix evaluateDerivativeWRTCoeffsANDSpatialSingleImpl(
    const DomainTypeOrIntegrationPointIndex& ipIndexOrPosition, int coeffsIndex, int spatialIndex,
    const On<TransformArgs>& transArgs) const;  


auto evaluateThirdDerivativeWRTCoeffsTwoTimesAndSpatialImpl(
    const DomainTypeOrIntegrationPointIndex& ipIndexOrPosition, const std::array<size_t, 2>& coeffsIndex,
    const Along<AlongArgs...>& alongArgs, const On<TransformArgs>& transArgs) const; 

CoeffDerivMatrix evaluateThirdDerivativeWRTCoeffsTwoTimesAndSpatialSingleImpl(
    const DomainTypeOrIntegrationPointIndex& ipIndexOrPosition, const std::array<size_t, 2>& coeffsIndex,
    const int spatialIndex, const Along<AlongArgs...>& alongArgs, const On<TransformArgs>& transArgs) const; 

Expressions

Expression templates are usually used in linear algebra libraries, e.g., Eigen or Blaze. The syntax is similar to the one provided by UML but only acts on local functions. Expression templates are used here to combine existing local functions in order to obtain new nested ones.

For example, consider the following code:

...
auto f = Ikarus::StandardLocalFunction(localBasis, coeffVectors0, sharedGeometry);
auto g = Ikarus::StandardLocalFunction(localBasis, coeffVectors1, sharedGeometry);

Two local functions that satisfy the interface described above are created. Now it is possible to combine these functions and get an object that also satisfies the concept above. Thus, the following is possible:

...
auto f = Ikarus::StandardLocalFunction(localBasis, coeffVectors0, sharedGeometry);
auto g = Ikarus::StandardLocalFunction(localBasis, coeffVectors1, sharedGeometry);
auto k = f + g;
k.evaluateDerivative(ipIndex, wrt(coeff(i), spatial(d)));

Currently, binary and unary expressions are supported. The following expressions are defined:

Name Mathematical formula Code Note
Sum f+g f+g f and g needs to be of the same size.
DotProduct fg=figi dot(f,g) f and g needs to be of the same size.
normSquared ff=fifi normSquared(f)
Negate f -f
sqrt f sqrt(f) The function f needs a scalar return type.
log logf log(f) The function f needs a scalar return type. Here, log is the natural logarithm.
pow fn pow<n>(f) The function f needs a scalar return type. n is an integer given during compile-time.
Scale af,aR a*f and f/a a has to satisfy std::is_arithmetic<..>
LinearStrains 12(H+HT),H=grad(f) linearStrains(f) If you call the formula on the left with on(gridElement), it will assume the transformed derivatives.
GreenLagrangianStrains 12(H+HT+HTH),H=grad(f) greenLagrangianStrains(f) If you call the formula on the left with on(gridElement), it will assume the transformed derivatives.

These expressions can also be nested. As a result, it is valid to write

auto f = Ikarus::StandardLocalFunction(localBasis, coeffVectors0);
auto g = Ikarus::StandardLocalFunction(localBasis, coeffVectors1);
auto k = -sqrt(dot(2*f+f,5*g));
k.evaluateDerivative(ipIndex, wrt(coeff(i), spatial(d)));

To use these expressions, there are additional exported static types for all expressions.

constexpr bool isLeaf; 
constexpr bool children; 

Note

To use these expression, simply include the header #include <dune/localfefunctions/expressions.hh>.

Tagging leaf local functions

In the context of mixed finite elements, there are usually several local functions that contribute to the energy. These stem from different local bases. For example, consider the Q1P0 element, where displacements are interpolated by using the four bilinear ansatz functions and the element-wise constant pressure field.

Then, to obtain gradients and Hessians, we need to differentiate w.r.t. different coefficients. This can be done by tagging the local function during construction.

using namespace Dune::Indices;
auto f = Ikarus::StandardLocalFunction(localBasis0, coeffVectors0, sharedGeometry, _0);
auto g = Ikarus::StandardLocalFunction(localBasis1, coeffVectors1, sharedGeometry, _1);
auto k = dot(f,g);
k.evaluateDerivative(ipIndex, wrt(coeff(_0,i,_1,j))); // Second derivative w.r.t. the nodal coefficients

To explain the last line above, let us consider that the function f is constructed as f=I=0nNLfL and similarly g=I=0mMKgK, where N and M are ansatz functions and fL and gK are nodal coefficients.

Thus, the above call translates to

2(fg)figj,

where the correct sizes of the result are derived at compile-time. The complete Hessian of fg can be calculated by the following:

using namespace Dune::Indices;
auto hessianDirichletEnergy(Matrix& h) {
  //... bind localBasis to some integration rule
  using namespace Dune::Indices;
  auto f = Ikarus::StandardLocalFunction(localBasis0, coeffVectors0, sharedGeometry, _0);
  auto g = Ikarus::StandardLocalFunction(localBasis1, coeffVectors1, sharedGeometry, _1);
  auto k = dot(f,g);
  constexpr int sizef = f.correctionSize; // spatial size of the correction of the coefficients of f
  constexpr int sizeg = g.correctionSize; // spatial size of the correction of the coefficients of g
  constexpr int coeffSizef = coeffVectors0.size();
  constexpr int coeffSizeg = coeffVectors1.size();

  Dune::MultiTypeBlockMatrix<Dune::MultiTypeBlockVector<MatrixBlock00,MatrixBlock01>,
                                       Dune::MultiTypeBlockVector<MatrixBlock10,MatrixBlock11> > KBlocked; 


  for (const auto& [gpIndex, gp] : k.viewOverIntegrationPoints()) {
    for (size_t I = 0; I < coeffSizef; ++I)
      for (size_t J = 0; J < coeffSizef; ++J)
        KBlocked[_0,_0].block<sizef, sizef>(I * sizef, J * sizef)
          += k.evaluateDerivative(gpIndex, wrt(coeff(_0,I,_0,J))) * gp.weight() * sharedGeometry->integrationElement(gp.position());

    for (size_t I = 0; I < coeffSizef; ++I)
      for (size_t J = 0; J < coeffSizeg; ++J)
        KBlocked[_0,_1].block<sizef, sizeg>(I * sizef, J * sizeg)
          += k.evaluateDerivative(gpIndex, wrt(coeff(_0,I,_1,J))) * gp.weight() * sharedGeometry->integrationElement(gp.position());

    for (size_t I = 0; I < coeffSizeg; ++I)
      for (size_t J = 0; J < coeffSizeg; ++J)
        KBlocked[_1,_1].block<sizeg, sizeg>(I * sizeg, J * sizeg)
          += k.evaluateDerivative(gpIndex, wrt(coeff(_1,I,_1,J))) * gp.weight() * sharedGeometry->integrationElement(gp.position());

    for (size_t I = 0; I < coeffSizeg; ++I)
      for (size_t J = 0; J < coeffSizef; ++J)
        KBlocked[_1,_0].block<sizef, sizeg>(I * sizeg, J * sizef)
          += k.evaluateDerivative(gpIndex, wrt(coeff(_1,I,_0,J))) * gp.weight() * sharedGeometry->integrationElement(gp.position());
    }
}

Writing your own expressions

It is also possible to write your own expressions. To do so, please take a look at the existing expressions. The sqrt and normSquared expressions are the most general unary and binary expressions implemented.

Implementing the return value

Implementing a return value is the first step in implementing an expression. This is done by using the following function:

template <typename LFArgs>
auto evaluateValueOfExpression(const LFArgs &lfArgs) const;

Warning

The interface dictates that the return value needs to be an Eigen type. Thus, even if a scalar double is to be returned, it is to

be wrapped in Eigen::Vector<double, 1>

Additionally, the evaluation of the derivative is to be implemented as shown below:

template <int DerivativeOrder, typename LFArgs>
auto evaluateDerivativeOfExpression(const LFArgs &lfArgs) const;

Evaluate the underlying functions

Expressions always act on existing expressions. Therefore, to have the correct return value for the expression, the underlying quantities are to be evaluated. this->m() is used to access unary functions, and this->l() and this->r() are used to access binary expressions.

To evaluate unary functions, the following syntax is used:

const auto mEvaluated = evaluateFunctionImpl(this->m(), lfArgs);

and for binary functions,

const auto l_Evaluated = evaluateFunctionImpl(this->l(), lfArgs);
const auto r_Evaluated = evaluateFunctionImpl(this->r(), lfArgs);

Because the expression conforms to the syntax of a local function, its derivative can also be evaluated.

In the function evaluateDerivativeOfExpression, the template argument DerivativeOrder contains the derivative order. Additionally, the derivative types can also be accessed using the static booleans, as shown below:

    static constexpr bool hasTwoCoeff;
    static constexpr bool hasSingleCoeff;
    static constexpr bool hasNoCoeff;
    static constexpr bool hasNoSpatial;
    static constexpr bool hasOneSpatialAll;
    static constexpr bool hasOneSpatialSingle;
    static constexpr bool hasOneSpatial;

Using the dot-product as a binary expression example, we have

template <int DerivativeOrder, typename LFArgs>
    auto evaluateDerivativeOfExpression(const LFArgs &lfArgs) const {
      const auto u = evaluateFunctionImpl(this->l(), lfArgs);
      const auto v = evaluateFunctionImpl(this->r(), lfArgs);
      if constexpr (DerivativeOrder == 1)  
      {
        const auto u_x = evaluateDerivativeImpl(this->l(), lfArgs); 
        const auto v_x = evaluateDerivativeImpl(this->r(), lfArgs);
        return Ikarus::eval(v.transpose() * u_x + u.transpose() * v_x); 
      } else if constexpr (DerivativeOrder == 2) {   
        const auto &[u_x, u_y] = evaluateFirstOrderDerivativesImpl(this->l(), lfArgs); 
        const auto &[v_x, v_y] = evaluateFirstOrderDerivativesImpl(this->r(), lfArgs);
        if constexpr (LFArgs::hasNoSpatial and LFArgs::hasTwoCoeff) { (6)
          const auto alonguArgs = replaceAlong(lfArgs, along(v)); (7)
          const auto alongvArgs = replaceAlong(lfArgs, along(u));

          const auto u_xyAlongv = evaluateDerivativeImpl(this->l(), alongvArgs); (8)
          const auto v_xyAlongu = evaluateDerivativeImpl(this->r(), alonguArgs);

          return Ikarus::eval(u_xyAlongv + transpose(u_x) * v_y + transpose(v_x) * u_y + v_xyAlongu);
        } else if constexpr (LFArgs::hasOneSpatial and LFArgs::hasSingleCoeff) { (9)
          const auto u_xy = evaluateDerivativeImpl(this->l(), lfArgs); (10)
          const auto v_xy = evaluateDerivativeImpl(this->r(), lfArgs);
          if constexpr (LFArgs::hasOneSpatialSingle and LFArgs::hasSingleCoeff) { (11)
            return Ikarus::eval(transpose(v) * u_xy + transpose(u_x) * v_y + transpose(v_x) * u_y
                                + transpose(u) * v_xy);
          } else if constexpr (LFArgs::hasOneSpatialAll and LFArgs::hasSingleCoeff) { (12)
            std::array<std::remove_cvref_t<decltype(Ikarus::eval(transpose(v) * u_xy[0]))>, gridDim> res; (13)
            for (int i = 0; i < gridDim; ++i)
              res[i] = Ikarus::eval(transpose(v) * u_xy[i] + transpose(u_x.col(i)) * v_y + transpose(v_x.col(i)) * u_y
                                    + transpose(u) * v_xy[i]);
            return res;
          }
        }
      } else if constexpr (DerivativeOrder == 3) { (14)
        if constexpr (LFArgs::hasOneSpatialSingle) {  (15)
          const auto argsForDyz = lfArgs.extractSecondWrtArgOrFirstNonSpatial(); (16)

          const auto &[u_x, u_y, u_z] = evaluateFirstOrderDerivativesImpl(this->l(), lfArgs); (17)
          const auto &[v_x, v_y, v_z] = evaluateFirstOrderDerivativesImpl(this->r(), lfArgs);
          const auto &[u_xy, u_xz]    = evaluateSecondOrderDerivativesImpl(this->l(), lfArgs); (18)
          const auto &[v_xy, v_xz]    = evaluateSecondOrderDerivativesImpl(this->r(), lfArgs);

          const auto alonguArgs             = replaceAlong(lfArgs, along(u));
          const auto alongvArgs             = replaceAlong(lfArgs, along(v));
          const auto argsForDyzalongv_xArgs = replaceAlong(argsForDyz, along(v_x)); (19)
          const auto argsForDyzalongu_xArgs = replaceAlong(argsForDyz, along(u_x));

          const auto u_xyzAlongv = evaluateDerivativeImpl(this->l(), alongvArgs); (20)
          const auto v_xyzAlongu = evaluateDerivativeImpl(this->r(), alonguArgs);
          const auto u_yzAlongvx = evaluateDerivativeImpl(this->l(), argsForDyzalongv_xArgs); (21)
          const auto v_yzAlongux = evaluateDerivativeImpl(this->r(), argsForDyzalongu_xArgs);

          return Ikarus::eval(u_xyzAlongv + transpose(u_xy) * v_z + transpose(u_xz) * v_y + v_yzAlongux + u_yzAlongvx
                              + transpose(v_xz) * u_y + transpose(v_xy) * u_z + v_xyzAlongu);
        } else if constexpr (LFArgs::hasOneSpatialAll) { (22)
          const auto &alongMatrix = std::get<0>(lfArgs.alongArgs.args); (23)

          const auto uTimesA = eval(u * alongMatrix);
          const auto vTimesA = eval(v * alongMatrix);

          const auto &[gradu, u_c0, u_c1]  = evaluateFirstOrderDerivativesImpl(this->l(), lfArgs); (24)
          const auto &[gradv, v_c0, v_c1]  = evaluateFirstOrderDerivativesImpl(this->r(), lfArgs);
          const auto &[gradu_c0, gradu_c1] = evaluateSecondOrderDerivativesImpl(this->l(), lfArgs); (25)
          const auto &[gradv_c0, gradv_c1] = evaluateSecondOrderDerivativesImpl(this->r(), lfArgs);

          const auto graduTimesA = (gradu * alongMatrix.transpose()).eval();
          const auto gradvTimesA = (gradv * alongMatrix.transpose()).eval();

          const auto argsForDyz = lfArgs.extractSecondWrtArgOrFirstNonSpatial();

          const auto alonguAArgs          = replaceAlong(lfArgs, along(uTimesA));
          const auto alongvAArgs          = replaceAlong(lfArgs, along(vTimesA));
          const auto alonggraduTimesAArgs = replaceAlong(argsForDyz, along(graduTimesA));
          const auto alonggradvTimesAArgs = replaceAlong(argsForDyz, along(gradvTimesA));

          const auto u_xyzAlongv            = evaluateDerivativeImpl(this->l(), alongvAArgs);
          const auto v_xyzAlongu            = evaluateDerivativeImpl(this->r(), alonguAArgs);
          const auto v_c0c1AlongGraduTimesA = evaluateDerivativeImpl(this->r(), alonggraduTimesAArgs);
          const auto u_c0c1AlongGradvTimesA = evaluateDerivativeImpl(this->l(), alonggradvTimesAArgs);
          decltype(eval(u_xyzAlongv)) res;

          res = u_xyzAlongv + v_xyzAlongu + v_c0c1AlongGraduTimesA + u_c0c1AlongGradvTimesA;
          for (int i = 0; i < gridDim; ++i)
            res += (transpose(u_c1) * gradv_c0[i] + transpose(v_c1) * gradu_c0[i] + transpose(v_c0) * gradu_c1[i]
                    + transpose(u_c0) * gradv_c1[i])
                   * alongMatrix(0, i);

          return res;

        }
      }
    }
    u.evaluateDerivative(gpIndex, wrt(spatialAll,coeff(i)));

and the second one if it is called as

    u.evaluateDerivative(gpIndex, wrt(spatial(0),coeff(i)));

and the third without any spatial derivative is returned using

    u.evaluateDerivative(gpIndex, wrt(coeff(i,j)));

Therefore, this function separates the two wrt arguments and returns the corresponding first order derivatives. 6. Compile-time branch for the case where no spatial derivatives are requested but only the derivatives w.r.t. coefficients are needed. 7. Creates a new argument variable where the along argument is replaced by v. 8. This function evaluates the derivatives of l w.r.t. to both wrt arguments. Furthermore, it takes the along argument since otherwise the returned object would be a 3-dimensional array. If we consider the left function as u(ξ,uI) and v of the same size as u this calls returns u_xyAlongv =2uiuIuJvi This is the same if the following is called:

    u.evaluateDerivative(gpIndex, wrt(coeff(i,j)),along(v));
  1. Compile-time branch for the case where one spatial derivative and one derivative w.r.t. the coefficients is needed.
  2. This function evaluates the derivatives of l w.r.t. to both wrt arguments. If we consider the left function as u(ξ,uI) this calls returns u_xy=2uξuIoru_xy=2uξ0uI The first one would be returned if it is called as

       u.evaluateDerivative(gpIndex, wrt(spatialAll,coeff(i)));
    

    and the second one if it is called as

    u.evaluateDerivative(gpIndex, wrt(spatial(0),coeff(i)));
    

    In the first case the result is stored in an array. Thus, in the first index, the derivative w.r.t. to the first spatial coordinate is stored. Therefore, we have

    spatialAllCoeffDeriv = u.evaluateDerivative(gpIndex, wrt(spatialAll,coeff(i)));
    spatialAllCoeffDeriv[0] // derivative as in u.evaluateDerivative(gpIndex, wrt(spatial(0),coeff(i)));
    spatialAllCoeffDeriv[1] // derivative as in u.evaluateDerivative(gpIndex, wrt(spatial(1),coeff(i)));
    
  3. Compile-time branch for the case where one single spatial derivatives and one derivative w.r.t. coefficients is needed.

  4. Compile-time branch for the case where all spatial derivatives and one derivative w.r.t. coefficients is needed.
  5. The return type here is an array of single spatial derivatives and each derived w.r.t. the coefficient. Thus, the type inside the array must be deduced here.
  6. Compile-time branch for third order derivatives
  7. Compile-time branch for single spatial derivatives
  8. To obtain derivatives w.r.t. to the second and third wrt argument, we extract the arguments here. E.g., for the following request:

    u.evaluateDerivative(gpIndex, wrt(spatialAll,coeff(i,j)),along(matrix));
    

    This call would extract the arguments as

    newArgs =  "wrt(coeff(i,j)),along(matrix))" //THIS IS NO VALID SYNTAX
    

    This can then be used as

    u.evaluateDerivative(gpIndex, newArgs);
    
  9. As in the second order derivative case, it returns all the three first order derivatives. E.g., for the case,

    u.evaluateDerivative(gpIndex, wrt(spatialAll,coeff(i,j)),along(matrix));
    

    the returned values would be u_x=uξu_y=uuIu_z=uuJ 18. This returns the derivatives w.r.t. the given spatial direction and w.r.t. the first and second coefficient. E.g., for the case,

    u.evaluateDerivative(gpIndex, wrt(spatialAll,coeff(i,j)),along(matrix));
    

    the returned values would be u_xy=2uξuIu_xz=2uξuJ 19. Creates a new argument variable where the along argument is replaced by v_x. 20. This return call would be

    u.evaluateDerivative(gpIndex, wrt(spatial(0),coeff(i,j),along(v));
    

    In mathematical notation, it returns u_xyzAlongv  =3uiξ0uIuJvi 21. This return would be

    v_x = v.evaluateDerivative(gpIndex, wrt(spatial(0));
    u_yzAlongvx = u.evaluateDerivative(gpIndex, wrt(coeff(i,j),along(v_x));
    

    In mathematical notation, it returns u_yzAlongvx=2uiuIuJ[vξ0]i 22. Compile-time branch for all spatial derivatives 23. Obtain the along argument defined, for example, in

    u.evaluateDerivative(gpIndex, wrt(spatialAll,coeff(i,j),along(matrix));
    
  10. Similar to the single spatial case

  11. Similar to the single spatial case

If your expression is working, it can be added to dune/localfefunctions/expressions.hh by submitting a PR to dune-localfefunctions.


  1. Oliver Sander. DUNE—The Distributed and Unified Numerics Environment. Volume 140. Springer Nature, 2020. doi:10.1007/978-3-030-59702-3

  2. Philipp Grohs, Hanne Hardering, Oliver Sander, and Markus Sprecher. Projection-based finite elements for nonlinear function spaces. SIAM J Numer Anal, 57(1):404–428, 2019. doi:10.1137/18M1176798