Physics-Informed Motion (Neo-Hookean)

A MeshGraphNet trained on displacement alone has no opinion about whether the motion it predicts is motion tissue could undergo. An element may inflate, thin past what myocardium allows, or invert outright, and an L2 loss notices only to the extent that the vertices land in the wrong place. This module adds a neo-Hookean strain energy to that loss, which prices exactly those deformations, and exposes the Cauchy stress the same constitutive law implies.

The energy is

\[W = \frac{\mu}{2}(I_1 - 3) - \mu \ln J + \frac{\lambda}{2} (\ln J)^2\]

with \(I_1 = \operatorname{tr}(F^T F)\) and \(J = \det F\), evaluated from the deformation gradient of each tetrahedron. Spatial derivatives come from PhysicsNeMo Sym’s least-squares gradient reconstruction, its method for unstructured meshes.

The two appearances of \(J\) are not the same quantity in code. An inverted element makes \(\det F\) non-positive, and \(\ln J\) would then poison the whole loss with a NaN, so the logarithmic terms use \(\max(J, 10^{-6})\): an inversion costs a large finite penalty and stays trainable rather than ending the run. The incompressibility penalty \((J - 1)^2\) uses the raw determinant, which is signed and therefore already prices an inversion correctly.

That clamp is also why the inversion count matters. It keeps an inverted element finite, which is exactly what would let one pass unnoticed, so PhysicsInformedMotion.inverted_element_count reports the unclamped determinant’s non-positive entries and is the only signal that the predicted motion turned tissue inside out.

Requirements

The shape model must be volumetric: a strain energy needs volume elements, and the template’s own cells are what supply them. A surface model has no interior and no deformation gradient can be formed on it. Tutorial 16 builds such a model; see Tutorials.

physicsnemo.sym supplies PhysicsInformer and ships inside nvidia-physicsnemo, so no separate install is needed. It is imported lazily.

Training method

class physiotwin4d.TrainPhysicsNeMoPhysicsInformedMotion(log_level=20)

Bases: TrainPhysicsNeMoMGN

Train a MeshGraphNet whose loss also prices the tissue’s strain energy.

Identical to physiotwin4d.TrainPhysicsNeMoMGN apart from the loss, which becomes data + lambda_physics * (energy + incompressibility). The data term is scored on normalized targets, as before, while the physics term is scored in millimeters and kilopascals, so predictions are returned to physical units before the residual sees them.

Call set_mechanics(), set_elements() and set_reference_meshes() before training, unless lambda_physics is zero – which reproduces the data-only MeshGraphNet exactly and is the ablation a physics-informed run is measured against.

model_tag = 'physics_informed_motion'

__init__(log_level=20)

Initialize the physics-informed MeshGraphNet training method.

Parameters:

log_level (int | str) – Logging level. Default: logging.INFO.

set_mechanics(residual, lambda_physics=0.1)

Set the constitutive residual and how heavily it is weighted.

Parameters:

• residual (Optional[PhysicsInformedMotion]) – Configured PhysicsInformedMotion, or None to train on displacement alone.

• lambda_physics (float) – Weight of the physics term. 0.0 skips it entirely, reproducing the data-only MeshGraphNet, which is the ablation a physics-informed run is measured against.

Raises:

ValueError – If lambda_physics is negative, or is positive without a residual to evaluate.

Return type:

None

property inverted_element_count: int

Elements whose Jacobian went non-positive during training.

Non-zero means the network predicted motion that turns tissue inside out somewhere, which the clamp keeps trainable but does not make physical.

set_reference_meshes(reference_meshes)

Set each subject’s fitted reference mesh, keyed by subject id.

These are the undeformed configurations the strain energy is measured against – the same meshes the manifests name as fitted_reference_mesh, whose points the targets are defined at.

Return type:

None

set_elements(tets)

Set the (n_tet, 4) template elements the residual is summed over.

Return type:

None

train(train_dataset, val_dataset, stats, context, epochs, output_dir, template_mesh, template_coords, resume_from=None)

Bind each sample to its subject’s reference geometry, then train.

Return type:

tuple[Module, list[float], list[dict]]

checkpoint_fields()

Return architecture-specific fields to store in the checkpoint.

Return type:

dict

Constitutive law and geometry

Physics-informed motion training for PhysicsNeMo mesh-stage models.

physiotwin4d.TrainPhysicsNeMoMGN scores predicted motion against measured displacement alone, so nothing in its loss rules out motion no myocardium could undergo: locally inverted elements, non-physical dilation, a wall that thins past what tissue allows. This module adds a neo-Hookean strain energy to that loss, which prices those deformations, and exposes the Cauchy stress the same constitutive law implies so predicted motion can be rendered as stress.

The energy needs volume elements, so the shape model must be a tetrahedral one: the template’s own cells become the elements, which is what makes one set of element node ids valid for every subject and phase at once.

Reference configuration. The residual is measured against each subject’s fitted reference geometry, never the shared template. The stored targets are phase.points - fitted_reference.points, so the fitted reference is the undeformed state; measuring against the population mean instead would charge every subject a strain energy for merely being shaped unlike the mean, which confuses variation between subjects with deformation within one.

Two formulations of the same energy live here on purpose. The symbolic one (neo_hookean_pde()) is what PhysicsNeMo Sym differentiates and evaluates during training; the tensor one (NeoHookeanResidual) is what computes the Cauchy stress for export, which the symbolic path does not hand back. The tests cross-check them against each other.

PhysicsNeMo Sym is an optional dependency imported lazily, so import physiotwin4d works without it. It ships inside nvidia-physicsnemo.

physiotwin4d.train_physicsnemo_physics_informed_motion.tet_volumes(points, tets)

Return each tetrahedron’s volume and the volume lumped onto each node.

The nodal volumes are the quadrature weights the residual is averaged with, so a region contributes in proportion to the tissue it holds rather than to how finely it happens to be meshed.

Parameters:

• points (ndarray) – (n_points, 3) node positions.

• tets (ndarray) – (n_tet, 4) node ids per tetrahedron.

Return type:

tuple[ndarray, ndarray]

Returns:

(element_volumes, nodal_volumes), shaped (n_tet,) and (n_points,).

Raises:

ValueError – If any element is inverted or degenerate. Templates come from physiotwin4d.ContourTools.trim_tetrahedra_to_surface(), which holds every cell above a scaled Jacobian of 0.1, so a violation here means the template is broken rather than merely tight.

physiotwin4d.train_physicsnemo_physics_informed_motion.tet_edges(tets)

Return the unique undirected (n_edge, 2) edges of a tetrahedral mesh.

This is the stencil the least-squares gradient reconstruction fits over.

Return type:

ndarray

physiotwin4d.train_physicsnemo_physics_informed_motion.edge_matrix(points, tets)

Return the (n_tet, 3, 3) matrix whose columns are an element’s edges.

Return type:

Tensor

physiotwin4d.train_physicsnemo_physics_informed_motion.compute_deformation_gradient(reference_points, displacement, tets, reference_inverse=None)

Return the per-element deformation gradient F.

F = Ds @ Dm^-1, with the columns of Dm the reference edge vectors of an element and the columns of Ds its deformed ones – exact for the linear shape functions a tetrahedron carries.

Parameters:

• reference_points (Tensor) – (n_points, 3) undeformed node positions.

• displacement (Tensor) – (n_points, 3) predicted displacement, same units.

• tets (Tensor) – (n_tet, 4) node ids.

• reference_inverse (Optional[Tensor]) – Cached (n_tet, 3, 3) inverse of Dm. It depends only on the reference configuration, so passing it back in avoids re-inverting it once per phase of the same subject.

Return type:

Tensor

Returns:

(n_tet, 3, 3) deformation gradients.

class physiotwin4d.train_physicsnemo_physics_informed_motion.NeoHookeanResidual(mu_kpa=10.0, lambda_lame_kpa=100.0, log_level=20)

Compressible neo-Hookean constitutive law, evaluated on tensors.

The strain energy density is

W = (mu / 2) (I1 - 3) - mu ln(J) + (lambda / 2) ln(J)^2

with I1 = tr(F^T F) and J = det(F). Its Cauchy stress is sigma = (mu / J)(B - I) + (lambda ln(J) / J) I for B = F F^T.

Parameters:

• mu_kpa (float) – Shear modulus, in kilopascals.

• lambda_lame_kpa (float) – First Lame parameter, in kilopascals.

• log_level (int | str) – Logging level. Default: logging.INFO.

__init__(mu_kpa=10.0, lambda_lame_kpa=100.0, log_level=20)

Initialize the base class with logging configuration.

Parameters:

• class_name – Name for the class (used in log messages). If None, uses the class name. Default: None

• log_level (int | str) – Logging level. Can be an integer (logging.DEBUG, logging.INFO, logging.WARNING, logging.ERROR, logging.CRITICAL) or a string (‘DEBUG’, ‘INFO’, ‘WARNING’, ‘ERROR’, ‘CRITICAL’). Default: logging.INFO

• log_to_file – Optional file path to write logs to in addition to console output. Default: None

property inverted_element_count: int

Elements whose Jacobian went non-positive since this was constructed.

Non-zero means the deformation turns tissue inside out somewhere. The clamp in jacobian() keeps the energy finite so training can continue, which is also what would let an inversion pass unnoticed.

jacobian(deformation_gradient)

Return det(F) clamped away from zero, counting any inversion.

Return type:

Tensor

strain_energy(deformation_gradient)

Return the per-element strain energy density, in kilopascals.

Return type:

Tensor

incompressibility(deformation_gradient)

Return (J - 1)^2, the soft penalty on volume change.

Return type:

Tensor

cauchy_stress(deformation_gradient)

Return the (..., 3, 3) Cauchy stress tensor, in kilopascals.

Return type:

Tensor

physiotwin4d.train_physicsnemo_physics_informed_motion.neo_hookean_pde(mu_kpa, lambda_lame_kpa)

Return the neo-Hookean energy as a PhysicsNeMo Sym PDE.

The symbolic form is what PhysicsInformedMotion differentiates: it is written in terms of the displacement fields u, v and w, so PhysicsNeMo Sym supplies their spatial derivatives and assembles F = I + grad(u) itself.

Parameters:

• mu_kpa (float) – Shear modulus, in kilopascals.

• lambda_lame_kpa (float) – First Lame parameter, in kilopascals.

Return type:

Any

Returns:

A PDE exposing neo_hookean_energy and incompressibility.

class physiotwin4d.train_physicsnemo_physics_informed_motion.PhysicsInformedMotion(tets, n_points, mu_kpa=10.0, lambda_lame_kpa=100.0, device=None, log_level=20)

Evaluate the neo-Hookean residual of a predicted displacement field.

Spatial derivatives come from PhysicsNeMo Sym’s least-squares gradient reconstruction, the method built for unstructured meshes: it fits a gradient at every node over that node’s edge neighborhood. The neighborhood is fixed by the template’s topology, so the connectivity is built once here rather than per batch.

Parameters:

• tets (ndarray) – (n_tet, 4) template element node ids.

• n_points (int) – Node count of the template.

• mu_kpa (float) – Shear modulus, in kilopascals.

• lambda_lame_kpa (float) – First Lame parameter, in kilopascals.

• device (Optional[device]) – Device the residual is evaluated on. Defaults to the GPU when there is one, since this is the most expensive part of the loss and has to sit where the predictions are. Pass it explicitly when training somewhere other than the default device.

• log_level (int | str) – Logging level. Default: logging.INFO.

__init__(tets, n_points, mu_kpa=10.0, lambda_lame_kpa=100.0, device=None, log_level=20)

Initialize the base class with logging configuration.

Parameters:

• class_name – Name for the class (used in log messages). If None, uses the class name. Default: None

• log_level (int | str) – Logging level. Can be an integer (logging.DEBUG, logging.INFO, logging.WARNING, logging.ERROR, logging.CRITICAL) or a string (‘DEBUG’, ‘INFO’, ‘WARNING’, ‘ERROR’, ‘CRITICAL’). Default: logging.INFO

• log_to_file – Optional file path to write logs to in addition to console output. Default: None

property device: torch.device

Device this residual’s connectivity and symbolic graph were built on.

Fixed at construction: PhysicsInformer is given the device when its graph is compiled, so the residual cannot be moved afterwards. The trainer checks this against the device it predicts on.

Read from a tensor the residual actually owns, so the CUDA index is concrete even when the caller passed an index-less torch.device("cuda").

property inverted_element_count: int

Nodes whose Jacobian went non-positive since this was built.

Non-zero means the network predicted motion that turns tissue inside out somewhere. The energy clamps J to stay finite and trainable, so without this count an inversion would leave no trace.

__call__(reference_points, displacement_mm, nodal_volumes)

Return the volume-weighted (strain energy, incompressibility).

Parameters:

• reference_points (Tensor) – (n_points, 3) undeformed positions of this subject, in millimeters.

• displacement_mm (Tensor) – (n_points, 3) predicted displacement, in millimeters rather than in the normalized units the data loss is scored on.

• nodal_volumes (Tensor) – (n_points,) quadrature weights from tet_volumes().

Return type:

tuple[Tensor, Tensor]

Returns:

Two scalars, each a nodal-volume-weighted mean over the mesh.

class physiotwin4d.train_physicsnemo_physics_informed_motion.TrainPhysicsNeMoPhysicsInformedMotion(log_level=20)

Train a MeshGraphNet whose loss also prices the tissue’s strain energy.

Identical to physiotwin4d.TrainPhysicsNeMoMGN apart from the loss, which becomes data + lambda_physics * (energy + incompressibility). The data term is scored on normalized targets, as before, while the physics term is scored in millimeters and kilopascals, so predictions are returned to physical units before the residual sees them.

Call set_mechanics(), set_elements() and set_reference_meshes() before training, unless lambda_physics is zero – which reproduces the data-only MeshGraphNet exactly and is the ablation a physics-informed run is measured against.

model_tag: str = 'physics_informed_motion'

__init__(log_level=20)

Initialize the physics-informed MeshGraphNet training method.

Parameters:

log_level (int | str) – Logging level. Default: logging.INFO.

lambda_physics: float

set_mechanics(residual, lambda_physics=0.1)

Set the constitutive residual and how heavily it is weighted.

Parameters:

• residual (Optional[PhysicsInformedMotion]) – Configured PhysicsInformedMotion, or None to train on displacement alone.

• lambda_physics (float) – Weight of the physics term. 0.0 skips it entirely, reproducing the data-only MeshGraphNet, which is the ablation a physics-informed run is measured against.

Raises:

ValueError – If lambda_physics is negative, or is positive without a residual to evaluate.

Return type:

None

property inverted_element_count: int

Elements whose Jacobian went non-positive during training.

Non-zero means the network predicted motion that turns tissue inside out somewhere, which the clamp keeps trainable but does not make physical.

set_reference_meshes(reference_meshes)

Set each subject’s fitted reference mesh, keyed by subject id.

These are the undeformed configurations the strain energy is measured against – the same meshes the manifests name as fitted_reference_mesh, whose points the targets are defined at.

Return type:

None

set_elements(tets)

Set the (n_tet, 4) template elements the residual is summed over.

Return type:

None

train(train_dataset, val_dataset, stats, context, epochs, output_dir, template_mesh, template_coords, resume_from=None)

Bind each sample to its subject’s reference geometry, then train.

Return type:

tuple[Module, list[float], list[dict]]

checkpoint_fields()

Return architecture-specific fields to store in the checkpoint.

Return type:

dict

Notes

The reference configuration is the subject’s own fit, not the mean. The stored targets are phase.points - fitted_reference.points, so the fitted reference is the undeformed state. Measuring the residual against the population mean instead would charge every subject a strain energy for merely being shaped unlike the mean, confusing variation between subjects with deformation within one.

Two formulations of one law. The symbolic energy (neo_hookean_pde()) is what PhysicsNeMo Sym differentiates during training; the tensor one (NeoHookeanResidual) computes the Cauchy stress for export, which the symbolic path does not hand back. tests/test_physics_informed_motion.py cross-checks them against each other on the same field.

Loss scale. The data term is scored on normalized displacement and the physics term in millimeters and kilopascals, so lambda_physics is a value to sweep rather than one to trust. The two components are accumulated separately so they can be reported apart.

See Also

• Training a Mesh-Stage Model

• Scoring a Mesh-Stage Model Against Images