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)
F = I + grad(u)
# 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)
# a += Variation((0.5/factor * InnerProduct(grad(u), grad(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")
mesh = meshio.read("netgenMesh.msh")
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()
XDMFFile("netGenMesh.xdmf").read(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)
]
F = Identity(len(u)) + grad(u)
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.