Thank you very much for the reply @joachim. I was able to get the code to run. But I am confused about the output. In my understanding, I applied pressure on the inner surface of the sphere but the output shows that the highest magnitude of displacement is on the outer surface of the sphere?! Please help me understand whether I am applying pressure on the inner surface correctly.
a += Variation((-InnerProduct(pressure, u)).Compile()*ds('inner'))
The full working code is below:
from netgen.csg import CSGeometry, Plane, Pnt, Sphere, Vec
from ngsolve import (BilinearForm, Det, Grad, GridFunction, Id,
InnerProduct, Mesh, Norm, Parameter, SetVisualization, Trace, Variation, VectorH1,
ds, dx, specialcf)
from ngsolve.webgui import Draw
# Prepare mesh
yz_plane = Plane (Pnt(0, 0, 0), Vec(-1, 0, 0)).bc("xface")
zx_plane = Plane (Pnt(0, 0, 0), Vec(0, -1, 0)).bc("yface")
xy_plane = Plane (Pnt(0, 0, 0), Vec(0, 0, -1)).bc("zface")
inner_sphere = Sphere(Pnt(0, 0, 0), 5).maxh(0.5).bc("inner")
outer_sphere = Sphere(Pnt(0, 0, 0), 25).maxh(5).bc("outer")
one_eight_sphere = ((outer_sphere - inner_sphere) * yz_plane * zx_plane * xy_plane).mat('rubber')
geo = CSGeometry()
geo.Add(one_eight_sphere)
mesh = geo.GenerateMesh()
mesh.Refine()
mesh.SecondOrder()
#mesh.Export('TestSphere.msh', 'Gmsh2 Format')
mesh = Mesh(mesh)
# Prepare constitutive equation
E, nu = 1.794, 0.3
mu = 0.5*E/(1 + nu)
lam = E * nu / ((1 + nu)*(1 - 2*nu))
def C(u):
F = Id(3) + Grad(u)
return F.trans * F
def NeoHooke(C):
return 0.5 * mu * (Trace(C - Id(3)) + 2 * mu/lam * Det(C)**(-0.5*lam/nu) - 1)
# Prepare incremental load
factor = Parameter(0)
factor.Set(1.5)
pressure = factor * specialcf.normal(3)
# Define the FESpace
fes = VectorH1(mesh, order=2, dirichletx='xface', dirichlety='yface', dirichletz='zface')
u = fes.TrialFunction()
# Set up the governing weak-form
a = BilinearForm(fes, symmetric=False)
a += Variation(NeoHooke(C(u)).Compile()*dx)
a += Variation((-InnerProduct(pressure, u)).Compile()*ds('inner'))
# Set up the solver
u = GridFunction(fes)
u.vec[:] = 0
def SolveNewton():
res = u.vec.CreateVector()
for it in range(5):
print ("it", it, "energy = ", a.Energy(u.vec))
a.Apply(u.vec, res)
a.AssembleLinearization(u.vec)
inv = a.mat.Inverse(fes.FreeDofs() )
u.vec.data -= inv*res
# Solve!
SolveNewton()
Draw(Norm(u), mesh, "Displacement")
SetVisualization(deformation=True)