3D nonlinear elasticity set one component of displacement as zero on a face

Hi
I am trying to learn NGSPy. I am looking at one-eight of a thick spherical with internal pressure. I want to restrict the displacements on the y-z plane to have u_x = 0, in the z-x plane to have u_y = 0, and in the x-y plane to have u_z = 0.

Can I set the face to have Dirichlet Boundary condition? Does that mean all the three components of the displacement have to be specified on the boundaries marked as `dirichlet`? I was not able to find an example in the documentation that sets just one component of  displacement to 0. Can you kindly give a sample how to do it?

I have attached the code for the geometry, the load and the FESpace:

from ngsolve import Draw, Redraw, Id, Grad, Trace, Det, Parameter, specialcf, H1
from netgen.csg import Plane, Sphere, Pnt, Vec, CSGeometry

# 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')

# 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(0.05)
pressure = factor * specialcf.normal(3)

# Define the FESpace
fes = H1(mesh, order=2, dirichlet='xface|yface|zface', dim=3)

Hi,

if you replacer H1(... dim=3) by VectorH1 (and you don't give the dim-flag), you can use "dirichletx", "dirichlety" and "dirichletz" flags to set Dirichlet-constraints on the individual Cartesian coordinates.

Joachim

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)

The setup is fine, but I think its a problem with this version of Neo-Hooke material law. 
It looks good if I use the following one (but then we need loadsteps for convergence).

def NeoHooke(C):
    # return 0.5 * mu * (Trace(C - Id(3)) + 2 * mu/lam * Det(C)**(-0.5*lam/nu) - 1)
    return mu/2 * (Trace(C)-3) - mu/2*log(Det(C)) + lam/2*(sqrt(Det(C))-1)**2
 

Thank you very much for your reply.