Hi all,

I started to use NGSolve for solving some problems in modeling hyperelastic materials. I followed the examples here closely.

The example I want to calculate is the inflation of an ellipsoidal structure.

I generated my mesh as follows (units are in mm)

```
from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw
from ngsolve.solvers import NewtonMinimization
ngsglobals.msg_level = 1
import matplotlib.pyplot as plt
from scipy.stats import qmc
import scipy as sp
import numpy as np
np.random.seed(137)
rs_endo = 2.5e-2 * 10 ** 3
rl_endo = 9.0e-2 * 10 ** 3
rs_epi = 3.5e-2 * 10 ** 3
rl_epi = 9.7e-2 * 10 ** 3
ax = Axes((1,0,0), n=X, h=Y)
ax2 = Axes((1,0,0), n=X, h=Y)
maxh=3.0
ell_endo = Ellipsoid (Axis ((0,0,0), X), rl_endo, rs_endo)
ell_endo.faces.name = "endocardium"
ell_endo.faces.col = (0.5,0.5,0)
ell_epi = Ellipsoid (Axis ((0,0,0), X), rl_epi, rs_epi)
ell_epi.faces.name = "epicardium"
ell_epi.faces.col = (1,0,0)
cutplane=HalfSpace((-25,0,0), X)
cutplane.faces.name="top"
cutplane.faces.col=(0,0,1)
lv = (ell_epi-ell_endo)-cutplane
geo = OCCGeometry(lv)
mesh=Mesh(geo.GenerateMesh(maxh=maxh))
```

This works fine. Now I can replicate some of the features laid out in the elasticity examples of the NGS24 workshop. I adapted things for almost incompressible Elasticity

```
fes = VectorH1(mesh, order=1, dirichlet="top")
u, deltau = fes.TnT()
gfu = GridFunction(fes)
loadpar = Parameter(0)
normal = specialcf.normal(mesh.dim)
kPa2MPa = 1e-3
mu = 33.3 * kPa2MPa
kappa = 1000 * kPa2MPa
appl_pressure = Parameter(0.0) # in kPa
def psi(F):
C = F.trans * F
J = Det(F)
Cbar = J**(-2/3) * C
Wiso = mu/2*(Trace(Cbar)-3)
Wvol = kappa/2 * log(J)**2
return Wiso + Wvol
a = BilinearForm(fes, symmetric=False)
defgrad = Id(mesh.dim) + Grad(u)
# Energy volumetric
a += Variation(psi(defgrad)*dx)
# dead load (pressure acts against normal direction this swallows the minus sign
a += appl_pressure *kPa2MPa*InnerProduct(normal, deltau)*ds("endocardium")
def SolveNonlinearElast(fes, pressure, numloadsteps=10):
print("Solve Nonlinear Elast for p=", pressure)
# Solution
gfu = ng.GridFunction(fes)
gfu.vec[:] = 0
# Set params
incpr = pressure / numloadsteps
newval = 0.0
conv = 0
converged = True
with TaskManager():
for loadstep in range(numloadsteps):
print ("loadstep", loadstep+1, end='')
newval += incpr
print(" | val", newval)
appl_pressure .Set(newval)
i1, _ = NewtonMinimization(a, gfu, printing=True, linesearch=False, inverse="pardiso", maxerr=1e-6, maxit=10)
conv = conv + i1
if conv < 0:
converged = False
return gfu, converged
```

However, what I would actually like to model is a followerload which is defined as

P\boldsymbol{N} = -p_\mathrm{appl}JF^{-\top}\boldsymbol{N} \quad\quad\mathrm{on}~\Gamma_\mathrm{follow}

with

F = I + \nabla u,~J=\mathrm{det}(F),~P=\frac{\partial \Psi(F)}{\partial F} and N being the normal vector in the reference configuration (aka specialcf.normal).

The above example actually models

P\boldsymbol{N} = -p_\mathrm{appl}\boldsymbol{N} \quad\quad\mathrm{on}~\Gamma_\mathrm{follow}

Since follower loads don’t allow for a potential formulation I wanted to ask how I could introduce this here in my example?