Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
<tt>field_lazy</tt>

field expressions: concept, grammar and class hierarchy

Description

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.

Expression features

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:

  1. lvalues: they could be located on left- and the right-hand-sides of an assignment. A lvalue is thus writable and its interface provides a loop on the degrees-of-freedom (dof). The base class for a lvalue is thus called 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.
  2. assembled rvalues: they could be only located on left–hand-side of an assignment. They are read-only and provide a loop on the degrees-of-freedom for accessing the field values. The base class for an assembled rvalue is thus called field_rdof_base. It could be any linear expression involving field, possibly sliced or masked.
  3. unassembled rvalues: they also could be only located on left–hand-side of an assignment. They do not provide a loop on the degrees-of-freedom: instead, for obtaining the field values, a summation of all local contributions on the finite element 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

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.

Duality by integration

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.

Detailed specifications

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.


Implementation

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.