I am trying is to occlude a cylinder with internal suction pressure (the aim is to occlude it as much as possible). I am using purely incompressible hyperelastic material with the 2-field formulation in Fenics. However, I can not occlude the cylinder after a certain limit. The figures below are the initial configuration and the deformed configuration of the cylinder.

I can not occlude the cylinder any further. However, if I use other FE packages, e.g., FEBio, I can almost fully occlude the cylinder. I am wondering what could be the possible reason that I can not occlude it further.

There are 2-3 things that might be wrong. First, although I am applying pressure normal to the inner face using the FacetNormal function, after a certain stage of deformation the pressure is not acting in the normal direction. If that is the case somehow we need to calculate the FacetNormal function in every step of deformation (do not know how to do that). Second, I am using the 2-field formulation whereas FeBio uses the 3-field formulation. I really do not know whether that would be helpful or not. But, I think the 2-field formulation should be fine.

Could you please help me to find the issue?

```
from dolfin import *
import matplotlib.pyplot as plt
import os
import dolfin as dlf
import numpy as np
import math
from mshr import *
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
# Kinematics
def pk1Stress(u,c1,pressure):
I = Identity(V.mesh().geometry().dim()) # Identity tensor
F = I + grad(u) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
Ic = tr(C) # Invariants of deformation tensors
J = det(F)
E = 0.5*(C-I)
pk1 = 2*c1*F-pressure*inv(F).T
return pk1, (J-1)
## geometry
def geometry_3d():
c_length = 2.0; # L
c_inner_radius = 1.0; # A
c_outer_radius = 2.0; # B
mesh = Mesh('./QuarterCylinder.xml')
boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
x0 = AutoSubDomain(lambda x: np.absolute(x[0]) < 0.01)
y0 = AutoSubDomain(lambda x: np.absolute(x[1]) < 0.01)
z0 = AutoSubDomain(lambda x: near(x[2], 0))
innerface = AutoSubDomain(lambda x: (np.sqrt(x[0]*x[0]+x[1]*x[1])- c_inner_radius) <= 0.1 ) #..
x0 .mark(boundary_parts, 1)
y0 .mark(boundary_parts, 2)
z0 .mark(boundary_parts, 3)
innerface.mark(boundary_parts, 4)
return boundary_parts
# Create mesh and define function space
facet_function = geometry_3d()
mesh = facet_function.mesh()
gdim = mesh.geometry().dim()
dx = Measure("dx")
ds = Measure("ds", subdomain_data=facet_function, subdomain_id=4)
print('Number of nodes: ',mesh.num_vertices())
print('Number of cells: ',mesh.num_cells())
# Limit quadrature degree
dx = dx(degree=4)
ds = ds(degree=4)
# Create function space (in mixed space approach, have to define elements)
element_u = VectorElement("CG", mesh.ufl_cell(), 2)
element_p = FiniteElement("CG", mesh.ufl_cell(), 1)
mixed_element = MixedElement([element_u,element_p])
V = FunctionSpace(mesh, mixed_element)
# Define test and trial functions
dup = TrialFunction(V)
_u_p = Function(V)
u, p = split(_u_p)
_u, _p = TestFunctions(V)
# Define Dirichlet boundary
bc0 = DirichletBC(V.sub(0).sub(0), Constant(0.), facet_function, 1)
bc1 = DirichletBC(V.sub(0).sub(1), Constant(0.), facet_function, 2)
bc2 = DirichletBC(V.sub(0).sub(2), Constant(0.), facet_function, 3)
bcs = [bc0,bc1,bc2]
# material parameters
c1 = 1000.0
# pressure on the inner face of the cylinder
facenormal = FacetNormal(mesh)
innerpress = Constant(-0.0)
I = Identity(V.mesh().geometry().dim())
F = I + grad(u)
press_trac = det(F)*dot(inv(F).T, innerpress*facenormal) # pressure
# # Total potential energy
I = Identity(V.mesh().geometry().dim())
pk1, p_bar = pk1Stress(u,c1,p)
dgF = I + grad(u)
F1 = inner(pk1, grad(_u))*dx - dot(press_trac, _u)*ds
F2 = p_bar*_p*dx
F = F1+F2
DJ = derivative(F, _u_p,dup)
# Create nonlinear variational problem and solve
problem = NonlinearVariationalProblem(F, _u_p, bcs=bcs, J=DJ)
solver = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['relative_tolerance'] = 1e-6
solver.parameters['newton_solver']['linear_solver'] = 'lu' # 'mumps', 'lu'
# Extract solution components
u, p = _u_p.split() # seperate functions
u.rename("u", "displacement")
p.rename("p", "pressure")
# Save solution in VTK format
file_results = XDMFFile("./Results/TestTriaxialLoading/Triaxial_3field.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
# Time stepping parameters
Tend = 1
Nsteps = 10
time = np.linspace(0, Tend, Nsteps+1)
dt = time[1]-time[0]
for i in range(len(time)):
tt = time[i]
print('time: ', time[i])
# Increase traction
innerpress.assign(2900.0*time[i]) #Pa
# solve and save disp
solver.solve()
# save xdmf file
file_results.write(u, time[i])
print('done')
```