# Clamped semi-cylindrical Naghdi shell under point load

Authors: Tian Yang (FEniCSx-Shells), Matteo Brunetti (FEniCS-Shells)

+++

This demo program solves the nonlinear Naghdi shell equations for a
semi-cylindrical shell loaded by a point force.

This problem is a standard reference for testing shell finite element
formulations, see [1]. The numerical locking issue is cured using enriched
finite element including cubic bubble shape functions and Partial Selective
Reduced Integration (PRSI) [2].

It is assumed the reader understands most of the basic functionality of the
new FEniCSx Project.

This demo then illustrates how to:

- Define and solve a nonlinear Naghdi shell problem with a curved stress-free
configuration given as analytical expression in terms of two curvilinear
coordinates.
- Use the PSRI approach to simultaneously cure shear- and membrane-locking
issues.

We begin by importing the necessary functionality from DOLFINx, UFL and
PETSc.

```python
import typing
from pathlib import Path
```

```python
from mpi4py import MPI
from petsc4py import PETSc

import matplotlib.pyplot as plt
import numpy as np

import dolfinx
import ufl
```

```python
from basix.ufl import blocked_element, element, enriched_element, mixed_element
from dolfinx.fem import Expression, Function, dirichletbc, functionspace, locate_dofs_topological
from dolfinx.fem.bcs import DirichletBC
from dolfinx.fem.function import Function as _Function
from dolfinx.fem.petsc import NonlinearProblem, apply_lifting, assemble_vector, set_bc
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from dolfinx.nls.petsc import NewtonSolver
from ufl import grad, inner, split
```

We consider a semi-cylindrical shell of radius $r$ and axis length $L$. The
shell is made of a linear elastic isotropic homogeneous material with Young
modulus $E$ and Poisson ratio $\nu$. The (uniform) shell thickness is denoted
by $t$. The Lamé moduli $\lambda$, $\mu$ are introduced to write later the 2D
constitutive equation in plane-stress:

```python
r = 1.016
L = 3.048
E, nu = 2.0685e7, 0.3
mu = E / (2.0 * (1.0 + nu))
lmbda = 2.0 * mu * nu / (1.0 - 2.0 * nu)
t = 0.03
```

The midplane of the initial (stress-free) configuration $\vec{\phi_0}$ of the
shell is given in the form of an analytical expression:

$$
\vec{\phi}_0(\xi_1, \xi_2) \subset \mathbb{R}³
$$

where $\xi_1 \in [-\pi/2, \pi/2]$ and $\xi_2 \in [0, L]$ are the curvilinear
coordinates. In this case, they represent the angular and axial coordinates,
respectively.

+++

We generate a mesh in the $(\xi_1, \xi_2)$ space with quadrilateral cells

```python
mesh = create_rectangle(
    MPI.COMM_WORLD, np.array([[-np.pi / 2, 0], [np.pi / 2, L]]), [20, 20], CellType.triangle
)
tdim = mesh.topology.dim  # = 2
```

We provide the analytical expression of the initial shape as a `ufl`
expression

```python
x = ufl.SpatialCoordinate(mesh)
phi0_ufl = ufl.as_vector([r * ufl.sin(x[0]), x[1], r * ufl.cos(x[0])])
```

Given the analytical expression of midplane, we define the unit normal as
below:

$$
\vec{n}  = \frac{\partial_1 \phi_0 \times \partial_2 \phi_0}{\| \partial_1 \phi_0 \times
\partial_2 \phi_0 \|}
$$

```python
def unit_normal(phi):
    n = ufl.cross(phi.dx(0), phi.dx(1))
    return n / ufl.sqrt(inner(n, n))


n0_ufl = unit_normal(phi0_ufl)
```

We define a local orthonormal frame $\{\vec{t}_{01}, \vec{t}_{02}, \vec{n}\}$
of the initial configuration $\phi_0$ by rotating the global Cartesian basis
$\vec{e}_i$ with a rotation matrix $\mathbf{R}_0$:

$$
\vec{t}_{0i} = \mathbf{R}_0 \vec{e}_i , \quad \vec{n} = \vec{t}_{03},
$$

A convienient choice of $\vec{t}_{01}$ and $\vec{t}_{02}$ (when $\vec{n}
\nparallel \vec{e}_2 $) could be:

$$
\vec{t}_{01} = \frac{\vec{e}_2 \times \vec{n}}{\| \vec{e}_2 \times \vec{n}\|} \\
\vec{t}_{02} =   \vec{n} \times \vec{t}_{01}
$$

The corresponding rotation matrix $\mathbf{R}_0$:

$$
\mathbf{R}_0 = [\vec{t}_{01}; \vec{t}_{02}; \vec{n}]
$$

```python
def tangent_1(n):
    e2 = ufl.as_vector([0, 1, 0])
    t1 = ufl.cross(e2, n)
    t1 = t1 / ufl.sqrt(inner(t1, t1))
    return t1


def tangent_2(n, t1):
    t2 = ufl.cross(n, t1)
    t2 = t2 / ufl.sqrt(inner(t2, t2))
    return t2


# The analytical expression of t1 and t2
t1_ufl = tangent_1(n0_ufl)
t2_ufl = tangent_2(n0_ufl, t1_ufl)


# The analytical expression of R0
def rotation_matrix(t1, t2, n):
    R = ufl.as_matrix([[t1[0], t2[0], n[0]], [t1[1], t2[1], n[1]], [t1[2], t2[2], n[2]]])
    return R


R0_ufl = rotation_matrix(t1_ufl, t2_ufl, n0_ufl)
```

The kinematics of the Nadghi shell model is defined by the following vector
fields:
- $\vec{\phi}$: the position of the midplane in the deformed configuration,
or equivalently, the displacement $\vec{u} = \vec{\phi} - \vec{\phi}_0$
- $\vec{d}$: the director, a unit vector giving the orientation of fiber at
the midplane. (not necessarily normal to the midsplane because of shears)

According to [3], the director $\vec{d}$ in the deformed configuration can be
parameterized with two successive rotation angles $\theta_1, \theta_2$

$$
\vec{t}_i = \mathbf{R} \vec{e}_i, \quad \mathbf{R}  = \text{exp}[\theta_1
\hat{\mathbf{t}}_1] \text{exp}[\theta_2 \hat{\mathbf{t}}_{02}] \mathbf{R}_0
$$

The rotation matrix $\mathbf{R}$ represents three successive rotations:
- First one: the initial rotation matrix $\mathbf{R}_0$
- Second one :$\text{exp}[\theta_2 \hat{\mathbf{t}}_{02}]$ rotates a vector
about the axis $\vec{t}_{02}$ of $\theta_2$ angle;
- Third one : $\text{exp}[\theta_1 \hat{\mathbf{t}}_1]$ rotates a vector
about the axis $\vec{t}_{1}$ of $\theta_1$ angle, and $\vec{t}_1 = \text{exp}
[\theta_2\hat{\mathbf{t}}_{02}] \vec{t}_{01}$

The rotation matrix $\mathbf{R}$ on the other hand it is equivalent to rotate
around the fixed axis $\vec{e}_1$ and $\vec{e}_2$ (Proof see [3]):

$$
\mathbf{R} = \mathbf{R}_0 \text{exp}[\theta_2 \hat{\mathbf{e}}_{2}]
\text{exp}[\theta_1 \hat{\mathbf{e}}_1]
$$

Therefore, the director $\vec{d}$ is updated with $(\theta_1, \theta_2)$ by:

$$
\vec{d} =\mathbf{R} \vec{e}_3 = \mathbf{R}_0 \vec{\Lambda}_3, \quad \vec{\Lambda}_3 =
[\sin(\theta_2)\cos(\theta_1), -\sin(\theta_1), \cos(\theta_2)\cos(\theta_1)]^\text{T}
$$

Note: the above formular becomes singular when $\theta_1 = \pm \pi/2, ...$, (See Chapter 4.2.1 in
[3] for details)

```python
def director(R0, theta):
    """Updates the director with two successive elementary rotations"""
    Lm3 = ufl.as_vector(
        [
            ufl.sin(theta[1]) * ufl.cos(theta[0]),
            -ufl.sin(theta[0]),
            ufl.cos(theta[1]) * ufl.cos(theta[0]),
        ]
    )
    d = ufl.dot(R0, Lm3)
    return d
```

In our 5-parameter Naghdi shell model the configuration of the shell is
assigned by:
- the 3-component vector field $\vec{u}$ representing the
displacement with respect to the initial configuration $\vec{\phi}_0$
- the 2-component vector field $\vec{\theta}$ representing the angle
variation of the director $\vec{d}$ with respect to initial unit normal
$\vec{n}$

+++

Following [1], we use a $[P_2 + B_3]^3$ element for $\vec{u}$ and a $[P_2]^2$
element for $\vec{\theta}$ and collect them in the state vector $\vec{q} =
[\vec{u}, \vec{\theta}]$:

```python
cell = mesh.basix_cell()
P2 = element("Lagrange", cell, degree=2)
B3 = element("Bubble", cell, degree=3)
P2B3 = enriched_element([P2, B3])

naghdi_shell_element = mixed_element(
    [blocked_element(P2B3, shape=(3,)), blocked_element(P2, shape=(2,))]
)
naghdi_shell_FS = functionspace(mesh, naghdi_shell_element)
```

Then, we define `Function`, `TrialFunction` and `TestFunction` objects to
express the variational forms and we split the mixed function into two
subfunctions for displacement and rotation.

```python
q_func = Function(naghdi_shell_FS)  # current configuration
q_trial = ufl.TrialFunction(naghdi_shell_FS)
q_test = ufl.TestFunction(naghdi_shell_FS)

u_func, theta_func = split(q_func)  # current displacement and rotation
```

We calculate the deformation gradient and the first and second fundamental
forms:

- Deformation gradient $\mathbf{F}$

$$
\mathbf{F} = \nabla \vec{\phi} \quad  (F_{ij} = \frac{\partial \phi_i}{\partial \xi_j});
\quad \vec{\phi} =
\vec{\phi}_0 +
\vec{u} \quad i = 1,2,3; j = 1,2
$$

- Metric tensor $\mathbf{a} \in \mathbb{S}^2_+$ and curvature tensor
$\mathbf{b} \in \mathbb{S}^2$ (First and second fundamental form)

$$
\begin{aligned}
\mathbf{a} &= {\nabla \vec{\phi}} ^{T} \nabla \vec{\phi} \\
\mathbf{b} &= -\frac{1}{2}({\nabla \vec{\phi}} ^{T} \nabla \vec{d} + {\nabla \vec{d}} ^{T}
\nabla \vec{\phi})
\end{aligned}
$$

In the initial configuration, $\vec{d} = \vec{n}$, $\vec{\phi} =
\vec{\phi}_0$, the conresponding initial tensors are $\mathbf{a}_0$,
$\mathbf{b}_0$

```python
# Current deformation gradient
F = grad(u_func) + grad(phi0_ufl)

# current director
d = director(R0_ufl, theta_func)

# initial metric and curvature tensor a0 and b0
a0_ufl = grad(phi0_ufl).T * grad(phi0_ufl)
b0_ufl = -0.5 * (grad(phi0_ufl).T * grad(n0_ufl) + grad(n0_ufl).T * grad(phi0_ufl))
```

We define strain measures of the Naghdi shell model:
- Membrane strain tensor $\boldsymbol{\varepsilon}(\vec{u})$

$$
\boldsymbol{\varepsilon} (\vec{u})= \frac{1}{2} \left ( \mathbf{a}(\vec{u}) - \mathbf{a}_0 \right)
$$

- Bending strain tensor $\boldsymbol{\kappa}(\vec{u}, \vec{\theta})$

$$
\boldsymbol{\kappa}(\vec{u}, \vec{\theta}) = \mathbf{b}(\vec{u}, \vec{\theta}) - \mathbf{b}_0
$$

- transverse shear strain vector $\vec{\gamma}(\vec{u}, \vec{\theta})$

$$
\begin{aligned}
\vec{\gamma}(\vec{u}, \vec{\theta}) & = {\nabla \vec{\phi}(\vec{u})}^T \vec{d}(\vec{\theta})
- {\nabla\vec{\phi}_0}^T
\vec{n} \\
& = {\nabla \vec{\phi}(\vec{u})}^T \vec{d}(\vec{\theta}) \quad \text{if zero initial shears}
\end{aligned}
$$

```python
def epsilon(F):
    """Membrane strain"""
    return 0.5 * (F.T * F - a0_ufl)


def kappa(F, d):
    """Bending strain"""
    return -0.5 * (F.T * grad(d) + grad(d).T * F) - b0_ufl


def gamma(F, d):
    """Transverse shear strain"""
    return F.T * d
```

In curvilinear coordinates, the stiffness modulus of linear isotropic
material is defined as:

- Membrane and bending stiffness modulus $A^{\alpha\beta\sigma\tau}$,
$D^{\alpha\beta\sigma\tau}$
(contravariant components)

$$
\frac{A^{\alpha\beta\sigma\tau}}t=12\frac{D^{\alpha\beta\sigma\tau}}{t^3}=
\frac{2\lambda\mu}{\lambda+2\mu}
a_0^{\alpha\beta}a_0^{\sigma\tau}+\mu(a_0^{\alpha\sigma}a_0^{\beta\tau}+a_0^{\alpha\tau}
a_0^{\beta\sigma})
$$

- Shear stiffness modulus $S^{\alpha\beta}$ (contravariant components)

$$
\frac{S^{\alpha\beta}}t = \alpha_s \mu a_0^{\alpha\beta} , \quad \alpha_s =
\frac{5}{6}: \text{shear factor}
$$

where $a_0^{\alpha\beta}$ is the contravariant components of the initial
metric tensor $\mathbf{a}_0$

```python
a0_contra_ufl = ufl.inv(a0_ufl)
j0_ufl = ufl.det(a0_ufl)

i, j, l, m = ufl.indices(4)  # noqa: E741
A_contra_ufl = ufl.as_tensor(
    (
        ((2.0 * lmbda * mu) / (lmbda + 2.0 * mu)) * a0_contra_ufl[i, j] * a0_contra_ufl[l, m]
        + 1.0
        * mu
        * (a0_contra_ufl[i, l] * a0_contra_ufl[j, m] + a0_contra_ufl[i, m] * a0_contra_ufl[j, l])
    ),
    [i, j, l, m],
)
```

We define the resultant stress measures:

- Membrane stress tensor $\mathbf{N}$

$$
\mathbf{N} = \mathbf{A} : \boldsymbol{\varepsilon}
$$

- Bending stress tensor $\mathbf{M}$

$$
\mathbf{M} = \mathbf{D} : \boldsymbol{\kappa}
$$

- Shear stress vector $\vec{T}$

$$
\vec{T} = \mathbf{S} \cdot \vec{\gamma}
$$

```python
N = ufl.as_tensor(t * A_contra_ufl[i, j, l, m] * epsilon(F)[l, m], [i, j])
M = ufl.as_tensor((t**3 / 12.0) * A_contra_ufl[i, j, l, m] * kappa(F, d)[l, m], [i, j])
T = ufl.as_tensor((t * mu * 5.0 / 6.0) * a0_contra_ufl[i, j] * gamma(F, d)[j], [i])
```

We define elastic strain energy density $\psi_{m}$, $\psi_{b}$, $\psi_{s}$ for membrane, bending
and shear,
respectively.

$$
\psi_m = \frac{1}{2} \mathbf{N} : \boldsymbol{\varepsilon}; \quad
\psi_b = \frac{1}{2} \mathbf{M} : \boldsymbol{\kappa}; \quad
\psi_s = \frac{1}{2} \vec{T} \cdot \vec{\gamma}
$$

They are per unit surface in the initial configuration:

```python
psi_m = 0.5 * inner(N, epsilon(F))
psi_b = 0.5 * inner(M, kappa(F, d))
psi_s = 0.5 * inner(T, gamma(F, d))
```

Shear and membrane locking is treated using the partial reduced selective
integration proposed in Arnold and Brezzi [2].

We introduce a parameter $\alpha \in \mathbb{R}$ that splits the membrane and
shear energy in the energy functional into a weighted sum of two parts:

$$
\begin{aligned}\Pi_{N}(u,\theta)&=\Pi^b(u_h,\theta_h)+\alpha\Pi^m(u_h)+(1-\alpha)\Pi^m(u_h)\\&+
\alpha\Pi^s(u_h,\theta_h)
+(1-\alpha)\Pi^s(u_h,\theta_h)-W_{\mathrm{ext}},\end{aligned}
$$

We apply reduced integration to the parts weighted by the factor $(1-\alpha)$

More details:
- Optimal choice $\alpha = \frac{t^2}{h^2}$, $h$ is the diameter of the cell
- Full integration : Gauss quadrature of degree 4 (6 integral points for triangle)
- Reduced integration : Gauss quadrature of degree 2 (3 integral points for triangle).
- While [1] suggests a 1-point reduced integration, we observed that this
  leads to spurious modes in the present case.

```python
# Full integration of order 4
dx_f = ufl.Measure("dx", domain=mesh, metadata={"quadrature_degree": 4})

# Reduced integration of order 2
dx_r = ufl.Measure("dx", domain=mesh, metadata={"quadrature_degree": 2})

# Calculate the factor alpha as a function of the mesh size h
h = ufl.CellDiameter(mesh)
alpha_FS = functionspace(mesh, element("DG", cell, 0))
alpha_expr = Expression(t**2 / h**2, alpha_FS.element.interpolation_points)
alpha = Function(alpha_FS)
alpha.interpolate(alpha_expr)

# Full integration part of the total elastic energy
Pi_PSRI = psi_b * ufl.sqrt(j0_ufl) * dx_f
Pi_PSRI += alpha * psi_m * ufl.sqrt(j0_ufl) * dx_f
Pi_PSRI += alpha * psi_s * ufl.sqrt(j0_ufl) * dx_f

# Reduced integration part of the total elastic energy
Pi_PSRI += (1.0 - alpha) * psi_m * ufl.sqrt(j0_ufl) * dx_r
Pi_PSRI += (1.0 - alpha) * psi_s * ufl.sqrt(j0_ufl) * dx_r

# External work part (zero in this case)
W_ext = 0.0
Pi_PSRI -= W_ext
```

The residual and jacobian are the first and second order derivatives of the
total potential energy, respectively

```python
F = ufl.derivative(Pi_PSRI, q_func, q_test)
J = ufl.derivative(F, q_func, q_trial)
```

Next, we prescribe the dirichlet boundary conditions:
- fully clamped boundary conditions on the top boundary ($\xi_2 = 0$):
- $u_{1,2,3} = \theta_{1,2} = 0$

```python
def clamped_boundary(x):
    return np.isclose(x[1], 0.0)


fdim = tdim - 1
clamped_facets = locate_entities_boundary(mesh, fdim, clamped_boundary)

u_FS, _ = naghdi_shell_FS.sub(0).collapse()
theta_FS, _ = naghdi_shell_FS.sub(1).collapse()

# u1, u2, u3 = 0 on the clamped boundary
u_clamped = Function(u_FS)
clamped_dofs_u = locate_dofs_topological((naghdi_shell_FS.sub(0), u_FS), fdim, clamped_facets)
bc_clamped_u = dirichletbc(u_clamped, clamped_dofs_u, naghdi_shell_FS.sub(0))

# theta1, theta2 = 0 on the clamped boundary
theta_clamped = Function(theta_FS)
clamped_dofs_theta = locate_dofs_topological(
    (naghdi_shell_FS.sub(1), theta_FS), fdim, clamped_facets
)
bc_clamped_theta = dirichletbc(theta_clamped, clamped_dofs_theta, naghdi_shell_FS.sub(1))
```

- symmetry boundary conditions on the left and right side ($\xi_1 = \pm
\pi/2$):
- $u_3 = \theta_2 = 0$

```python
def symm_boundary(x):
    return np.isclose(abs(x[0]), np.pi / 2)


symm_facets = locate_entities_boundary(mesh, fdim, symm_boundary)

symm_dofs_u = locate_dofs_topological(
    (naghdi_shell_FS.sub(0).sub(2), u_FS.sub(2)), fdim, symm_facets
)
bc_symm_u = dirichletbc(u_clamped, symm_dofs_u, naghdi_shell_FS.sub(0).sub(2))

symm_dofs_theta = locate_dofs_topological(
    (naghdi_shell_FS.sub(1).sub(1), theta_FS.sub(1)), fdim, symm_facets
)
bc_symm_theta = dirichletbc(theta_clamped, symm_dofs_theta, naghdi_shell_FS.sub(1).sub(1))

bcs = [bc_clamped_u, bc_clamped_theta, bc_symm_u, bc_symm_theta]
```

The loading is exerted by a point force along the $z$ direction applied at
the midpoint of the bottom boundary.
Since `PointSource` function is not available by far in new FEniCSx, we
achieve the same functionality according to the method detailed in [4].

```python
def compute_cell_contributions(V, points):
    """Returns the cell containing points and the values of the basis functions
    at that point"""
    # Determine what process owns a point and what cells it lies within
    mesh = V.mesh
    point_ownership_data = dolfinx.cpp.geometry.determine_point_ownership(
        mesh._cpp_object, points, 1e-6
    )

    owning_points = np.asarray(point_ownership_data.dest_points).reshape(-1, 3)
    cells = point_ownership_data.dest_cells

    # Pull owning points back to reference cell
    mesh_nodes = mesh.geometry.x
    cmap = mesh.geometry.cmap
    ref_x = np.zeros((len(cells), mesh.geometry.dim), dtype=mesh.geometry.x.dtype)
    for i, (point, cell) in enumerate(zip(owning_points, cells)):
        geom_dofs = mesh.geometry.dofmap[cell]
        ref_x[i] = cmap.pull_back(point.reshape(-1, 3), mesh_nodes[geom_dofs])

    # Create expression evaluating a trial function (i.e. just the basis function)
    u = ufl.TrialFunction(V.sub(0).sub(2))
    num_dofs = V.sub(0).sub(2).dofmap.dof_layout.num_dofs * V.sub(0).sub(2).dofmap.bs
    if len(cells) > 0:
        # NOTE: Expression lives on only this communicator rank
        expr = dolfinx.fem.Expression(u, ref_x, comm=MPI.COMM_SELF)
        values = expr.eval(mesh, np.asarray(cells, dtype=np.int32))

        # Strip out basis function values per cell
        basis_values = values[: num_dofs : num_dofs * len(cells)]
    else:
        basis_values = np.zeros((0, num_dofs), dtype=dolfinx.default_scalar_type)
    return cells, basis_values
```

```python
# Point source
if mesh.comm.rank == 0:
    points = np.array([[0.0, L, 0.0]], dtype=mesh.geometry.x.dtype)
else:
    points = np.zeros((0, 3), dtype=mesh.geometry.x.dtype)

cells, basis_values = compute_cell_contributions(naghdi_shell_FS, points)
```

We define a custom `NonlinearProblem` which is able to include the point
force.

```python
class NonlinearProblemPointSource(NonlinearProblem):
    def __init__(
        self,
        F: ufl.form.Form,
        u: _Function,
        bcs: typing.List[DirichletBC] = [],  # noqa: UP006
        J: ufl.form.Form = None,
        cells=[],
        basis_values=[],
        PS: float = 0.0,
    ):
        super().__init__(F, u, bcs, J)
        self.PS = PS
        self.cells = cells
        self.basis_values = basis_values
        self.function_space = u.function_space

    def F(self, x: PETSc.Vec, b: PETSc.Vec) -> None:
        with b.localForm() as b_local:
            b_local.set(0.0)
        assemble_vector(b, self._L)

        # Add point source
        if len(self.cells) > 0:
            for cell, basis_value in zip(self.cells, self.basis_values):
                dofs = self.function_space.sub(0).sub(2).dofmap.cell_dofs(cell)
                with b.localForm() as b_local:
                    b_local.setValuesLocal(
                        dofs, basis_value * self.PS, addv=PETSc.InsertMode.ADD_VALUES
                    )

        apply_lifting(b, [self._a], bcs=[self.bcs], x0=[x], alpha=-1.0)
        b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        set_bc(b, self.bcs, x, -1.0)
```

We use the standard Newton iteration.

```python
problem = NonlinearProblemPointSource(F, q_func, bcs, J, cells, basis_values)

solver = NewtonSolver(mesh.comm, problem)

# Set Newton solver options
solver.rtol = 1e-6
solver.atol = 1e-6
solver.max_it = 30
solver.convergence_criterion = "incremental"
solver.report = True

# Modify the linear solver in each Newton iteration
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
```

Finally, we can solve the quasi-static problem, incrementally increasing the
loading from $0$N to $2000$N

```python
PS_diff = 50.0
n_step = 40

# Store the displacement at the point load
if mesh.comm.rank == 0:
    u3_list = np.zeros(n_step + 1)
    PS_list = np.arange(0, PS_diff * (n_step + 1), PS_diff)

q_func.x.array[:] = 0.0

bb_point = np.array([[0.0, L, 0.0]], dtype=np.float64)

for i in range(1, n_step + 1):
    problem.PS = PS_diff * i
    n, converged = solver.solve(q_func)
    assert converged
    q_func.x.scatter_forward()
    if mesh.comm.rank == 0:
        print(f"Load step {i:d}, Number of iterations: {n:d}, Load: {problem.PS:.2f}", flush=True)
    # Calculate u3 at the point load
    u3_bb = None
    u3_func = q_func.sub(0).sub(2).collapse()
    if len(cells) > 0:
        u3_bb = u3_func.eval(bb_point, cells[0])[0]
    u3_bb = mesh.comm.gather(u3_bb, root=0)
    if mesh.comm.rank == 0:
        for u3 in u3_bb:
            if u3 is not None:
                u3_list[i] = u3
                break
```

We write the outputs of $\vec{u}$, $\vec{\theta}$, and $\vec{\phi}$ in the
second order Lagrange space.

```python
# Interpolate phi_ufl into CG2 Space
u_P2B3 = q_func.sub(0).collapse()
theta_P2 = q_func.sub(1).collapse()

# Interpolate phi in the [P2]³ Space
phi_FS = functionspace(mesh, blocked_element(P2, shape=(3,)))
phi_expr = Expression(phi0_ufl + u_P2B3, phi_FS.element.interpolation_points)
phi_func = Function(phi_FS)
phi_func.interpolate(phi_expr)

# Interpolate u in the [P2]³ Space
u_P2 = Function(phi_FS)
u_P2.interpolate(u_P2B3)

results_folder = Path("results/nonlinear_naghdi/semi_cylinder")
results_folder.mkdir(exist_ok=True, parents=True)

with dolfinx.io.VTXWriter(mesh.comm, results_folder / "u_naghdi.bp", [u_P2]) as vtx:
    vtx.write(0)

with dolfinx.io.VTXWriter(mesh.comm, results_folder / "theta_naghdi.bp", [theta_P2]) as vtx:
    vtx.write(0)

with dolfinx.io.VTXWriter(mesh.comm, results_folder / "phi_naghdi.bp", [phi_func]) as vtx:
    vtx.write(0)
```

The results for the transverse displacement at the point of application of
the force are validated against a standard reference from the literature,
obtained using Abaqus S4R element and a structured mesh of 40 times 40
elements, see [1]:

```python
if mesh.comm.rank == 0:
    fig = plt.figure()
    reference_u3 = 1.0e-2 * np.array(
        [
            0.0,
            5.421,
            16.1,
            22.195,
            27.657,
            32.7,
            37.582,
            42.633,
            48.537,
            56.355,
            66.410,
            79.810,
            94.669,
            113.704,
            124.751,
            132.653,
            138.920,
            144.185,
            148.770,
            152.863,
            156.584,
            160.015,
            163.211,
            166.200,
            168.973,
            171.505,
        ]
    )
    reference_P = 2000.0 * np.array(
        [
            0.0,
            0.05,
            0.1,
            0.125,
            0.15,
            0.175,
            0.2,
            0.225,
            0.25,
            0.275,
            0.3,
            0.325,
            0.35,
            0.4,
            0.45,
            0.5,
            0.55,
            0.6,
            0.65,
            0.7,
            0.75,
            0.8,
            0.85,
            0.9,
            0.95,
            1.0,
        ]
    )
    plt.plot(-u3_list, PS_list, label="FEniCSx-Shells 20 x 20")
    plt.plot(reference_u3, reference_P, "or", label="Sze (Abaqus S4R)")
    plt.xlabel("Displacement (mm)")
    plt.ylabel("Load (N)")
    plt.legend()
    plt.grid()
    plt.tight_layout()
    plt.savefig(results_folder / "comparisons.pdf")
```

References:

[1] K. Sze, X. Liu, and S. Lo. Popular benchmark problems for geometric
nonlinear analysis of shells. Finite Elements in Analysis and Design,
40(11):1551 - 1569, 2004.

[2] D. Arnold and F.Brezzi, Mathematics of Computation, 66(217): 1-14, 1997.
https://www.ima.umn.edu/~arnold//papers/shellelt.pdf

[3] P. Betsch, A. Menzel, and E. Stein. On the parametrization of finite
rotations in computational mechanics: A classification of concepts with
application to smooth shells. Computer Methods in Applied Mechanics and
Engineering, 155(3):273 - 305, 1998.

[4] https://fenicsproject.discourse.group/t/point-sources-redux/13496/4