# Solving Hyperleasticity: Correctly defining boundaries

Hi everyone,

I was trying to translate the hyperelasticity example from FEniCS to ngsolve, following the nonlinear elasticity tutorial, and Navier Stokes tutorial to implement homogenous and non-homogenous dirichlet boundary conditions on different parts of the boundary.

However, it seems that the boundary conditions are not getting assigned properly (looking at the GUI). Would anyone be able to help me pin down the error? You can see from the norm of the error, that it does not converge.

``````import netgen.gui
from math import pi as myPi
from ngsolve import *
from netgen.csg import unit_cube

mesh = Mesh(unit_cube.GenerateMesh(maxh=0.1))
V = VectorH1(mesh, order=1, dirichlet="left|right")
u = GridFunction(V, str="u")
u.vec[:] = 0.
du, v = V.TnT()

# material parameters and forces
E, nu = 10.0, 0.3
mu = E / 2./(1+nu)
lmbda = 2 * mu * nu/(1-2*nu)
bodyForce = CoefficientFunction((0, -0.5, 0.))

# definition of the boundaries
scale = 0.5
yo, zo = 0.5, 0.5
thta = myPi/3.
uLeft = CoefficientFunction((0., 0., 0.))
uRight = CoefficientFunction((
0.,
scale * (yo + (y - yo) * cos(thta) - (z - zo) * sin(thta) - y),
scale * (zo + (y - yo) * sin(thta) + (z-zo) * cos(thta) - z)
))
u.Set(uLeft, definedon=mesh.Boundaries("left"))
u.Set(uRight, definedon=mesh.Boundaries("right"))

# Define strain measures
I = Id(mesh.dim)
I1 = Trace(F.trans * F)
J = Det(F)
psi = mu * (I1 - 3) - mu * log(J) + lmbda/2. * (J-1)**2

# definition of bilinear and linear forms
a = BilinearForm(V, symmetric=False)
a += Variation(psi.Compile() * dx)
a += Variation((-1 * InnerProduct(bodyForce, du)).Compile() * dx)

# create the residual and iteration vectors
res = u.vec.CreateVector()
w = u.vec.CreateVector()

resNorm = 1.0

while resNorm > 1.e-8:
a.Apply(u.vec, res)
a.AssembleLinearization(u.vec)
inv = a.mat.Inverse(V.FreeDofs())
w.data = inv * res
resNorm = InnerProduct(w, res)
print(f"Norm: {resNorm}")
u.vec.data -= w

Draw(u, mesh, "displacement")
SetVisualization(deformation=True)``````

gives

``````Norm: 6.686688205442815
Norm: 0.49337768606525656
Norm: 0.34295300814873925
Norm: 0.15001895765763176
Norm: 0.06517533253371934
Norm: 0.019404532237654374
Norm: 0.18895838479738908
Norm: 2.024789522475844
Norm: 5815.697445306668
Norm: 49695.17791867306
Norm: 58282.3479021152
Norm: 11584203088583.24
Norm: 2.3448080388890567e+22
Norm: 3.154397165784633e+26
Norm: 8.26906306359288e+25
Norm: 2.167685693089545e+25
Norm: 5.682457983539969e+24``````

u.Set(cf, definedon=…)
sets u to cf on the definedon part and to 0 everywhere else, set the boundary values like this:

``u.Set(mesh.BoundaryCF({"left" : uLeft, "right" : uRight}),  definedon=mesh.Boundaries("left|right"))``

Best
Christopher

Thanks Christopher!

Is this a newer feature? NGSolve-6.2.1909 does not seem to have it. I will try updating NGSolve to the latest version and test it again…

``AttributeError: 'ngsolve.comp.Mesh' object has no attribute 'BoundaryCF'``

Ok this works with the newer version. But even with the newer version (NGSolve-6.2.2009), I get a diverging Newton.

``````Norm: 2.6808208867801797
Norm: 577.2816648847577
Norm: 152009.87301735353
Norm: 7366462.124902223
Norm: 65416518034032.02
Norm: 5.554793589592823e+16
Norm: 1.4678499972939826e+16
Norm: 3848066201629095.0
Norm: 1008756719615877.0
Norm: 264477504711945.2``````

I will take a look at implementing non-homogenous BC’s more carefully.

After inspecting a bit more the location of the boundaries “left”, “right”, … I noticed that the boundaries that I need (x=0 and x=1) correspond to “back” and “front”. Indeed after replacing them accordingly and applying the boundary condition in steps, rather than at once, the solution seems to be reasonable till a point, and then it diverges:

``````for i in range(1, 11):
factor.Set(i/10)
u.Set(mesh.BoundaryCF({"back" : uLeft, "front" : uRight}),  definedon=mesh.Boundaries("back|front"))
print(i)
resNorm = 1.0
iter = 0
while resNorm > 1.e-8 and iter < maxiters:
a.Apply(u.vec, res)
a.AssembleLinearization(u.vec)
inv = a.mat.Inverse(V.FreeDofs())
w.data = inv * res
resNorm = sqrt(abs(InnerProduct(w, w)))
print(f"Norm: {resNorm}")
u.vec.data -= w
iter += 1
Draw(u, mesh, "displacement")
SetVisualization(deformation=True)``````

gives

``````1
Norm: 0.840274935028221
Norm: 0.010631338197297788
Norm: 2.840043781994656e-05
Norm: 5.625360578002826e-10
2
Norm: 1.3145809712371617
Norm: 0.03864670201207834
Norm: 0.0003690373463827511
Norm: 8.934665769419485e-08
Norm: 1.54388135479357e-14
3
Norm: 1.8523594043744418
Norm: 0.08405462602954474
Norm: 0.0018541634265310004
Norm: 1.9279392913312044e-06
Norm: 8.193550778213933e-12
4
Norm: 2.4123885921399872
Norm: 0.14726391288538263
Norm: 0.0060619644879045765
Norm: 1.9854480113520826e-05
Norm: 1.1251011825818185e-09
5
Norm: 2.816514170226768
Norm: 0.4758279175877349
Norm: 0.07470240485416736
Norm: 0.14433503658998675
Norm: 0.2576930420365274
Norm: 0.15947077282736843
Norm: 0.04193568718574943
Norm: 0.04151801466555532
Norm: 0.0044976089152604015
Norm: 5.306089515763066e-05
Norm: 8.324583871537812e-09
6
Norm: 4.623385880239417
Norm: 2.913358091896825
Norm: 2.617097160078447
Norm: 2.9116418166392988
Norm: 6.296120155081865
Norm: 12.355324294951478
Norm: 19.201276563870323
Norm: 58.62691253475249
Norm: 230.766474679662``````

I think, this is due to inaccurate approximation of the jacobian in automatic differentiation of the residual. I don’t see any other reason on why Newton wouldn’t converge after a point for this problem. I checked the deformation and it seems reasonable till the 6th loading step.

Hi bhaveshshrimali,

are you sure your material law is correct? For zero external forces and zero Dirichlet boundary condition the identity function is not the trivial solution. For the Neo-Hookean material law mu/2 should be the correct scaling for the first term:

``psi = mu/2 * (I1 - 3) - mu * log(J) + lmbda/2. * (J-1)**2``

By setting the boundary conditions for every loadstep you reset your solution rather than using the solution of the previous loadstep as starting point for the next one. I would suggest to set the boundary conditions only at the beginning:

``````u.Set(mesh.BoundaryCF({"back" : uLeft, "front" : uRight}),  definedon=mesh.Boundaries("back|front"))
for i in range(1, 11):
factor.Set(i/10)
...``````

The boundary condition at the right boundary is quite large. As on the opposite boundary you have zero Dirichlet condition you can set them by

``u.Set(x*uRight)``

which will provide a better initial guess helping Newton to converge.

Best,
Michael

Thanks Michael,

Indeed it’s mu/2., I somehow forgot to update it here.

``````import netgen.gui
from math import pi as myPi
from ngsolve import *
from netgen.csg import unit_cube

mesh = Mesh(unit_cube.GenerateMesh(maxh=0.06))
V = VectorH1(mesh, order=1, dirichlet="back|front")
print(mesh.GetBoundaries())
u, v = V.TnT()

# material parameters and forces
E, nu = 10.0, 0.3
mu = E / 2./(1+nu)
lmbda = 2 * mu * nu/(1-2*nu)
bodyForce = CoefficientFunction((0., -0.5, 0.))

# Define strain measures
I = Id(mesh.dim)
# help(F)
I1 = Trace(F.trans * F)
J = Det(F)
psi = mu/2. * (I1 - 3.) - mu * log(J) + lmbda/2. * (J - 1.0)**2
# S = mu * F - mu * (F.trans).Inv + lmbda * J * (J - 1.0) * (F.trans).Inv

# definition of bilinear and linear forms
a = BilinearForm(V, symmetric=False)
a += Variation(psi.Compile() * dx)
a += Variation((-1 * InnerProduct(bodyForce, u)).Compile() * dx)

u = GridFunction(V)
u.vec[:] = 0.

# definition of the boundaries
scale = 0.5
yo, zo = 0.5, 0.5
thta = myPi/3.
factor=Parameter(1.0)
uLeft = CoefficientFunction((0.0, 0., 0.))
uRight = CoefficientFunction((
0.,
scale * (yo + (y - yo) * cos(thta) - (z - zo) * sin(thta) - y),
scale * (zo + (y - yo) * sin(thta) + (z - zo) * cos(thta) - z)
))
u.Set(mesh.BoundaryCF({"back" : uLeft, "front" : x*uRight}),  definedon=mesh.Boundaries("back|front"))

res = u.vec.CreateVector()
w = u.vec.CreateVector()

maxiters = 20

a.Apply(u.vec, res)
a.AssembleLinearization(u.vec)
inv = a.mat.Inverse(V.FreeDofs(), inverse="sparsecholesky")
w.data = inv * res
resNorm = sqrt(abs(InnerProduct(w, w)))
u.vec.data -= w
iter = 0
while resNorm > 1.e-8 and iter < maxiters:
a.Apply(u.vec, res)
a.AssembleLinearization(u.vec)
inv = a.mat.Inverse(V.FreeDofs(), inverse="sparsecholesky")
w.data = inv * res
resNorm = sqrt(abs(InnerProduct(w, w)))
print(f"Norm: {resNorm}")
u.vec.data -= w
iter += 1
Draw(u, mesh, "displacement")
SetVisualization(deformation=True)``````

I will take a careful look today at the problem definition and bc’s more carefully today.

The source of the above exercise is the FEniCS demo on hyperelasticity:
https://fenicsproject.org/docs/dolfin/latest/python/demos/hyperelasticity/demo_hyperelasticity.py.html

I have implemented the above problem by myself using pure-NumPy too, so I will check further why this has problems converging…

So Newton indeed converges for a smaller mesh/problem, but fails on a reasonable mesh (65k elements). Unless I am doing something fundamentally wrong in NGSolve, I would expect that this should converge…

And to add, the exact same mesh when run with dolfin, converges in 5 iterations. I would try and assemble the linearized system instead of using the NewtonSolver class in Dolfin, but I expect that it would converge then too:

Here is the code:

``````import netgen.gui
from math import pi as myPi
from ngsolve import *
from netgen.csg import unit_cube

mesh = Mesh(unit_cube.GenerateMesh(maxh=0.05))
# for i in range(2):
#     mesh.Refine()

V = VectorH1(mesh, order=2, dirichlet="back|front")
print(mesh.GetBoundaries())
u, v = V.TnT()

# material parameters and forces
E, nu = 10.0, 0.3
mu = E / 2./(1+nu)
lmbda = 2 * mu * nu/(1-2*nu)
bodyForce = CoefficientFunction((0., -0.5, 0.))

# Define strain measures
I = Id(mesh.dim)
# help(F)
I1 = Trace(F.trans * F)
J = Det(F)
psi = mu/2. * (I1 - 3.) - mu * log(J) + lmbda/2. * (J - 1.0)**2
# S = mu * F - mu * (F.trans).Inv + lmbda * J * (J - 1.0) * (F.trans).Inv
factor=Parameter(1000.0)
# definition of bilinear and linear forms
a = BilinearForm(V, symmetric=False)
a += Variation(psi.Compile() * dx)
a += Variation((-1 * InnerProduct(bodyForce, u)).Compile() * dx)

u = GridFunction(V)
u.vec[:] = 0.

# definition of the boundaries
scale = 0.5
yo, zo = 0.5, 0.5
thta = myPi/3.

uLeft = CoefficientFunction((0.0, 0., 0.))
uRight = CoefficientFunction((
0.,
scale * (yo + (y - yo) * cos(thta) - (z - zo) * sin(thta) - y),
scale * (zo + (y - yo) * sin(thta) + (z - zo) * cos(thta) - z)
))
u.Set(mesh.BoundaryCF({"back" : uLeft, "front" : x*uRight}),  definedon=mesh.Boundaries("back|front"))

res = u.vec.CreateVector()
# resNew = u.vec.CreateVector()
w = u.vec.CreateVector()

maxiters = 50

a.Apply(u.vec, res)
a.AssembleLinearization(u.vec)
inv = a.mat.Inverse(V.FreeDofs(), inverse="pardiso")
w.data = inv * res
resNorm = sqrt(abs(InnerProduct(w, w)))
u.vec.data -= w
iter = 0
while resNorm > 1.e-8 and iter < maxiters:
factor.Set(1.e6)
a.Apply(u.vec, res)
a.AssembleLinearization(u.vec)

inv = a.mat.Inverse(V.FreeDofs(), inverse="pardiso")
w.data = inv * res
resNorm = sqrt(abs(InnerProduct(w, w)))
print(f"Norm: {resNorm}")
u.vec.data -= w
iter += 1
Draw(u, mesh, "displacement", deformation=True)
SetVisualization(deformation=True)``````

Dolfin code:

``````import numpy as np, meshio
from netgen.csg import unit_cube
from dolfin import *

def generateMesh():
msh = unit_cube.GenerateMesh(maxh=0.05)
msh.Export("netgenMesh.msh", "Gmsh2 Format")
points, cells = mesh.points, np.vstack([cell.data for cell in mesh.cells if cell.type == "tetra"])
meshio.xdmf.write("netGenMesh.xdmf", meshio.Mesh(points, cells = {"tetra": cells}))
return 1

generateMesh()
mesh = Mesh()
# mesh = UnitCubeMesh(24, 16, 16)

class Problem(NonlinearProblem):
def __init__(self, J, F, bcs):
self.bilinear_form = J
self.linear_form = F
self.bcs = bcs
NonlinearProblem.__init__(self)

def F(self, b, x):
assemble(self.linear_form, tensor=b)
for bc in self.bcs:
bc.apply(b, x)

def J(self, A, x):
assemble(self.bilinear_form, tensor=A)
for bc in self.bcs:
bc.apply(A)

E, nu = Constant(10), Constant(0.3)
mu = E/2./(1+nu)
lmbda = 2*mu*nu/(1-2*nu)

V = VectorFunctionSpace(mesh, "CG", 1)
u = Function(V)
u_ = TrialFunction(V)
v = TestFunction(V)
bodyForce = as_vector([0, -0.5, 0])
surFaceTraction = as_vector([0.1, 0., 0.])
left = "near(x[0], 0) && on_boundary"
right = "near(x[0], 1) && on_boundary"
rightBoundary = Expression((
"0.0",
"(0.5 + (x[1] - 0.5) * cos(pi/3) - (x[2] - 0.5) * sin(pi/3) - x[1])/2.",
"(0.5 + (x[1] - 0.5) * sin(pi/3) + (x[2] - 0.5) * cos(pi/3) - x[2])/2."
), pi=np.pi, degree=2)
bcs = [
DirichletBC(V, Constant((0,0,0)), left),
DirichletBC(V, rightBoundary, right)
]
I1 = inner(F, F)
J = det(F)
psi = mu/2. * (I1-3) - mu * ln(J) + lmbda/2 * (J-1)**2
pi = derivative(psi, u, v) * dx - inner(bodyForce, v) * dx - inner(surFaceTraction, v) * ds
Jac = derivative(pi, u, u_)

problem = Problem(Jac, pi, bcs)
solver = NewtonSolver()
solver.solve(problem, u.vector())
from vedo.dolfin import plot as vplot
vplot(u, mode="displaced mesh")``````

Newton iterations: NGSolve

``````Norm: 6298.188667502683
Norm: 16417474.934260502
Norm: 5484205427.2328205
Norm: 216770135728.68326
Traceback (most recent call last):
File ".\hyperelasticityDolfin.py", line 74, in <module>
inv = a.mat.Inverse(V.FreeDofs(), inverse="pardiso")
KeyboardInterrupt``````

Newton Iterations: Dolfin

``````Newton iteration 0: r (abs) = 4.895e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 1.295e-01 (tol = 1.000e-10) r (rel) = 2.646e-02 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 1.470e-02 (tol = 1.000e-10) r (rel) = 3.004e-03 (tol = 1.000e-09)
Newton iteration 3: r (abs) = 1.394e-03 (tol = 1.000e-10) r (rel) = 2.848e-04 (tol = 1.000e-09)
Newton iteration 4: r (abs) = 7.730e-06 (tol = 1.000e-10) r (rel) = 1.579e-06 (tol = 1.000e-09)
Newton iteration 5: r (abs) = 1.394e-09 (tol = 1.000e-10) r (rel) = 2.847e-10 (tol = 1.000e-09)
Newton solver finished in 5 iterations and 5 linear solver iterations.``````

Hi bhaveshshrimali,

when looking at the very first update the increment seems to overshoot such that Newton is not able any more to find the correct solution. I guess this is induced by the boundary conditions. Maybe Dolfin handles boundary conditions differently than NGSolve and thus converges easier.

The most simple adaption is to use a damped Newton to avoid the overshooting. You can use the build in Newton solver to play around with damping factors:

``solvers.Newton(a, u, dampfactor=0.25, maxit=maxiters, maxerr=1.e-8, printing=True)``

Another possibility is to solve a linear elasticity problem with the given boundary data and hope that the result gives an initial guess without inducing overshooting behaviour.

Attached you’ll find your code with the Newton damping (and some little changes).

Best
Michael

https://ngsolve.org/media/kunena/attachments/889/hyperelastic.py

Attachment: hyperelastic.py

Thanks a lot Michael! It seems a rather peculiar example to try. I will try another simpler hyperelasticity problem (simpler in terms of the boundary conditions)

I had a question: does NGSolve use all threads by default? Can I control the number of threads in command line usage? I am working on a Windows laptop.

Bhavesh

To assemble bilinearforms or invert matrices in parallel you need to activate the TaskManager. Then NGSolve uses all Threads per default. To limit the maximal threads you can use the SetNumThreads function or the NGS_NUM_THREADS environment variable, see also this tutorial.

If you use umfpack/pardiso you might need additional commands. For Linux I set the environment variables MKL_NUM_THREADS and OMP_NUM_THREADS.

Best,
Michael