Subdomains with different materials, solver not converging

dolfinx
1

Hello,

I am trying to solve a problem with different materials. I have simple 2D beam with tensile loading on the right end and a square in the middle which consists of a stiffer material. My code is largely based on Defining subdomains for different materials — FEniCSx tutorial and Hyperelasticity — FEniCSx tutorial. However, in my case the solver does not converge (RuntimeError: Newton solver did not converge because maximum number of iterations reached).

The following code produces the error message:

import ufl
import numpy as np
from dolfinx import fem, mesh, plot, cpp, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI


# structural parameters
Lx = 2
Ly = 1
N = 100
# lower left and upper right corner of inner square with second material
ll = [0.75, 0.25]
ur = [1.25, 0.75]
# force on right end
f_right = (0.5, 0.0)

# material parameters
E1 = default_scalar_type(1.)
E2 = default_scalar_type(10.)
nu = default_scalar_type(0.3)

# mesh
domain = mesh.create_rectangle(MPI.COMM_WORLD, [[0.0, 0.0], [Lx, Ly]], [Lx*N, Ly*N], mesh.CellType.quadrilateral)
Q = fem.functionspace(domain, ("DG", 0))
V = fem.functionspace(domain, ("Lagrange", 1 , (domain.geometry.dim, )))

# left end of beam
def left(x):
    return np.isclose(x[0], 0)

# right end of beam
def right(x):
    return np.isclose(x[0], Lx)

# first material
def domain1(x):
    return np.logical_and(np.logical_or(x[0] <= ll[0], x[0] >= ur[0]), np.logical_or(x[1] <= ll[1], x[1] >= ur[1]))

# second material
def domain2(x):
    return np.logical_and(np.logical_and(ll[0] <= x[0], x[0] <= ur[0]), np.logical_and(ll[1] <= x[1], x[1] <= ur[1]))

# mark cells for bcs
fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, left)
right_facets = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, right)
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag_bc = mesh.meshtags(domain, domain.topology.dim - 1, marked_facets[sorted_facets], marked_values[sorted_facets])

# define lame constants on both materials
cells1 = mesh.locate_entities(domain, domain.topology.dim, domain1)
cells2 = mesh.locate_entities(domain, domain.topology.dim, domain2)
mu = fem.Function(Q)
mu.x.array[cells1] = np.full_like(cells1, E1 / (2 * (1 + nu)), dtype=default_scalar_type)
mu.x.array[cells2] = np.full_like(cells2, E2 / (2 * (1 + nu)), dtype=default_scalar_type)
lam = fem.Function(Q)
lam.x.array[cells1] = np.full_like(cells1, E1 * nu / ((1 + nu) * (1 - 2 * nu)), dtype=default_scalar_type)
lam.x.array[cells2] = np.full_like(cells2, E2 * nu / ((1 + nu) * (1 - 2 * nu)), dtype=default_scalar_type)

# dirchlet bc
u_zero = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)
left = fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
dirichlet_bc = [fem.dirichletbc(u_zero, left, V)]

# test and solution functions
v = ufl.TestFunction(V)
u = fem.Function(V)
# traction 
t = fem.Constant(domain, default_scalar_type(f_right))

# spatial dimension
d = len(u)
# identity tensor
I = ufl.variable(ufl.Identity(d))
# lagrange strain
eps = ufl.sym(ufl.grad(u))

# strain energy
sig = 2 * mu * eps + lam * ufl.tr(eps) * I
psi = 0.5 * ufl.inner(eps, sig)

# measure inner domain
dx = ufl.Measure('dx', domain=domain)
# measure right end
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag_bc)

# potential energy (integrate strain energy on inner and outer material, traction energy on right boundary)
Pi = psi*dx - ufl.dot(t, u)*ds
# variation of Pi 
F = ufl.derivative(Pi, u, v)

# define problem F(u) = 0
problem = NonlinearProblem(F, u, dirichlet_bc)
# set solver and options
solver = NewtonSolver(domain.comm, problem)
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"

# solve
solver.solve(u)

In an attempt to avoid the error, I have used the following approach:

import ufl
import numpy as np
from dolfinx import fem, mesh, plot, cpp, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI


# structural parameters
Lx = 2
Ly = 1
N = 100
# lower left and upper right corner of inner square with second material
ll = [0.75, 0.25]
ur = [1.25, 0.75]
# force on right end
f_right = (0.5, 0.0)

# material parameters
E1 = default_scalar_type(1.)
E2 = default_scalar_type(10.)
nu = default_scalar_type(0.3)

# mesh
domain = mesh.create_rectangle(MPI.COMM_WORLD, [[0.0, 0.0], [Lx, Ly]], [Lx*N, Ly*N], mesh.CellType.quadrilateral)
Q = fem.functionspace(domain, ("DG", 0))
V = fem.functionspace(domain, ("Lagrange", 1 , (domain.geometry.dim, )))

# left end of beam
def left(x):
    return np.isclose(x[0], 0)

# right end of beam
def right(x):
    return np.isclose(x[0], Lx)

# inner square for second material
def square(x):
    return np.logical_and(np.logical_and(0.75 <= x[0], x[0] <= 1.25), np.logical_and(0.25 <= x[1], x[1] <= 0.75)) 

# mark cells for bcs
fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, left)
right_facets = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, right)
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag_bc = mesh.meshtags(domain, domain.topology.dim - 1, marked_facets[sorted_facets], marked_values[sorted_facets])

# mark cells for inner material
marked_facets = mesh.locate_entities(domain, domain.topology.dim, square)
sorted_facets = np.argsort(marked_facets)
marked_values = np.full_like(marked_facets, 1)
facet_tag_square = mesh.meshtags(domain, domain.topology.dim, marked_facets[sorted_facets], marked_values[sorted_facets])

# dirchlet bc
u_zero = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)
left = fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
dirichlet_bc = [fem.dirichletbc(u_zero, left, V)]

# test and solution functions
v = ufl.TestFunction(V)
u = fem.Function(V)
# traction 
t = fem.Constant(domain, default_scalar_type(f_right))

# spatial dimension
d = len(u)
# identity tensor
I = ufl.variable(ufl.Identity(d))
# lagrange strain
eps = ufl.sym(ufl.grad(u))

# lame coefficients
mu = fem.Constant(domain, 1 / (2 * (1 + nu)))
lam = fem.Constant(domain, 1 * nu / ((1 + nu) * (1 - 2 * nu)))
# strain energy
sig = 2 * mu * eps + lam * ufl.tr(eps) * I
psi = 0.5 * ufl.inner(eps, sig)

# measure inner domain
dx = ufl.Measure('dx', domain=domain, subdomain_data=facet_tag_square)
# measure right end
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag_bc)

# potential energy (integrate strain energy on inner and outer material, traction energy on right boundary)
Pi = E1*psi*dx - E1*psi*dx(1) + E2*psi*dx(1) - ufl.dot(t, u)*ds(2)
# variation of Pi 
F = ufl.derivative(Pi, u, v)

# define problem F(u) = 0
problem = NonlinearProblem(F, u, dirichlet_bc)
# set solver and options
solver = NewtonSolver(domain.comm, problem)
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"

# solve
solver.solve(u)

The second code does converge and yields the correct solution. However, it is rather inelegant and I would prefer to get the first approach to work for various reasons. As far as the first approach goes, assigning Lamé constants to elements and even computing the potential energy seems to work just fine. I am at a loss as to where the problem lies and appreciate any help or guidance.

I am using FEniCSx v0.7.0 on Docker.

2

the diff between the two codes is

@@ -34,13 +34,9 @@
 def right(x):
     return np.isclose(x[0], Lx)
 
-# first material
-def domain1(x):
-    return np.logical_and(np.logical_or(x[0] <= ll[0], x[0] >= ur[0]), np.logical_or(x[1] <= ll[1], x[1] >= ur[1]))
-
-# second material
-def domain2(x):
-    return np.logical_and(np.logical_and(ll[0] <= x[0], x[0] <= ur[0]), np.logical_and(ll[1] <= x[1], x[1] <= ur[1]))
+# inner square for second material
+def square(x):
+    return np.logical_and(np.logical_and(0.75 <= x[0], x[0] <= 1.25), np.logical_and(0.25 <= x[1], x[1] <= 0.75)) 
 
 # mark cells for bcs
 fdim = domain.topology.dim - 1
@@ -51,15 +47,11 @@
 sorted_facets = np.argsort(marked_facets)
 facet_tag_bc = mesh.meshtags(domain, domain.topology.dim - 1, marked_facets[sorted_facets], marked_values[sorted_facets])
 
-# define lame constants on both materials
-cells1 = mesh.locate_entities(domain, domain.topology.dim, domain1)
-cells2 = mesh.locate_entities(domain, domain.topology.dim, domain2)
-mu = fem.Function(Q)
-mu.x.array[cells1] = np.full_like(cells1, E1 / (2 * (1 + nu)), dtype=default_scalar_type)
-mu.x.array[cells2] = np.full_like(cells2, E2 / (2 * (1 + nu)), dtype=default_scalar_type)
-lam = fem.Function(Q)
-lam.x.array[cells1] = np.full_like(cells1, E1 * nu / ((1 + nu) * (1 - 2 * nu)), dtype=default_scalar_type)
-lam.x.array[cells2] = np.full_like(cells2, E2 * nu / ((1 + nu) * (1 - 2 * nu)), dtype=default_scalar_type)
+# mark cells for inner material
+marked_facets = mesh.locate_entities(domain, domain.topology.dim, square)
+sorted_facets = np.argsort(marked_facets)
+marked_values = np.full_like(marked_facets, 1)
+facet_tag_square = mesh.meshtags(domain, domain.topology.dim, marked_facets[sorted_facets], marked_values[sorted_facets])
 
 # dirchlet bc
 u_zero = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)
@@ -79,17 +71,20 @@
 # lagrange strain
 eps = ufl.sym(ufl.grad(u))
 
+# lame coefficients
+mu = fem.Constant(domain, 1 / (2 * (1 + nu)))
+lam = fem.Constant(domain, 1 * nu / ((1 + nu) * (1 - 2 * nu)))
 # strain energy
 sig = 2 * mu * eps + lam * ufl.tr(eps) * I
 psi = 0.5 * ufl.inner(eps, sig)
 
 # measure inner domain
-dx = ufl.Measure('dx', domain=domain)
+dx = ufl.Measure('dx', domain=domain, subdomain_data=facet_tag_square)
 # measure right end
 ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag_bc)
 
 # potential energy (integrate strain energy on inner and outer material, traction energy on right boundary)
-Pi = psi*dx - ufl.dot(t, u)*ds
+Pi = E1*psi*dx - E1*psi*dx(1) + E2*psi*dx(1) - ufl.dot(t, u)*ds(2)
 # variation of Pi 
 F = ufl.derivative(Pi, u, v)

It would be appreciated if you explained in your own words what these differences are, as it may help people answer your question.

3

Hi,
the difference is that the first code defines the Lamé parameters mu and lam as functions which assume different values depending on which part of the domain they are evaluated on (inner or outer material). Using these, the stored energy density is integrated over the whole domain in order to compute the potential energy. As indicated, this is the approach presented in Defining subdomains for different materials — FEniCSx tutorial.
In the second code, the Lamé parameters are defined as constants. Generally, the parameters, and thus the stored energy density itself, scale linearly with the Young’s modulus E (E1 for the outer material, E2 for the inner). Hence, one can set E=1 when computing the stored energy density, integrate on the whole domain multiplied with E1, then subtract the density integrated over the inner material mutliplied with E1, and then add the density integrated over the inner material but this time multiplied with E2. Thus, the stored energy for the whole domain is obtained and the potential energy can be computed as before. However, this approach is more of a workaround and comes with some problems, most notably that the stress field for the whole domain cannot be obtained as easily.

4

I have found a mistake in the function assigning coordinates to the outer material in the first code. The correct function is:

def domain1(x):
    return np.logical_or(np.logical_or(x[0] <= ll[0], x[0] >= ur[0]), np.logical_or(x[1] <= ll[1], x[1] >= ur[1]))

The code works fine now.