Draft: Unified transformed local finite elements

Simon Praetoriusrequested to merge
feature/transformed-finite-elements into master

Motivation

dune-Functions currently needs several different kinds of local finite-element transformations:

  • Scalar bases such as Lagrange keep their values unchanged, but derivatives must be transformed from reference to physical coordinates.
  • Raviart-Thomas, Brezzi-Douglas-Marini, and Nedelec elements transform their vector-valued ranges with contravariant or covariant Piola maps. Their natural derivatives are divergence and curl.
  • Cubic Hermite, Morley, and Argyris elements are not affine-equivalent and need an element-dependent transformation of the complete basis-function set.
  • More complex elements may need several transformations in sequence and may change scalar, vector, or tensor range types.

These cases were previously implemented through unrelated wrappers, mixins, and element-specific finite elements. This merge request introduces one experimental framework for describing them while retaining element-specific transformation details in the corresponding basis headers.

The framework also replaces the fixed value/Jacobian interface with derivative tags. A finite element only implements the quantities that are mathematically meaningful for it, for example Derivatives::Divergence for an H(div) element or Derivatives::Curl for an H(curl) element.

Central abstractions

TransformedLocalFiniteElement

TransformedLocalFiniteElement is the common finite-element adapter. It binds:

  1. A reference local finite element.
  2. An element-dependent context.
  3. A transformation policy for physical-basis evaluation.
  4. Optionally, the physical-function pullback provided by the same policy.

It provides:

  • referenceBasis() for the untransformed reference basis.
  • physicalBasis() for values and derivatives on the physical element.
  • localBasis(), selected through LocalBasisMode, for compatibility with existing basis-node conventions.
  • The reference local coefficients.
  • localInterpolation(), selected independently through InterpolationMode.
  • precomputeIdentity() for caches that must distinguish reference finite elements, polynomial orders, and orientation variants.

LocalBasisMode::reference is used where existing callers expect localBasis() to retain the classical reference-element meaning, notably for scalar Lagrange-like bases. LocalBasisMode::physical is used for elements whose traditional local basis is already the mapped physical basis.

InterpolationMode::reference exposes the reference interpolation and is the default. InterpolationMode::transformed uses Transformation::localFunctionPullback() before invoking the reference interpolation. Basis and interpolation selection are independent: there is no separate interpolation-transformation type. In transformed mode, one policy constructor argument initializes independent basis and interpolation adapters, which are then bound to the same context.

The default precomputation identity uses the address of the bound reference finite element. A reference finite element can provide its own precomputeIdentity() returning LocalBasisPrecomputeIdentity, including a generation counter when the same storage can be rebuilt with different contents.

TransformedLocalBasis

The transformed basis exposes derivative-dependent types and operations:

template<class Derivative>
using PrecomputeBuffer = ...;

template<class Derivative>
using DerivativeRange = ...;

precompute(derivative, x, buffer);
finalize(derivative, x, buffer, values);
evaluate(derivative, x, values);

The transformation policy determines which derivative tags are available and the physical result type of each derivative. The current tag vocabulary is:

  • Value
  • Jacobian
  • Gradient
  • Divergence
  • Curl
  • Hessian
  • Laplacian

The combined evaluate() operation is the convenient user interface. precompute() and finalize() form the cache-oriented interface.

Bind contexts

A transformation is bound to a context rather than directly to a geometry. Reusable contexts include:

  • ElementBindContext
  • GeometryBindContext
  • SimplexVertexMeshSizeContext
  • SimplexEdgeOrientationContext
  • SimplexVertexMeshSizeAndEdgeOrientationContext

The basic contexts provide a stable geometry snapshot. Richer contexts add element information required by non-affine-equivalent elements, such as edge orientations and per-vertex mesh sizes.

Transformation classes

Reference evaluation

ReferenceLocalBasisEvaluator evaluates the mathematical operator named by a derivative tag on the reference element:

  • values use evaluateFunction();
  • Jacobians use evaluateJacobian();
  • gradient, divergence, and curl use dedicated local-basis methods when available and otherwise derive the requested operator from reference Jacobians;
  • Hessians use evaluateHessian() when available and otherwise use second partial derivatives;
  • Laplacians use evaluateLaplacian() when available and otherwise take the trace of reference Hessians.

ReferenceEvaluation separates the requested physical derivative from the reference operator required by a transformation. Its operator policy normally keeps the derivative unchanged, but may remap it. For example, a physical Laplacian requires the complete reference Hessian. This also allows an element author to provide a custom reference evaluator or operator-selection policy without modifying the pipeline.

Piola transformations request reference values and their natural reference operators directly: divergence for contravariant Piola and curl for covariant Piola. They do not evaluate a complete reference Jacobian merely to extract the natural derivative.

Ordered evaluation pipelines

BasisEvaluationPipeline applies an ordered sequence of typed transformation stages after reference evaluation. Each stage declares its output type as a function of the derivative tag, local basis, context, and input range. Stages may therefore change scalar, vector, or tensor dimensions.

Two reusable stage adapters distinguish the main transformation classes:

  • BasisSetTransformationStage transforms the complete vector of basis functions and may mix basis-function indices.
  • RangeTransformationStage transforms each basis-function result independently and may change its range type.

The selected reference operator is evaluated in precompute(). Intermediate stage vectors are stored in the caller-owned precompute buffer and are reused between finalizations. The cached reference evaluation remains immutable.

Geometry derivative transformation

GeometryDerivativeTransformation is a pointwise range transformation for scalar physical derivatives. GeometryDerivativeStage adapts it for use in a pipeline.

It implements:

  • unchanged scalar values;
  • physical Jacobians;
  • ambient tangential gradients on embedded geometries;
  • affine, codimension-zero Hessians and Laplacians.

Non-affine Hessians and Laplacians are deliberately rejected. Their complete chain rule needs higher derivatives of the geometry that are not readily available through the current geometry interface.

Piola transformations

ContravariantPiolaTransformation implements H(div) values, divergence, and the corresponding interpolation pullback.

CovariantPiolaTransformation implements H(curl) values, curl, and the corresponding interpolation pullback.

Both support rectangular geometry Jacobians and separate reference and physical range types. Piola transformations remain complete transformation policies rather than pipeline stages because reference operator selection, natural derivative mapping, and interpolation pullback are tightly coupled.

Element-specific basis-set transformations

The transformation formulas for cubic Hermite, Morley, and Argyris remain in their respective basis headers. They are wrapped in BasisSetTransformationStage and composed with GeometryDerivativeStage.

This keeps specialized finite-element mathematics close to the element while reusing the common binding, evaluation, derivative, and caching infrastructure.

Implementing another transformed finite element

A finite-element author provides:

  1. A reference local finite element with basis, coefficients, and interpolation.
  2. A bind context containing the geometry and any additional element data, such as orientations or vertex mesh sizes.
  3. A basis transformation policy with derivative-dependent PrecomputeBuffer and DerivativeRange types, bind(), precompute(), and finalize(). For staged transformations, this is normally a BasisEvaluationPipeline.
  4. If physical functions must be pulled back for interpolation, the same transformation policy additionally provides localFunctionPullback() and the adapter selects InterpolationMode::transformed. Otherwise the default InterpolationMode::reference exposes the reference interpolation.
  5. A LocalBasisMode selecting whether compatibility localBasis() denotes the reference or physical basis.

Pointwise transformations should normally be implemented as a policy wrapped in RangeTransformationStage. Transformations mixing basis-function indices should use BasisSetTransformationStage. Piola-like mappings whose natural reference operator and interpolation pullback are inseparable can remain complete transformation policies.

For example, a Piola finite element uses one transformation type for both operations:

using FiniteElement = TransformedLocalFiniteElement<
  ReferenceFiniteElement,
  Context,
  PiolaTransformation,
  LocalBasisMode::physical,
  InterpolationMode::transformed>;

Scalar derivative and basis-set transformations normally keep the default reference interpolation.

If several reference finite elements can coexist, for example orientation or hp variants, their precomputation identities must differ. The default reference finite-element address provides this distinction while the owning finite-element map remains unchanged. Rebuilt storage can provide an explicit generation-aware identity.

Migrated finite elements

The following bases use the new framework:

Transformation family Finite elements
Scalar physical derivatives Lagrange, refined Lagrange, Rannacher-Turek
Contravariant Piola Raviart-Thomas, Brezzi-Douglas-Marini
Covariant Piola Nedelec
Double-Contravariant Piola Hellan-Hermann-Johnson !596
Basis-set transformation followed by geometry derivative transformation Cubic Hermite, Morley, Argyris

The discrete global-basis derivative implementation now detects a physical basis and evaluates its physical Jacobian directly. Conventional finite elements retain the previous reference-Jacobian chain-rule path.

Caching

Reference-basis evaluation is independent of the current physical geometry. The API therefore separates evaluation into:

  1. precompute(): evaluate and retain the required reference quantity in a buffer that also owns reusable pipeline scratch storage.
  2. finalize(): apply the transformations that depend on the bound context.

QuadratureBasisFunctionCache demonstrates this protocol:

  • It stores precompute buffers independently for every basis-tree leaf and derivative tag.
  • A cache key combines LocalBasisPrecomputeIdentity with the address of the quadrature-rule object. Equality is checked on the complete key, so hash collisions cannot produce incorrect evaluations.
  • Rules returned by QuadratureRules::rule() reside in global storage and have stable process-lifetime addresses. Other quadrature-rule objects must outlive the cache.
  • The default basis identity is the bound reference finite-element address, distinguishing RT, BDM, and Nedelec orientation variants as well as independently stored hp variants.
  • Precomputations are stored in an unordered_map per derivative tag.
  • It rebinds cached leaf entries to the current local-view nodes.
  • It reuses reference evaluations only when both the reference finite-element identity and quadrature-rule identity match.
  • It repeats finalization for every bound physical element.
  • clear() explicitly invalidates precomputations and bound-node references. It must be called after updating a basis unless custom generation-aware identities are provided.
  • It exposes basis(), size(), and order() as a small assembly-facing facade.

For multi-stage transformations, intermediate vectors are part of the precompute buffer. This avoids allocating temporary vectors in every finalize() call while preserving the reference part of the buffer as immutable.

Tests and compatibility

The transformed-basis tests use derivative tags and compare supported derivatives with finite-difference approximations. Additional coverage includes:

  • scalar tangential gradients on embedded geometries;
  • covariant and contravariant surface Piola maps;
  • divergence and curl transformations;
  • Piola interpolation duality;
  • composed basis-set and geometry transformations;
  • reusable pipeline scratch storage;
  • explicit rejection of non-affine Hessians;
  • hybrid triangle/quadrilateral quadrature caching;
  • heterogeneous derivative result types in composite bases;
  • distinct cache entries for Piola orientation variants;
  • explicit cache invalidation and reuse after clear();
  • rebinding after local-view destruction;
  • physical derivatives of discrete Argyris, Morley, and cubic Hermite functions.

The transformed-basis-benchmark example assembles an RT mass plus divergence-divergence form and compares manual reference evaluation and Piola mapping, combined evaluate(), split precompute()/finalize(), and quadrature-cache evaluation. All paths are checked against the same assembly checksum. It also contains an experimental vector-backed cache lookup for comparison with the installed unordered_map implementation.

The classical localBasis() and localInterpolation() behavior is preserved per basis through LocalBasisMode and InterpolationMode. A constrained compatibility evaluateJacobian() remains available for affine, codimension-zero Piola elements without adding a generic Piola Jacobian derivative tag.

Edited by Simon Praetorius

Merge request reports