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.