Variational formulations for elliptic problems require the specification of the functions \( F \) and \( G \) as outlined in the basics. This is done by writing a class that adheres to the NonlinearVariationalFunctional concept, see the doxygen documentation. In the following, we will go through the essential components of this class.
The first ingredient is to specify whether this is a symmetric problem with a variational functional (like the Poisson equation), or just a weak formulation (like advection-diffusion equations). This is most conveniently done by deriving from the base class FunctionalBase, that accepts a template parameter for that: The enum values VariationalFunctional and WeakFormulation can be used.
In addition, the scalar field type (usually double, but could also be complex, float, or some automatic-differentiation-augmented type) must be announced.
class Functional: public Kaskade::FunctionalBase<Kaskade::VariationalFunctional>
{
public:
using Scalar = double;
Next the number and type of variables and test functions. For systems such as Stokes flow or KKT systems in optimal control, there are multiple variables coming from different function spaces, which need to be specified statically by providing a VariableSetDescription<...> type - containing a static list of variables, their number of components, and the index of the space they belong to. This is often passed as a template parameter to the class we define, such that the finite element spaces to be used can easily be configured from outside. For symmetric problems, ansatz and test variables are the same, but for nonsymmetric weak formulations they might be different.
Moreover, for nonlinear problems one may want to provide a linearization point of the equations that is not contained in the ansatz space (e.g., an analytically given approximation), and then compute a correction (Newton step) in the ansatz space. Therefore, the type of variables used to describe this linearization point has to be given as well - most often, it is just the type of ansatz variables itself.
using AnsatzVars = VariableSetDescription<...>;
using TestVars = AnsatzVars;
using OriginVars = AnsatzVars;
Then the functions \( F \) and \( G \) need to be defined, providing also their derivatives. This is done by defining a member class for each of these two functions, the DomainCache for \( F \) and the BoundaryCache for \( G \).
These mainly define the function values (zeroth order derivative, method d0, only for variational problems), and first and second order directional derivatives (methods d1 and d2, respectively). These classes are instantiated multiply for multi-threaded assembly of the stiffness matrix and right hand side, and shall therefore be thread-safe.
class DomainCache
{
public:
They are constructed by the assembler (or by the convenience helper function linearizationAt) using a particular constructor convention. The variational functional object F is provided, as well as the linearization point \( u \) and a bitfield of flags announcing which subset of the three values shall be computed - the matrix (d2), the right hand side (d1), or the value of the functional (d0).
DomainCache(Functional const& F,
typename OriginVars::VariableSet const& u,
int flags);
For evaluating the integrals outlined in the basics, the function \( F \) must be evaluated at a position \( x \in \Omega \), which is provided in two steps - as the grid cell \( T \) in which the point lies, and its cell-local coordinate \( \xi \in T \). This two-step procedure is done for performance reasons, since some computations such as index look-ups need be done only once for each cell, and not again for each quadrature point. In addition, again for performance reasons, an object eval is provided, that already contains shape functions evaluated at the given local position and can be used for efficient evaluation of the linearization point variables.
void moveTo(Cell const& cell);
void evaluateAt(LocalPosition const& xi, Evaluators const& eval);
Then the actual function and derivative evaluations need to be provided. For the variational functional value \( J(u) \) this is simple, as it only depends on the linearization point \( u \) given at DomainCache construction.
Scalar d0() const;
The derivatives are directional derivatives in the given direction, for which value and spatial derivative are provided in a VariationalArg object. The direction is specified for a particular variable identified by its number, i.e. the row number for test variables and column number for ansatz variables. row and column are specified statically as template parameters.
template <int row>
D1Result<TestVars,row> // a vector type
d1(VariationalArg<Scalar,dim> const& argT) const;
template <int row, int col>
D2Result<TestVars,row,
AnsatzVars,col> // a matrix type
d2(VariationalArg<Scalar,dim> const& argT,
VariationalArg<Scalar,dim> const& argA) const;
};
The return type deserves some explanation. For multiple-component variables such as displacement in solid mechanics or velocity in fluid mechanics, the finite element functions are defined in terms of scalar FE function spaces such as the ubiquitous Lagrange spaces. Consequently, the test function \( \phi \) is scalar, with value and derivative provided in argT, whereas conceptually the test function would be vectorial with dimension \( m \). Now the \( i \)-th entry of the return vector of d1 shall contain the directional derivative in direction of the (vectorial) test function \( \phi e_i \), where \( e_i \) is the \( i \)-th unit vector.
For the second derivative, of course, two directions are given, related to the test space and the ansatz space, respectively. The same return value construction as for d1 is used, such that for multi-component variables a matrix is returned.
The BoundaryCache does exactly the same for the function \( G \), with the only exception that instead of a grid cell, here a FaceIterator pointing to the current boundary face is provided.
PDEs come in a huge variety, and thus we need to specify some structural information in addition. This is done via the member class templates D1 and D2, parametrized by row (i.e. equation number) and column (i.e. variable number). They express properties of individual blocks of the stiffness matrix in particular for the case of PDE systems - if a block is present in the first place, if it is symmetric, and if it shall be lumped.
template <int row, int col>
struct D2
{
constexpr static bool present = true;
constexpr static bool symmetric = true;
constexpr static bool lumped = false;
};
Non-present blocks will not be assembled at all, and the assembly of symmetric blocks is somewhat more efficient than of nonsymmetric blocks.
For the corresponding D1, parametrized only by row, in particular the presence of right hand side blocks is of interest.
Finally, we need to specify which quadrature order shall be used on any grid cell or boundary face, given the polynomial order of the shape functions.
int integrationOrder(Cell const& cell,
int shapeFunctionOrder,
bool boundary);
};
Have a look at the tutorial section for more insights into the implementation of different types of problems.
template <class FaceIterator>
void moveTo(FaceIterator const& fi)
{
useNitsche = fi->centerUnitOuterNormal()[0] > 0.5;
if (useNitsche)
nitsche.moveTo(fi);
else
dirichlet.moveTo(fi);
}
template <class Position, class Evaluators>
void evaluateAt(Position const& xi, Evaluators const& evaluators)
{
auto u = vars.template value<uIdx>(evaluators);
if (useNitsche)
{
auto ux = vars.template derivative<uIdx>(evaluators);
nitsche.setBoundaryData(xi,penalty,u,uDirichletBoundaryValue,ux,unitMatrix<Scalar,dim>());
}
else
dirichlet.setBoundaryData(penalty,u,uDirichletBoundaryValue);
}
Scalar d0() const
{
if (useNitsche)
return nitsche.d0();
else
return dirichlet.d0();
}
template<int row>
Dune::FieldVector<Scalar,1> d1(Kaskade::VariationalArg<Scalar,dim> const& argT) const
{
if (useNitsche)
return nitsche.d1(argT);
else
return dirichlet.d1(argT);
}
template<int row, int col>
Dune::FieldMatrix<Scalar,1,1> d2(Kaskade::VariationalArg<Scalar,dim> const &argT,
Kaskade::VariationalArg<Scalar,dim> const &argA) const
{
if (useNitsche)
return nitsche.d2(argT,argA);
else
return dirichlet.d2(argT,argA);
}
Types of penalty approaches
DirichletPenaltyBoundary<typename AnsatzVars::GridView,1> dirichlet;
DirichletNitscheBoundary<typename AnsatzVars::GridView,1> nitsche