---
title: Curvilinear elements and metric terms
summary: Elements are curved with the same polynomial technology as the solution. Metric terms are computed so that a uniform flow stays exactly uniform on any mesh, and they are evaluated inside the differentiated code so that shape derivatives come for free.
section: theory
order: 4
level: 2
keywords: [metric identities, free-stream preservation, geometric conservation law, isoparametric, shape derivative]
---
<!-- gyre docs: expanded Markdown of theory/curvilinear-metrics.md (generated tables, verification numbers and chart data are written out). Base level of this page: L2 Practice; passages at another level are bracketed by 'begin level N' / 'end level N' comments and start with [LN Name]. Built from revision ceb1f4b0be. -->

High-order accuracy near curved walls needs curved elements: a straight-sided mesh limits the
whole solution to second order regardless of $p$, and its kinks generate spurious entropy. Each Gyre element carries $(p_g+1)^d$ geometry
nodes on an LGL lattice (degree $p_g$, default 4) that define its polynomial mapping
$\mathbf x(\boldsymbol\xi)$. Geometry nodes on shared faces are shared, so any smooth displacement
of the *global* node array keeps the mesh watertight — which is how surface design variables
enter.

<!-- begin level 2 (Practice) -->
**[L2 Practice] Free-stream preservation**

On a curved element the divergence of a constant flux is discretised as
$\sum_a \mathbf D_a\,(J\mathbf a^a)\cdot\mathbf F_\text{const}$. It vanishes only if the *discrete*
metric identities

$$
\sum_{a=1}^{d} \mathbf D_a\,(J a^a_b) = 0 \qquad (b = 1,\dots,d)
$$

hold. If they do not, a uniform flow generates spurious waves at every curved element. Gyre
satisfies the identities to round-off in 2-D and 3-D, including with viscous terms:

**Table: Uniform-flow residual on warped meshes** (Free-stream preservation on curved meshes)

| dim | elements | p | geometry degree | viscous | min J / max J | max \|Rw\| | surface closure |
|---|---|---|---|---|---|---|---|
| 2-D | 6x5 | 3 | 4 | no | 0.451047 | 1.3492e-15 | 3.0531e-16 |
| 2-D | 6x5 | 3 | 4 | yes | 0.451047 | 1.3368e-15 | 3.0531e-16 |
| 2-D | 4x4 | 5 | 4 | no | 0.453556 | 2.8560e-15 | 4.9960e-16 |
| 2-D | 4x4 | 5 | 4 | yes | 0.453556 | 2.9091e-15 | 4.9960e-16 |
| 3-D | 3x3x3 | 2 | 2 | no | 0.276159 | 2.6316e-16 | 3.4694e-17 |
| 3-D | 3x3x3 | 2 | 2 | yes | 0.276159 | 2.6316e-16 | 3.4694e-17 |
| 3-D | 3x3x3 | 3 | 3 | no | 0.417089 | 3.2526e-16 | 6.9823e-17 |
| 3-D | 3x3x3 | 3 | 3 | yes | 0.417089 | 3.2823e-16 | 6.9823e-17 |

Provenance: gyre 0.4.0; revision `39979fcf90` (commit `39979fcf902ad6116eb550b881b4278d4599c8e2`); 2026-11-12T10:22:33-0500; device cpu; quick suite. **Revision check: Same code as this build.** Produced by 39979fcf902ad6116eb550b881b4278d4599c8e2. Every commit between it and the commit this site was built from (ceb1f4b0bed41de53e586df10858079f5a2bcdd8) changes only documentation text, evidence records or bookkeeping: the code is the same. Full record: [free_stream.json](../data/verification/free_stream.json).

<!-- end level 2 -->

<!-- begin level 3 (Mechanics) -->
**[L3 Mechanics] How the metrics are computed**

The geometry is first interpolated from the geometry nodes to the solution nodes (degree $p$),
then differentiated with the same $\mathbf D$ used for the solution.

*Two dimensions.* The cross-product form is already divergence-free discretely:

$$
J\mathbf a^1 = (y_\eta,\,-x_\eta), \qquad J\mathbf a^2 = (-y_\xi,\; x_\xi), \qquad J = x_\xi y_\eta - x_\eta y_\xi .
$$

*Three dimensions.* Products such as $y_\eta z_\zeta$ are of degree $2p$ and their derivatives no
longer commute after interpolation, so the cross-product form violates the identities. Gyre uses
Kopriva's conservative curl form,

$$
J a^i_n = -\hat{\mathbf x}_i\cdot\nabla_\xi \times \left( I^p\!\left( X_l\,\nabla_\xi X_m \right) \right),
\qquad (n, m, l)\ \text{cyclic},
$$

where $I^p$ is nodal interpolation at degree $p$. Being the discrete curl of a polynomial, each
contravariant vector is discretely divergence-free by construction.

Outward face normals and surface Jacobians are the face values of the same arrays:
$\mathbf n\, s_J = \pm J\mathbf a^a$ on the two faces normal to $\xi_a$, which makes the normals of two
elements sharing a face exactly opposite — a requirement for conservation.

<!-- end level 3 -->

<!-- begin level 3 (Mechanics) -->
**[L3 Mechanics] Geometry is data, not configuration**

Metric terms are not precomputed constants: `Discretization.geometry(X)` is an ordinary JAX
function of the global node coordinates $\mathbf X$. The residual is $\mathbf R(\mathbf U, \mathbf X,
\theta)$ and the chain

$$
\text{surface displacement } \mathbf d \;\xrightarrow{\text{deformation}}\; \mathbf X(\mathbf d)
\;\xrightarrow{\text{metrics}}\; (J, J\mathbf a, \mathbf n, s_J) \;\xrightarrow{\text{operator}}\; \mathbf R,\, J
$$

is differentiated end-to-end by reverse-mode AD. There is no separate "mesh sensitivity" code to
keep consistent with the solver. The verification study on mesh deformation checks the
derivative of a geometric integral through this chain against finite differences.

**Verification status: Differentiable mesh deformation** (study `mesh_deformation`): **PASSED**, 9 of 9 checks passed, 2.45 s.

Surface displacements are propagated to all high-order geometry nodes by inverse-distance weighting. The study deforms a NACA 0012 O-mesh with a thickness bump of increasing amplitude and checks that elements stay valid, that wall nodes follow the prescription exactly while far-field nodes stay fixed, and that the derivative of a geometric quantity through deformation and metric evaluation matches finite differences.

| Check | Value | Expected | Result |
|---|---|---|---|
| scale 0.0: mesh valid (min J > 0) | 4.3328e-05 | > 0 | pass |
| scale 0.0: surface/far-field prescription defect | 0 | < 1e-13 | pass |
| scale 0.5: mesh valid (min J > 0) | 4.3327e-05 | > 0 | pass |
| scale 0.5: surface/far-field prescription defect | 0 | < 1e-13 | pass |
| scale 1.0: mesh valid (min J > 0) | 4.3327e-05 | > 0 | pass |
| scale 1.0: surface/far-field prescription defect | 0 | < 1e-13 | pass |
| scale 2.0: mesh valid (min J > 0) | 4.3325e-05 | > 0 | pass |
| scale 2.0: surface/far-field prescription defect | 0 | < 1e-13 | pass |
| d(weighted area)/d(surface displacement): AD vs FD | 1.0828e-08 | < 1e-7 | pass |

Provenance: gyre 0.4.0; revision `39979fcf90` (commit `39979fcf902ad6116eb550b881b4278d4599c8e2`); 2026-11-12T10:22:42-0500; device cpu; quick suite. **Revision check: Same code as this build.** Produced by 39979fcf902ad6116eb550b881b4278d4599c8e2. Every commit between it and the commit this site was built from (ceb1f4b0bed41de53e586df10858079f5a2bcdd8) changes only documentation text, evidence records or bookkeeping: the code is the same. Full record: [mesh_deformation.json](../data/verification/mesh_deformation.json).

<!-- end level 3 -->

<!-- begin level 4 (Theory) -->
**[L4 Theory] The finite-volume level $p = 0$**

The coarsest p-multigrid level uses one node per element. Differentiating a one-point
"polynomial" gives no metrics, so this level takes its geometry from the degree-$p_g$
description by exact reduction: the cell Jacobian is the mean $\sum w J / 2^d$, face normals are
the face integrals $\int \mathbf n\,\mathrm dS$ divided by the reference face area. For the
inviscid terms the resulting scheme is a first-order finite-volume method on the *curved* cell
whose faces close exactly (geometric conservation), verified in
`test_finite_volume_geometry_reduction`. The viscous terms keep only the face-normal two-point
difference at this level ([penalty parameter](viscous-br2.md#penalty-parameter)).

<!-- end level 4 -->

<!-- begin level 3 (Mechanics) -->
**[L3 Mechanics] Mesh generators**

All built-in generators evaluate a smooth mapping on the LGL lattice, so elements are curved to
degree $p_g$:

| generator | topology | notes |
|---|---|---|
| `airfoil_omesh` | O-mesh around a NACA 4-digit section or a spline fit of coordinate data | sharp trailing edge and leading edge at element vertices; wall-normal lines from a smoothed normal field; geometric wall-normal growth from a first cell of the requested height |
| `annulus` | O-mesh between two circles | cylinder flows |
| `bump_channel` | channel with a Gaussian bump | smooth-wall verification |
| `box` | Cartesian block, 2-D or 3-D | optional smooth warp that preserves periodicity |

A detail worth knowing: classical cosine clustering along the chord has zero parametric speed
at the trailing edge, which makes the Jacobian vanish at the trailing-edge vertex. The aerofoil
generator uses $x = \sin^2(\tfrac{c\pi}{2} v)/\sin^2(\tfrac{c\pi}{2})$ with $c < 1$, which keeps
the smooth $\sqrt{x}$ behaviour at the nose but a finite tangent at the trailing edge.

<!-- end level 3 -->

<!-- begin level 3 (Mechanics) -->
**[L3 Mechanics] Meshes from files**

A mesh read from a file (`mesh.generator: file`) becomes the same global node array and
element lattice as a generated one, so everything on this page applies unchanged. For a Gmsh
file the reader maps Gmsh's node ordering of each element onto the tensor-product lattice,
re-orients left-handed elements, and interpolates the equispaced Gmsh nodes to the
Legendre–Gauss–Lobatto distribution of the same degree (the identity up to degree 2). Lattice
nodes are numbered topologically (a vertex by its file node, an edge or face node by the
vertices of that edge or face and its position on it), so neighbouring elements share nodes
exactly and nodes that merely coincide in space, at a slit or a periodic pair, stay distinct.
Elements may meet in any relative orientation. A mesh is refused if the Jacobian determinant
is not positive at every geometry node. The
[imported-mesh study](../verification/imported-meshes.md) tests these steps.

<!-- end level 3 -->