diff --git a/demos/jax/elastoplasticity/Hosford_yield_surface.md b/demos/jax/elastoplasticity/Hosford_yield_surface.md
index 70125c3..2e2aa37 100644
--- a/demos/jax/elastoplasticity/Hosford_yield_surface.md
+++ b/demos/jax/elastoplasticity/Hosford_yield_surface.md
@@ -49,7 +49,7 @@ elastic_model = LinearElasticIsotropic(E, nu)
 a = 10.0
 sig0 = 500.0
 material = GeneralIsotropicHardening(
-    elastic_model, lambda sig: sig0, lambda sig: Hosford_stress(sig, a=a)
+    elastic_model, lambda p: sig0, lambda sig: Hosford_stress(sig, a=a)
 )
 ```
 
diff --git a/demos/mfront/finite_strain_elastoplasticity/finite_strain_elastoplasticity.md b/demos/mfront/finite_strain_elastoplasticity/finite_strain_elastoplasticity.md
index 06ba9d8..f28a69e 100644
--- a/demos/mfront/finite_strain_elastoplasticity/finite_strain_elastoplasticity.md
+++ b/demos/mfront/finite_strain_elastoplasticity/finite_strain_elastoplasticity.md
@@ -62,11 +62,14 @@ H0.setEntryName("HardeningSlope");
 };
 ```
 
-## `FEniCS` implementation
+## `FEniCSx` implementation
 
 We define a box mesh representing half of a beam oriented along the $x$-direction. The beam will be fully clamped on its left side and symmetry conditions will be imposed on its right extremity. The loading consists of a uniform self-weight.
 
-<img src="finite_strain_plasticity_solution.png" width="500">
+```{image} finite_strain_plasticity_solution.png
+:align: center
+:width: 500px
+```
 
 ```{code-cell} ipython3
 import numpy as np
@@ -84,12 +87,16 @@ from dolfinx_materials.utils import (
     nonsymmetric_tensor_to_vector,
 )
 
+
+comm = MPI.COMM_WORLD
+rank = comm.rank
+
 current_path = os.getcwd()
 
 length, width, height = 1.0, 0.04, 0.1
-nx, ny, nz = 30, 5, 8
+nx, ny, nz = 20, 4, 6
 domain = mesh.create_box(
-    MPI.COMM_WORLD,
+    comm,
     [(0, -width / 2, -height / 2.0), (length, width / 2, height / 2.0)],
     [nx, ny, nz],
     cell_type=mesh.CellType.tetrahedron,
@@ -134,10 +141,11 @@ material = MFrontMaterial(
     "LogarithmicStrainPlasticity",
     material_properties={"YieldStrength": 250e6, "HardeningSlope": 1e6},
 )
-print(material.behaviour.getBehaviourType())
-print(material.behaviour.getKinematic())
-print(material.gradient_names, material.gradient_sizes)
-print(material.flux_names, material.flux_sizes)
+if rank == 0:
+    print(material.behaviour.getBehaviourType())
+    print(material.behaviour.getKinematic())
+    print(material.gradient_names, material.gradient_sizes)
+    print(material.flux_names, material.flux_sizes)
 ```
 
 In this large-strain setting, the `QuadratureMapping` acts from the deformation gradient $\boldsymbol{F}=\boldsymbol{I}+\nabla\boldsymbol{u}$ to the first Piola-Kirchhoff stress $\boldsymbol{P}$. We must therefore register the deformation gradient as `Identity(3)+grad(u)`. 
@@ -178,7 +186,7 @@ Finally, we setup the nonlinear problem, the corresponding Newton solver and sol
 
 problem = NonlinearMaterialProblem(qmap, Res, Jac, u, bcs)
 
-newton = NewtonSolver(MPI.COMM_WORLD)
+newton = NewtonSolver(comm)
 newton.rtol = 1e-4
 newton.atol = 1e-4
 newton.convergence_criterion = "residual"
@@ -203,7 +211,8 @@ for i, t in enumerate(load_steps[1:]):
 
     converged, it = problem.solve(newton, print_solution=False)
 
-    print(f"Increment {i+1} converged in {it} iterations.")
+    if rank == 0:
+        print(f"Increment {i+1} converged in {it} iterations.")
 
     p0 = qmap.project_on("EquivalentPlasticStrain", ("DG", 0))
 
@@ -211,7 +220,10 @@ for i, t in enumerate(load_steps[1:]):
     vtk.write_function(p0, t)
 
     w = u.sub(2).collapse()
-    results[i + 1, 0] = max(np.abs(w.vector.array))
+    local_max = max(np.abs(w.vector.array))
+    # Perform the reduction to get the global maximum on rank 0
+    global_max = comm.reduce(local_max, op=MPI.MAX, root=0)
+    results[i + 1, 0] = global_max
     results[i + 1, 1] = t
 vtk.close()
 ```
@@ -220,11 +232,12 @@ During the load incrementation, we monitor the evolution of the maximum vertical
 The load-displacement curve exhibits a classical elastoplastic behavior rapidly followed by a stiffening behavior due to membrane catenary effects. 
 
 ```{code-cell} ipython3
-plt.figure()
-plt.plot(results[:, 0], results[:, 1], "-oC3")
-plt.xlabel("Displacement")
-plt.ylabel("Load")
-plt.show()
+if rank==0:
+    plt.figure()
+    plt.plot(results[:, 0], results[:, 1], "-oC3")
+    plt.xlabel("Displacement")
+    plt.ylabel("Load")
+    plt.show()
 ```
 
 ## References
diff --git a/demos/mfront/finite_strain_elastoplasticity/finite_strain_elastoplasticity.py b/demos/mfront/finite_strain_elastoplasticity/finite_strain_elastoplasticity.py
index 32783d9..95009c8 100644
--- a/demos/mfront/finite_strain_elastoplasticity/finite_strain_elastoplasticity.py
+++ b/demos/mfront/finite_strain_elastoplasticity/finite_strain_elastoplasticity.py
@@ -1,21 +1,38 @@
-# %% [markdown]
+# ---
+# jupyter:
+#   jupytext:
+#     text_representation:
+#       extension: .py
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.16.1
+#   kernelspec:
+#     display_name: Python 3
+#     language: python
+#     name: python3
+# ---
+
 # # Finite-strain elastoplasticity within the logarithmic strain framework
 #
-# This demo is dedicated to the resolution of a finite-strain elastoplastic problem using the logarithmic strain framework proposed in [@miehe_anisotropic_2002].
+# This demo is dedicated to the resolution of a finite-strain elastoplastic problem using the logarithmic strain framework proposed in {cite:p}`miehe2002anisotropic`.
+#
+# ```{tip}
+# This simulation is a bit heavy to run so we suggest running it in parallel.
+# ```
 #
 # ## Logarithmic strains
 #
-# This framework expresses constitutive relations between the Hencky strain measure $\boldsymbol{H} = \dfrac{1}{2}\log (\boldsymbol{F}^T\cdot\boldsymbol{F})$ and its dual stress measure $\boldsymbol{T}$. This approach makes it possible to extend classical small strain constitutive relations to a finite-strain setting. In particular, the total (Hencky) strain can be split additively into many contributions (elastic, plastic, thermal, swelling, etc.) e.g. $\boldsymbol{H}=\boldsymbol{H}^e+\boldsymbol{H}^p$. Its trace is also linked with the volume change $J=\exp(\operatorname{tr}(\boldsymbol{H}))$. As a result, the deformation gradient $\boldsymbol{F}$ is used for expressing the Hencky strain $\boldsymbol{H}$, a small-strain constitutive law is then written for the $(\boldsymbol{H},\boldsymbol{T})$-pair and the dual stress $\boldsymbol{T}$ is then post-processed to an appropriate stress measure such as the Cauchy stress $\boldsymbol{\sigma}$ or Piola-Kirchhoff stresses.
+# This framework expresses constitutive relations between the Hencky strain measure $\boldsymbol{H} = \dfrac{1}{2}\log (\boldsymbol{F}^\text{T}\cdot\boldsymbol{F})$ and its dual stress measure $\boldsymbol{T}$. This approach makes it possible to extend classical small strain constitutive relations to a finite-strain setting. In particular, the total (Hencky) strain can be split additively into many contributions (elastic, plastic, thermal, swelling, etc.) e.g. $\boldsymbol{H}=\boldsymbol{H}^e+\boldsymbol{H}^p$. Its trace is also linked with the volume change $J=\exp(\operatorname{tr}(\boldsymbol{H}))$. As a result, the deformation gradient $\boldsymbol{F}$ is used for expressing the Hencky strain $\boldsymbol{H}$, a small-strain constitutive law is then written for the $(\boldsymbol{H},\boldsymbol{T})$-pair and the dual stress $\boldsymbol{T}$ is then post-processed to an appropriate stress measure such as the Cauchy stress $\boldsymbol{\sigma}$ or Piola-Kirchhoff stresses.
 #
 # ## `MFront` implementation
 #
-# The logarithmic strain framework discussed in the previous paragraph consists merely as a pre-processing and a post-processing stages of the behaviour integration. The pre-processing stage compute the logarithmic strain and its increment and the post-processing stage inteprets the stress resulting from the behaviour integration as the dual stress $\boldsymbol{T}$ and convert it to the Cauchy stress.
+# The logarithmic strain framework discussed in the previous paragraph consists merely as a pre-processing and a post-processing stages of the behavior integration. The pre-processing stage compute the logarithmic strain and its increment and the post-processing stage inteprets the stress resulting from the behavior integration as the dual stress $\boldsymbol{T}$ and convert it to the Cauchy stress.
 #
-# `MFront` provides the `@StrainMeasure` keyword that allows to specify which strain measure is used by the behaviour. When choosing the `Hencky` strain measure, `MFront` automatically generates those pre- and post-processing stages, allowing the user to focus on the behaviour integration.
+# `MFront` provides the `@StrainMeasure` keyword that allows to specify which strain measure is used by the behavior. When choosing the `Hencky` strain measure, `MFront` automatically generates those pre- and post-processing stages, allowing the user to focus on the behavior integration.
 #
 # This leads to the following implementation (see the [small-strain elastoplasticity example](https://thelfer.github.io/mgis/web/mgis_fenics_small_strain_elastoplasticity.html) for details about the various implementations available):
 #
-# ```cpp
+# ```
 # @DSL Implicit;
 #
 # @Behaviour LogarithmicStrainPlasticity;
@@ -45,15 +62,16 @@
 # };
 # ```
 #
-# ## `FEniCS` implementation
+# ## `FEniCSx` implementation
 #
 # We define a box mesh representing half of a beam oriented along the $x$-direction. The beam will be fully clamped on its left side and symmetry conditions will be imposed on its right extremity. The loading consists of a uniform self-weight.
 #
-#
-# <img src="finite_strain_plasticity_solution.png" width="500">
-#
+# ```{image} finite_strain_plasticity_solution.png
+# :align: center
+# :width: 500px
+# ```
 
-# %%
+# +
 import numpy as np
 import matplotlib.pyplot as plt
 import os
@@ -69,12 +87,16 @@
     nonsymmetric_tensor_to_vector,
 )
 
+
+comm = MPI.COMM_WORLD
+rank = comm.rank
+
 current_path = os.getcwd()
 
 length, width, height = 1.0, 0.04, 0.1
-nx, ny, nz = 30, 5, 8
+nx, ny, nz = 20, 4, 6
 domain = mesh.create_box(
-    MPI.COMM_WORLD,
+    comm,
     [(0, -width / 2, -height / 2.0), (length, width / 2, height / 2.0)],
     [nx, ny, nz],
     cell_type=mesh.CellType.tetrahedron,
@@ -109,26 +131,25 @@ def right(x):
 du = ufl.TrialFunction(V)
 v = ufl.TestFunction(V)
 u = fem.Function(V, name="Displacement")
+# -
 
-# %% [markdown]
-# The `MFrontMaterial` instance is loaded from the `MFront` `LogarithmicStrainPlasticity` behaviour. This behaviour is a finite-strain behaviour (`material.is_finite_strain=True`) which relies on a kinematic description using the total deformation gradient $\boldsymbol{F}$. By default, a `MFront` behaviour always returns the Cauchy stress as the stress measure after integration. However, the stress variable dual to the deformation gradient is the first Piola-Kirchhoff (PK1) stress. An internal option of the MGIS interface is therefore used in the finite-strain context to return the PK1 stress as the "flux" associated to the "gradient" $\boldsymbol{F}$. Both quantities are non-symmetric tensors, aranged as a 9-dimensional vector in 3D following [`MFront` conventions on tensors](http://tfel.sourceforge.net/tensors.html).
+# The `MFrontMaterial` instance is loaded from the `MFront` `LogarithmicStrainPlasticity` behavior. This behavior is a finite-strain behavior (`material.is_finite_strain=True`) which relies on a kinematic description using the total deformation gradient $\boldsymbol{F}$. By default, a `MFront` behavior always returns the Cauchy stress as the stress measure after integration. However, the stress variable dual to the deformation gradient is the first Piola-Kirchhoff (PK1) stress. An internal option of the MGIS interface is therefore used in the finite-strain context to return the PK1 stress as the "flux" associated to the "gradient" $\boldsymbol{F}$. Both quantities are non-symmetric tensors, aranged as a 9-dimensional vector in 3D following [`MFront` conventions on tensors](https://thelfer.github.io/tfel/web/tensors.html).
 
-# %%
 material = MFrontMaterial(
     os.path.join(current_path, "src/libBehaviour.so"),
     "LogarithmicStrainPlasticity",
     material_properties={"YieldStrength": 250e6, "HardeningSlope": 1e6},
 )
-print(material.behaviour.getBehaviourType())
-print(material.behaviour.getKinematic())
-print(material.gradient_names, material.gradient_sizes)
-print(material.flux_names, material.flux_sizes)
+if rank == 0:
+    print(material.behaviour.getBehaviourType())
+    print(material.behaviour.getKinematic())
+    print(material.gradient_names, material.gradient_sizes)
+    print(material.flux_names, material.flux_sizes)
 
-# %% [markdown]
-# The `QuadratureMap` object must therefore register the deformation gradient as `Identity(3)+grad(u)`.
 
+# In this large-strain setting, the `QuadratureMapping` acts from the deformation gradient $\boldsymbol{F}=\boldsymbol{I}+\nabla\boldsymbol{u}$ to the first Piola-Kirchhoff stress $\boldsymbol{P}$. We must therefore register the deformation gradient as `Identity(3)+grad(u)`.
 
-# %%
+# +
 def F(u):
     return nonsymmetric_tensor_to_vector(ufl.Identity(gdim) + ufl.grad(u))
 
@@ -139,75 +160,81 @@ def dF(u):
 
 qmap = QuadratureMap(domain, 2, material)
 qmap.register_gradient("DeformationGradient", F(u))
+# -
 
-# %%
-sig = qmap.fluxes["FirstPiolaKirchhoffStress"]
-Res = (ufl.dot(sig, dF(v)) - ufl.dot(selfweight, v)) * qmap.dx
+# We will work in a Total Lagrangian formulation, writing the weak form of equilibrium on the reference configuration $\Omega_0$, thereby defining the nonlinear residual weak form as:
+# Find $\boldsymbol{u}\in V$ such that:
+#
+# $$
+# \int_{\Omega_0} \boldsymbol{P}(\boldsymbol{F}(\boldsymbol{u})):\nabla \boldsymbol{v} \,\text{d}\Omega - \int_{\Omega_0} \boldsymbol{f}\cdot\boldsymbol{v}\,\text{d}\Omega = 0 \quad \forall \boldsymbol{v}\in V
+# $$
+# where $\boldsymbol{f}$ is the self-weight.
+#
+# The corresponding Jacobian form is computed via automatic differentiation. As for the [small-strain elastoplasticity example](https://thelfer.github.io/mgis/web/mgis_fenics_small_strain_elastoplasticity.html), state variables include the `ElasticStrain` and `EquivalentPlasticStrain` since the same behavior is used as in the small-strain case with the only difference that the total strain is now given by the Hencky strain measure. In particular, the `ElasticStrain` is still a symmetric tensor (vector of dimension 6). Note that it has not been explicitly defined as a state variable in the `MFront` behavior file since this is done automatically when using the `IsotropicPlasticMisesFlow` parser.
+
+PK1 = qmap.fluxes["FirstPiolaKirchhoffStress"]
+Res = (ufl.dot(PK1, dF(v)) - ufl.dot(selfweight, v)) * qmap.dx
 Jac = qmap.derivative(Res, u, du)
 
-# %% [markdown]
-# The loading is then defined and, as for the [small-strain elastoplasticity example](https://thelfer.github.io/mgis/web/mgis_fenics_small_strain_elastoplasticity.html), state variables include the `ElasticStrain` and `EquivalentPlasticStrain` since the same behaviour is used as in the small-strain case with the only difference that the total strain is now given by the Hencky strain measure. In particular, the `ElasticStrain` is still a symmetric tensor (vector of dimension 6). Note that it has not been explicitly defined as a state variable in the `MFront` behaviour file since this is done automatically when using the `IsotropicPlasticMisesFlow` parser.
-#
-# Finally, we setup a few parameters of the Newton non-linear solver.
+# Finally, we setup the nonlinear problem, the corresponding Newton solver and solve the load-stepping problem.
 
-# %%
+# + tags=["hide-output"]
 problem = NonlinearMaterialProblem(qmap, Res, Jac, u, bcs)
 
-newton = NewtonSolver(MPI.COMM_WORLD)
+newton = NewtonSolver(comm)
 newton.rtol = 1e-4
 newton.atol = 1e-4
-newton.convergence_criterion = "incremental"
+newton.convergence_criterion = "residual"
 newton.report = True
 
 # Set solver options
 ksp = newton.krylov_solver
 opts = PETSc.Options()
 option_prefix = ksp.getOptionsPrefix()
-opts[f"{option_prefix}ksp_type"] = "cg"
-opts[f"{option_prefix}pc_type"] = "gamg"
+opts[f"{option_prefix}ksp_type"] = "preonly"
+opts[f"{option_prefix}pc_type"] = "lu"
 opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
 ksp.setFromOptions()
 
 Nincr = 30
 load_steps = np.linspace(0.0, 1.0, Nincr + 1)
 
-vtk_u = io.VTKFile(domain.comm, "results_u.pvd", "w")
-vtk_p = io.VTKFile(domain.comm, "results_p.pvd", "w")
+vtk = io.VTKFile(domain.comm, f"results/{material.name}.pvd", "w")
 results = np.zeros((Nincr + 1, 2))
 for i, t in enumerate(load_steps[1:]):
     selfweight.value[-1] = -50e6 * t
 
     converged, it = problem.solve(newton, print_solution=False)
 
-    print(f"Increment {i+1} converged in {it} iterations.")
+    if rank == 0:
+        print(f"Increment {i+1} converged in {it} iterations.")
 
     p0 = qmap.project_on("EquivalentPlasticStrain", ("DG", 0))
 
-    vtk_u.write_function(u, t)
-    vtk_p.write_function(p0, t)
+    vtk.write_function(u, t)
+    vtk.write_function(p0, t)
 
     w = u.sub(2).collapse()
-    results[i + 1, 0] = max(np.abs(w.vector.array))
+    local_max = max(np.abs(w.vector.array))
+    # Perform the reduction to get the global maximum on rank 0
+    global_max = MPI.COMM_WORLD.reduce(local_max, op=MPI.MAX, root=0)
+    results[i + 1, 0] = global_max
     results[i + 1, 1] = t
-vtk_u.close()
-vtk_p.close()
+vtk.close()
+# -
 
-# %% [markdown]
 # During the load incrementation, we monitor the evolution of the maximum vertical downwards displacement.
-#
-# This simulation is a bit heavy to run so we suggest running it in parallel:
-# ```bash
-# mpirun -np 4 python3 finite_strain_elastoplasticity.py
-# ```
-
-# %% [markdown]
-# The load-displacement curve exhibits a classical elastoplastic behaviour rapidly followed by a stiffening behaviour due to membrane catenary effects.
+# The load-displacement curve exhibits a classical elastoplastic behavior rapidly followed by a stiffening behavior due to membrane catenary effects.
 
-# %%
-plt.figure()
-plt.plot(results[:, 0], results[:, 1], "-o")
-plt.xlabel("Displacement")
-plt.ylabel("Load")
-plt.show()
+if rank==0:
+    plt.figure()
+    plt.plot(results[:, 0], results[:, 1], "-oC3")
+    plt.xlabel("Displacement")
+    plt.ylabel("Load")
+    plt.show()
 
-# %%
+# ## References
+#
+# ```{bibliography}
+# :filter: docname in docnames
+# ```