Rheolef
7.2
an efficient C++ finite element environment
|
field expressions: concept, grammar and class hierarchy
This page explains the design of the class hierarchy used for the representation of field
expressions and how they fit together. Casual users probably need not concern themselves with these details, but it may be useful for both advanced users and Rheolef's developers.
Expressions involving fields are computed in a lazy way: nothing is done until a conversion to the field
class is explicitly required. This approach allows many memory and computing time optimizations. Complex expressions are evaluated in only one loop instead of several ones with temporary variables.
Such an efficient implementation is obtained by using the expression template technique in C++. In consequence, if uh
is a field
, then, for instance, uh+uh
is no longer a field
: it is a field expression. It leads to represent at compile time an expression tree by a hierarchy of classes. Each class represent a node of the expression tree while leaves, also called terminals, are the field themselves.
Three kind of field expressions are identified:
field_wdof_base
. It could be a field
itself or any sliced or masked sub-indexing operation on it, such as field_wdof_sliced
or field_wdof_indirect
.field_rdof_base
. It could be any linear expression involving field
, possibly sliced or masked.geo
is required. The base class for an unassembled rvalue is thus called field_lazy_base
. It could be the result of the interpolate
or the integrate
functions, and either linear or nonlinear expressions involving them.This corresponds to the following base class inheritance diagram:
field_lazy_base <-- field_rdof_base <-- field_wdof_base <-- field
Any class that derives from a field_xxxx_base
is called here a field_xxx
and is represented by a C++ concept. Since concepts in C++ would be available only with the forthcoming 2020 standard, these concepts are actually represented by traits in Rheolef.
Nonlinear expressions, such as exp(uh)
are no more piecewise polynomial functions, but instead piecewise exponential-of-polynomial functions. Thus, they should be re-interpolated in a suitable finite element space
in order to recover a piecewise polynomial function, i.e. a true field
.
Note that degrees-of-freedom are in general not directly the field values at some mesh node locations, except for nodal elements, such as Lagrange elements. For instance, the Raviart-Thomas element involves non-nodal degrees-of-freedom. Thus, exp(uh)
could not be evaluated in general by applying the exp
function to all degrees-of-freedom. Instead it requires an element-by-element loop for converting the exp(uh)
local field values to degrees-of-freedom. In consequence, the result of interpolate
could not be a field_rdof
, but instead a field_lazy
.
Finally, a field_lazy
could be or not a piecewise polynomial function, depending upon the linearity of the considered field expression. Since the considered field expression is known at compile-time, this property is also identified at compile-time: it is available inside all the field_lazy
classes as a static boolean flag named is_field_piecewise_polynomial
. When this flag is true, a field_lazy
object could be directly assigned to a field_wdof
. Otherwise, an invocation to interpolate
is required. If such an invocation is missing, a compile-time error message is issued.
Instead of using an invocation to interpolate
, a general nonlinear field_lazy
expression could be converted to a piecewise polynomial function by an invocation to integrate
after a multiplication by a suitable test
function. From a mathematical point of view, the result of integrate
is a linear application on a finite element space
associated to the test
function. The linear application belongs to the dual of the finite element space
. Since, in finite dimension, the dual space could be identified to the space itself, this linear application could be represented by a field
, the so-called Riesz representative. This is precisely the result returned by the integrate
function.
The argument to the integrate
function is a variational expression. The base class for a variational expression is thus called field_test_base
and its associated concept is field_test
. It could be a test
itself, any multiplication of a field_test
by a field_lazy
and any linear combination of field_test
expressions.
The field_wdof
specification could be expressed recursively as a grammar:
field_wdof: field | field_wdof [integer] # field_wdof_sliced | field_wdof [domain] # field_wdof_indirect
where domain
stands for any string
or geo
variable. Here, we used a conventional yacc(1) like grammar notations and, when a class proxy is used by the implementation, its name is indicated on the right, after the '#' symbol.
Conversely, the field_rdof
concept writes:
field_rdof: field_wdof | field_rdof [integer] # field_rdof_sliced_const | field_rdof [domain] # field_rdof_indirect_const | +- field_rdof # field_rdof_unary | scalar +- field_rdof # field_rdof_unary | field_rdof +- scalar # field_rdof_unary | scalar * field_rdof # field_rdof_unary | field_rdof * scalar # field_rdof_unary | field_rdof / scalar # field_rdof_unary | field_rdof +- field_rdof
where +-
stands for either the +
or the -
operator and *
could be replaced by dot
or ddot
calls, for respectively vector- or tensor-valued field expressions. The constant
is any scalar, vector or tensor valued constant expression, represented respectively by Float
, point
and tensor
. The const
specifier behind the field
means that this is a read-only variable.
The field_lazy
concept is defined by:
field_lazy: field_rdof | constant | interpolate (space, field_lazy) | integrate (opt_domain field_test opts) | +- field_lazy | field_lazy +- field_lazy | field_lazy * / field_lazy | grad (field_lazy) | f (field_lazy...) | compose (g, field_lazy...) | compose (field_lazy, X)
The opt_domain
stands for an optional domain specification and opts for an optional integrate_option
variable. The grad
operator could be any other classical diffferential, e.g. div
or curl
. Also, f
stands for any standard mathematical function, e.g. exp
, and g
for any used-defined one, while the ...
means an argument list of any length. The X
argument is any characteristic
variable. Finally, the field_test
concept involved by the integrate
call is defined by:
field_test: test | field_test*field_lazy | field_lazy*field_test
Note the crossed recursive definitions of the field_lazy
and field_test
concepts.
In Rheolef, the previous field's expression class hierarchy is designed so that virtual functions are avoided where their overhead would significantly impair performance: this approach is referred in C++ as static polymorphism. Rheolef achieves static polymorphism with the curiously recurring template pattern (CRTP). In this pattern, the base class, for instance, field_wdof_base
, is in fact a template class, and the derived class, or instance, field_wdof_sliced
, inherits the base class with the derived class itself as a template argument. In this case, field_wdof_sliced
inherits from field_wdof_base<field_wdof_sliced>
. This allows Rheolef to resolve the polymorphic function calls at compile time.
Finally, the implementation of field
expression heavily relies on C++ template techniques and the interested reader is invited to consult the following very nice book:
D. Vandevoorde and N. M. Josuttis and D. Gregor, C++ templates: the complete guide, 2nd ed., Addison Wesley, 2017.