Ferrite-FEM / Ferrite-FEM/FerriteViz.jl

Tabulate reference shape values in transfer_solution instead of PointValues (7-15x on 3D grids)

Open
#153 0 comments 0 reactions 0 assignees View on GitHub
Dominant language
Julia
Stars
36
Forks
11
Avg merge
14d 16h
Merged PRs (30d)
1

Description

**Aimed at after #150 is merged** — this touches `transfer_solution`, which that PR already rewrites, so doing it now would only make the diff harder to review.

Follow-up to [this review comment](https://github.com/Ferrite-FEM/FerriteViz.jl/pull/150#discussion_r3652736897) on `src/dataset.jl`:

> TODO for later: Can we replace PointValues with a plain function values? We do not need the geometry here, which is dead weight.

The premise is right, but the geometry turns out to be the smaller half of the prize.

## The actual win: tabulate, don't just drop the geometry

`transfer_solution` currently does, per tessellation vertex per cell:

```julia
Ferrite.reinit!(pv, local_coords, ξ)
val = Ferrite.function_value(pv, 1, @views(u[celldofs_field]))
```

Every call re-evaluates every basis function at `ξ` and computes a Jacobian that is then discarded (for an `IdentityMapping` with `update_gradients=false`, `Nx` *is* `Nξ` — Ferrite aliases them at `FunctionValues.jl:70` and `apply_mapping!` is a literal no-op at `:185`).

But the tessellation puts the **same reference coordinates on every cell of a given reference shape**. So the shape values can be tabulated once into an `nbf × nverts_per_cell` table, and the inner loop becomes a small matrix–vector product. `PointValues` cannot express this — its API is `reinit!` per point.

Simply swapping `PointValues` for per-vertex `reference_shape_value` calls is roughly a wash. Tabulating is what pays.

## Measurements

Both variants given the same `Val(refdim)` function barrier that `_transfer_solution!` already uses, so neither is penalised by type instability (an earlier measurement without it was misleading):

| case | `PointValues` | tabulated | speedup |
|---|---|---|---|
| 2D quad, vector Lagrange 1 | 0.19 ms / 0.3 MiB | 0.06 ms / 0.3 MiB | **3.0×** |
| 2D quad, vector Lagrange 2 | 0.27 ms / 0.2 MiB | 0.06 ms / 0.2 MiB | **4.3×** |
| 3D hex, scalar Lagrange 1 | 0.84 ms / 0.2 MiB | 0.06 ms / 0.2 MiB | **15.2×** |
| 3D hex, vector Lagrange 2 | 1.84 ms / 0.3 MiB | 0.23 ms / 0.3 MiB | **8.0×** |
| 3D tet, scalar Lagrange 2 | 1.04 ms / 0.7 MiB | 0.14 ms / 0.6 MiB | **7.2×** |

This is in the interactive path — `transfer_solution` re-runs on every `FerriteViz.update!` — so it is worth more than the raw numbers suggest.

Correctness: bit-for-bit identical (`max |diff| == 0.0`, not merely within tolerance) against the current `PointValues` implementation for scalar and vectorised `Lagrange` orders 1–2, `DiscontinuousLagrange`, `Serendipity` and `CrouzeixRaviart`, on quadrilateral, triangle, hexahedron, tetrahedron and wedge grids.

## Design concern: keep the mapping a dispatch point, not an `if`

This optimisation is **only valid for `Ferrite.IdentityMapping`**. Ferrite currently defines three mapping types (`FunctionValues.jl:146-148`): `IdentityMapping`, `CovariantPiolaMapping`, `ContravariantPiolaMapping`. The Piola ones need the Jacobian *and* the cell (for `get_direction`), so they cannot share the tabulated path — see #151, where H(curl)/H(div) fields currently fail outright.

The two should therefore land as one dispatch on `Ferrite.mapping_type(ip_field)` rather than as two unrelated patches:

- `IdentityMapping` → tabulated reference shape values, no geometry at all.
- `CovariantPiolaMapping` / `ContravariantPiolaMapping` → `CellValues` over a quadrature rule whose points are the tessellation reference coordinates, with `reinit!(cv, cell, coords)`. Public API, and already correct.

Structuring it that way means a fourth mapping type added upstream is a new method, not an edit to a growing `if`. Doing #151 and this issue as separate ad-hoc branches would be the thing to avoid.

Note the shared piece: *both* paths want "the tessellation reference coordinates for this reference shape", and the Piola path wants them as a `QuadratureRule`. Worth extracting once, alongside the existing `tess_cache`, since a mixed grid needs one per reference shape either way.

## Concern: ansatz spaces with derivative dofs (Hermite)

`transfer_solution` already carries this caveat:

```julia
# NOTE this does not work for ansatz spaces where derivatives are mixed in (e.g. Hermite)
ncomps = length(Ferrite.reference_shape_value(ip_field, ξ0, 1))
```

Hermite-type spaces mix value dofs and derivative dofs, so the reference-to-physical relation is not a plain per-basis-function identity: the derivative dofs have to be rescaled by the Jacobian. Ferrite has no Hermite interpolation today (no hits in `Ferrite/src`), so this is latent rather than broken — but it is exactly the kind of space that would report `IdentityMapping` while *not* being safe for a tabulated, geometry-free path, and it would fail silently with plausible-looking numbers rather than erroring.

Two things to decide when implementing:

1. Whether dispatching on `mapping_type` is a strong enough guard, or whether the identity path should additionally assert something about the interpolation (e.g. that it is one of the known nodal families) and fall back to the `CellValues` path otherwise. Falling back is cheap insurance: the slow path is still correct.
2. Whether to keep the NOTE where it is or promote it to a documented limitation, since after this change there would be two places relying on the same assumption.

## Scope

- `src/dataset.jl`: `transfer_solution` / `_transfer_solution!`.
- Tests should pin the equivalence rather than just the values: run both paths and compare, for at least one mixed-refshape grid.
- Worth a benchmark in the suite if one is added later, since this is an interactive-path cost.

Contributor guide

No contributing guide indexed for this repository

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.