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.
gmsh.model.add("temp")
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)
l5 = gmsh.model.geo.addLine(p6, p7, 5)
l6 = gmsh.model.geo.addLine(p7, p8, 6)
l7 = gmsh.model.geo.addLine(p8, p9, 7)
l8 = gmsh.model.geo.addLine(p9, p6, 8)
cloop1 = gmsh.model.geo.addCurveLoop([l1,l2,l3,l4], 9)
cloop2 = gmsh.model.geo.addCurveLoop([l5,l6,l7,l8], 10)
s1 = gmsh.model.geo.addPlaneSurface([cloop1], 11)
s2 = gmsh.model.geo.addPlaneSurface([cloop2,cloop1], 12)
gmsh.model.geo.synchronize()
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.model.mesh.setRecombine(2, s1)
gmsh.model.mesh.setRecombine(2, s2)
bottom = gmsh.model.addPhysicalGroup(1, [l5], 1)
right = gmsh.model.addPhysicalGroup(1, [l6], 2)
top = gmsh.model.addPhysicalGroup(1, [l7], 3)
left = gmsh.model.addPhysicalGroup(1, [l8], 4)
interface = gmsh.model.addPhysicalGroup(1, [l1,l2,l3,l4], 5)
inclusion = gmsh.model.addPhysicalGroup(2, [s1], 1)
matrix = gmsh.model.addPhysicalGroup(2, [s2], 2)
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]
metadata = {"quadrature_degree": 4}
dx = ufl.Measure("dx", domain=msh, subdomain_data=cell_markers, metadata=metadata)
ds = ufl.Measure('ds', domain=msh, subdomain_data=facet_markers, metadata=metadata)
dS = ufl.Measure("dS", domain=msh, subdomain_data=facet_markers, metadata=metadata)
# Define elements
celltype = ufl.Cell("quadrilateral", 2)
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)
# Deformation gradient
F = ufl.variable(I + ufl.grad(u))
# 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):
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
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