Isoparametric mapping#

Every solid element in femorph-solver is isoparametric: the same shape functions \(N_i(\boldsymbol{\xi})\) interpolate both the geometry and the displacement field. This chapter develops the geometric mapping between reference and physical space, the resulting Jacobian transform, and the conditions under which the mapping is invertible.

Reference and physical coordinates#

Each element is the image of a fixed reference element \(\hat{\Omega} \subset \mathbb{R}^{d}\) under a smooth map

(1)#\[\mathbf{x}(\boldsymbol{\xi}) = \sum_{i = 1}^{n_\mathrm{nodes}} N_i(\boldsymbol{\xi})\, \mathbf{x}_i, \qquad \boldsymbol{\xi} \in \hat{\Omega},\]

with \(\mathbf{x}_i\) the global coordinates of nodal point \(i\) and \(N_i(\boldsymbol{\xi})\) the same shape function used for the displacement field

\[\mathbf{u}(\boldsymbol{\xi}) = \sum_{i = 1}^{n_\mathrm{nodes}} N_i(\boldsymbol{\xi})\, \mathbf{u}_i.\]

Reference elements in femorph-solver:

Kernel

Reference \(\hat{\Omega}\)

Shape function family

HEX8

\([-1, +1]^{3}\) (cube)

Trilinear Lagrange

HEX20

\([-1, +1]^{3}\) (cube)

Serendipity

TET10

Unit tetrahedron, volume coords \(L_0 + L_1 + L_2 + L_3 = 1\)

Quadratic in volume coordinates

WEDGE15

Triangle \(\times [-1, +1]\)

Quadratic triangle ⊗ Hermite line

PYR13

Collapsed cube (Duffy 1982)

Rational pyramid (Bedrosian 1992)

QUAD4_PLANE / QUAD4_SHELL

\([-1, +1]^{2}\) (square)

Bilinear Lagrange

BEAM2

\([-1, +1]\) (line)

Linear axial / Hermite cubic transverse

The Jacobian#

Differentiating (1) with respect to \(\boldsymbol{\xi}\) defines the Jacobian of the geometric map,

(2)#\[\mathbf{J}(\boldsymbol{\xi}) = \frac{\partial \mathbf{x}}{\partial \boldsymbol{\xi}} = \sum_{i = 1}^{n_\mathrm{nodes}} \frac{\partial N_i}{\partial \boldsymbol{\xi}}\, \mathbf{x}_i^{\!\top} \in \mathbb{R}^{d \times d}.\]

Two derived quantities matter throughout the FE assembly:

  • Determinant: \(\det \mathbf{J}\) is the local volumetric distortion factor. The volume integral \(\int_{\Omega_e} f\,\mathrm{d}V\) transforms to

    \[\int_{\Omega_e} f\,\mathrm{d}V = \int_{\hat\Omega} f(\mathbf{x}(\boldsymbol{\xi}))\, \det \mathbf{J}(\boldsymbol{\xi}) \,\mathrm{d}\boldsymbol{\xi}.\]

    Gauss quadrature on the reference cube / tet / triangle combines the per-rule weights on \(\hat\Omega\) with the per-quadrature-point \(\det \mathbf{J}\).

  • Inverse: physical-frame derivatives of the shape functions follow from the chain rule,

    (3)#\[\frac{\partial N_i}{\partial \mathbf{x}} = \mathbf{J}^{-\!\top}\, \frac{\partial N_i}{\partial \boldsymbol{\xi}},\]

    which is the row of the strain-displacement matrix \(\mathbf{B}\) corresponding to node \(i\).

Invertibility and element distortion#

The isoparametric map is well-posed at a quadrature point \(\boldsymbol{\xi}\) only if \(\det \mathbf{J}( \boldsymbol{\xi}) > 0\). A negative or zero Jacobian indicates either a “tangled” element (nodes ordered the wrong way) or one distorted past the convex limit — both produce non-physical strain fields and break the FE solution.

femorph-solver’s element kernels do not currently raise on near-singular Jacobians; downstream callers can use Model.element_quality() (TA-roadmap) to flag elements whose minimum \(\det \mathbf{J}\) falls below a chosen threshold.

A common quality measure is the Jacobian ratio

\[r_J = \frac{\min_{\boldsymbol{\xi} \in \mathrm{Gauss}} \det \mathbf{J}(\boldsymbol{\xi})} {\max_{\boldsymbol{\xi} \in \mathrm{Gauss}} \det \mathbf{J}(\boldsymbol{\xi})},\]

with \(r_J = 1\) for an undistorted element and \(r_J \to 0\) for a degenerate one. Bathe §5 and Cook §6.13 give detailed treatments and recommended thresholds (\(r_J > 0.05\) is the de-facto Ansys / Abaqus warning level).

How femorph-solver uses the Jacobian#

  • \(\mathbf{B}\) and \(\det \mathbf{J}\) are computed per Gauss point inside the element kernels’ ke / me / eel_batch calls. The Python-level kernels (HEX8 etc. in femorph_solver.elements) call into C++ batch routines through nanobind so the loop runs in native code; the equations are still those above.

  • Stress recovery (Model.eel / compute_nodal_stress()) evaluates \(\boldsymbol{\varepsilon} = \mathbf{B}\, \mathbf{u}_e\) at the element nodes rather than the Gauss points, so the per-node strain is directly available — see Cook §6.14.

References#

  • Zienkiewicz, O. C. and Taylor, R. L. (2013) The Finite Element Method: Its Basis and Fundamentals, 7th ed., §4 (mapping and integration), §5 (curved isoparametric elements).

  • Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002) Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, §6.6 (isoparametric formulation), §6.13 (element distortion).

  • Hughes, T. J. R. (2000) The Finite Element Method — Linear Static and Dynamic Finite Element Analysis, Dover, §3.6 (isoparametric mapping for the bilinear quad).

  • Bathe, K.-J. (2014) Finite Element Procedures, 2nd ed., §5 (3D isoparametric formulation, distortion measures).

  • Ergatoudis, J. G., Irons, B. M. and Zienkiewicz, O. C. (1968) “Curved isoparametric quadrilateral elements for finite element analysis,” International Journal of Solids and Structures 4 (1), 31–42 (original isoparametric serendipity treatment).