Add nodal gradient/Hessian recovery utilities

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, ...)

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

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",
+]

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")