Discussion: Tensor-like output arguments for local basis evaluation

Summary

This issue proposes evolving local basis evaluation from fixed std::vector<RangeType> output containers toward tensor-like output views.

The immediate motivation is to support:

  • direct evaluation into caller-owned storage;
  • one contiguous allocation for all quadrature points;
  • scalar values and derivatives without artificial FieldVector<T,1> and FieldMatrix<T,1,dim> wrappers;
  • strided and sliced output storage;
  • reuse of the same evaluation interface in transformed local finite elements;
  • fewer allocations in quadrature caches and transformation pipelines.

An experimental implementation exists in dune-localfunctions/dune/localfunctions/lagrange/lagrangesimplex.hh, together with TensorView and TensorWrapper utilities.

This is intended as a discussion of the interface and migration strategy, not as a request to immediately replace the established local basis API.

Current interface

The classical local basis interface chooses both the mathematical result type and its storage:

void evaluateFunction(
  DomainType const& x,
  std::vector<RangeType>& out) const;

void evaluateJacobian(
  DomainType const& x,
  std::vector<JacobianType>& out) const;

For a scalar basis this commonly produces:

std::vector<FieldVector<T,1>>
std::vector<FieldMatrix<T,1,dim>>

This representation is compatible with existing DUNE interfaces, but has some limitations:

  • the local basis controls the output container;
  • evaluating at several quadrature points naturally creates nested vectors;
  • scalar ranges retain artificial dimensions of extent one;
  • alternative layouts require copying after evaluation;
  • views into larger assembly buffers cannot be passed directly;
  • transformations often need adapters to reshape scalar, vector, and matrix results.

Experimental interface

The experimental overload accepts a writable tensor-like object:

template<class Out>
  requires (Out::rank() == 1)
void evaluateFunction(DomainType const& x, Out out) const;

template<class Out>
  requires (Out::rank() == 2)
void evaluateJacobian(DomainType const& x, Out out) const;

The minimal output interface is currently:

Out::rank();
out.extent(i);
out(i, j, k, ...);

Unlike std::mdspan, contiguity is not required. The accessor may refer to a contiguous array, strided storage, or a nested legacy container.

For a scalar basis, the natural shapes are:

value:              [basis function]
gradient/Jacobian:  [basis function, coordinate]
hessian:            [basis function, coordinate, coordinate]

For a vector-valued basis they become:

value:              [basis function, range component]
Jacobian:           [basis function, range component, coordinate]

The existing overloads remain available and forward through an accessor-based view:

void evaluateFunction(DomainType const& x,
                      std::vector<RangeType>& out) const
{
  out.resize(size());
  evaluateFunction(x, TensorView(
    [&out](auto i) -> auto& { return out[i][0]; },
    extents(size())));
}

This keeps the old interface source-compatible while making the tensor interface the implementation primitive.

Tensor views and wrappers

Two related adapters are being explored:

  • TensorView stores an accessor and extents. It can represent arbitrary writable storage and supports slicing by fixing its first index.
  • TensorWrapper recursively adapts nested containers such as std::vector<FieldMatrix<...>> to the tensor interface.

Slicing is relevant for quadrature-wide storage. Given a tensor with shape

[quadrature point, basis function, components...]

the expression

auto pointOutput = output[iq];

produces the output view expected by a pointwise local basis evaluation.

Before these utilities become shared API, several details need resolution:

  • index_type should use the extents' index type, not its rank type.
  • Inferring nested extents from container[0] is invalid for empty containers.
  • Nested containers must be rectangular, or their extents must be supplied explicitly.
  • The required interface should not imply mdspan properties such as contiguity, uniqueness, or striding unless they are actually required.
  • The ownership and lifetime rules of accessor-based slices need clear documentation.

Application to transformed local finite elements

The experimental transformed finite-element framework in dune-functions uses:

precompute(derivative, x, precomputeBuffer);
finalize(derivative, x, precomputeBuffer, output);
evaluate(derivative, x, output);

Currently, output and intermediate storage are based on:

std::vector<DerivativeRange>

The pipeline therefore still materializes vectors of FieldVector, FieldMatrix, or tensor-like aggregate objects. A tensor-output overload only in the reference basis would retain most of these intermediate allocations and shape conversions.

A complete tensor-oriented pipeline would additionally need derivative shape metadata, for example:

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

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

The complete basis evaluation then has the shape:

[basis function, DerivativeExtents...]

Pipeline stages can use this structure as follows:

  • A pointwise range transformation fixes the basis-function index and transforms one component tensor.
  • A basis-set transformation receives the complete basis-function dimension and may mix entries along that dimension.
  • Piola transformations operate directly on tensor components instead of adapting aggregate values through vector-view helpers.

The tensor representation also gives scalar and vector-valued bases one common shape model. It avoids treating FieldVector<T,1> as both a scalar and a one-component vector depending on context.

Application to quadrature caching

The current quadrature cache stores evaluations approximately as:

std::vector<std::vector<DerivativeRange>>

The outer vector indexes quadrature points and each inner vector owns the evaluation of all local basis functions at one point.

A tensor-oriented cache could instead own one allocation with shape:

[quadrature point, basis function, derivative components...]

Potential benefits include:

  • one allocation instead of one allocation per quadrature point;
  • less pointer chasing;
  • predictable memory layout;
  • improved cache locality;
  • direct evaluation into the final cache storage;
  • easier transfer to matrix-free and batched assembly kernels;
  • potential compiler vectorization over basis functions or components.

With a right-major layout and quadrature as the first index, a complete pointwise basis evaluation remains contiguous:

basis.evaluate(tag, x, output[iq]);

Other layouts should remain possible when an application wants a component-major or basis-major organization.

Reference data and transformation workspace

The transformed pipeline currently stores immutable reference evaluations and mutable intermediate vectors together in each point's precompute buffer.

For a quadrature rule this can retain one complete intermediate workspace per quadrature point, although finalization normally processes the points sequentially.

A tensor-oriented cache suggests separating:

Cached reference data:
  [quadrature point, basis function, reference components...]

Reusable transformation workspace:
  one set of intermediate tensors

Physical result:
  [quadrature point, basis function, physical components...]

This would preserve immutable cached reference data while reusing a single workspace during finalization. It may reduce both allocation count and retained memory, especially for multi-stage transformations.

This separation is independent of the public tensor interface and could be implemented incrementally.

Expected performance characteristics

The expected improvements are primarily related to storage and allocation:

  • fewer dynamic allocations;
  • fewer nested container indirections;
  • better spatial locality;
  • reusable caller-owned buffers;
  • direct evaluation into assembly-specific layouts.

There is no guarantee that an accessor-based TensorView is always as fast as a concrete mdspan. A callable accessor should often inline, but a standard layout mapping exposes stronger information about strides and contiguity.

Benchmarks should compare at least:

  1. Existing std::vector<RangeType> evaluation.
  2. A direct contiguous mdspan or mdarray.
  3. An accessor-based TensorView.
  4. A legacy nested container adapted through TensorWrapper.
  5. Pointwise quadrature evaluation into nested vectors.
  6. Evaluation into one rule-wide contiguous tensor.
  7. Transformed evaluation with per-point and shared intermediate workspaces.

Both optimized and debug builds are relevant: optimized builds measure inlining and vectorization, while debug builds expose abstraction overhead that may affect development workflows.

User-interface implications

The tensor interface gives users control over storage:

basis.evaluate(tag, x, outputView);

This is useful for:

  • normal owning containers;
  • mdspan and mdarray;
  • slices of larger element or quadrature buffers;
  • strided views;
  • device or accelerator storage, provided the accessor is suitable;
  • compatibility wrappers around established containers.

However, callers must know the required extents. The interface therefore needs a reliable shape query. Possible forms include:

auto extents = basis.extents(tag);

or traits combining static component extents with the runtime number of basis functions.

The return type of a slice is another compatibility question. Existing code expects aggregate operations such as:

values[i] * values[j]
jacobians[i][0][j]

A generic tensor slice may not provide FieldVector arithmetic. We need to decide whether:

  • tensor views are only writable evaluation targets;
  • cache access returns tensor views plus tensor algorithms;
  • common rank-one and rank-two slices adapt to existing vector/matrix types;
  • compatibility accessors remain available alongside the tensor-native API.

Possible migration path

A staged migration would reduce the risk of combining storage, derivative typing, and transformation changes in one step.

  1. Define and document a minimal writable tensor-output concept.
  2. Harden TensorView and TensorWrapper, including empty-container and rectangular-storage behavior.
  3. Add tensor-output overloads to selected dune-localfunctions bases while retaining the classical overloads.
  4. Add explicit derivative field and shape traits.
  5. Prototype contiguous rule-wide storage in the dune-functions quadrature cache.
  6. Separate immutable reference data from mutable transformation workspace.
  7. Add tensor-output overloads to TransformedLocalBasis.
  8. Migrate reusable pipeline stages and Piola transformations.
  9. Evaluate whether the tensor interface should become the primary local basis API in a later release.

Open questions

Feedback is requested in particular on the following points:

  1. Is the minimal interface rank(), extent(i), and operator() sufficient, or should it use a standard mdspan concept more directly?
  2. Should arbitrary non-contiguous accessors be supported, or should evaluation require a strided layout?
  3. Where should TensorView, TensorWrapper, and the tensor-output concept live: dune-common, dune-localfunctions, or another module?
  4. How should derivative shapes be represented for scalar, vector, matrix, and higher-rank finite elements?
  5. Should the basis-function index always be the leading tensor dimension?
  6. Should local basis evaluation resize owning outputs, or require correctly sized views and expose a separate allocation helper?
  7. Should quadrature-point batching be part of the local basis interface, or remain an application/cache concern implemented through slicing?
  8. How should tensor slices interoperate with existing FieldVector and FieldMatrix algorithms?
  9. Is a shared transformation workspace acceptable for sequential finalization, and how should parallel finalization request independent workspaces?
  10. Which existing local finite elements should be used as representative prototypes beyond scalar Lagrange bases?

Initial recommendation

The proposed direction appears useful, especially for quadrature caching and transformed finite elements. The first common API should remain deliberately small and should not require contiguous storage.

The next useful prototype would combine:

  • tensor-output evaluation for a scalar Lagrange basis;
  • legacy compatibility through TensorView;
  • one rule-wide contiguous quadrature buffer;
  • a benchmark against the current nested-vector cache.

The transformed pipeline should be migrated only after shape traits and the basic output concept have proven practical in dune-localfunctions.