Solving Hyperleasticity: Correctly defining boundaries

bhaveshshrimali
Topic Author
Offline
Junior Member

4 years 1 week ago #3333 by bhaveshshrimali

Replied by bhaveshshrimali on topic Solving Hyperleasticity: Correctly defining boundaries

Thanks Michael,

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

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.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)
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

# 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:
fenicsproject.org/docs/dolfin/latest/pyt...erelasticity.py.html

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

bhaveshshrimali
Topic Author
Offline
Junior Member

4 years 1 week ago #3337 by bhaveshshrimali

Replied by bhaveshshrimali on topic Solving Hyperleasticity: Correctly defining boundaries

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:

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:

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

Code:

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

Code:

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.

mneunteufel
Offline
Premium Member

4 years 6 days ago #3338 by mneunteufel

Replied by mneunteufel on topic Solving Hyperleasticity: Correctly defining boundaries

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:

Code:

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

File Attachment:

File Name: hyperelastic.py
File Size:2 KB

Attachments:

hyperelastic.py

bhaveshshrimali
Topic Author
Offline
Junior Member

4 years 6 days ago #3341 by bhaveshshrimali

Replied by bhaveshshrimali on topic Solving Hyperleasticity: Correctly defining boundaries

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

mneunteufel
Offline
Premium Member

4 years 6 days ago #3343 by mneunteufel

Replied by mneunteufel on topic Solving Hyperleasticity: Correctly defining boundaries

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

The following user(s) said Thank You: bhaveshshrimali

Time to create page: 0.136 seconds