# Nonlinear incompressible inclusion problem with SNES solver

Dear community,

I am now implementing a nonlinear elastic problem with two subdomains, matrix and inclusion. Both the materials are Neo-Hookean, while the matrix is compressible and the inclusion is incompressible. To account for the incompressibility of the inlcuison, a pressure field is defined within the inclusion. I would like to restrict the degree of freedom of the pressure field within the inclusion only, this is achieved by `multiphenicsx` package. My code is as follows.

# Problem definitions

``````import gmsh
import numpy as np
import ufl
from typing import List, Tuple, Optional

from mpi4py import MPI
from petsc4py import PETSc

import dolfinx
from dolfinx import fem, io, mesh, plot, nls, log, la
import meshio
import dolfinx.io.gmshio as gmshio

from matplotlib import pyplot as plt

import multiphenicsx.fem as mfem

def generate_2D_full_comp(L, W, a, b, sz_in, sz_out, filename):
'''
This function generates a rectangular composite with an ellipsoidal inclusion
at its centre

Parameters
----------
L : TYPE float
DESCRIPTION Length of the composite.
W : TYPE float
DESCRIPTION Width of the composite.
a : TYPE float
DESCRIPTION Major axis of the composite.
b : TYPE float
DESCRIPTION Minor axis of the composite .
filename : TYPE
DESCRIPTION.

Returns
-------
None.

'''
gmsh.initialize()
gmsh.option.setNumber("General.Terminal",0) # Uncomment this line if testing the meshing function for the first time.

p1 = gmsh.model.geo.addPoint(L/2, W/2, 0, sz_in, 1)  # Center
p2 = gmsh.model.geo.addPoint(L/2 + a, W/2, 0, sz_in, 2)  # Right ellipse
p3 = gmsh.model.geo.addPoint(L/2, b + W/2, 0, sz_in, 3)  # Top ellipse
p4 = gmsh.model.geo.addPoint(L/2 - a, W/2, 0, sz_in, 4) # Left ellipse
p5 = gmsh.model.geo.addPoint(L/2, -b + W/2, 0, sz_in, 5) # Bottom ellipse

p6 = gmsh.model.geo.addPoint(0.0, 0.0, 0, sz_out, 6) # Left bottom corner
p7 = gmsh.model.geo.addPoint(L, 0.0, 0, sz_out, 7) # Right bottom corner
p8 = gmsh.model.geo.addPoint(L, W, 0, sz_out, 8) # Right top corner
p9 = gmsh.model.geo.addPoint(0.0, W, 0, sz_out, 9) # Left top corner

l1 = gmsh.model.geo.addEllipseArc(p2, p1, p2, p3, 1)
l2 = gmsh.model.geo.addEllipseArc(p3, p1, p3, p4, 2)
l3 = gmsh.model.geo.addEllipseArc(p4, p1, p4, p5, 3)
l4 = gmsh.model.geo.addEllipseArc(p5, p1, p5, p2, 4)

gmsh.model.geo.synchronize()
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.model.mesh.setRecombine(2, s1)
gmsh.model.mesh.setRecombine(2, s2)

gmsh.model.mesh.generate(2)
gmsh.write(filename)

filename = "Composite_block"
L_0 = 5.0
R_0 = 0.1

sz_in = R_0 / 10
sz_out = L_0 / 5
# Generate mesh
msh_filename =  filename + ".msh"
generate_2D_full_comp(L_0, L_0, R_0, R_0, sz_in, sz_out, msh_filename)

msh, cell_markers, facet_markers = gmshio.read_from_msh(msh_filename, MPI.COMM_WORLD, gdim=2)
tdim = msh.topology.dim
fdim = tdim - 1
# Find cells(subdomains)
cells_inclusion = cell_markers.indices[cell_markers.values == 1]  # Subdomain: Inclusion
cells_matrix = cell_markers.indices[cell_markers.values == 2]  # Subdomain: Matrix
# Find facets(boundaries)
facets_bottom = facet_markers.indices[facet_markers.values == 1]
facets_right = facet_markers.indices[facet_markers.values == 2]
facets_top = facet_markers.indices[facet_markers.values == 3]
facets_left = facet_markers.indices[facet_markers.values == 4]
facets_interface = facet_markers.indices[facet_markers.values == 5]

# Define elements
u_el = ufl.VectorElement("CG", celltype, 3)
p_el = ufl.FiniteElement("CG", celltype, 2)

# Define function spaces
u_space = fem.FunctionSpace(msh, u_el)
p_space = fem.FunctionSpace(msh, p_el)

# Find dofs in each subdomain
dofs_inc_p = fem.locate_dofs_topological(p_space, tdim, cells_inclusion)
dofs_u = np.arange(u_space.dofmap.index_map.size_local)

# Define restrictions
restriction_u = mfem.DofMapRestriction(u_space.dofmap, dofs_u)
restriction_p = mfem.DofMapRestriction(p_space.dofmap, dofs_inc_p)
restrictions = [restriction_u, restriction_p]

# Define dirichlet boundary conditions
# Find subspaces
ux_space, _ = u_space.sub(0).collapse()
uy_space, _ = u_space.sub(1).collapse()
# Find dofs at boundaries
dofs_left_ux = fem.locate_dofs_topological((u_space.sub(0), ux_space), fdim, facets_left)
dofs_bottom_uy = fem.locate_dofs_topological((u_space.sub(1), uy_space), fdim, facets_bottom)
dofs_right_ux = fem.locate_dofs_topological((u_space.sub(0), ux_space), fdim, facets_right)
dofs_top_uy = fem.locate_dofs_topological((u_space.sub(1), uy_space), fdim, facets_top)

# Define the displacement boundary condition interpolation expression
class strain_expr:
def __init__(self):
self.m = 0.0
def eval_x(self, x):
return np.full(x.shape[1], self.m*x[0])
def eval_y(self, x):
return np.full(x.shape[1], self.m*x[1])
def eval_z(self, x):
return np.full(x.shape[1], self.m*x[2])
def fix_m(self, x):
return np.full(x.shape[1], self.m)

no_disp = 0.0
u_left_x = fem.Function(ux_space)
u_left_x_expr = strain_expr()
u_left_x_expr.m = no_disp
u_left_x.interpolate(u_left_x_expr.eval_x)

u_bottom_y = fem.Function(uy_space)
u_bottom_y_expr = strain_expr()
u_bottom_y_expr.m = no_disp
u_bottom_y.interpolate(u_bottom_y_expr.eval_y)

u_right_x = fem.Function(ux_space)
u_right_x_expr = strain_expr()
u_right_x_expr.m = 0.01
u_right_x.interpolate(u_right_x_expr.eval_x)

bc_left_ux = fem.dirichletbc(u_left_x, dofs_left_ux, u_space.sub(0))
bc_bottom_uy = fem.dirichletbc(u_bottom_y, dofs_bottom_uy, u_space.sub(1))
bc_right_ux = fem.dirichletbc(u_right_x, dofs_right_ux, u_space.sub(0))

bcs = [bc_left_ux, bc_bottom_uy, bc_right_ux]

# Define the weak form
# Test function and solution functions
(u, p) = (fem.Function(u_space), fem.Function(p_space))
(u_trial, p_trial) = (ufl.TrialFunction(u_space), ufl.TrialFunction(p_space))
(du, dp) = (ufl.TestFunction(u_space), ufl.TestFunction(p_space))

G_i = fem.Constant(msh, PETSc.ScalarType(200.0))
G_m = fem.Constant(msh, PETSc.ScalarType(200.0))
lambda_i = fem.Constant(msh, PETSc.ScalarType(1000.0))
lambda_m = fem.Constant(msh, PETSc.ScalarType(1000.0))

# Kinematics
# Identity matrix
I = ufl.Identity(tdim)
# Jacobian
J  = ufl.det(F)
# Inverse of the deformation gradient
f = ufl.inv(F)

# Free energies
# Strain energy
V_G_field = fem.FunctionSpace(msh, ("DG", 0))
G_field = fem.Function(V_G_field)
G_field.x.array[cells_matrix] = np.full_like(cells_matrix, G_m.value, dtype=PETSc.ScalarType)
G_field.x.array[cells_inclusion] = np.full_like(cells_inclusion, G_i.value, dtype=PETSc.ScalarType)
G_field.name = "Shear_modulus"

V_lambda_ = fem.FunctionSpace(msh, ("DG", 0))
lambda_ = fem.Function(V_lambda_)
lambda_.x.array[cells_matrix] = np.full_like(cells_matrix, lambda_m.value, dtype=PETSc.ScalarType)
lambda_.x.array[cells_inclusion] = np.full_like(cells_inclusion, lambda_i.value, dtype=PETSc.ScalarType)
lambda_.name = "Bulk_modulus"

U_strain = (G_field / 2.0) * (Ic - tdim) - G_field * ufl.ln(J) + lambda_ / 2 * (J - 1)**2
P_el = ufl.diff(U_strain, F)

P_isochoric = p * J * f.T
``````

# Definition of residual function and the solver

``````R_incompressible = [ufl.inner(P_el, ufl.grad(du)) * dx + ufl.inner(P_isochoric, ufl.grad(du)) * dx(1), dp * (J-1) * dx(1)]

dR_incompressible = [[ufl.derivative(R_incompressible[0], u, u_trial), ufl.derivative(R_incompressible[0], p, p_trial)],
[ufl.derivative(R_incompressible[1], u, u_trial), ufl.derivative(R_incompressible[1], p, p_trial)]]

restrictions_incompressible = [restriction_u, restriction_p]

class IncompressibleInclusionProblem:
'''
Nonlinear problem class compatible with PETSC.SNES solver.
'''

def __init__(self, R: List[ufl.Form], dR: List[List[ufl.Form]],
solutions: Tuple[fem.Function, fem.Function],
bcs: List[fem.DirichletBCMetaClass],
P: Optional[List[List[ufl.Form]]]=None):
'''

Parameters
----------
R : List[ufl.Form]
DESCRIPTION Residual.
dR : List[List[ufl.Form]]
DESCRIPTION Jacobian.
solutions : Tuple[fem.Function, fem.Function]
DESCRIPTION Solution functions.
bcs : List[dolfinx.DirichletBC]
DESCRIPTION Dirichlet boundary conditions.
Returns
-------
None.

'''

self._R = fem.form(R)
self._dR = fem.form(dR)
self._obj_vec = mfem.petsc.create_vector_block(self._R, restrictions_incompressible)
self._solutions = solutions
self._bcs = bcs
self._P = P

def create_snes_solution(self):
x = mfem.petsc.create_vector_block(self._R, restriction=restrictions_incompressible)
with mfem.petsc.BlockVecSubVectorWrapper(x, [u_space.dofmap, p_space.dofmap], restrictions_incompressible) as x_wrapper:
for x_wrapper_local, sub_solution in zip(x_wrapper, self._solutions):
with sub_solution.vector.localForm() as sub_solution_local:
x_wrapper_local[:] = sub_solution_local
return x

def update_solutions(self, x: PETSc.Vec):
with mfem.petsc.BlockVecSubVectorWrapper(x, [u_space.dofmap, p_space.dofmap], restrictions_incompressible) as x_wrapper:
for x_wrapper_local, sub_solution in zip(x_wrapper, self._solutions):
with sub_solution.vector.localForm() as sub_solution_local:
x_wrapper_local[:] = sub_solution_local

def obj(self, snes: PETSc.SNES, x: PETSc.Vec):
self.R(snes, x, self._obj_vec)
return self._obj_vec.norm()

def R(self, snes: PETSc.SNES, x: PETSc.Vec, R_vec: PETSc.Vec):
with R_vec.localForm() as R_vec_local:
R_vec_local.set(0.0)
mfem.petsc.assemble_vector_block(
R_vec, self._R, self._dR, self._bcs, x0=x, scale=-1.0,
restriction = restrictions_incompressible, restriction_x0=restrictions_incompressible)
def dR(self, snes: PETSc.SNES, x: PETSc.Vec, dR_mat: PETSc.Mat, P_mat: PETSc.Mat):
dR_mat.zeroEntries()
mfem.petsc.assemble_matrix_block(
dR_mat, self._dR, self._bcs, diagonal=1.0,
restriction = (restrictions_incompressible, restrictions_incompressible))
dR_mat.assemble()
if self._P is not None:
P_mat.zeroEntries()
mfem.petsc.assemble_matrix_block(
P_mat, self._P, self._bcs, diagonal=1.0,
restriction=(restrictions_incompressible, restrictions_incompressible))
P_mat.assemble()

problem_incompressible =IncompressibleInclusionProblem(R_incompressible, dR_incompressible, (u, p), bcs)
R_vec_incompressible = mfem.petsc.create_vector_block(problem_incompressible._R, restriction=restrictions_incompressible)
dR_mat_incompressible = mfem.petsc.create_matrix_block(problem_incompressible._dR, restriction=(restrictions_incompressible, restrictions_incompressible))

# Solve
snes = PETSc.SNES().create(msh.comm)
snes.setType('newtonls')
snes.setTolerances(max_it=10)
snes.getKSP().setType("preonly")
snes.getKSP().getPC().setType("lu")
snes.getKSP().getPC().setFactorSolverType("mumps")
snes.getKSP().setTolerances(atol=1e-9, rtol=1.0e-9)
PETSc.Options()['snes_monitor'] = None
PETSc.Options()['snes_linesearch_monitor'] = None
snes.setFromOptions()
snes.setObjective(problem_incompressible.obj)
snes.setFunction(problem_incompressible.R, R_vec_incompressible)
snes.setJacobian(problem_incompressible.dR, J=dR_mat_incompressible, P=None)

solution_incompressible = problem_incompressible.create_snes_solution()
snes.solve(None, solution_incompressible)
problem_incompressible.update_solutions(solution_incompressible)
solution_incompressible.destroy()
snes.destroy()
``````

But the problem does not converge and the residual energy during iterations does not change, the convergence details are:

``````  0 SNES Function norm 8.876329999380e+02
Line search: Using full step: obj 8.876329999380e+02 obj 3.400456520180e-04
1 SNES Function norm 3.400456520180e-04
Line search: Using full step: obj 3.400456520180e-04 obj 3.400456520180e-04
2 SNES Function norm 3.400456520180e-04
Line search: Using full step: obj 3.400456520180e-04 obj 3.400456520180e-04
3 SNES Function norm 3.400456520180e-04
Line search: Using full step: obj 3.400456520180e-04 obj 3.400456520180e-04
4 SNES Function norm 3.400456520180e-04
Line search: Using full step: obj 3.400456520180e-04 obj 3.400456520180e-04
5 SNES Function norm 3.400456520180e-04
Line search: Using full step: obj 3.400456520180e-04 obj 3.400456520180e-04
6 SNES Function norm 3.400456520180e-04
Line search: Using full step: obj 3.400456520180e-04 obj 3.400456520180e-04
7 SNES Function norm 3.400456520180e-04
Line search: Using full step: obj 3.400456520180e-04 obj 3.400456520180e-04
8 SNES Function norm 3.400456520180e-04
Line search: Using full step: obj 3.400456520180e-04 obj 3.400456520180e-04
9 SNES Function norm 3.400456520180e-04
Line search: Using full step: obj 3.400456520180e-04 obj 3.400456520180e-04
10 SNES Function norm 3.400456520180e-04
``````