diff --git a/examples/adaptive/hessian_recovery.py b/examples/adaptive/hessian_recovery.py new file mode 100644 index 00000000..86b96f75 --- /dev/null +++ b/examples/adaptive/hessian_recovery.py @@ -0,0 +1,69 @@ +""" +Nodal Hessian recovery on a fedoo mesh +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Compute the recovered nodal gradient and Hessian of an analytic scalar +field on a 2D mesh using ``fedoo.recover_gradient`` and +``fedoo.recover_hessian``. The output is written to a vtk file so you can +inspect it in ParaView. + +Pair with ``mmgpy.metrics.create_metric_from_hessian`` (and +``fedoo.to_upper_diagonal`` to get mmg's row-major upper-triangular +packing) to drive anisotropic adaptive remeshing. +""" + +import numpy as np +import fedoo as fd + +############################################################################### +# Build a 2D mesh and an analytic scalar field with a known closed-form Hessian. +mesh = fd.mesh.rectangle_mesh(nx=51, ny=51, elm_type="tri3") +x, y = mesh.nodes.T + +# f(x, y) = sin(2*pi*x) * sin(2*pi*y) +# d2f/dx2 = -4*pi^2 * sin(2*pi*x) * sin(2*pi*y) +# d2f/dy2 = -4*pi^2 * sin(2*pi*x) * sin(2*pi*y) +# d2f/dxdy = 4*pi^2 * cos(2*pi*x) * cos(2*pi*y) +field = np.sin(2 * np.pi * x) * np.sin(2 * np.pi * y) + +############################################################################### +# Recover gradient and Hessian. Both calls reuse fedoo's cached +# GP-to-Node projection matrix and shape-derivative tables - everything is +# vectorized, no per-vertex Python loop. +g = fd.recover_gradient(mesh, field) # (n_nodes, 2) +H = fd.recover_hessian(mesh, field) # (n_nodes, 2, 2) + +# Compare against the analytic Hessian on interior nodes. +k = (2 * np.pi) ** 2 +H_xx_exact = -k * np.sin(2 * np.pi * x) * np.sin(2 * np.pi * y) +H_xy_exact = +k * np.cos(2 * np.pi * x) * np.cos(2 * np.pi * y) + +interior = (x > 0.1) & (x < 0.9) & (y > 0.1) & (y < 0.9) +err_xx = np.abs(H[interior, 0, 0] - H_xx_exact[interior]).max() +err_xy = np.abs(H[interior, 0, 1] - H_xy_exact[interior]).max() +print(f"max |H_xx - exact| (interior) = {err_xx:.3e}") +print(f"max |H_xy - exact| (interior) = {err_xy:.3e}") + +############################################################################### +# Save a vtk file with the recovered fields for visualisation in ParaView. +ds = fd.DataSet( + mesh, + data={ + "field": field, + "grad": g.T, # (2, n_nodes) for fedoo's vector layout + "H_xx": H[:, 0, 0], + "H_yy": H[:, 1, 1], + "H_xy": H[:, 0, 1], + }, + data_type="node", +) +ds.save("hessian_recovery.vtk") +print("Wrote hessian_recovery.vtk") + +############################################################################### +# To drive an mmg adaptive remesh, pack the Hessian into mmg's +# row-major upper-triangular layout and feed it to +# ``mmgpy.metrics.create_metric_from_hessian``:: +# +# metric = fd.to_upper_diagonal(H) # (n_nodes, 3) for 2D +# # mmgpy.metrics.create_metric_from_hessian(metric, target_error=1e-3, ...) diff --git a/fedoo/__init__.py b/fedoo/__init__.py index adb69d05..7a565713 100644 --- a/fedoo/__init__.py +++ b/fedoo/__init__.py @@ -24,6 +24,12 @@ WeakFormSum, read_data, ) +from .util.recovery import ( + recover_gradient, + recover_hessian, + to_upper_diagonal, + to_voigt, +) try: from .util.viewer import viewer diff --git a/fedoo/util/__init__.py b/fedoo/util/__init__.py index fbb8b573..c0f3880d 100644 --- a/fedoo/util/__init__.py +++ b/fedoo/util/__init__.py @@ -8,3 +8,12 @@ # from .voigt_tensors import StrainTensorList, StressTensorList # from .simple_plot import mesh_plot_2d, field_plot_2d # from .abaqus_inp import ReadINP + +from .recovery import recover_gradient, recover_hessian, to_voigt, to_upper_diagonal + +__all__ = [ + "recover_gradient", + "recover_hessian", + "to_voigt", + "to_upper_diagonal", +] diff --git a/fedoo/util/recovery.py b/fedoo/util/recovery.py new file mode 100644 index 00000000..fcb23e13 --- /dev/null +++ b/fedoo/util/recovery.py @@ -0,0 +1,212 @@ +"""Vectorized FE-aware nodal gradient and Hessian recovery on fedoo meshes. + +Useful as a fast input for metric-based adaptive remeshing (e.g. +``mmgpy.metrics.create_metric_from_hessian``). The recovery uses a double +Galerkin L2-projection: per-element gradient at Gauss points -> lumped +GP-to-Node averaging via fedoo's cached projection matrix -> repeat for the +Hessian. All work is a few einsums plus a handful of sparse matvecs; no +per-vertex Python loop. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import numpy as np + +if TYPE_CHECKING: + from fedoo.core.mesh import Mesh + + +def _physical_shape_derivatives(mesh: Mesh, n_elm_gp: int | None = None): + """Return (dNdx, proj, n_gp) for the requested integration order. + + dNdx[el, gp, n, i] is dN_n/dx_i at gauss point gp of element el. + proj is the cached sparse GP-to-Node lumped averaging matrix. + """ + if n_elm_gp is None: + from fedoo.lib_elements.element_list import get_default_n_gp + + n_elm_gp = get_default_n_gp(mesh.elm_type) + + # Triggers init_interpolation + compute_jacobian_with_inverse, populating + # inv_jacobian_matrix on the cached element interpolation object. + mesh._get_gaussian_quadrature_mat(n_elm_gp) + + elm = mesh._elm_interpolation[n_elm_gp] + inv_J = elm.inv_jacobian_matrix # (n_el, n_gp, ndim_phys, ndim_xi) + + if inv_J.shape[-1] != inv_J.shape[-2]: + raise NotImplementedError( + "recover_gradient/recover_hessian require a square Jacobian " + "(volumetric or planar elements). Got inv_J shape " + f"{inv_J.shape} on element type '{mesh.elm_type}'." + ) + + # Reuse the per-element-class cache populated in __init__ + # (e.g. ElementTriangle.__init__ at fedoo/lib_elements/triangle.py:15). + # dN_xi[g, x, n] = dN_n/dxi_x at gp g + dN_xi = np.stack(elm.shape_function_derivative_gp, axis=0) + + # dNdx[e, g, n, i] = sum_x dN_xi[g, x, n] * inv_J[e, g, i, x] + dNdx = np.einsum("gxn, egix -> egni", dN_xi, inv_J) + + proj = mesh._get_gausspoint2node_mat(n_elm_gp) + return dNdx, proj, n_elm_gp + + +def _project_gp_field(field_gp: np.ndarray, proj) -> np.ndarray: + """Project a GP field of shape (n_el, n_gp, ...) to nodes via proj.""" + n_el, n_gp = field_gp.shape[:2] + tail = field_gp.shape[2:] + # GP layout in proj: column index = el + gp*n_el (mesh.py:1138-1142), + # i.e. elements vary fastest. Transposing axes 0 and 1 then ravelling + # in C order produces exactly that layout. + flat = field_gp.transpose(1, 0, *range(2, field_gp.ndim)).reshape(n_el * n_gp, -1) + out = proj @ flat + return out.reshape(out.shape[0], *tail) + + +def _validate_scalar_field(mesh: Mesh, field: np.ndarray) -> np.ndarray: + field = np.asarray(field, dtype=np.float64) + if field.shape != (mesh.n_nodes,): + raise ValueError( + f"field must have shape (n_nodes={mesh.n_nodes},), got {field.shape}" + ) + return field + + +def _recover_gradient(field: np.ndarray, mesh: Mesh, dNdx, proj) -> np.ndarray: + f_elem = field[mesh.elements] + g_gp = np.einsum("en, egni -> egi", f_elem, dNdx) + return _project_gp_field(g_gp, proj) + + +def recover_gradient( + mesh: Mesh, + field: np.ndarray, + n_elm_gp: int | None = None, +) -> np.ndarray: + """Recover the nodal gradient of a scalar field on a fedoo mesh. + + Computes the gradient at every Gauss point via the FE shape function + derivatives, then averages back to nodes through fedoo's cached lumped + GP-to-Node projection (``mesh._get_gausspoint2node_mat``). + + Restricted to volumetric or planar elements (square Jacobian); raises + ``NotImplementedError`` for shell, beam, or surface-in-3D elements. + Boundary nodes use one-sided element averaging and are therefore less + accurate than interior nodes. + + Parameters + ---------- + mesh : fedoo.Mesh + The mesh on which the field is defined. + field : (n_nodes,) ndarray + Scalar nodal values. + n_elm_gp : int, optional + Number of Gauss points per element. Defaults to the element default. + + Returns + ------- + (n_nodes, ndim) ndarray + Recovered nodal gradient. + """ + field = _validate_scalar_field(mesh, field) + dNdx, proj, _ = _physical_shape_derivatives(mesh, n_elm_gp) + return _recover_gradient(field, mesh, dNdx, proj) + + +def recover_hessian( + mesh: Mesh, + field: np.ndarray, + n_elm_gp: int | None = None, +) -> np.ndarray: + """Recover the nodal Hessian of a scalar field via double L2-projection. + + Step 1 recovers the gradient as a continuous P1 nodal vector field. + Step 2 takes the gradient of each component and projects back to nodes, + giving a full (possibly non-symmetric) tensor; the result is then + symmetrized. + + Restricted to volumetric or planar elements (square Jacobian); raises + ``NotImplementedError`` for shell, beam, or surface-in-3D elements. + Boundary nodes use one-sided element averaging and are therefore less + accurate than interior nodes - acceptable when feeding an mmg metric, + which clips eigenvalues by ``hmin``/``hmax``. + + Pack the result with :func:`to_voigt` for the fedoo convention + (``[XX, YY, ZZ, XY, XZ, YZ]``) or :func:`to_upper_diagonal` for the + row-major upper triangle (``[m11, m12, m13, m22, m23, m33]``), the layout + consumed by mmg's ``MMG{2,3}D_Set_tensorSols``. + + Parameters + ---------- + mesh : fedoo.Mesh + field : (n_nodes,) ndarray + Scalar nodal values. + n_elm_gp : int, optional + Number of Gauss points per element. Defaults to the element default. + + Returns + ------- + (n_nodes, ndim, ndim) ndarray + Symmetric Hessian tensor at every node. + """ + field = _validate_scalar_field(mesh, field) + dNdx, proj, _ = _physical_shape_derivatives(mesh, n_elm_gp) + g_node = _recover_gradient(field, mesh, dNdx, proj) + + g_elem = g_node[mesh.elements] + # H_gp[e, g, j, i] = sum_n g_elem[e, n, j] * dNdx[e, g, n, i] + H_gp = np.einsum("enj, egni -> egji", g_elem, dNdx) + H = _project_gp_field(H_gp, proj) + + return 0.5 * (H + H.swapaxes(-1, -2)) + + +def to_voigt(T: np.ndarray) -> np.ndarray: + """Pack a symmetric tensor field to fedoo Voigt order. + + 3D ``(n, 3, 3)`` -> ``(6, n)`` in order ``[XX, YY, ZZ, XY, XZ, YZ]``. + 2D ``(n, 2, 2)`` -> ``(3, n)`` in order ``[XX, YY, XY]``. + + Matches the storage convention of + :class:`fedoo.util.voigt_tensors._SymmetricTensorList`. + """ + T = np.asarray(T) + if T.ndim != 3 or T.shape[-1] != T.shape[-2]: + raise ValueError(f"expected (n, d, d) tensor field, got {T.shape}") + d = T.shape[-1] + if d == 3: + return np.stack( + [T[:, 0, 0], T[:, 1, 1], T[:, 2, 2], T[:, 0, 1], T[:, 0, 2], T[:, 1, 2]], + axis=0, + ) + if d == 2: + return np.stack([T[:, 0, 0], T[:, 1, 1], T[:, 0, 1]], axis=0) + raise ValueError(f"unsupported tensor dimension {d}; expected 2 or 3") + + +def to_upper_diagonal(T: np.ndarray) -> np.ndarray: + """Pack a symmetric tensor field to its row-major upper triangle. + + 3D ``(n, 3, 3)`` -> ``(n, 6)`` in order ``[m11, m12, m13, m22, m23, m33]``. + 2D ``(n, 2, 2)`` -> ``(n, 3)`` in order ``[m11, m12, m22]``. + + This is the layout consumed by mmg's ``MMG3D_Set_tensorSols`` / + ``MMG2D_Set_tensorSols`` (and therefore by mmgpy), but the helper is + consumer-agnostic and applies to any symmetric tensor field. + """ + T = np.asarray(T) + if T.ndim != 3 or T.shape[-1] != T.shape[-2]: + raise ValueError(f"expected (n, d, d) tensor field, got {T.shape}") + d = T.shape[-1] + if d == 3: + return np.stack( + [T[:, 0, 0], T[:, 0, 1], T[:, 0, 2], T[:, 1, 1], T[:, 1, 2], T[:, 2, 2]], + axis=1, + ) + if d == 2: + return np.stack([T[:, 0, 0], T[:, 0, 1], T[:, 1, 1]], axis=1) + raise ValueError(f"unsupported tensor dimension {d}; expected 2 or 3") diff --git a/tests/test_recovery.py b/tests/test_recovery.py new file mode 100644 index 00000000..cb17ff1b --- /dev/null +++ b/tests/test_recovery.py @@ -0,0 +1,140 @@ +"""Tests for fedoo.util.recovery (nodal gradient and Hessian recovery).""" + +import numpy as np +import pytest + +import fedoo as fd + + +def _interior_mask(mesh, margin=0.15): + """Boolean mask of nodes that are at least `margin` (relative bbox) + away from the bounding box, i.e. away from the mesh boundary where the + lumped GP-to-Node averaging is one-sided. + """ + crd = mesh.nodes + lo = crd.min(axis=0) + hi = crd.max(axis=0) + span = hi - lo + eps = margin * span + return np.all((crd > lo + eps) & (crd < hi - eps), axis=1) + + +def _eval_field(mesh, f): + return f(*mesh.nodes.T) + + +@pytest.mark.parametrize( + "elm_type, nx, ny", + [("quad4", 21, 21), ("tri3", 21, 21), ("tri6", 21, 21)], +) +def test_recover_hessian_2d_quadratic(elm_type, nx, ny): + """Quadratic 2D field: H is constant; recovery must hit it on interior nodes.""" + mesh = fd.mesh.rectangle_mesh(nx=nx, ny=ny, elm_type=elm_type) + + # f(x,y) = x^2 + x*y + 2*y^2 => H = [[2, 1], [1, 4]] + field = _eval_field(mesh, lambda x, y: x * x + x * y + 2 * y * y) + H = fd.recover_hessian(mesh, field) + + assert H.shape == (mesh.n_nodes, 2, 2) + # Symmetric by construction (post-symmetrize step). + np.testing.assert_allclose(H[:, 0, 1], H[:, 1, 0], atol=1e-12) + + H_expected = np.array([[2.0, 1.0], [1.0, 4.0]]) + interior = _interior_mask(mesh) + np.testing.assert_allclose( + H[interior], + np.broadcast_to(H_expected, (interior.sum(), 2, 2)), + atol=1e-9, + ) + + +def test_recover_hessian_2d_linear_field_is_zero(): + """Linear field has zero Hessian everywhere (boundary too).""" + mesh = fd.mesh.rectangle_mesh(nx=11, ny=11, elm_type="tri3") + field = _eval_field(mesh, lambda x, y: 2.0 * x + 3.0 * y - 0.5) + H = fd.recover_hessian(mesh, field) + np.testing.assert_allclose(H, 0.0, atol=1e-10) + + +def test_recover_hessian_3d_quadratic_hex8(): + """Quadratic 3D field on hex8: H constant on interior nodes.""" + mesh = fd.mesh.box_mesh(nx=11, ny=11, nz=11, elm_type="hex8") + # f(x,y,z) = x^2 + 2*y^2 + 3*z^2 + x*y + # H = [[2, 1, 0], [1, 4, 0], [0, 0, 6]] + field = _eval_field(mesh, lambda x, y, z: x * x + 2 * y * y + 3 * z * z + x * y) + H = fd.recover_hessian(mesh, field) + + assert H.shape == (mesh.n_nodes, 3, 3) + np.testing.assert_allclose(H, H.swapaxes(-1, -2), atol=1e-12) + + H_expected = np.array([[2.0, 1.0, 0.0], [1.0, 4.0, 0.0], [0.0, 0.0, 6.0]]) + interior = _interior_mask(mesh) + np.testing.assert_allclose( + H[interior], + np.broadcast_to(H_expected, (interior.sum(), 3, 3)), + atol=1e-9, + ) + + +def test_recover_gradient_2d(): + """Gradient of f(x,y) = 3x + 2y is the constant vector [3, 2].""" + mesh = fd.mesh.rectangle_mesh(nx=11, ny=11, elm_type="quad4") + field = _eval_field(mesh, lambda x, y: 3.0 * x + 2.0 * y) + g = fd.recover_gradient(mesh, field) + + assert g.shape == (mesh.n_nodes, 2) + np.testing.assert_allclose(g, np.broadcast_to([3.0, 2.0], g.shape), atol=1e-10) + + +def test_to_upper_diagonal_3d_round_trip(): + """to_upper_diagonal must recover the (n, 6) packing in mmg's order.""" + rng = np.random.default_rng(0) + n = 7 + H = rng.normal(size=(n, 3, 3)) + H = 0.5 * (H + H.swapaxes(-1, -2)) + + packed = fd.to_upper_diagonal(H) + assert packed.shape == (n, 6) + # Order [m11, m12, m13, m22, m23, m33]. + np.testing.assert_array_equal(packed[:, 0], H[:, 0, 0]) + np.testing.assert_array_equal(packed[:, 1], H[:, 0, 1]) + np.testing.assert_array_equal(packed[:, 2], H[:, 0, 2]) + np.testing.assert_array_equal(packed[:, 3], H[:, 1, 1]) + np.testing.assert_array_equal(packed[:, 4], H[:, 1, 2]) + np.testing.assert_array_equal(packed[:, 5], H[:, 2, 2]) + + +def test_to_upper_diagonal_2d(): + rng = np.random.default_rng(1) + n = 5 + H = rng.normal(size=(n, 2, 2)) + H = 0.5 * (H + H.swapaxes(-1, -2)) + + packed = fd.to_upper_diagonal(H) + assert packed.shape == (n, 3) + np.testing.assert_array_equal(packed[:, 0], H[:, 0, 0]) + np.testing.assert_array_equal(packed[:, 1], H[:, 0, 1]) + np.testing.assert_array_equal(packed[:, 2], H[:, 1, 1]) + + +def test_to_voigt_3d(): + """to_voigt must produce (6, n) in fedoo Voigt order [XX,YY,ZZ,XY,XZ,YZ].""" + rng = np.random.default_rng(2) + n = 4 + H = rng.normal(size=(n, 3, 3)) + H = 0.5 * (H + H.swapaxes(-1, -2)) + + voigt = fd.to_voigt(H) + assert voigt.shape == (6, n) + np.testing.assert_array_equal(voigt[0], H[:, 0, 0]) + np.testing.assert_array_equal(voigt[1], H[:, 1, 1]) + np.testing.assert_array_equal(voigt[2], H[:, 2, 2]) + np.testing.assert_array_equal(voigt[3], H[:, 0, 1]) + np.testing.assert_array_equal(voigt[4], H[:, 0, 2]) + np.testing.assert_array_equal(voigt[5], H[:, 1, 2]) + + +def test_recover_hessian_validates_shape(): + mesh = fd.mesh.rectangle_mesh(nx=5, ny=5, elm_type="quad4") + with pytest.raises(ValueError, match="shape"): + fd.recover_hessian(mesh, np.zeros(mesh.n_nodes - 1))