Nonlinear elasticity with blocked Newton solver

Author: Henrik N. T. Finsberg

SPDX-License-Identifier: MIT

In this example we will solve a nonlinear elasticity problem using a blocked Newton solver. We consider a unit cube domain \(\Omega = [0, 1]^3\) with Dirichlet boundary conditions on the left face and traction force on the right face, and we seek a displacement field \(\mathbf{u}: \Omega \to \mathbb{R}^3\) that solves the momentum balance equation

\[\begin{split} \begin{align} \nabla \cdot \mathbf{P} = 0, \quad \mathbf{X} \in \Omega, \\ \end{align} \end{split}\]

where \(\mathbf{P}\) is the first Piola-Kirchhoff stress tensor, and \(\mathbf{X}\) is the reference configuration. We consider a Neo-Hookean material model, where the strain energy density is given by

\[ \begin{align} \psi = \frac{\mu}{2}(\text{tr}(\mathbf{C}) - 3), \end{align} \]

and the first Piola-Kirchhoff stress tensor is given by

\[ \begin{align} \mathbf{P} = \frac{\partial \psi}{\partial \mathbf{F}} \end{align} \]

where \(\mathbf{F} = \nabla \mathbf{u} + \mathbf{I}\) is the deformation gradient, \(\mathbf{C} = \mathbf{F}^T \mathbf{F}\) is the right Cauchy-Green tensor, \(\mu\) is the shear modulus, and \(p\) is the pressure. We also enforce the incompressibility constraint

\[ \begin{align} J = \det(\mathbf{F}) = 1, \end{align} \]

so that the total Lagrangian is given by

\[ \begin{align} \mathcal{L}(\mathbf{u}, p) = \int_{\Omega} \psi \, dx - \int_{\partial \Omega} t \cdot \mathbf{u} \, ds + \int_{\Omega} p(J - 1) \, dx. \end{align} \]

Here \(t\) is the traction force which is set to \(10\) on the right face of the cube and \(0\) elsewhere. The Euler-Lagrange equations for this problem are given by: Find \(\mathbf{u} \in V\) and \(p \in Q\) such that

\[\begin{split} \begin{align} D_{\delta \mathbf{u} } \mathcal{L}(\mathbf{u}, p) = 0, \quad \forall \delta \mathbf{u} \in V, \\ D_{\delta p} \mathcal{L}(\mathbf{u}, p) = 0, \quad \forall \delta p \in Q, \end{align} \end{split}\]

where \(V\) is the displacement space and \(Q\) is the pressure space. For this we select \(Q_2/P_1\) elements i.e second order Lagrange elements for \(\mathbf{u}\) and first order discontinuous polynomial cubical elements for \(p\), which is a stable element for incompressible elasticity [ABeiraodVL+13]. Note also that the Euler-Lagrange equations can be derived automatically using ufl.

import logging
from mpi4py import MPI
import numpy as np
import ufl
import dolfinx
import scifem

Initialize logging and set log level to info

logging.basicConfig(level=logging.INFO)

We create the mesh and the function spaces

mesh = dolfinx.mesh.create_unit_cube(
    MPI.COMM_WORLD, 3, 3, 3, dolfinx.mesh.CellType.hexahedron, dtype=np.float64
)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2, (3,)))
Q = dolfinx.fem.functionspace(mesh, ("DPC", 1,))

And the test and trial functions

v = ufl.TestFunction(V)
q = ufl.TestFunction(Q)
du = ufl.TrialFunction(V)
dp = ufl.TrialFunction(Q)
u = dolfinx.fem.Function(V)
p = dolfinx.fem.Function(Q)

Next we create the facet tags for the left and right faces

def left(x):
    return np.isclose(x[0], 0)

def right(x):
    return np.isclose(x[0], 1)

facet_tags = scifem.create_entity_markers(
    mesh, mesh.topology.dim - 1, [(1, left, True), (2, right, True)]
)

We create the Dirichlet boundary conditions on the left face

facets_left = facet_tags.find(1)
dofs_left = dolfinx.fem.locate_dofs_topological(V, 2, facets_left)
u_bc_left = dolfinx.fem.Function(V)
u_bc_left.x.array[:] = 0
bc = dolfinx.fem.dirichletbc(u_bc_left, dofs_left)

Define the deformation gradient, right Cauchy-Green tensor, and invariants of the deformation tensors

d = len(u)
I = ufl.Identity(d)             # Identity tensor
F = I + ufl.grad(u)             # Deformation gradient
C = F.T*F                       # Right Cauchy-Green tensor
I1 = ufl.tr(C)                  # First invariant of C
J  = ufl.det(F)                 # Jacobian of F

Traction for to be applied on the right face

t = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(10.0))

N = ufl.FacetNormal(mesh)

# Material parameters and strain energy density
mu = dolfinx.fem.Constant(mesh, 10.0)
psi = (mu / 2)*(I1 - 3)

We for the total Lagrangian

L = psi*ufl.dx - ufl.inner(t * N, u)*ufl.ds(subdomain_data=facet_tags, subdomain_id=2)  + p * (J - 1) * ufl.dx

and take the first variation of the total Lagrangian to obtain the residual

r_u = ufl.derivative(L, u, v)
r_p = ufl.derivative(L, p, q)
R = [r_u, r_p]

We do the same for the second variation to obtain the Jacobian

K = [
    [ufl.derivative(r_u, u, du), ufl.derivative(r_u, p, dp)],
    [ufl.derivative(r_p, u, du), ufl.derivative(r_p, p, dp)],
]

Now we can create the Newton solver and solve the problem

petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
solver = scifem.NewtonSolver(R, K, [u, p], max_iterations=25, bcs=[bc], petsc_options=petsc_options)

/tmp/ipykernel_1500/2054379347.py:2: DeprecationWarning: NewtonSolver is deprecated in favor of `dolfinx.fem.petsc.NonlinearProblem`.
  solver = scifem.NewtonSolver(R, K, [u, p], max_iterations=25, bcs=[bc], petsc_options=petsc_options)

We can also set a callback function that is called before and after the solve, which takes the solver object as input

def pre_solve(solver: scifem.NewtonSolver):
    print(f"Starting solve with {solver.max_iterations} iterations")

def post_solve(solver: scifem.NewtonSolver):
    print(f"Solve completed in with correction norm {solver.dx.norm(0)}")

solver.set_pre_solve_callback(pre_solve)
solver.set_post_solve_callback(post_solve)

solver.solve()

INFO:scifem.solvers:Newton iteration 1: r (abs) = 60.799721691168955 (tol=1e-06), r (rel) = 1.0 (tol=1e-08)

INFO:scifem.solvers:Newton iteration 2: r (abs) = 9.289281588000586 (tol=1e-06), r (rel) = 0.15278493600982776 (tol=1e-08)

INFO:scifem.solvers:Newton iteration 3: r (abs) = 1.9945156715633932 (tol=1e-06), r (rel) = 0.0328046842334986 (tol=1e-08)

Starting solve with 25 iterations
Solve completed in with correction norm 728.8606816927028
Starting solve with 25 iterations
Solve completed in with correction norm 114.65409622584185
Starting solve with 25 iterations
Solve completed in with correction norm 19.979899461096647
Starting solve with 25 iterations

INFO:scifem.solvers:Newton iteration 4: r (abs) = 0.013328841102457887 (tol=1e-06), r (rel) = 0.00021922536373047044 (tol=1e-08)

INFO:scifem.solvers:Newton iteration 5: r (abs) = 1.4961407138852281e-06 (tol=1e-06), r (rel) = 2.4607690171426225e-08 (tol=1e-08)

INFO:scifem.solvers:Newton iteration 6: r (abs) = 1.0791974183168362e-13 (tol=1e-06), r (rel) = 1.7750038787983262e-15 (tol=1e-08)

Solve completed in with correction norm 0.12246721729906325
Starting solve with 25 iterations
Solve completed in with correction norm 1.2223032723898574e-05
Starting solve with 25 iterations
Solve completed in with correction norm 7.377106113520911e-13

6

Finally, we can visualize the solution using pyvista

import pyvista
p = pyvista.Plotter()
topology, cell_types, geometry = dolfinx.plot.vtk_mesh(V)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
linear_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh))
grid["u"] = u.x.array.reshape((geometry.shape[0], 3))
p.add_mesh(linear_grid, style="wireframe", color="k")
warped = grid.warp_by_vector("u", factor=1.5)
p.add_mesh(warped, show_edges=False)
p.show_axes()
if not pyvista.OFF_SCREEN:
    p.show()
else:
    figure_as_array = p.screenshot("displacement.png")

INFO:matplotlib.font_manager:generated new fontManager

2026-08-05 06:37:23.682 (   0.624s) [    7F5BF08C2140]vtkXOpenGLRenderWindow.:1460  WARN| bad X server connection. DISPLAY=:99.0
INFO:trame_server.utils.namespace:Translator(prefix=None)

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
File /dolfinx-env/lib/python3.12/site-packages/aiohttp/web_urldispatcher.py:546, in StaticResource.__init__(self, prefix, directory, name, expect_handler, chunk_size, show_index, follow_symlinks, append_version)
    545 try:
--> 546     directory = Path(directory).expanduser().resolve(strict=True)
    547 except FileNotFoundError as error:

File /usr/lib/python3.12/pathlib.py:1242, in Path.resolve(self, strict)
   1241 try:
-> 1242     s = self._flavour.realpath(self, strict=strict)
   1243 except OSError as e:

File <frozen posixpath>:432, in realpath(filename, strict)

File <frozen posixpath>:477, in _joinrealpath(path, rest, strict, seen)

FileNotFoundError: [Errno 2] No such file or directory: '/dolfinx-env/lib/python3.12/site-packages/trame_client/module/vue3-www'

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
Cell In[21], line 12
      8 warped = grid.warp_by_vector("u", factor=1.5)
      9 p.add_mesh(warped, show_edges=False)
     10 p.show_axes()
     11 if not pyvista.OFF_SCREEN:
---> 12     p.show()
     13 else:
     14     figure_as_array = p.screenshot("displacement.png")

File /dolfinx-env/lib/python3.12/site-packages/pyvista/_deprecate_positional_args.py:243, in _deprecate_positional_args.<locals>._inner_deprecate_positional_args.<locals>.inner_f(*args, **kwargs)
    239             warn_external(msg, PyVistaDeprecationWarning)
    241         warn_positional_args()
--> 243 return f(*args, **kwargs)

File /dolfinx-env/lib/python3.12/site-packages/pyvista/plotting/plotter.py:8231, in Plotter.show(self, title, window_size, interactive, auto_close, interactive_update, full_screen, screenshot, return_img, cpos, jupyter_backend, return_viewer, return_cpos, before_close_callback, store_image_depth, **kwargs)
   8228         jupyter_backend = self._theme.jupyter_backend
   8230     if jupyter_backend is None or jupyter_backend.lower() != 'none':
-> 8231         jupyter_disp = handle_plotter(self, backend=jupyter_backend, **jupyter_kwargs)
   8233 self.render()
   8235 # initial double render needed for certain passes when offscreen

File /dolfinx-env/lib/python3.12/site-packages/pyvista/jupyter/notebook.py:74, in handle_plotter(plotter, backend, screenshot, **kwargs)
     71     if backend in ['server', 'client', 'trame', 'html']:
     72         from pyvista.trame.jupyter import show_trame  # noqa: PLC0415
---> 74         return show_trame(plotter, mode=backend, **kwargs)
     76 except ImportError as e:
     77     # Trame was explicitly requested but not available
     78     _ensure_entry_points()

File /dolfinx-env/lib/python3.12/site-packages/pyvista/trame/jupyter.py:394, in show_trame(plotter, mode, name, server_proxy_enabled, server_proxy_prefix, jupyter_extension_enabled, collapse_menu, add_menu, add_menu_items, default_server_rendering, handler, **kwargs)
    391 kwargs.setdefault('height', dh)
    393 if mode == 'html':
--> 394     return EmbeddableWidget(plotter, **kwargs)
    396 if jupyter_extension_enabled is None:
    397     jupyter_extension_enabled = pv.global_theme.trame.jupyter_extension_enabled

File /dolfinx-env/lib/python3.12/site-packages/pyvista/trame/jupyter.py:137, in EmbeddableWidget.__init__(self, plotter, width, height, **kwargs)
    135     msg = 'Please install `ipywidgets`.'
    136     raise ImportError(msg)
--> 137 scene = plotter.export_html(filename=None)
    138 src = scene.getvalue().replace('"', '&quot;')
    139 # eventually we could maybe expose this, but for now make sure we're at least
    140 # consistent with matplotlib's color (light gray)

File /dolfinx-env/lib/python3.12/site-packages/pyvista/plotting/plotter.py:784, in BasePlotter.export_html(self, filename)
    781     msg = 'Please install trame dependencies: pip install "pyvista[jupyter]"'
    782     raise ImportError(msg)
--> 784 data = self.export_vtksz(filename=None)
    785 buffer = io.StringIO()
    786 write_html(data, buffer)

File /dolfinx-env/lib/python3.12/site-packages/pyvista/plotting/plotter.py:837, in BasePlotter.export_vtksz(self, filename, format)
    835 server = get_server(pv.global_theme.trame.jupyter_server_name)
    836 if not server.running:
--> 837     elegantly_launch(pv.global_theme.trame.jupyter_server_name)
    839 view = PyVistaLocalView(self, trame_server=server)
    841 content = view.export(format=format)

File /dolfinx-env/lib/python3.12/site-packages/pyvista/trame/jupyter.py:487, in elegantly_launch(*args, **kwargs)
    484 # Basically monkey patches asyncio to support this
    485 nest_asyncio2.apply()
--> 487 return asyncio.run(launch_it())

File /dolfinx-env/lib/python3.12/site-packages/nest_asyncio2.py:115, in run(main, debug, loop_factory)
    113 task = asyncio.ensure_future(main, loop=loop)
    114 try:
--> 115     return loop.run_until_complete(task)
    116 finally:
    117     if not task.done():

File /dolfinx-env/lib/python3.12/site-packages/nest_asyncio2.py:230, in _patch_loop.<locals>.run_until_complete(self, future)
    227 if self._ready:
    228     self._run_once()
--> 230 return f.result()

File /usr/lib/python3.12/asyncio/futures.py:203, in Future.result(self)
    201 self.__log_traceback = False
    202 if self._exception is not None:
--> 203     raise self._exception.with_traceback(self._exception_tb)
    204 return self._result

File /usr/lib/python3.12/asyncio/tasks.py:314, in Task.__step_run_and_handle_result(***failed resolving arguments***)
    310 try:
    311     if exc is None:
    312         # We use the `send` method directly, because coroutines
    313         # don't have `__iter__` and `__next__` methods.
--> 314         result = coro.send(None)
    315     else:
    316         result = coro.throw(exc)

File /dolfinx-env/lib/python3.12/site-packages/pyvista/trame/jupyter.py:482, in elegantly_launch.<locals>.launch_it()
    481 async def launch_it():
--> 482     await launch_server(*args, **kwargs).ready

File /dolfinx-env/lib/python3.12/site-packages/pyvista/trame/jupyter.py:215, in launch_server(server, port, host, wslink_backend, **kwargs)
    213 if server._running_stage == 0:
    214     server.controller.on_server_ready.add(on_ready)
--> 215     server.start(
    216         exec_mode='task',
    217         host=host,
    218         port=port,
    219         open_browser=False,
    220         show_connection_info=False,
    221         disable_logging=True,
    222         timeout=0,
    223         backend=wslink_backend,
    224     )
    225 # else, server is already running or launching
    226 return server

File /dolfinx-env/lib/python3.12/site-packages/trame_server/core.py:725, in Server.start(self, port, thread, open_browser, show_connection_info, disable_logging, backend, follow_symlinks, exec_mode, timeout, host, **kwargs)
    722 CoreServer.configure(options)
    724 self._running_stage = 1
--> 725 task = CoreServer.server_start(
    726     options,
    727     **{  # Do a proper merging/override
    728         **kwargs,
    729         "disableLogging": disable_logging,
    730         "backend": backend,
    731         "exec_mode": exec_mode,
    732     },
    733 )
    735 # Manage exit life cycle unless coroutine
    736 if exec_mode == "main":

File /dolfinx-env/lib/python3.12/site-packages/trame_server/protocol.py:52, in CoreServer.server_start(options, disableLogging, backend, exec_mode, **kwargs)
     43 @staticmethod
     44 def server_start(
     45     options,
   (...)     50 ):
     51     # NOTE: **kwargs to wslink's start_webserver are currently unused
---> 52     return server.start_webserver(
     53         options=options,
     54         protocol=CoreServer,
     55         disableLogging=disableLogging,
     56         backend=backend,
     57         exec_mode=exec_mode,
     58         **kwargs,
     59     )

File /dolfinx-env/lib/python3.12/site-packages/wslink/server.py:257, in start_webserver(options, protocol, disableLogging, backend, exec_mode, **_)
    254     server_config["handle_signals"] = not options.nosignalhandlers
    256 # Create the webserver and start it
--> 257 ws_server = create_webserver(server_config, backend=backend)
    259 # Once we have python 3.7 minimum, we can start the server with asyncio.run()
    260 # asyncio.run(ws_server.start())
    261 
    262 # Until then, we can start the server this way
    263 try:

File /dolfinx-env/lib/python3.12/site-packages/wslink/server.py:168, in create_webserver(server_config, backend)
    167 def create_webserver(server_config, backend="aiohttp"):
--> 168     return backends.create_webserver(server_config, backend=backend)

File /dolfinx-env/lib/python3.12/site-packages/wslink/backends/__init__.py:5, in create_webserver(server_config, backend)
      2 if backend == "aiohttp":
      3     from .aiohttp import create_webserver  # noqa: PLC0415
----> 5     return create_webserver(server_config)
      7 if backend == "generic":
      8     from .generic import create_webserver  # noqa: PLC0415

File /dolfinx-env/lib/python3.12/site-packages/wslink/backends/aiohttp/__init__.py:218, in create_webserver(server_config)
    215     return ReverseWebAppServer(server_config)
    217 # Normal web server
--> 218 return WebAppServer(server_config)

File /dolfinx-env/lib/python3.12/site-packages/wslink/backends/aiohttp/__init__.py:118, in WebAppServer.__init__(self, server_config)
    116     # Resolve / => index.html
    117     self.app.router.add_route("GET", "/", _root_handler)
--> 118     self.app.add_routes(routes)
    120 self.app[STATE_KEY] = {}

File /dolfinx-env/lib/python3.12/site-packages/aiohttp/web_app.py:379, in Application.add_routes(self, routes)
    378 def add_routes(self, routes: Iterable[AbstractRouteDef]) -> list[AbstractRoute]:
--> 379     return self.router.add_routes(routes)

File /dolfinx-env/lib/python3.12/site-packages/aiohttp/web_urldispatcher.py:1259, in UrlDispatcher.add_routes(self, routes)
   1257 registered_routes = []
   1258 for route_def in routes:
-> 1259     registered_routes.extend(route_def.register(self))
   1260 return registered_routes

File /dolfinx-env/lib/python3.12/site-packages/aiohttp/web_routedef.py:87, in StaticDef.register(self, router)
     86 def register(self, router: UrlDispatcher) -> list[AbstractRoute]:
---> 87     resource = router.add_static(self.prefix, self.path, **self.kwargs)
     88     routes = resource.get_info().get("routes", {})
     89     return list(routes.values())

File /dolfinx-env/lib/python3.12/site-packages/aiohttp/web_urldispatcher.py:1183, in UrlDispatcher.add_static(self, prefix, path, name, expect_handler, chunk_size, show_index, follow_symlinks, append_version)
   1181 if prefix.endswith("/"):
   1182     prefix = prefix[:-1]
-> 1183 resource = StaticResource(
   1184     prefix,
   1185     path,
   1186     name=name,
   1187     expect_handler=expect_handler,
   1188     chunk_size=chunk_size,
   1189     show_index=show_index,
   1190     follow_symlinks=follow_symlinks,
   1191     append_version=append_version,
   1192 )
   1193 self.register_resource(resource)
   1194 return resource

File /dolfinx-env/lib/python3.12/site-packages/aiohttp/web_urldispatcher.py:548, in StaticResource.__init__(self, prefix, directory, name, expect_handler, chunk_size, show_index, follow_symlinks, append_version)
    546     directory = Path(directory).expanduser().resolve(strict=True)
    547 except FileNotFoundError as error:
--> 548     raise ValueError(f"'{directory}' does not exist") from error
    549 if not directory.is_dir():
    550     raise ValueError(f"'{directory}' is not a directory")

ValueError: '/dolfinx-env/lib/python3.12/site-packages/trame_client/module/vue3-www' does not exist

References

[ABeiraodVL+13]

Ferdinando Auricchio, Lourenço Beirão da Veiga, Carlo Lovadina, Alessandro Reali, Robert L Taylor, and Peter Wriggers. Approximation of incompressible large deformation elastic problems: some unresolved issues. Computational Mechanics, 52:1153–1167, 2013.