[Preview] Hyperelasticity via Input-Convex Neural Networks (ICNN) (PyTorch)

[Preview] Hyperelasticity via Input-Convex Neural Networks (ICNN) (PyTorch)#

This tutorial demonstrates how to define a complex hyperelastic constitutive model using PyTorch automatic differentiation (AD) and integrate it with FEniCSx. In solid mechanics, standard models (like Neo-Hookean) have simple analytical derivatives. However, data-driven constitutive equations (like neural network surrogates) require automatic differentiation to evaluate stresses and tangent operators (Jacobians) at the quadrature points of the finite element mesh.

Here, we stick to the Isihara hyperelasticity model [Isihara et al., 1951]. The tutorial is based on the work of Thakolkaran et al. [2022] on unsupervised deep-learning hyperelasticity without stress data (NN-EUCLID): EUCLID-code/EUCLID-hyperelasticity-NN.

To get familiarized with hyperelasticity problems, we advise to take a look first at the following basic tutorials:

Numerical Tours of Computational Mechanics with FEniCSx by Jérémy Bleyer.
The FEniCSx tutorial by Jørgen S. Dokken and Garth N. Wells.

Problem Formulation#

We solve a hyperelasticity boundary value problem on a 2D square domain \(\Omega\) containing elliptic holes, subjected to displacement-controlled tensile loading.

Let \(V\) be the functional space of admissible displacement fields. Under the assumption of quasi-static displacement-controlled loading and in the absence of body forces, the weak formulation of the equilibrium equation is expressed as follows.

Find \(\mathbf{u} \in V\) satisfying the Dirichlet boundary conditions such that:

\[ F(\mathbf{u}; \mathbf{v}) = \int\limits_\Omega \mathbf{P}(\mathbf{F}) : \nabla\mathbf{v} \, \mathrm{d}\mathbf{x} = 0, \quad \forall \mathbf{v} \in V_0 \]

where:

\(\mathbf{F} = \mathbf{I} + \nabla\mathbf{u}\) is the deformation gradient.
\(\mathbf{P}(\mathbf{F}) = \frac{\partial W}{\partial \mathbf{F}}\) is the first Piola-Kirchhoff stress tensor derived from the strain energy density \(W(\mathbf{F})\).
\(\mathbf{v} \in V_0\) is a test function vanishing on the Dirichlet boundaries.

Implementation#

Preamble#

from mpi4py import MPI

from utilities import build_square_with_elliptic_holes

import basix
import dolfinx.plot
import ufl
from dolfinx import default_scalar_type, fem, mesh
from dolfinx.fem.petsc import NonlinearProblem

# Define external function for Piola-Kirchhoff stress and its derivative
from dolfinx_external_operator import (
    FEMExternalOperator,
    evaluate_external_operators,
    evaluate_operands,
    replace_external_operators,
)
from dolfinx_external_operator.petsc import assemble_residual_with_callback

# isort: split

from functools import partial

import matplotlib.pyplot as plt
import numpy as np
import pyvista
import torch
import torch.func

Geometry and boundary conditions#

We solve a 2D tension problem of a square specimen with elliptic holes. The boundary conditions consist of zero displacement at the bottom, and a prescribed vertical displacement at the top. We use quadratic Lagrange elements for the displacement field \(\mathbf{u}\).

L = 1.0  # Length
W = 1.0  # Height

domain, facet_tags, facet_tags_labels = build_square_with_elliptic_holes(L=L, lc=0.02)
V = fem.functionspace(domain, ("Lagrange", 2, (domain.geometry.dim,)))


def bottom(x):
    return np.isclose(x[1], 0.0)


def top(x):
    return np.isclose(x[1], W)


fdim = domain.topology.dim - 1
bottom_facets = mesh.locate_entities_boundary(domain, fdim, bottom)
top_facets = mesh.locate_entities_boundary(domain, fdim, top)

marked_facets = np.hstack([bottom_facets, top_facets])
marked_values = np.hstack([np.full_like(bottom_facets, 1), np.full_like(top_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

bottom_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
top_dofs_uy = fem.locate_dofs_topological(V.sub(1), facet_tag.dim, facet_tag.find(2))
top_dofs_ux = fem.locate_dofs_topological(V.sub(0), facet_tag.dim, facet_tag.find(2))
u_D_bottom = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)
u_D_top = fem.Constant(domain, 0.0)

bcs_u = [
    fem.dirichletbc(u_D_top, top_dofs_uy, V.sub(1)),
    fem.dirichletbc(fem.Constant(domain, 0.0), top_dofs_ux, V.sub(0)),
    fem.dirichletbc(u_D_bottom, bottom_dofs, V),
]

Constitutive Formulation#

We model the hyperelastic response using an Input-Convex Neural Network (ICNN) designed to output a strain energy density \(W(\mathbf{F})\) that is convex with respect to the right Cauchy-Green strain invariants. Objectivity is satisfied by expressing \(W\) as a function of the right Cauchy-Green deformation tensor \(\mathbf{C} = \mathbf{F}^T\mathbf{F}\).

Following [Thakolkaran et al., 2022], the strain energy density is formulated as:

\[ W(\mathbf{F}) = W_{\mathbf{Q}, \mathbf{\mathcal{A}}}^{\text{NN}}(\mathbf{E}(\mathbf{F})) + W^0 + \mathbf{H}:\mathbf{E}, \]

where:

\(\mathbf{E} = \frac{1}{2}(\mathbf{C} - \mathbf{I})\) is the Green-Lagrange strain tensor, representing the deformation.
\(W_{\mathbf{Q}, \mathbf{\mathcal{A}}}^{\text{NN}}(\mathbf{E})\) is the neural network mapping the strain tensor \(\mathbf{E}\) (via its invariants) to a scalar energy.
\(W^0\) is a scalar correction offsetting the energy density to zero in the reference configuration:

\[ W^0 = -\left. W_{\mathbf{Q}, \mathbf{\mathcal{A}}}^{\text{NN}}(\mathbf{E}(\mathbf{F})) \right|_{\mathbf{F} = \mathbf{I}} \]

\(\mathbf{H}\) is a stress correction tensor ensuring that the stress vanishes in the undeformed state (\(\mathbf{F}=\mathbf{I}\)):

\[ \mathbf{H} = -\left. \frac{\partial W_{\mathbf{Q}, \mathbf{\mathcal{A}}}^{\text{NN}}(\mathbf{E})}{\partial \mathbf{F}} \right|_{\mathbf{F}=\mathbf{I}} \]

The first Piola-Kirchhoff stress \(\mathbf{P}\) is computed as:

\[ \mathbf{P}(\mathbf{F}) = \frac{\partial W(\mathbf{F})}{\partial \mathbf{F}} = \mathbf{P}^\text{NN} + \mathbf{F}\mathbf{H}, \]

where \(\mathbf{P}^\text{NN} := \frac{\partial W^{\text{NN}}}{\partial \mathbf{F}}\) and the tangent modulus is:

\[ \mathbb{C}_{ijkl} = \frac{\partial P_{ij}(\mathbf{F})}{\partial F_{kl}} = \mathbb{C}^\text{NN}_{ijkl} + \delta_{ik} H_{lj} ,\]

where \(\mathbb{C}^\text{NN}_{ijkl} := \frac{\partial^2 W^{\text{NN}}}{\partial F_{ij} \partial F_{kl}}\).

We will load a pre-trained network modeling the Isihara material behaviour to evaluate the strain energy density \(W^{\text{NN}}\). Then, using PyTorch automatic differentiation (AD) tools, we will compute the values of the first Piola-Kirchhoff stress \(\mathbf{P}\) and the tangent \(\mathbb{C}\). The programs resulting from AD will then be used to define the external operators associated with \(\mathbf{P}(\mathbf{F})\) and \(\mathbb{C}\) respectively.

Input-Convex Neural Network (ICNN) in PyTorch#

We preserve the same structure of convexLinear and ICNN as defined in the original implementation using custom layers that enforce positive weights in the hidden layers to maintain input-convexity.

Here below, we represent the \(2 \times 2\) tensors \(\mathbf{F}, \mathbf{P}\), and \(\mathbf{H}\) in flattened vector forms, e.g.

\[ \mathbf{P} = [P_{11}, P_{12}, P_{21}, P_{22}]^T .\]

class convexLinear(torch.nn.Module):
    """Custom linear layer with positive weights and no bias"""

    def __init__(
        self,
        size_in,
        size_out,
    ):
        super().__init__()
        self.size_in, self.size_out = size_in, size_out
        weights = torch.Tensor(size_out, size_in)
        self.weights = torch.nn.Parameter(weights)

        # initialize weights
        torch.nn.init.kaiming_uniform_(self.weights, a=np.sqrt(5))

    def forward(self, x):
        w_times_x = torch.mm(x, torch.nn.functional.softplus(self.weights.t()))
        return w_times_x


class ICNN(torch.nn.Module):
    def __init__(self, n_input, n_hidden, n_output, dropout):
        super().__init__()
        # Create Module dicts for the hidden and skip-connection layers
        self.layers = torch.nn.ModuleDict()
        self.skip_layers = torch.nn.ModuleDict()
        self.depth = len(n_hidden)
        self.dropout = dropout[0]
        self.p_dropout = dropout[1]

        self.layers[str(0)] = torch.nn.Linear(n_input, n_hidden[0]).float()
        # Create create NN with number of elements in n_hidden as depth
        for i in range(1, self.depth):
            self.layers[str(i)] = convexLinear(n_hidden[i - 1], n_hidden[i]).float()
            self.skip_layers[str(i)] = torch.nn.Linear(n_input, n_hidden[i]).float()

        self.layers[str(self.depth)] = convexLinear(n_hidden[self.depth - 1], n_output).float()
        self.skip_layers[str(self.depth)] = convexLinear(n_input, n_output).float()

    def forward(self, x):
        # Get F components
        F11 = x[:, 0:1]
        F12 = x[:, 1:2]
        F21 = x[:, 2:3]
        F22 = x[:, 3:4]

        # Compute right Cauchy green strain Tensor
        C11 = F11**2 + F21**2
        C12 = F11 * F12 + F21 * F22
        C21 = F11 * F12 + F21 * F22
        C22 = F12**2 + F22**2

        # Compute strain invariants
        I1 = C11 + C22 + 1.0
        I2 = C11 + C22 - C12 * C21 + C11 * C22
        I3 = C11 * C22 - C12 * C21

        # Apply transformation to invariants
        K1 = I1 * torch.pow(I3, -1 / 3) - 3.0
        K2 = I2 * torch.pow(I3, -2 / 3) - 3.0
        J = torch.sqrt(I3)
        K3 = (J - 1) ** 2

        # Concatenate feature
        x_input = torch.cat((K1, K2, K3), 1).float()

        z = x_input.clone()
        z = self.layers[str(0)](z)
        for layer in range(1, self.depth):
            skip = self.skip_layers[str(layer)](x_input)
            z = self.layers[str(layer)](z)
            z += skip
            z = torch.nn.functional.softplus(z)
            z = 1 / 12.0 * torch.square(z)
            if self.training:
                if self.dropout:
                    z = torch.nn.functional.dropout(z, p=self.p_dropout)
        y = self.layers[str(self.depth)](z) + self.skip_layers[str(self.depth)](x_input)
        return y


n_input = 3
n_output = 1
n_hidden = [64, 64, 64]
dropout = [True, 0.2]
model = ICNN(n_input=n_input, n_hidden=n_hidden, n_output=n_output, dropout=dropout)

Once the structure of the ICNN model is defined we can load a pretrained one. In particular, we explore the Isihara model [Isihara et al., 1951] with high noise.

model.load_state_dict(torch.load("Isihara_noise=high.pth"))
model.eval()

ICNN(
  (layers): ModuleDict(
    (0): Linear(in_features=3, out_features=64, bias=True)
    (1): convexLinear()
    (2): convexLinear()
    (3): convexLinear()
  )
  (skip_layers): ModuleDict(
    (1): Linear(in_features=3, out_features=64, bias=True)
    (2): Linear(in_features=3, out_features=64, bias=True)
    (3): convexLinear()
  )
)

Train your own ICNN model!

Here we are not studying how to train the ICNN model. Instead, we load the pretrained one Isihara_noise=high.pth. If you wish to explore other hyperelastic models we encourage you to follow the instructions in the [original repository of EUCLID] (EUCLID-code/EUCLID-hyperelasticity-NN), which we outline here

In the folder with demo_hyperelasticity.py, clone EUCLID-hyperelasticity-NN

git clone https://github.com/EUCLID-code/EUCLID-hyperelasticity-NN

Install pandas and [other dependencies] (EUCLID-code/EUCLID-hyperelasticity-NN), if needed

pip install pandas

Train a new model

cd drivers/
python main.py Isihara high

Note: Although we rely here on the CPU-based PyTorch installation, it is motivated by keeping the external operators demos light. For training of new models, we suggest to install the normal GPU-based PyTorch package.

Then in the code above try

model.load_state_dict(torch.load("Isihara_noise=high.pth"))
model.eval()

Note: there is a [bug] (EUCLID-code/EUCLID-hyperelasticity-NN#2) related to NumPy>=2.0. If the error persists, try [this fork] (a-latyshev/EUCLID-hyperelasticity-NN) with a fix.

PyTorch-based Tangent and Stress Evaluation#

We first precompute the constant stress correction tensor \(\mathbf{H}\) and the constant tangent correction \(\mathbb{C}^\text{cor}_{ijkl} = \frac{\partial \mathbf{P}^{\text{cor}}_{ij}}{\partial \mathbf{F}_{kl}} = \delta_{ik} H_{lj}\), since they do not change throughout the loading.

F_0 = torch.zeros((1, 4))
F_0[:, 0] = 1
F_0[:, 3] = 1

F_0.requires_grad = True
W_NN_0 = model(F_0)
P_NN_0 = torch.autograd.grad(W_NN_0, F_0, torch.ones_like(W_NN_0), create_graph=True)[0].squeeze()

# Precompute correction matrix H for stress correction
H_flat = -P_NN_0.detach()
H = torch.tensor(
    [
        [H_flat[0], H_flat[1], 0, 0],
        [H_flat[2], H_flat[3], 0, 0],
        [0, 0, H_flat[0], H_flat[1]],
        [0, 0, H_flat[2], H_flat[3]],
    ],
    dtype=P_NN_0.dtype,
    device=P_NN_0.device,
)

Secondly, we define the behaviour of the constitutive model locally, at a single Gauss point, and then globally by leveraging the modern (JAX-inspired) PyTorch features such that

automatic differentiation (torch.func.jacfwd) to compute \(\mathbb{C}^\text{NN}\),
vectorization (torch.func.vmap) to extrapolate the local behaviour onto all Gauss points,
and JIT compilation (torch.compile) to speed up computations.

Overall, as you will see later, the workflow is very similar to the previous tutorial using JAX.

Install PyTorch 2.0 or later

To be able to use torch.func.jacfwd, torch.func.vmap and torch.compile, make sure that you use torch=>2.0.

Currently, there is an [issue] (pytorch/pytorch#160508) with combining all three together. Compile just model instead.

Since PyTorch compiler Inductor currently has tracing limitations when directly compiling nested higher-order functional AD operators (vmap / jacfwd), we compile the underlying neural network model itself instead. This avoids compiler-level limitations while capturing the bulk of the compilation speedup, since the neural network evaluations dominate the computational cost.

# JIT-compile the model to optimize forward and backward passes
model = torch.compile(model)

Now we apply AD and vectorize the functions that compute the stress state. The final dP_dF_impl computes both the stress state and the tangent at all Gauss points at once.

# Define vectorized stress and tangent computation using torch.func


def compute_stress_local(F_single):
    # F_single: shape (4,)
    # Define local energy function for the neural network
    def energy_fn(f):
        return model(f.unsqueeze(0)).squeeze()

    # Stress from neural network (first derivative)
    P_NN = torch.func.grad(energy_fn)(F_single)

    # Stress correction P_cor = F_single @ H (using flat matrix multiplication for efficiency)
    P_cor = F_single @ H.to(F_single)

    P = P_NN + P_cor
    # Return stress as primary output and as auxiliary data
    return P, P


# Vectorize the tangent evaluation over the batch dimension using jacfwd
# with has_aux=True returns a tuple of (tangent, auxiliary_stress)
vectorized_stress_and_tangent = torch.func.vmap(torch.func.jacfwd(compute_stress_local, has_aux=True), in_dims=(0,))


def dP_dF_impl(Fvals):
    F = torch.from_numpy(np.ascontiguousarray(Fvals)).reshape(-1, 4)  # NumPy-compatible without copy
    # Evaluate stress and tangent simultaneously in a single vectorized pass
    dP, P = vectorized_stress_and_tangent(F)

    return dP.reshape(-1).detach().numpy(), P.reshape(-1).detach().numpy()

Similarly to the previous tutorials with plasticity, we do not implement explicitly the evaluation of the stress. Instead, we obtain its values during the evaluation of its derivative and then update the values of the operator in the main Newton loop. This happens because automatic differentiation may evaluate the value of a function while computing its derivative.

FEniCSx Integration of PyTorch model via External Operators#

The stress state has one operand, which is the strain tensor \(\mathbf{F}\) (defined as gradU) that depends on the displacement field \(\mathbf{u}\). Here we define both of them and the quadrature space Q, where the external operator will live.

u = fem.Function(V)
v = ufl.TestFunction(V)
d = len(u)
gradU = ufl.variable(ufl.Identity(d) + ufl.grad(u))  # \mathb{F} tensor

# Create a quadrature element and function space for tensor-valued P
quadrature_degree = 2
cell_name = domain.topology.cell_name()
Qe = basix.ufl.quadrature_element(cell_name, value_shape=(d, d), degree=quadrature_degree, scheme="default")
Q = fem.functionspace(domain, Qe)

Finally, we wrap the PyTorch execution code into a Python callable, P_external, and define a FEMExternalOperator that can be used to define the variational formulation.

def P_external(derivatives):
    if derivatives == (1,):
        return dP_dF_impl
    else:
        raise NotImplementedError(f"No external function is defined for the requested derivative {derivatives}.")


P = FEMExternalOperator(gradU, function_space=Q, external_function=P_external)

Note

We don’t need to cover the case if derivatives == (0,) for just calling evaluation of P because, as it was previously mentioned, thanks to AD we can compute the values of P while computing its derivative. Yet, we still need to properly update the values of P, which will be done in constitutive_update, a function defined later below while defining a nonlinear solver.

Variational Forms#

We define the weak form using UFL as usual, compute its directional derivative, and replace it with the replace_external_operators function to let DOLFINx assembling the forms.

metadata = {"quadrature_degree": 2}
dx = ufl.Measure("dx", domain=domain, metadata=metadata)

u_hat = ufl.TrialFunction(V)
F = ufl.inner(ufl.grad(v), P) * dx
J = ufl.derivative(F, u, u_hat)

F_replaced, F_external_operators = replace_external_operators(F)
J_replaced, J_external_operators = replace_external_operators(J)

F_form = fem.form(F_replaced)
J_form = fem.form(J_replaced)

Solving the Problem#

We set up the PETSc SNES solver using a Newton solver, wrapping the constitutive update in a callback function. This function affects the performance the most, since it will be called on every single iteration of the Newton solver to update the values of \(\mathbf{P}\) and \(\mathbb{C}\).

def constitutive_update(
    F_external_operators: list[FEMExternalOperator],
    J_external_operators: list[FEMExternalOperator],
):
    """Update the constitutive model by evaluating the external operators."""
    evaluated_operands = evaluate_operands(F_external_operators)
    # `dP_dF_impl` will be called here
    ((_, P_new),) = evaluate_external_operators(J_external_operators, evaluated_operands)
    # manual update the values of the external operator
    P.ref_coefficient.x.array[:] = P_new


petsc_options = {
    "snes_type": "vinewtonrsls",
    "snes_linesearch_type": "basic",
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
    "snes_atol": 1.0e-7,
    "snes_rtol": 1.0e-7,
    "snes_max_it": 50,
    "snes_monitor": "",
}

problem = NonlinearProblem(
    F_replaced, u, J=J_replaced, bcs=bcs_u, petsc_options_prefix="demo_hyperelasticity_", petsc_options=petsc_options
)

assemble_residual_with_callback_ = partial(
    assemble_residual_with_callback,
    problem.u,
    problem.F,
    problem.J,
    bcs_u,
    constitutive_update,  # external callback with respect to SNES
    [F_external_operators, J_external_operators],  # input arguments of the callback
)
problem.solver.setFunction(assemble_residual_with_callback_, problem.b)

We incrementally apply a displacement-controlled tensile load by moving the top boundary.

# Apply a tensile load by incrementally increasing displacement on the top edge
n_steps = 100
max_displacement = 0.5
u.name = "displacement"
u.x.array[:] = 0
for step in range(1, n_steps + 1):
    u_D_top.value = step * max_displacement / n_steps  # moving the top boundary
    num_its, converged = problem.solve()
    assert converged, f"Newton solver did not converge at step {step}"
    u.x.scatter_forward()
    if domain.comm.rank == 0:
        print(f"Step {step}: Displacement {u_D_top.value:.2f}, Newton its: {num_its}")

  0 SNES Function norm 5.749637490610e-01

  1 SNES Function norm 8.708253721394e-05

  2 SNES Function norm 1.266275510398e-08
Step 1: Displacement 0.01, Newton its: displacement[0]
  0 SNES Function norm 5.726081790327e-01

  1 SNES Function norm 8.061158348093e-05

  2 SNES Function norm 1.122751598604e-08
Step 2: Displacement 0.01, Newton its: displacement[0]
  0 SNES Function norm 5.703076949711e-01

  1 SNES Function norm 7.469177952053e-05

  2 SNES Function norm 1.034738254651e-08
Step 3: Displacement 0.01, Newton its: displacement[0]
  0 SNES Function norm 5.680606125732e-01

  1 SNES Function norm 6.934848232740e-05

  2 SNES Function norm 9.856346482217e-09
Step 4: Displacement 0.02, Newton its: displacement[0]
  0 SNES Function norm 5.658652875654e-01

  1 SNES Function norm 6.458788265383e-05

  2 SNES Function norm 9.737752479412e-09
Step 5: Displacement 0.03, Newton its: displacement[0]
  0 SNES Function norm 5.637201217829e-01

  1 SNES Function norm 6.039587675794e-05

  2 SNES Function norm 9.803453728520e-09
Step 6: Displacement 0.03, Newton its: displacement[0]
  0 SNES Function norm 5.616235526180e-01

  1 SNES Function norm 5.673959224695e-05

  2 SNES Function norm 9.969330723327e-09
Step 7: Displacement 0.04, Newton its: displacement[0]
  0 SNES Function norm 5.595740714625e-01

  1 SNES Function norm 5.357227324486e-05

  2 SNES Function norm 1.020242615242e-08
Step 8: Displacement 0.04, Newton its: displacement[0]
  0 SNES Function norm 5.575702285808e-01

  1 SNES Function norm 5.083682568909e-05

  2 SNES Function norm 1.115731771911e-08
Step 9: Displacement 0.04, Newton its: displacement[0]
  0 SNES Function norm 5.556106256545e-01

  1 SNES Function norm 4.847340542741e-05

  2 SNES Function norm 1.186594950781e-08
Step 10: Displacement 0.05, Newton its: displacement[0]
  0 SNES Function norm 5.536939254464e-01

  1 SNES Function norm 4.642173435271e-05

  2 SNES Function norm 1.247642144818e-08
Step 11: Displacement 0.06, Newton its: displacement[0]
  0 SNES Function norm 5.518188467564e-01

  1 SNES Function norm 4.462676766693e-05

  2 SNES Function norm 1.358269996852e-08
Step 12: Displacement 0.06, Newton its: displacement[0]
  0 SNES Function norm 5.499841692034e-01

  1 SNES Function norm 4.304142262640e-05

  2 SNES Function norm 1.429572131714e-08
Step 13: Displacement 0.07, Newton its: displacement[0]
  0 SNES Function norm 5.481887306740e-01

  1 SNES Function norm 4.162442938956e-05

  2 SNES Function norm 1.506695146781e-08
Step 14: Displacement 0.07, Newton its: displacement[0]
  0 SNES Function norm 5.464314201560e-01

  1 SNES Function norm 4.034325270596e-05

  2 SNES Function norm 1.586988911980e-08
Step 15: Displacement 0.07, Newton its: displacement[0]
  0 SNES Function norm 5.447111829316e-01

  1 SNES Function norm 3.917262918224e-05

  2 SNES Function norm 1.703358017868e-08
Step 16: Displacement 0.08, Newton its: displacement[0]
  0 SNES Function norm 5.430270129764e-01

  1 SNES Function norm 3.809236380355e-05

  2 SNES Function norm 1.760074908600e-08
Step 17: Displacement 0.09, Newton its: displacement[0]
  0 SNES Function norm 5.413779603026e-01

  1 SNES Function norm 3.708844863445e-05

  2 SNES Function norm 1.873519642411e-08
Step 18: Displacement 0.09, Newton its: displacement[0]
  0 SNES Function norm 5.397631128382e-01

  1 SNES Function norm 3.615028088598e-05

  2 SNES Function norm 1.985845485204e-08
Step 19: Displacement 0.10, Newton its: displacement[0]
  0 SNES Function norm 5.381816132941e-01

  1 SNES Function norm 3.526834334837e-05

  2 SNES Function norm 2.060868166882e-08
Step 20: Displacement 0.10, Newton its: displacement[0]
  0 SNES Function norm 5.366326399336e-01

  1 SNES Function norm 3.443876468910e-05

  2 SNES Function norm 2.173976047771e-08
Step 21: Displacement 0.10, Newton its: displacement[0]
  0 SNES Function norm 5.351154144416e-01

  1 SNES Function norm 3.365772682532e-05

  2 SNES Function norm 2.235524583456e-08
Step 22: Displacement 0.11, Newton its: displacement[0]
  0 SNES Function norm 5.336291938032e-01

  1 SNES Function norm 3.292006364727e-05

  2 SNES Function norm 2.325744673506e-08
Step 23: Displacement 0.12, Newton its: displacement[0]
  0 SNES Function norm 5.321732723048e-01

  1 SNES Function norm 3.222373785084e-05

  2 SNES Function norm 2.439185345636e-08
Step 24: Displacement 0.12, Newton its: displacement[0]
  0 SNES Function norm 5.307469787278e-01

  1 SNES Function norm 3.156743638318e-05

  2 SNES Function norm 2.514735934925e-08
Step 25: Displacement 0.12, Newton its: displacement[0]
  0 SNES Function norm 5.293496750129e-01

  1 SNES Function norm 3.095028918161e-05

  2 SNES Function norm 2.582038434253e-08
Step 26: Displacement 0.13, Newton its: displacement[0]
  0 SNES Function norm 5.279807503953e-01

  1 SNES Function norm 3.036949636924e-05

  2 SNES Function norm 2.668997525020e-08
Step 27: Displacement 0.14, Newton its: displacement[0]
  0 SNES Function norm 5.266396202286e-01

  1 SNES Function norm 2.982312675107e-05

  2 SNES Function norm 2.735463172589e-08
Step 28: Displacement 0.14, Newton its: displacement[0]
  0 SNES Function norm 5.253257282215e-01

  1 SNES Function norm 2.930958958159e-05

  2 SNES Function norm 2.819046814021e-08
Step 29: Displacement 0.14, Newton its: displacement[0]
  0 SNES Function norm 5.240385419426e-01

  1 SNES Function norm 2.882860649492e-05

  2 SNES Function norm 2.967076303667e-08
Step 30: Displacement 0.15, Newton its: displacement[0]
  0 SNES Function norm 5.227775554376e-01

  1 SNES Function norm 2.837765204951e-05

  2 SNES Function norm 3.009598853790e-08
Step 31: Displacement 0.15, Newton its: displacement[0]
  0 SNES Function norm 5.215422795598e-01

  1 SNES Function norm 2.795676048829e-05

  2 SNES Function norm 3.024578677680e-08
Step 32: Displacement 0.16, Newton its: displacement[0]
  0 SNES Function norm 5.203322468389e-01

  1 SNES Function norm 2.756067094600e-05

  2 SNES Function norm 3.186801484963e-08
Step 33: Displacement 0.17, Newton its: displacement[0]
  0 SNES Function norm 5.191470083138e-01

  1 SNES Function norm 2.719304846752e-05

  2 SNES Function norm 3.210209025216e-08
Step 34: Displacement 0.17, Newton its: displacement[0]
  0 SNES Function norm 5.179861382037e-01

  1 SNES Function norm 2.684693256986e-05

  2 SNES Function norm 3.275713990921e-08
Step 35: Displacement 0.17, Newton its: displacement[0]
  0 SNES Function norm 5.168492142155e-01

  1 SNES Function norm 2.652462389253e-05

  2 SNES Function norm 3.398779734971e-08
Step 36: Displacement 0.18, Newton its: displacement[0]
  0 SNES Function norm 5.157358437801e-01

  1 SNES Function norm 2.622448668387e-05

  2 SNES Function norm 3.471718583914e-08
Step 37: Displacement 0.18, Newton its: displacement[0]
  0 SNES Function norm 5.146456426460e-01

  1 SNES Function norm 2.594157450638e-05

  2 SNES Function norm 3.536956483730e-08
Step 38: Displacement 0.19, Newton its: displacement[0]
  0 SNES Function norm 5.135782361278e-01

  1 SNES Function norm 2.567708388737e-05

  2 SNES Function norm 3.703993116434e-08
Step 39: Displacement 0.20, Newton its: displacement[0]
  0 SNES Function norm 5.125332722662e-01

  1 SNES Function norm 2.543139118806e-05

  2 SNES Function norm 3.769301642020e-08
Step 40: Displacement 0.20, Newton its: displacement[0]
  0 SNES Function norm 5.115104005239e-01

  1 SNES Function norm 2.519905826327e-05

  2 SNES Function norm 3.852991070897e-08
Step 41: Displacement 0.20, Newton its: displacement[0]
  0 SNES Function norm 5.105092867323e-01

  1 SNES Function norm 2.498064715394e-05

  2 SNES Function norm 3.889925186871e-08
Step 42: Displacement 0.21, Newton its: displacement[0]
  0 SNES Function norm 5.095296078018e-01

  1 SNES Function norm 2.477335818158e-05

  2 SNES Function norm 4.096078214321e-08
Step 43: Displacement 0.21, Newton its: displacement[0]
  0 SNES Function norm 5.085710497564e-01

  1 SNES Function norm 2.458020893740e-05

  2 SNES Function norm 4.205717780190e-08
Step 44: Displacement 0.22, Newton its: displacement[0]
  0 SNES Function norm 5.076333049493e-01

  1 SNES Function norm 2.439340468114e-05

  2 SNES Function norm 4.174961541144e-08
Step 45: Displacement 0.23, Newton its: displacement[0]
  0 SNES Function norm 5.067160795858e-01

  1 SNES Function norm 2.421830430955e-05

  2 SNES Function norm 4.343845476343e-08
Step 46: Displacement 0.23, Newton its: displacement[0]
  0 SNES Function norm 5.058190829578e-01

  1 SNES Function norm 2.405191487266e-05

  2 SNES Function norm 4.365511928447e-08
Step 47: Displacement 0.23, Newton its: displacement[0]
  0 SNES Function norm 5.049420363119e-01

  1 SNES Function norm 2.388877683851e-05

  2 SNES Function norm 4.617608332470e-08
Step 48: Displacement 0.24, Newton its: displacement[0]
  0 SNES Function norm 5.040846634167e-01

  1 SNES Function norm 2.373713382620e-05

  2 SNES Function norm 4.575080199777e-08
Step 49: Displacement 0.24, Newton its: displacement[0]
  0 SNES Function norm 5.032466990636e-01

  1 SNES Function norm 2.358859512848e-05

  2 SNES Function norm 4.586479606545e-08
Step 50: Displacement 0.25, Newton its: displacement[0]
  0 SNES Function norm 5.024278821889e-01

  1 SNES Function norm 2.344396133356e-05

  2 SNES Function norm 4.762134262714e-08
Step 51: Displacement 0.26, Newton its: displacement[0]
  0 SNES Function norm 5.016279587940e-01

  1 SNES Function norm 2.330466640200e-05

  2 SNES Function norm 4.808103163538e-08
Step 52: Displacement 0.26, Newton its: displacement[0]
  0 SNES Function norm 5.008466800308e-01

  1 SNES Function norm 2.317143847779e-05

  2 SNES Function norm 4.942227684810e-08
Step 53: Displacement 0.27, Newton its: displacement[0]
  0 SNES Function norm 5.000837981192e-01

  1 SNES Function norm 2.303351475547e-05

  2 SNES Function norm 5.122917771753e-08
Step 54: Displacement 0.27, Newton its: displacement[0]
  0 SNES Function norm 4.993390831064e-01

  1 SNES Function norm 2.290552029698e-05

  2 SNES Function norm 5.195606934457e-08
Step 55: Displacement 0.28, Newton its: displacement[0]
  0 SNES Function norm 4.986122970536e-01

  1 SNES Function norm 2.277766478094e-05

  2 SNES Function norm 5.387677345812e-08
Step 56: Displacement 0.28, Newton its: displacement[0]
  0 SNES Function norm 4.979032121551e-01

  1 SNES Function norm 2.265305619007e-05

  2 SNES Function norm 5.526741775715e-08
Step 57: Displacement 0.28, Newton its: displacement[0]
  0 SNES Function norm 4.972116017662e-01

  1 SNES Function norm 2.252556793758e-05

  2 SNES Function norm 5.490686753222e-08
Step 58: Displacement 0.29, Newton its: displacement[0]
  0 SNES Function norm 4.965372498101e-01

  1 SNES Function norm 2.240245001431e-05

  2 SNES Function norm 5.744259270703e-08
Step 59: Displacement 0.29, Newton its: displacement[0]
  0 SNES Function norm 4.958799382803e-01

  1 SNES Function norm 2.228081168600e-05

  2 SNES Function norm 5.875156023072e-08
Step 60: Displacement 0.30, Newton its: displacement[0]
  0 SNES Function norm 4.952394536037e-01

  1 SNES Function norm 2.215599770495e-05

  2 SNES Function norm 5.909674118022e-08
Step 61: Displacement 0.30, Newton its: displacement[0]
  0 SNES Function norm 4.946155906868e-01

  1 SNES Function norm 2.203533173254e-05

  2 SNES Function norm 6.026441944635e-08
Step 62: Displacement 0.31, Newton its: displacement[0]
  0 SNES Function norm 4.940081423013e-01

  1 SNES Function norm 2.191027148086e-05

  2 SNES Function norm 6.084178201312e-08
Step 63: Displacement 0.32, Newton its: displacement[0]
  0 SNES Function norm 4.934169057680e-01

  1 SNES Function norm 2.179125834204e-05

  2 SNES Function norm 6.172072132423e-08
Step 64: Displacement 0.32, Newton its: displacement[0]
  0 SNES Function norm 4.928416855408e-01

  1 SNES Function norm 2.166613289667e-05

  2 SNES Function norm 6.412154442028e-08
Step 65: Displacement 0.33, Newton its: displacement[0]
  0 SNES Function norm 4.922822884833e-01

  1 SNES Function norm 2.154374600775e-05

  2 SNES Function norm 6.548336816734e-08
Step 66: Displacement 0.33, Newton its: displacement[0]
  0 SNES Function norm 4.917385166696e-01

  1 SNES Function norm 2.142043590524e-05

  2 SNES Function norm 6.530064941200e-08
Step 67: Displacement 0.34, Newton its: displacement[0]
  0 SNES Function norm 4.912101857876e-01

  1 SNES Function norm 2.130029724617e-05

  2 SNES Function norm 6.845768404147e-08
Step 68: Displacement 0.34, Newton its: displacement[0]
  0 SNES Function norm 4.906971121226e-01

  1 SNES Function norm 2.117443545147e-05

  2 SNES Function norm 7.005243697507e-08
Step 69: Displacement 0.34, Newton its: displacement[0]
  0 SNES Function norm 4.901991063740e-01

  1 SNES Function norm 2.105513213085e-05

  2 SNES Function norm 7.049532130868e-08
Step 70: Displacement 0.35, Newton its: displacement[0]
  0 SNES Function norm 4.897159902977e-01

  1 SNES Function norm 2.093022243995e-05

  2 SNES Function norm 7.168811381574e-08
Step 71: Displacement 0.35, Newton its: displacement[0]
  0 SNES Function norm 4.892475909654e-01

  1 SNES Function norm 2.080500989116e-05

  2 SNES Function norm 7.252751647621e-08
Step 72: Displacement 0.36, Newton its: displacement[0]
  0 SNES Function norm 4.887937286887e-01

  1 SNES Function norm 2.068221824813e-05

  2 SNES Function norm 7.633401377864e-08
Step 73: Displacement 0.36, Newton its: displacement[0]
  0 SNES Function norm 4.883542351548e-01

  1 SNES Function norm 2.055890118121e-05

  2 SNES Function norm 7.644995016555e-08
Step 74: Displacement 0.37, Newton its: displacement[0]
  0 SNES Function norm 4.879289336557e-01

  1 SNES Function norm 2.043613084504e-05

  2 SNES Function norm 8.030945498345e-08
Step 75: Displacement 0.38, Newton its: displacement[0]
  0 SNES Function norm 4.875176632725e-01

  1 SNES Function norm 2.030298206242e-05

  2 SNES Function norm 7.916427411705e-08
Step 76: Displacement 0.38, Newton its: displacement[0]
  0 SNES Function norm 4.871202592864e-01

  1 SNES Function norm 2.018656546926e-05

  2 SNES Function norm 8.278124554019e-08
Step 77: Displacement 0.39, Newton its: displacement[0]
  0 SNES Function norm 4.867365545887e-01

  1 SNES Function norm 2.006973476185e-05

  2 SNES Function norm 8.507299668319e-08
Step 78: Displacement 0.39, Newton its: displacement[0]
  0 SNES Function norm 4.863663940436e-01

  1 SNES Function norm 1.993656507259e-05

  2 SNES Function norm 8.539452541966e-08
Step 79: Displacement 0.40, Newton its: displacement[0]
  0 SNES Function norm 4.860096187883e-01

  1 SNES Function norm 1.981255505934e-05

  2 SNES Function norm 8.759381094377e-08
Step 80: Displacement 0.40, Newton its: displacement[0]
  0 SNES Function norm 4.856660706378e-01

  1 SNES Function norm 1.968671462515e-05

  2 SNES Function norm 9.012769185717e-08
Step 81: Displacement 0.41, Newton its: displacement[0]
  0 SNES Function norm 4.853356014943e-01

  1 SNES Function norm 1.956346772465e-05

  2 SNES Function norm 8.656907273266e-08
Step 82: Displacement 0.41, Newton its: displacement[0]
  0 SNES Function norm 4.850180571370e-01

  1 SNES Function norm 1.944624282003e-05

  2 SNES Function norm 9.065909014071e-08
Step 83: Displacement 0.41, Newton its: displacement[0]
  0 SNES Function norm 4.847132899364e-01

  1 SNES Function norm 1.932334896408e-05

  2 SNES Function norm 9.420391353122e-08
Step 84: Displacement 0.42, Newton its: displacement[0]
  0 SNES Function norm 4.844211556210e-01

  1 SNES Function norm 1.919269994432e-05

  2 SNES Function norm 9.593043888466e-08
Step 85: Displacement 0.42, Newton its: displacement[0]
  0 SNES Function norm 4.841415074398e-01

  1 SNES Function norm 1.908029549640e-05

  2 SNES Function norm 1.004046245423e-07

  3 SNES Function norm 4.643424577816e-08
Step 86: Displacement 0.43, Newton its: displacement[0]
  0 SNES Function norm 4.838742052376e-01

  1 SNES Function norm 1.895911070480e-05

  2 SNES Function norm 1.009439594583e-07

  3 SNES Function norm 4.850885635247e-08
Step 87: Displacement 0.43, Newton its: displacement[0]
  0 SNES Function norm 4.836191089626e-01

  1 SNES Function norm 1.883813919295e-05

  2 SNES Function norm 1.035329018750e-07

  3 SNES Function norm 5.322646091498e-08
Step 88: Displacement 0.44, Newton its: displacement[0]
  0 SNES Function norm 4.833760822513e-01

  1 SNES Function norm 1.871530837060e-05

  2 SNES Function norm 1.024686116612e-07

  3 SNES Function norm 5.097578955372e-08
Step 89: Displacement 0.45, Newton its: displacement[0]
  0 SNES Function norm 4.831449903722e-01

  1 SNES Function norm 1.860627754488e-05

  2 SNES Function norm 1.019659585402e-07

  3 SNES Function norm 5.035976186749e-08
Step 90: Displacement 0.45, Newton its: displacement[0]
  0 SNES Function norm 4.829256960207e-01

  1 SNES Function norm 1.848226921550e-05

  2 SNES Function norm 1.076677242440e-07

  3 SNES Function norm 4.964468227774e-08
Step 91: Displacement 0.46, Newton its: displacement[0]
  0 SNES Function norm 4.827180722962e-01

  1 SNES Function norm 1.836576542857e-05

  2 SNES Function norm 1.115530515312e-07

  3 SNES Function norm 5.607081571056e-08
Step 92: Displacement 0.46, Newton its: displacement[0]
  0 SNES Function norm 4.825219878627e-01

  1 SNES Function norm 1.825308637961e-05

  2 SNES Function norm 1.118724654336e-07

  3 SNES Function norm 5.641347301063e-08
Step 93: Displacement 0.47, Newton its: displacement[0]
  0 SNES Function norm 4.823373213241e-01

  1 SNES Function norm 1.813196400999e-05

  2 SNES Function norm 1.150667923027e-07

  3 SNES Function norm 5.933208728863e-08
Step 94: Displacement 0.47, Newton its: displacement[0]
  0 SNES Function norm 4.821639399076e-01

  1 SNES Function norm 1.802285006013e-05

  2 SNES Function norm 1.143248288768e-07

  3 SNES Function norm 5.323371818596e-08
Step 95: Displacement 0.47, Newton its: displacement[0]
  0 SNES Function norm 4.820017260682e-01

  1 SNES Function norm 1.790232623818e-05

  2 SNES Function norm 1.161375198641e-07

  3 SNES Function norm 5.815432044340e-08
Step 96: Displacement 0.48, Newton its: displacement[0]
  0 SNES Function norm 4.818505599508e-01

  1 SNES Function norm 1.779532788228e-05

  2 SNES Function norm 1.163541952879e-07

  3 SNES Function norm 5.885077660921e-08
Step 97: Displacement 0.48, Newton its: displacement[0]
  0 SNES Function norm 4.817103216757e-01

  1 SNES Function norm 1.768857662133e-05

  2 SNES Function norm 1.251853487961e-07

  3 SNES Function norm 5.988215822513e-08
Step 98: Displacement 0.49, Newton its: displacement[0]
  0 SNES Function norm 4.815808926209e-01

  1 SNES Function norm 1.757595263388e-05

  2 SNES Function norm 1.247764993860e-07

  3 SNES Function norm 6.139597851943e-08
Step 99: Displacement 0.49, Newton its: displacement[0]
  0 SNES Function norm 4.814621617310e-01

  1 SNES Function norm 1.747564505452e-05

  2 SNES Function norm 1.212384187038e-07

  3 SNES Function norm 6.255341944064e-08
Step 100: Displacement 0.50, Newton its: displacement[0]

Visualize the deformed configuration of the specimen.

try:
    import matplotlib.cm as mcm
    import matplotlib.colors as mcolors
    import pyvista

    import dolfinx.plot

    plotter = pyvista.Plotter(window_size=[600, 400], off_screen=True)
    topology, cell_types, x = dolfinx.plot.vtk_mesh(V)
    grid = pyvista.UnstructuredGrid(topology, cell_types, x)
    vals = np.zeros((x.shape[0], 3))
    vals[:, : len(u)] = u.x.array.reshape((x.shape[0], len(u)))
    grid["u"] = vals

    # Calculate displacement magnitude
    mag = np.linalg.norm(vals[:, : len(u)], axis=1)
    grid["mag"] = mag

    warped = grid.warp_by_vector("u", factor=1)
    plotter.add_mesh(warped, scalars="mag", cmap="viridis", show_edges=False, show_scalar_bar=False)
    plotter.view_xy()
    plotter.camera.tight()
    image = plotter.screenshot(None, transparent_background=True, return_img=True)

    fig, ax = plt.subplots(figsize=(7, 5))
    ax.imshow(image)
    ax.axis("off")

    # Create matching colorbar using matplotlib
    norm = mcolors.Normalize(vmin=mag.min(), vmax=mag.max())
    smap = mcm.ScalarMappable(cmap="viridis", norm=norm)
    smap.set_array([])
    fig.colorbar(smap, ax=ax, label=r"Displacement magnitude, $\sqrt{u_x^2 + u_y^2}$")

    fig.savefig("displacement_nn.png", bbox_inches="tight", dpi=200, transparent=True)
    plt.show()

except ImportError:
    print("pyvista required for this plot")

2026-07-07 12:15:30.976 (  84.410s) [    7F65EBCBA140]vtkXOpenGLRenderWindow.:1458  WARN| bad X server connection. DISPLAY=:99.0

../_images/fc963e111374d23678e49f7662161f459ad639ea28a0f0a1a35626e214402bdc.png

Verification against Analytical UFL Baseline#

To verify the results obtained with the ICNN model wrapped via external operators, we implement the pure UFL implementation of the Isihara model.

The variational formulation of the baseline problem is:

Find \(\mathbf{u}_{\text{UFL}} \in V\) satisfying the Dirichlet boundary conditions such that:

\[ F_{\text{UFL}}(\mathbf{u}_{\text{UFL}}; \mathbf{v}) = \int\limits_\Omega \mathbf{P}(\mathbf{F}) : \nabla\mathbf{v} \, \mathrm{d}\mathbf{x} = 0, \quad \forall \mathbf{v} \in V_0 \]

where the stress tensor is evaluated symbolically via \(\mathbf{P}(\mathbf{F}) = \frac{\partial W_{\text{Isihara}}}{\partial \mathbf{F}}\) using UFL’s automatic differentiation tool ufl.diff.

The strain energy density function \(W_{\text{Isihara}}(\mathbf{F})\) is defined explicitly as:

\[ W_{\text{Isihara}}(\mathbf{F}) = 0.5(\bar{I}_1 - 3) + (\bar{I}_2 - 3) + (\bar{I}_1 - 3)^2 + 1.5(J - 1)^2 \]

where the deviatoric invariants \(\bar{I}_1\) and \(\bar{I}_2\) are:

\(\bar{I}_1 = J^{-2/3} I_1\)
\(\bar{I}_2 = J^{-4/3} I_2\)
\(J = \det(\mathbf{F})\) is the volume ratio (determinant of the deformation gradient tensor \(\mathbf{F}\)).
\(I_1 = \text{tr}(\mathbf{C}) = \text{tr}(\mathbf{F}^T\mathbf{F}) + 1.0\) and \(I_2 = I_1 + J^2 - 1.0\) under the 2D plane strain assumption.

The UFL formulation is very straightforward

u_UFL = fem.Function(V)
v = ufl.TestFunction(V)
d = len(u_UFL)
# Deformation gradient
F_ = ufl.variable(ufl.Identity(d) + ufl.grad(u_UFL))

C = F_.T * F_
J_ = ufl.det(F_)
I1 = ufl.tr(C) + 1.0
I2 = I1 + J_**2 - 1.0

# Isihara model parameters
I1_bar = (J_ ** (-2.0 / 3.0)) * I1
I2_bar = (J_ ** (-4.0 / 3.0)) * I2

W_Isihara = 0.5 * (I1_bar - 3.0) + (I2_bar - 3.0) + (I1_bar - 3.0) ** 2 + 1.5 * (J_ - 1.0) ** 2

P = ufl.diff(W_Isihara, F_)

F_UFL = ufl.inner(ufl.grad(v), P) * dx

We could then explicitly apply ufl.derivative to F_UFL to obtain the Jacobian J_UFL but this will be done automatically in NonlinearProblem under the hood.

Now we simply define the nonlinear problem with the same petsc_options as previously and solve the UFL-based problem.

problem_UFL = NonlinearProblem(
    F_UFL, u_UFL, bcs=bcs_u, petsc_options=petsc_options, petsc_options_prefix="UFL_hyperelasticity_"
)

# Apply a tensile load by incrementally increasing displacement on the top edge
u_UFL.name = "UFL_displacement"
u_UFL.x.array[:] = 0
for step in range(1, n_steps + 1):
    u_D_top.value = step * max_displacement / n_steps
    num_its, converged = problem_UFL.solve()
    assert converged, f"Newton solver did not converge at step {step}"
    u_UFL.x.scatter_forward()
    if domain.comm.rank == 0:
        print(f"Step {step}: Displacement {u_D_top.value:.3f}, Newton its: {num_its}")

  0 SNES Function norm 5.730767073990e-01
  1 SNES Function norm 8.865300274904e-05

SNES Function norm 1.312771416962e-08
Step 1: Displacement 0.005, Newton its: UFL_displacement[0]
SNES Function norm 5.705566827913e-01
SNES Function norm 8.237825585370e-05
SNES Function norm 1.137315352558e-08
Step 2: Displacement 0.010, Newton its: UFL_displacement[0]
SNES Function norm 5.680931751910e-01
SNES Function norm 7.655929870430e-05

SNES Function norm 9.991684648725e-09
Step 3: Displacement 0.015, Newton its: UFL_displacement[0]
SNES Function norm 5.656847430312e-01
SNES Function norm 7.122527694131e-05
SNES Function norm 8.913796479586e-09
Step 4: Displacement 0.020, Newton its: UFL_displacement[0]
SNES Function norm 5.633299363685e-01
SNES Function norm 6.639373225685e-05

SNES Function norm 8.047983053948e-09
Step 5: Displacement 0.025, Newton its: UFL_displacement[0]
SNES Function norm 5.610273099902e-01
SNES Function norm 6.206845884518e-05
SNES Function norm 7.309914147031e-09
Step 6: Displacement 0.030, Newton its: UFL_displacement[0]
SNES Function norm 5.587754329979e-01

SNES Function norm 5.823892550296e-05
SNES Function norm 6.644401631866e-09
Step 7: Displacement 0.035, Newton its: UFL_displacement[0]
SNES Function norm 5.565728970219e-01
SNES Function norm 5.488113814197e-05

SNES Function norm 6.031754704934e-09
Step 8: Displacement 0.040, Newton its: UFL_displacement[0]
SNES Function norm 5.544183231203e-01
SNES Function norm 5.195980858306e-05
SNES Function norm 5.481423989748e-09
Step 9: Displacement 0.045, Newton its: UFL_displacement[0]
SNES Function norm 5.523103674304e-01

SNES Function norm 4.943147519131e-05
SNES Function norm 5.019115526424e-09
Step 10: Displacement 0.050, Newton its: UFL_displacement[0]
SNES Function norm 5.502477256600e-01
SNES Function norm 4.724807109121e-05
SNES Function norm 4.671620728081e-09
Step 11: Displacement 0.055, Newton its: UFL_displacement[0]
SNES Function norm 5.482291365175e-01

SNES Function norm 4.536040448602e-05
SNES Function norm 4.452654974798e-09
Step 12: Displacement 0.060, Newton its: UFL_displacement[0]
SNES Function norm 5.462533841923e-01
SNES Function norm 4.372110950292e-05

SNES Function norm 4.354659719729e-09
Step 13: Displacement 0.065, Newton its: UFL_displacement[0]
SNES Function norm 5.443192999992e-01
SNES Function norm 4.228680402388e-05
SNES Function norm 4.350899914610e-09
Step 14: Displacement 0.070, Newton its: UFL_displacement[0]
SNES Function norm 5.424257633042e-01
SNES Function norm 4.101938696434e-05

SNES Function norm 4.405680641047e-09
Step 15: Displacement 0.075, Newton its: UFL_displacement[0]
SNES Function norm 5.405717018409e-01
SNES Function norm 3.988656369521e-05
SNES Function norm 4.485090702418e-09
Step 16: Displacement 0.080, Newton its: UFL_displacement[0]
SNES Function norm 5.387560915209e-01
SNES Function norm 3.886177695379e-05

SNES Function norm 4.563008588571e-09
Step 17: Displacement 0.085, Newton its: UFL_displacement[0]
SNES Function norm 5.369779558327e-01
SNES Function norm 3.792374439528e-05
SNES Function norm 4.622335693794e-09
Step 18: Displacement 0.090, Newton its: UFL_displacement[0]
SNES Function norm 5.352363649117e-01
SNES Function norm 3.705578309123e-05

SNES Function norm 4.653682811962e-09
Step 19: Displacement 0.095, Newton its: UFL_displacement[0]
SNES Function norm 5.335304343525e-01
SNES Function norm 3.624505911080e-05
SNES Function norm 4.653335312188e-09
Step 20: Displacement 0.100, Newton its: UFL_displacement[0]
SNES Function norm 5.318593238273e-01
SNES Function norm 3.548185462229e-05

SNES Function norm 4.621367509458e-09
Step 21: Displacement 0.105, Newton its: UFL_displacement[0]
SNES Function norm 5.302222355598e-01
SNES Function norm 3.475890576448e-05
SNES Function norm 4.560167800968e-09
Step 22: Displacement 0.110, Newton its: UFL_displacement[0]
SNES Function norm 5.286184126998e-01

SNES Function norm 3.407083546540e-05
SNES Function norm 4.473380013845e-09
Step 23: Displacement 0.115, Newton its: UFL_displacement[0]
SNES Function norm 5.270471376328e-01
SNES Function norm 3.341368627101e-05

SNES Function norm 4.365183463122e-09
Step 24: Displacement 0.120, Newton its: UFL_displacement[0]
SNES Function norm 5.255077302534e-01
SNES Function norm 3.278454708416e-05
SNES Function norm 4.239825740138e-09
Step 25: Displacement 0.125, Newton its: UFL_displacement[0]
SNES Function norm 5.239995462260e-01

SNES Function norm 3.218126225553e-05
SNES Function norm 4.101338265697e-09
Step 26: Displacement 0.130, Newton its: UFL_displacement[0]
SNES Function norm 5.225219752503e-01
SNES Function norm 3.160220963515e-05
SNES Function norm 3.953376035248e-09
Step 27: Displacement 0.135, Newton its: UFL_displacement[0]
SNES Function norm 5.210744393475e-01

SNES Function norm 3.104613447678e-05
SNES Function norm 3.799145068336e-09
Step 28: Displacement 0.140, Newton its: UFL_displacement[0]
SNES Function norm 5.196563911755e-01
SNES Function norm 3.051202744211e-05
SNES Function norm 3.641382271083e-09
Step 29: Displacement 0.145, Newton its: UFL_displacement[0]
SNES Function norm 5.182673123851e-01
SNES Function norm 2.999903672320e-05
SNES Function norm 3.482371535577e-09
Step 30: Displacement 0.150, Newton its: UFL_displacement[0]
SNES Function norm 5.169067120203e-01

SNES Function norm 2.950640611343e-05
SNES Function norm 3.323978657829e-09
Step 31: Displacement 0.155, Newton its: UFL_displacement[0]
SNES Function norm 5.155741249689e-01
SNES Function norm 2.903343251372e-05

SNES Function norm 3.167696524891e-09
Step 32: Displacement 0.160, Newton its: UFL_displacement[0]
SNES Function norm 5.142691104666e-01
SNES Function norm 2.857943780176e-05

SNES Function norm 3.014692791210e-09
Step 33: Displacement 0.165, Newton its: UFL_displacement[0]
SNES Function norm 5.129912506562e-01
SNES Function norm 2.814375116277e-05
SNES Function norm 2.865856807878e-09
Step 34: Displacement 0.170, Newton its: UFL_displacement[0]
SNES Function norm 5.117401492031e-01

SNES Function norm 2.772569892860e-05
SNES Function norm 2.721845643990e-09
Step 35: Displacement 0.175, Newton its: UFL_displacement[0]
SNES Function norm 5.105154299680e-01
SNES Function norm 2.732459971445e-05

SNES Function norm 2.583121862840e-09
Step 36: Displacement 0.180, Newton its: UFL_displacement[0]
SNES Function norm 5.093167357367e-01
SNES Function norm 2.693976320938e-05
SNES Function norm 2.449989295947e-09
Step 37: Displacement 0.185, Newton its: UFL_displacement[0]
SNES Function norm 5.081437270066e-01
SNES Function norm 2.657049140633e-05

SNES Function norm 2.322624432056e-09
Step 38: Displacement 0.190, Newton its: UFL_displacement[0]
SNES Function norm 5.069960808289e-01
SNES Function norm 2.621608138739e-05

SNES Function norm 2.201100712296e-09
Step 39: Displacement 0.195, Newton its: UFL_displacement[0]
SNES Function norm 5.058734897058e-01
SNES Function norm 2.587582900861e-05
SNES Function norm 2.085410134130e-09
Step 40: Displacement 0.200, Newton its: UFL_displacement[0]
SNES Function norm 5.047756605418e-01

SNES Function norm 2.554903302494e-05
SNES Function norm 1.975484630074e-09
Step 41: Displacement 0.205, Newton its: UFL_displacement[0]
SNES Function norm 5.037023136472e-01
SNES Function norm 2.523499930804e-05

SNES Function norm 1.871205350285e-09
Step 42: Displacement 0.210, Newton its: UFL_displacement[0]
SNES Function norm 5.026531817924e-01
SNES Function norm 2.493304493441e-05
SNES Function norm 1.772419027760e-09
Step 43: Displacement 0.215, Newton its: UFL_displacement[0]
SNES Function norm 5.016280093121e-01

SNES Function norm 2.464250196249e-05
SNES Function norm 1.678946316507e-09
Step 44: Displacement 0.220, Newton its: UFL_displacement[0]
SNES Function norm 5.006265512574e-01
SNES Function norm 2.436272081249e-05
SNES Function norm 1.590590452221e-09
Step 45: Displacement 0.225, Newton its: UFL_displacement[0]
SNES Function norm 4.996485725940e-01
SNES Function norm 2.409307317310e-05
SNES Function norm 1.507141791755e-09
Step 46: Displacement 0.230, Newton its: UFL_displacement[0]
SNES Function norm 4.986938474452e-01

SNES Function norm 2.383295441117e-05
SNES Function norm 1.428386112252e-09
Step 47: Displacement 0.235, Newton its: UFL_displacement[0]
SNES Function norm 4.977621583781e-01
SNES Function norm 2.358178548219e-05

SNES Function norm 1.354103368816e-09
Step 48: Displacement 0.240, Newton its: UFL_displacement[0]
SNES Function norm 4.968532957308e-01
SNES Function norm 2.333901435123e-05
SNES Function norm 1.284076170981e-09
Step 49: Displacement 0.245, Newton its: UFL_displacement[0]
SNES Function norm 4.959670569796e-01
SNES Function norm 2.310411695380e-05
SNES Function norm 1.218088422648e-09
Step 50: Displacement 0.250, Newton its: UFL_displacement[0]
SNES Function norm 4.951032461443e-01

SNES Function norm 2.287659773894e-05
SNES Function norm 1.155929089471e-09
Step 51: Displacement 0.255, Newton its: UFL_displacement[0]
SNES Function norm 4.942616732297e-01
SNES Function norm 2.265598983280e-05
SNES Function norm 1.097391641750e-09
Step 52: Displacement 0.260, Newton its: UFL_displacement[0]
SNES Function norm 4.934421537024e-01

SNES Function norm 2.244185486478e-05
SNES Function norm 1.042276300173e-09
Step 53: Displacement 0.265, Newton its: UFL_displacement[0]
SNES Function norm 4.926445080001e-01

SNES Function norm 2.223378251718e-05
SNES Function norm 9.903913258103e-10
Step 54: Displacement 0.270, Newton its: UFL_displacement[0]
SNES Function norm 4.918685610735e-01
SNES Function norm 2.203138983017e-05
SNES Function norm 9.415513121462e-10
Step 55: Displacement 0.275, Newton its: UFL_displacement[0]
SNES Function norm 4.911141419580e-01

SNES Function norm 2.183432031027e-05
SNES Function norm 8.955793111592e-10
Step 56: Displacement 0.280, Newton its: UFL_displacement[0]
SNES Function norm 4.903810833740e-01
SNES Function norm 2.164224289054e-05

SNES Function norm 8.523047910202e-10
Step 57: Displacement 0.285, Newton its: UFL_displacement[0]
SNES Function norm 4.896692213553e-01
SNES Function norm 2.145485076498e-05
SNES Function norm 8.115676954868e-10
Step 58: Displacement 0.290, Newton its: UFL_displacement[0]
SNES Function norm 4.889783949024e-01
SNES Function norm 2.127186014742e-05

SNES Function norm 7.732140799215e-10
Step 59: Displacement 0.295, Newton its: UFL_displacement[0]
SNES Function norm 4.883084456614e-01
SNES Function norm 2.109300896756e-05
SNES Function norm 7.370972403016e-10
Step 60: Displacement 0.300, Newton its: UFL_displacement[0]
SNES Function norm 4.876592176254e-01
SNES Function norm 2.091805554374e-05

SNES Function norm 7.030801771492e-10
Step 61: Displacement 0.305, Newton its: UFL_displacement[0]
SNES Function norm 4.870305568590e-01
SNES Function norm 2.074677724531e-05
SNES Function norm 6.710312890855e-10
Step 62: Displacement 0.310, Newton its: UFL_displacement[0]
SNES Function norm 4.864223112427e-01

SNES Function norm 2.057896916620e-05
SNES Function norm 6.408265659420e-10
Step 63: Displacement 0.315, Newton its: UFL_displacement[0]
SNES Function norm 4.858343302377e-01
SNES Function norm 2.041444282363e-05

SNES Function norm 6.123510696947e-10
Step 64: Displacement 0.320, Newton its: UFL_displacement[0]
SNES Function norm 4.852664646689e-01
SNES Function norm 2.025302489910e-05
SNES Function norm 5.854936225403e-10
Step 65: Displacement 0.325, Newton its: UFL_displacement[0]
SNES Function norm 4.847185665259e-01
SNES Function norm 2.009455602158e-05

SNES Function norm 5.601515640923e-10
Step 66: Displacement 0.330, Newton its: UFL_displacement[0]
SNES Function norm 4.841904887804e-01
SNES Function norm 1.993888961233e-05
SNES Function norm 5.362271497796e-10
Step 67: Displacement 0.335, Newton its: UFL_displacement[0]
SNES Function norm 4.836820852190e-01

SNES Function norm 1.978589078379e-05
SNES Function norm 5.136298705144e-10
Step 68: Displacement 0.340, Newton its: UFL_displacement[0]
SNES Function norm 4.831932102910e-01
SNES Function norm 1.963543531189e-05

SNES Function norm 4.922738120486e-10
Step 69: Displacement 0.345, Newton its: UFL_displacement[0]
SNES Function norm 4.827237189698e-01
SNES Function norm 1.948740866965e-05
SNES Function norm 4.720786675368e-10
Step 70: Displacement 0.350, Newton its: UFL_displacement[0]
SNES Function norm 4.822734666274e-01
SNES Function norm 1.934170512724e-05

SNES Function norm 4.529711446307e-10
Step 71: Displacement 0.355, Newton its: UFL_displacement[0]
SNES Function norm 4.818423089207e-01
SNES Function norm 1.919822692815e-05
SNES Function norm 4.348794028313e-10
Step 72: Displacement 0.360, Newton its: UFL_displacement[0]
SNES Function norm 4.814301016899e-01
SNES Function norm 1.905688352243e-05

SNES Function norm 4.177394001901e-10
Step 73: Displacement 0.365, Newton its: UFL_displacement[0]
SNES Function norm 4.810367008668e-01
SNES Function norm 1.891759087238e-05
SNES Function norm 4.014902098976e-10
Step 74: Displacement 0.370, Newton its: UFL_displacement[0]
SNES Function norm 4.806619623937e-01

SNES Function norm 1.878027081614e-05
SNES Function norm 3.860745051272e-10
Step 75: Displacement 0.375, Newton its: UFL_displacement[0]
SNES Function norm 4.803057421513e-01

SNES Function norm 1.864485049137e-05
SNES Function norm 3.714393887637e-10
Step 76: Displacement 0.380, Newton its: UFL_displacement[0]
SNES Function norm 4.799678958956e-01
SNES Function norm 1.851126181809e-05

SNES Function norm 3.575367595387e-10
Step 77: Displacement 0.385, Newton its: UFL_displacement[0]
SNES Function norm 4.796482792026e-01
SNES Function norm 1.837944102946e-05
SNES Function norm 3.443187508144e-10
Step 78: Displacement 0.390, Newton its: UFL_displacement[0]
SNES Function norm 4.793467474210e-01

SNES Function norm 1.824932825325e-05
SNES Function norm 3.317445126717e-10
Step 79: Displacement 0.395, Newton its: UFL_displacement[0]
SNES Function norm 4.790631556320e-01
SNES Function norm 1.812086714052e-05

SNES Function norm 3.197746662254e-10
Step 80: Displacement 0.400, Newton its: UFL_displacement[0]
SNES Function norm 4.787973586151e-01
SNES Function norm 1.799400452922e-05
SNES Function norm 3.083707269526e-10
Step 81: Displacement 0.405, Newton its: UFL_displacement[0]
SNES Function norm 4.785492108211e-01

SNES Function norm 1.786869015499e-05
SNES Function norm 2.975000879602e-10
Step 82: Displacement 0.410, Newton its: UFL_displacement[0]
SNES Function norm 4.783185663498e-01
SNES Function norm 1.774487638196e-05

SNES Function norm 2.871300306471e-10
Step 83: Displacement 0.415, Newton its: UFL_displacement[0]
SNES Function norm 4.781052789333e-01
SNES Function norm 1.762251797815e-05
SNES Function norm 2.772310120684e-10
Step 84: Displacement 0.420, Newton its: UFL_displacement[0]
SNES Function norm 4.779092019248e-01
SNES Function norm 1.750157190822e-05
SNES Function norm 2.677762471809e-10
Step 85: Displacement 0.425, Newton its: UFL_displacement[0]
SNES Function norm 4.777301882911e-01

SNES Function norm 1.738199715151e-05
SNES Function norm 2.587402337873e-10
Step 86: Displacement 0.430, Newton its: UFL_displacement[0]
SNES Function norm 4.775680906096e-01
SNES Function norm 1.726375454782e-05

SNES Function norm 2.500992131462e-10
Step 87: Displacement 0.435, Newton its: UFL_displacement[0]
SNES Function norm 4.774227610696e-01
SNES Function norm 1.714680665207e-05
SNES Function norm 2.418304046679e-10
Step 88: Displacement 0.440, Newton its: UFL_displacement[0]
SNES Function norm 4.772940514767e-01
SNES Function norm 1.703111761468e-05

SNES Function norm 2.339144929513e-10
Step 89: Displacement 0.445, Newton its: UFL_displacement[0]
SNES Function norm 4.771818132602e-01
SNES Function norm 1.691665306593e-05
SNES Function norm 2.263309587636e-10
Step 90: Displacement 0.450, Newton its: UFL_displacement[0]
SNES Function norm 4.770858974845e-01
SNES Function norm 1.680338002533e-05

SNES Function norm 2.190631658727e-10
Step 91: Displacement 0.455, Newton its: UFL_displacement[0]
SNES Function norm 4.770061548622e-01
SNES Function norm 1.669126680942e-05
SNES Function norm 2.120946647646e-10
Step 92: Displacement 0.460, Newton its: UFL_displacement[0]
SNES Function norm 4.769424357707e-01

SNES Function norm 1.658028295867e-05
SNES Function norm 2.054095979926e-10
Step 93: Displacement 0.465, Newton its: UFL_displacement[0]
SNES Function norm 4.768945902699e-01
SNES Function norm 1.647039916525e-05

SNES Function norm 1.989936370215e-10
Step 94: Displacement 0.470, Newton its: UFL_displacement[0]
SNES Function norm 4.768624681233e-01
SNES Function norm 1.636158721380e-05
SNES Function norm 1.928331503108e-10
Step 95: Displacement 0.475, Newton its: UFL_displacement[0]
SNES Function norm 4.768459188199e-01

SNES Function norm 1.625381992210e-05
SNES Function norm 1.869158009713e-10
Step 96: Displacement 0.480, Newton its: UFL_displacement[0]
SNES Function norm 4.768447915988e-01

SNES Function norm 1.614707109413e-05
SNES Function norm 1.812297913672e-10
Step 97: Displacement 0.485, Newton its: UFL_displacement[0]
SNES Function norm 4.768589354738e-01
SNES Function norm 1.604131547152e-05

SNES Function norm 1.757642031911e-10
Step 98: Displacement 0.490, Newton its: UFL_displacement[0]
SNES Function norm 4.768881992615e-01
SNES Function norm 1.593652869201e-05
SNES Function norm 1.705079995514e-10
Step 99: Displacement 0.495, Newton its: UFL_displacement[0]
SNES Function norm 4.769324316085e-01
SNES Function norm 1.583268725357e-05

  2 SNES Function norm 1.654515683162e-10
Step 100: Displacement 0.500, Newton its: UFL_displacement[0]

try:
    import matplotlib.cm as mcm
    import matplotlib.colors as mcolors
    import pyvista

    import dolfinx.plot

    plotter = pyvista.Plotter(window_size=[600, 400], off_screen=True)
    topology, cell_types, x = dolfinx.plot.vtk_mesh(V)
    grid = pyvista.UnstructuredGrid(topology, cell_types, x)
    vals = np.zeros((x.shape[0], 3))
    vals[:, : len(u_UFL)] = u_UFL.x.array.reshape((x.shape[0], len(u_UFL)))
    grid["u"] = vals

    # Calculate displacement magnitude
    mag = np.linalg.norm(vals[:, : len(u_UFL)], axis=1)
    grid["mag"] = mag

    warped = grid.warp_by_vector("u", factor=1)
    plotter.add_mesh(warped, scalars="mag", cmap="viridis", show_edges=False, show_scalar_bar=False)
    plotter.view_xy()
    plotter.camera.tight()
    image = plotter.screenshot(None, transparent_background=True, return_img=True)

    fig, ax = plt.subplots(figsize=(7, 5))
    ax.imshow(image)
    ax.axis("off")

    # Create matching colorbar using matplotlib
    norm = mcolors.Normalize(vmin=mag.min(), vmax=mag.max())
    smap = mcm.ScalarMappable(cmap="viridis", norm=norm)
    smap.set_array([])
    fig.colorbar(smap, ax=ax, label=r"Displacement magnitude, $\sqrt{u_x^2 + u_y^2}$")

    fig.savefig("displacement_ufl.png", bbox_inches="tight", dpi=200, transparent=True)
    plt.show()

except ImportError:
    print("pyvista required for this plot")

../_images/2c4a295c578d7167db607e597881efca446c12db976766ae4df09973ca07b9cc.png

Comparison of Deformed Configurations#

Below, we compare the deformed configuration predicted by the Neural Network external operator against the analytical UFL baseline.

Neural Network Operator

UFL Baseline (Analytical)

Error Analysis#

To quantitatively assess the difference between the Input-Convex Neural Network (ICNN) external operator and the analytical UFL baseline, we compute two error metrics at the end of the simulation:

Relative Maximum Nodal Displacement Error (relative \(L^\infty\) norm):

\[ e_{\infty, \text{rel}} = \frac{\max_i |u_i - u_{\text{UFL}, i}|}{\max_i |u_{\text{UFL}, i}|} \]

np.abs(u.x.array[:] - u_UFL.x.array[:]).max() / np.abs(u_UFL.x.array[:]).max()

np.float64(0.004955246427269522)

Absolute \(L^2\) Displacement Difference:

\[ e_{L^2} = \sqrt{\int\limits_\Omega (\mathbf{u} - \mathbf{u}_{\text{UFL}})^2 \, \mathrm{d}\mathbf{x}} \]

np.sqrt(MPI.COMM_WORLD.allreduce(fem.assemble_scalar(fem.form((u - u_UFL) ** 2 * dx)), op=MPI.SUM))

np.float64(0.0008536299984819733)

Errors look quite promising :)

References#

[IHT51] (1,2)

Akira Isihara, Natsuki Hashitsume, and Masao Tatibana. Statistical theory of rubber-like elasticity. v. (two-dimensional stretching). The Journal of Chemical Physics, 19(12):1508–1512, 1951. doi:10.1063/1.1748111.

[TJZ+22] (1,2)

Prakash Thakolkaran, Akshay Joshi, Yiwen Zheng, Moritz Flaschel, Laura De Lorenzis, and Siddhant Kumar. NN-EUCLID: Deep-learning hyperelasticity without stress data. Journal of the Mechanics and Physics of Solids, 169:1–29, 2022. doi:10.1016/j.jmps.2022.105076.