Nonlinear elasticity with TDNNS and MCS

18. Nonlinear elasticity with TDNNS and MCS#

This notebook shows how the TDNNS and MCS methods can be extended to nonlinear elasticity by introducing the deformation gradient as an independent unknown and enforce its relation to the displacement weakly. This yields a Hu-Washizu formulation in which the nonlinear material law is evaluated on a regular field, while the method-dependent weak pairing handles the low regularity of the displacement.

The two methods differ mainly in the chosen spaces for the involved fields. TDNNS uses tangential-continuous displacements and normal-normal continuous stresses, while MCS uses normal-continuous displacements and tangential-normal continuous stresses. By this, the nonlinear formulations share the same lifting principle but target different robustness regimes.

18.1. Hyperelasticity in reference configuration#

Let \(u:\Omega\to\mathbb R^d\) be the displacement and

\[ F = I + \nabla u \]

the deformation gradient. A hyperelastic material is described by an energy density \(\Psi(F)\), and the first Piola-Kirchhoff stress is

\[ P = \frac{\partial \Psi}{\partial F}(F). \]

For conforming \(H^1\) displacements, the primal problem formally minimizes

\[ \int_\Omega \Psi(I+\nabla u)\,dx - W_{\rm ext}(u). \]

For TDNNS- and MCS-type discretizations, however, \(\nabla u\) is only a distribution. Nonlinear expressions such as \(\Psi(I+\nabla u)\), \(\det(I+\nabla u)\), or \((I+\nabla u)^T(I+\nabla u)\) are then not meaningful in the classical sense. The deformation gradient therefore has to be lifted to an independent regular field.

18.2. Three-field Hu-Washizu lifting#

Introduce \(F\) as an independent unknown and enforce

\[ F = I + \nabla u \]

with a stress-like Lagrange multiplier \(P\). The common Hu-Washizu Lagrangian is

\[ \mathcal L(u,F,P) = \int_\Omega \Psi(F)\,dx - \langle F-I-\nabla u, P\rangle - W_{\rm ext}(u). \]

The roles of the three fields are:

unknown	role
\(u\)	displacement
\(F\)	lifted deformation gradient / strain representative
\(P\)	first Piola-Kirchhoff stress or stress-like Lagrange multiplier
pairing	method-dependent weak enforcement of \(F=I+\nabla u\)

The material law is evaluated on \(F\), which is chosen in an \(L^2\)-type finite element space. The singular part of the formulation is isolated in the weak pairing between \(\nabla u\) and \(P\).

18.3. TDNNS and MCS as dual choices#

The nonlinear methods inherit their structure from the corresponding linear mixed methods.

method	displacement trace	stress trace	typical target
TDNNS	tangential-continuous, \(H(\operatorname{curl})\)	normal-normal continuous	thin/anisotropic structures
MCS	normal-continuous, \(H(\operatorname{div})\)	tangential-normal continuous	incompressibility

Both methods avoid a fully \(H^1\)-conforming displacement field, and both need a weak interpretation of \(\nabla u\) in nonlinear elasticity.

18.4. Nonlinear TDNNS#

For TDNNS, the displacement is placed in an \(H(\operatorname{curl})\)-conforming space. The stress variable uses the TDNNS stress structure, and the distributional gradient is replaced by the TDNNS pairing

\[ \langle P,\nabla u\rangle = \sum_T\left( \int_T P:\nabla u\,dx - \int_{\partial T} P_{nn} u_n\,ds \right). \]

A full nonlinear TDNNS lifting seeks [Neunteufel, Pechstein, Schöberl. Three-field mixed finite element methods for nonlinear elasticity. Computer Methods in Applied Mechanics and Engineering, (2021).]

\[ (u,F,P)\in H(\operatorname{curl}) \times [L^2(\Omega)]^{d\times d} \times \widetilde H(\operatorname{div}\operatorname{div}), \]

where \(\widetilde H(\operatorname{div}\operatorname{div})=H(\operatorname{div}\operatorname{div}) \times [L^2]^{d\times d}_{\mathrm{skew}}\), for the Lagrangian

\[ \begin{align*} \mathcal{L}(u, F, P)=\int_{\Omega}W(F)\,dx - \langle F-I-\nabla u,P\rangle-\int_{\Omega}f\cdot u\, dx. \end{align*} \]

Here \(F\) is a regular unknown, so general nonlinear material laws can be used without applying them to a distributional displacement gradient.

18.4.1. Symmetric TDNNS lifting#

Since

\[ \operatorname{skw}\nabla u = \operatorname{skw}(\operatorname{curl}u)\in L^2, \]

there is no need to lift the complete deformation gradient. It is enough to lift its symmetric part. The deformation gradient is represented as

\[ F = F_{\rm sym} + \operatorname{skw}(\operatorname{curl}u), \]

and only the symmetric part of the stress multiplier is needed. With hybridization, one obtains the Lagrangian for \((u,F_{\mathrm{sym}},P_{\mathrm{sym}})\in H(\operatorname{curl})\times [L^2(\Omega)]^{d\times d}_{\mathrm{sym}}\times H(\operatorname{div}\operatorname{div})\)

\[\begin{split} \begin{aligned} \mathcal L(u,F_{\rm sym},P_{\rm sym},\alpha) =&\int_\Omega \Psi\left(F_{\rm sym}+\operatorname{skw}(\operatorname{curl}u)\right)\,dx + \langle \varepsilon(u),P_{\rm sym}\rangle \\ &-\int_\Omega (F_{\rm sym}-I):P_{\rm sym}\,dx +\sum_T\int_{\partial T} P_{{\rm sym},nn}\alpha_n\,ds - W_{\rm ext}(u,\alpha). \end{aligned} \end{split}\]

The next example uses this reduced lifting.

18.5. Cook’s membrane#

from ngsolve import *
from ngsolve.webgui import Draw
import ipywidgets as widgets
from ngsolve.meshes import MakeQuadMesh

L_1 = 48  # 48m
L_2 = 44  # 44m
L_3 = 16  # 16m


def GenerateMesh(num_el=2):
    mesh = MakeQuadMesh(
        nx=num_el,
        ny=num_el,
        mapping=lambda x, y: (L_1 * x, L_2 * x + L_2 * y - (L_2 - L_3) * x * y),
    )
    return mesh


mesh = GenerateMesh()

mu = Parameter(80.194)  # N/mm^2
lam = Parameter(400889.8)  # N/mm^2

par = Parameter(0)
force = CoefficientFunction((0, par * 32))


def NeoHooke(F):
    C = F.trans * F
    return mu / 2 * (Trace(C - Id(2)) - log(Det(C))) + lam / 2 * (Det(F) - 1) ** 2


def SkewT(v):
    if v.dim == 3:
        return 0.5 * CoefficientFunction(
            (0, -v[2], v[1], v[2], 0, -v[0], -v[1], v[0], 0), dims=(3, 3)
        )
    else:
        return 0.5 * CoefficientFunction((0, -v, v, 0), dims=(2, 2))


order = 3

# clamped at left, free elsewhere
V = HCurl(mesh, order=order, dirichlet="left")
Sigma = Discontinuous(HDivDiv(mesh, order=order))
Hyb = NormalFacetFESpace(mesh, order=order, dirichlet="left")
Gamma = MatrixValued(L2(mesh, order=order + 1), symmetric=True)
fes = V * Sigma * Gamma * Hyb
gf_solution = GridFunction(fes)

u, Psym, Fsym, uh = fes.TrialFunction()
Fmat = Fsym + SkewT(curl(u))

n = specialcf.normal(2)

a = BilinearForm(fes, symmetric=True, condense=True)
a += Variation((NeoHooke(Fmat)) * dx).Compile()
a += Variation(InnerProduct(Fsym - Grad(u) - Id(2), Psym) * dx).Compile()
a += Variation((u - uh) * n * Psym[n, n] * dx(element_boundary=True)).Compile()
a += Variation(-(force * uh.Trace() + force * u.Trace()) * ds("right")).Compile()

gf_u, gf_Psym, gf_Fsym, gf_uh = gf_solution.components
gf_F = gf_Fsym + SkewT(curl(gf_u))

# F=I
gf_Fsym.Set(Id(2))

scene = Draw(
    gf_u,
    mesh,
    deformation=True,
    settings={
        "camera": {
            "transformations": [{"type": "move", "dir": (0, 0, 1), "dist": -1.5}]
        }
    },
)
tw = widgets.Text(value="step = 0")
display(tw)
gf_history = GridFunction(fes, multidim=0)

num_steps = 40
with TaskManager():
    for step in range(num_steps):
        par.Set((step + 1) / num_steps)
        solvers.Newton(
            a,
            gf_solution,
            maxerr=1e-8,
            printing=False,
            inverse="sparsecholesky",
        )
        scene.Redraw()
        tw.value = f"step = {step+1}/{num_steps}"
        gf_history.AddMultiDimComponent(gf_solution.vec)

Draw(
    gf_history.components[0],
    mesh,
    animate=True,
    min=0,
    max=25,
    autoscale=True,
    deformation=True,
    settings={
        "camera": {
            "transformations": [{"type": "move", "dir": (0, 0, 1), "dist": -1.5}]
        }
    },
);

18.6. Thin cylindrical shell#

The next benchmark is a compressed cylindrical shell. The shell is a half cylinder of length \(L=15\) and thickness \(t=0.2\). The right lower edge is fixed, and a vertical load is applied on the right upper face. The mesh is intentionally coarse so that the example remains notebook-sized.

from ngsolve.meshes import MakeHexMesh
from netgen.occ import *
from math import pi, cos, sin


def GenerateCylindricalShellMesh(thickness=0.2, length=15, maxh=2, order=2):
    inner_radius = 9 - thickness / 2
    outer_radius = inner_radius + thickness

    outer_cylinder = Cylinder(
        (0, 0, 0), (0, 0, 1), r=outer_radius, h=length, mantle="outer"
    )
    inner_cylinder = Cylinder(
        (0, 0, 0), (0, 0, 1), r=inner_radius, h=length, mantle="inner"
    )

    shell = (outer_cylinder - inner_cylinder) * Box(
        (-length + thickness, -(length + thickness), 0), (0, length + thickness, length)
    )
    shell.faces.Min(Z).name = "back"
    shell.faces.Max(Z).name = "front"
    shell.faces.Max(X).name = "rightbot"
    shell.faces.Max(Y).name = "righttop"

    trf = gp_GTrsf(
        (
            outer_radius / inner_radius,
            0,
            0,
            0,
            outer_radius / inner_radius,
            0,
            0,
            0,
            1,
        ),
        (0, 0, 0),
    )
    shell.faces["inner"].Identify(
        shell.faces["outer"], "cs", IdentificationType.CLOSESURFACES, trf
    )

    return Mesh(
        OCCGeometry(shell).GenerateMesh(maxh=maxh, quad_dominated=True)
    ).Curve(order)


def CylindricalShellNeoHooke(mu, lam, F):
    return (
        mu / 2 * (Trace(F.trans * F) - 3)
        - mu * log(Det(F))
        + lam / 2 * log(Det(F)) ** 2
    )


def CylindricalShellSkew(v):
    return 0.5 * CoefficientFunction(
        (0, -v[2], v[1], v[2], 0, -v[0], -v[1], v[0], 0), dims=(3, 3)
    )


def CylindricalShellPoint(thickness=0.2):
    inner_radius = 9 - thickness / 2
    return (-0.5 * thickness, 2 * inner_radius + 0.5 * thickness, 15 - 0.5 * thickness)


def CylindricalShellDeflection(displacement, mesh, thickness=0.2):
    inner_radius = 9 - thickness / 2
    points = [
        CylindricalShellPoint(thickness),
        (-0.5 * thickness, inner_radius + 0.5 * thickness, 15 - 0.5 * thickness),
        (-0.5 * thickness, inner_radius + 0.5 * thickness, 7.5),
        (-1.0, inner_radius - 0.5 * thickness, 7.5),
    ]
    for point in points:
        try:
            return displacement(mesh(*point))[1]
        except Exception:
            pass
    raise RuntimeError("Could not find a cylindrical-shell evaluation point in the mesh.")


cyl_shell_mesh = GenerateCylindricalShellMesh()
Draw(cyl_shell_mesh);

def SolveCylindricalShellStandard(
    mesh,
    order=2,
    thickness=0.2,
    mu=6000,
    lam=24000,
    force_magnitude=-2.7,
    num_steps=15,
    printing=False,
):
    par = Parameter(0)
    force = CoefficientFunction((0, force_magnitude, 0))
    I3 = Id(3)

    fes = VectorH1(
        mesh,
        order=order,
        dirichletx="righttop|rightbot",
        dirichlety="rightbot",
        dirichletz="rightbot|back",
    )
    gf_solution = GridFunction(fes, name="standard shell")
    u = fes.TrialFunction()
    F = Grad(u) + I3

    a = BilinearForm(fes, symmetric=True, condense=True)
    a += Variation(
        CylindricalShellNeoHooke(mu, lam, F) * dx
        - par * force * u * ds("righttop")
    ).Compile()

    history = GridFunction(fes, multidim=0)
    for step in range(num_steps):
        print(f"Step {step+1}/{num_steps}")
        par.Set((step + 1) / num_steps)
        solvers.Newton(
            a,
            gf_solution,
            maxerr=1e-8,
            dampfactor=0.4,
            maxit=20,
            printing=printing,
            inverse="sparsecholesky",
        )
        history.AddMultiDimComponent(gf_solution.vec)

    return {
        "name": "standard",
        "u": gf_solution,
        "history": history,
        "deflection": CylindricalShellDeflection(gf_solution, mesh, thickness),
        "dofs": sum(fes.FreeDofs()),
    }


def SolveCylindricalShellTDNNS(
    mesh,
    order=2,
    thickness=0.2,
    mu=6000,
    lam=24000,
    force_magnitude=-2.7,
    num_steps=15,
    printing=False,
):
    par = Parameter(0)
    force = CoefficientFunction((0, force_magnitude, 0))
    I3 = Id(3)
    n = specialcf.normal(3)

    V = HCurl(mesh, order=order, dirichlet="rightbot")
    Sigma = Discontinuous(HDivDiv(mesh, order=order))
    Gamma = Discontinuous(HCurlCurl(mesh, order=order, orderinner=order + 1))
    Hyb = NormalFacetFESpace(mesh, order=order, dirichlet="righttop|rightbot|back")
    fes = V * Sigma * Gamma * Hyb
    gf_solution = GridFunction(fes, name="TDNNS shell")

    u, Psym, Fsym, uhat = fes.TrialFunction()
    F = Fsym + CylindricalShellSkew(curl(u))

    a = BilinearForm(fes, symmetric=True, symmetric_storage=True, condense=True)
    a += Variation(CylindricalShellNeoHooke(mu, lam, F) * dx).Compile()
    a += Variation(InnerProduct(Fsym - Grad(u) - I3, Psym) * dx).Compile()
    a += Variation(
        (u - uhat) * n * (Psym * n) * n * dx(element_boundary=True, bonus_intorder=0)
    ).Compile(False)
    a += Variation(-par * (force * uhat.Trace() + force * u.Trace()) * ds("righttop")).Compile()

    gf_u, gf_Psym, gf_Fsym, gf_uhat = gf_solution.components
    gf_Fsym.Set(I3)

    history = GridFunction(fes, multidim=0)
    for step in range(num_steps):
        print(f"Step {step+1}/{num_steps}")
        par.Set((step + 1) / num_steps)
        solvers.Newton(
            a,
            gf_solution,
            maxerr=1e-6,
            dampfactor=0.4,
            maxit=12,
            printing=printing,
            inverse="sparsecholesky",
        )
        history.AddMultiDimComponent(gf_solution.vec)

    return {
        "name": "TDNNS",
        "u": gf_u,
        "F": gf_Fsym + CylindricalShellSkew(curl(gf_u)),
        "history": history,
        "deflection": CylindricalShellDeflection(gf_u, mesh, thickness),
        "dofs": sum(fes.FreeDofs(True)),
    }

shell_steps = 15
SetHeapSize(1000 * 1000 * 1000)

with TaskManager():
    shell_standard = SolveCylindricalShellStandard(
        cyl_shell_mesh,
        num_steps=shell_steps,
    )
    shell_tdnns = SolveCylindricalShellTDNNS(
        cyl_shell_mesh,
        num_steps=shell_steps,
    )

for result in [shell_standard, shell_tdnns]:
    print(
        f"{result['name']:>8s}: "
        f"u_y(A) = {result['deflection']: .5e}"
    )

Step 1/15
Step 2/15
Step 3/15

Step 4/15
Step 5/15
Step 6/15

Step 7/15
Step 8/15
Step 9/15

Step 10/15
Step 11/15
Step 12/15

Step 13/15
Step 14/15
Step 15/15

WARNING: kwarg 'orderinner' is an undocumented flags option for class <class 'ngsolve.comp.HCurlCurl'>, maybe there is a typo?

Step 1/15

Step 2/15

Step 3/15

Step 4/15

Step 5/15

Step 6/15

Step 7/15

Step 8/15

Step 9/15

Step 10/15

Step 11/15

Step 12/15

Step 13/15

Step 14/15

Step 15/15

standard: u_y(A) = -9.75124e+00
   TDNNS: u_y(A) = -1.47853e+01

Draw(shell_standard["u"], cyl_shell_mesh, "standard shell", deformation=shell_standard["u"]);
Draw(BoundaryFromVolumeCF(shell_tdnns["u"]), cyl_shell_mesh, "TDNNS shell", deformation=BoundaryFromVolumeCF(shell_tdnns["u"]));

18.7. Nonlinear MCS#

The MCS extension follows the same Hu-Washizu idea but uses different spaces motivated by the linear MCS method. The displacement is normal-continuous, while the stress-like variable carries the tangential-normal continuity needed for the weak gradient pairing, [Fu, Neunteufel, Schöberl, Zdunek. A four-field mixed formulation for incompressible finite elasticity. Computer Methods in Applied Mechanics and Engineering, (2025).].

The nonlinear MCS formulation introduces a lifted deformation gradient and a pressure variable. A typical incompressible neo-Hookean density is

\[ \Psi(F,p) = \frac{\mu}{2}\left(\operatorname{tr}(F^TF)-d\right) -p\left(\det F-1\right), \]

where \(p\) enforces the incompressibility constraint \(\det F=1\). In the MCS formulation the lifted gradient is the regular representative used in the determinant, while the relation to \(u\) is enforced weakly through the MCS pairing.

We have the following fields:

unknown	possible space	meaning
\(u\)	`HDiv`	normal-continuous displacement
\(P\)	`HCurlDiv`-type stress space	stress-like multiplier
\(F\)	discontinuous `HCurlDiv`/matrix field	lifted deformation-gradient part
\(\widehat u\)	tangential facet space	hybrid tangential displacement trace
\(p\)	`L2`	pressure / incompressibility multiplier

The practical robustness message is different from TDNNS: the MCS formulation is aimed at nearly or fully incompressible nonlinear elasticity, i.e. \(J=\det F\approx 1\).

18.8. Compressed block#

The geometry is a unit cube with a smaller upper corner block attached on the region \(x,y\geq 0.5\). The load acts downward on the small top face. We compare the nonlinear MCS method with a standard displacement-pressure formulation.

from ngsolve import *
from ngsolve.webgui import Draw
import ipywidgets as widgets
from netgen.occ import *


def GenerateCompressedBlockMesh(h0=0.55, h1=0.35, h2=0.35, h3=0.35):
    box0 = Box(Pnt(0, 0, 0), Pnt(1, 1, 1))
    box1 = Box(Pnt(0.5, 0.5, 0), Pnt(1, 1, 1))
    left = box0 - box1
    left.faces.Max(Z).name = "top"
    left.faces.Min(Z).name = "bottom"
    left.faces.Min(X).name = "left"
    left.faces.Max(X).name = "right"
    left.faces.Min(Y).name = "back"
    left.faces.Max(Y).name = "front"

    box1.faces.Max(Z).name = "top1"
    box1.faces.Min(Z).name = "bottom"
    box1.faces.Min(X).name = "interface"
    box1.faces.Max(X).name = "right"
    box1.faces.Min(Y).name = "interface"
    box1.faces.Max(Y).name = "front"

    box1.solids.maxh = h1
    box1.faces.Max(Z).maxh = h2
    box1.faces.Max(X).maxh = h3
    box1.faces.Max(Y).maxh = h3

    return Mesh(OCCGeometry(Glue([left, box1])).GenerateMesh(maxh=h0))


def NeoHookeIncompressible(mu, F, p, mesh):
    return mu / 2 * (Trace(F.trans * F) - mesh.dim) - p * (Det(F) - 1)


def ProjectJacobian(jac_cf, mesh):
    gf_jac = GridFunction(L2(mesh, order=0), name="Jacobian")
    gf_jac.Set(jac_cf)
    return gf_jac, min(gf_jac.vec), max(gf_jac.vec)


def PointValue(cf, mesh, point=(1, 1, 1)):
    return cf(mesh(*point))


block_mesh = GenerateCompressedBlockMesh()
Draw(block_mesh);

def SolveCompressedBlockMCS(
    mesh,
    order=2,
    mu=8.0194,
    force0=-40,
    stabD=800,
    num_steps=8,
    maxerr=1e-7,
    printing=False,
):
    par = Parameter(0)
    force = CoefficientFunction((0, 0, force0))
    n = specialcf.normal(mesh.dim)
    h = specialcf.mesh_size

    V = HDiv(mesh, order=order, RT=True, dirichlet="bottom|right|front")
    Sigma = Discontinuous(HCurlDiv(mesh, order=order, ordertrace=-1))
    Gamma = Discontinuous(HCurlDiv(mesh, order=order, ordertrace=-1))
    Hyb = TangentialFacetFESpace(mesh, order=order, dirichlet="top|top1")
    Q = L2(mesh, order=order, lowest_order_wb=False)
    fes = V * Sigma * Gamma * Hyb * Q
    gf_solution = GridFunction(fes, name="MCS")

    u, P, F, uhat, p = fes.TrialFunction()

    def tang(v):
        return v - (v * n) * n

    Fmat = F.trans + (div(u) / mesh.dim + 1) * Id(mesh.dim)

    stab0 = GridFunction(L2(mesh, order=0))
    stiff_region = IfPos(x - 0.4, 1, 0) * IfPos(y - 0.4, 1, 0)
    stab0.Set(stiff_region * mu * stabD / h + (1 - stiff_region) * mu * stabD)
    
    a = BilinearForm(fes, symmetric=True, symmetric_storage=True, condense=True)
    a += Variation(
        (
            NeoHookeIncompressible(mu, Fmat, p, mesh)
            - 0.5 * mu * 1e-6 * p**2
            - InnerProduct(F.trans - Grad(u), P)
        )
        * dx
        + (0.5 * stab0 * tang(u - uhat) ** 2 - tang(u - uhat) * tang(P * n))
        * dx(element_boundary=True)
        - par * force * (u.Trace() + tang(uhat.Trace())) * ds("top1")
    ).Compile()

    gf_u, gf_P, gf_F, gf_uhat, gf_p = gf_solution.components
    gf_history = GridFunction(fes, multidim=0)

    for step in range(num_steps):
        par.Set((step + 1) / num_steps)
        solvers.Newton(
            a,
            gf_solution,
            maxerr=maxerr,
            maxit=30,
            printing=printing,
            inverse="sparsecholesky",
        )
        gf_history.AddMultiDimComponent(gf_solution.vec)

    Fh = gf_F.trans + (div(gf_u) / mesh.dim + 1) * Id(mesh.dim)
    gf_jac, jac_min, jac_max = ProjectJacobian(Det(Fh), mesh)
    return {
        "name": "MCS",
        "fes": fes,
        "solution": gf_solution,
        "u": gf_u,
        "p": gf_p,
        "F": Fh,
        "J": gf_jac,
        "J_min": jac_min,
        "J_max": jac_max,
        "history": gf_history,
        "point_u": PointValue(gf_u, mesh),
    }


def SolveCompressedBlockStandard(
    mesh,
    order=2,
    mu=8.0194,
    force0=-40,
    num_steps=8,
    maxerr=1e-7,
    printing=False,
):
    par = Parameter(0)
    force = CoefficientFunction((0, 0, force0))

    V = VectorH1(
        mesh,
        order=order,
        dirichletx="right|top|top1",
        dirichlety="front|top|top1",
        dirichletz="bottom",
    )
    Q = H1(mesh, order=max(order - 1, 1))
    fes = V * Q
    gf_solution = GridFunction(fes, name="standard")

    u, p = fes.TrialFunction()
    Fmat = Grad(u) + Id(mesh.dim)

    a = BilinearForm(fes, symmetric=True)
    a += Variation(
        (NeoHookeIncompressible(mu, Fmat, p, mesh) - 1e-8 * p**2) * dx
        - par * force * u * ds("top1")
    ).Compile()

    gf_u, gf_p = gf_solution.components
    gf_history = GridFunction(fes, multidim=0)

    for step in range(num_steps):
        par.Set((step + 1) / num_steps)
        solvers.Newton(
            a,
            gf_solution,
            maxerr=maxerr,
            maxit=30,
            printing=printing,
            inverse="sparsecholesky",
        )
        gf_history.AddMultiDimComponent(gf_solution.vec)

    Fh = Grad(gf_u) + Id(mesh.dim)
    gf_jac, jac_min, jac_max = ProjectJacobian(Det(Fh), mesh)
    return {
        "name": "standard",
        "fes": fes,
        "solution": gf_solution,
        "u": gf_u,
        "p": gf_p,
        "F": Fh,
        "J": gf_jac,
        "J_min": jac_min,
        "J_max": jac_max,
        "history": gf_history,
        "point_u": PointValue(gf_u, mesh),
    }

block_order = 2
block_force = -40
block_steps = 8

with TaskManager():
    block_mcs = SolveCompressedBlockMCS(
        block_mesh,
        order=block_order,
        force0=block_force,
        num_steps=block_steps,
    )
    block_std = SolveCompressedBlockStandard(
        block_mesh,
        order=block_order,
        force0=block_force,
        num_steps=block_steps,
    )

for result in [block_mcs, block_std]:
    print(
        f"{result['name']:>8s}: "
        f"u_z(1,1,1) = {result['point_u'][2]: .5e}, "
        f"min J = {result['J_min']: .5e}, max J = {result['J_max']: .5e}"
    )

     MCS: u_z(1,1,1) = -7.45066e-01, min J =  9.93900e-01, max J =  1.00414e+00
standard: u_z(1,1,1) = -6.71251e-01, min J = -7.99520e-02, max J =  1.59053e+00

block_mcs_u_bnd = BoundaryFromVolumeCF(block_mcs["u"])

Draw(block_mcs_u_bnd, block_mesh, "MCS displacement", deformation=block_mcs_u_bnd);
Draw(block_std["u"], block_mesh, "standard displacement", deformation=block_std["u"]);