Hello all,
I am trying to create an inflation simulation in NGSolve by minimizing the potential energy functional which is calculated by the stresses on the boundary - the work done by inflation, i.e. the functional and it’s derivative are given by
VEC = VectorH1(mesh, order=2, dim=mesh.dim)
p_cf = CoefficientFunction(p)
F = Id(mesh.dim) + Grad(u)
F_bnd = Id(mesh.dim) + Grad(u).Trace()
W = neo_hookean_energy_nonneg(F_bnd, mu, lmbda)
E_bnd = W*ds
E_vol = -p_cf*dx
S_bnd = E_bnd.DiffShape(X)
S_vol = E_vol.DiffShape(X)
S = S_bnd + S_vol
dJOmegaAuto = LinearForm(VEC)
dJOmegaAuto += S
This works without any errors.
However, when I try using this to compute the deformation on my test cube, I see ‘discontinuities’ in my mesh:
For the optimization, I am using the following code:
from math import isfinite
Qk = 1.0 # Zhang–Hager running weights
Ck = None # Zhang–Hager reference value (set on first call)
def zh_ls(Jold, g_vec, direction_gf, phi_gf, *, descent=True,
c1=1e-3, shrink=0.5, max_backtracks=30, move_frac=0.3,
alpha_hint=None, theta=0.85):
"""
Zhang-Hager non-monotone Armijo line search.
"""
global Qk, Ck
# --- boundary move-limit initial α ---
bm = float(Integrate(1, mesh, BND))
if bm <= 0:
return 0.0, Jold
rms_b = (Integrate(InnerProduct(direction_gf, direction_gf), mesh, BND)/bm + 1e-30)**0.5
if rms_b == 0.0:
return 0.0, Jold
alpha_ml = move_frac * hmax / rms_b
alpha = alpha_ml if alpha_hint is None else min(alpha_hint, alpha_ml)
# set reference Ck on first call
if Ck is None:
Ck = Jold
# directional derivative at 0
gTd = float(InnerProduct(g_vec, direction_gf.vec))
if not descent:
gTd = -gTd # reuse same inequality by flipping sign convention below
# --- snapshot & trial ---
phi_old = phi_gf.vec.CreateVector(); phi_old.data = phi_gf.vec
trial = phi_gf.vec.CreateVector()
def ok(Jnew, a):
rhs = Ck + c1 * a * ( -gTd if not descent else gTd )
return Jnew <= rhs
for _ in range(max_backtracks):
trial.data = phi_old
trial.data += alpha * direction_gf.vec
# commit trial deformation
phi_gf.vec.data = trial
mesh.SetDeformation(phi_gf)
# u := trial deformation for surface energy
gfu.vec.data = trial
Jnew = float(Integrate(energy_functional_nonneg(gfu), mesh))
if not isfinite(Jnew):
alpha *= shrink
continue
if ok(Jnew, alpha):
# update Zhang-Hager reference
Qk1 = theta*Qk + 1.0
Ck = (theta*Qk*Ck + Jnew) / Qk1
Qk = Qk1
return alpha, Jnew
alpha *= shrink
# revert on failure
phi_gf.vec.data = phi_old
mesh.SetDeformation(phi_gf)
gfu.vec.data = phi_old
return 0.0, Jold
The output rendered is generated by
fes = VectorH1(mesh, order=2, dirichlet=".*")
gfu = GridFunction(fes)
scene_u = Draw (gfu, mesh, "state")
and then updated as in the Zhang-Hager code snippet. Where could this error originate from? I am thinking that this might be related to the fact that I am using two different H1 spaces (VEC for the shape derivative and fes in the rendering snippet), but I wanted to make sure that I am not missing any crucial information regarding the inner workings of NGsolve (i.e. a rendering bug or some flag that I forgot to set).
Kind regards
