Dear community/developers,
I am using FEniCSx 0.9 with Docker on Microsoft Windows.
In the MWE posted below, I perform a 2-D displacement-controlled FEA and compute the displacement U_disp_ctrl and stress sigma_disp_ctrl. What I am interested in is: imposing tractions obtained by t_N = ufl.dot(sigma_disp_ctrl, N) with N = ufl.FacetNormal(mesh) in a subsequent load-controlled FEA, and having a displacement field that is consistent with the one from the displacement-controlled FEA:
import dolfinx
from dolfinx import fem
import ufl
import dolfinx.fem.petsc
import dolfinx.nls.petsc
from mpi4py import MPI
import numpy as np
l_x = 10.0
l_y = 20.0
n_x = 10
n_y = 20
gamma_max = 2.0 # Prescribed extension in mm
mesh = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, np.array([[0, 0], [l_x, l_y]]), [n_x, n_y], cell_type=dolfinx.mesh.CellType.quadrilateral)
rank = mesh.comm.Get_rank() # Get process rank (ID)
num_process = mesh.comm.Get_size()
tdim = mesh.topology.dim # Topological dimension
fdim = tdim - 1 # Facet dimension
gdim = mesh.geometry.dim
# BOUNDARY CONDITIONS
def boundary_bottom(x):
return np.isclose(x[1], 0.0)
def boundary_top(x):
return np.isclose(x[1], l_y)
def boundary_left(x):
return np.isclose(x[0], 0.0)
def boundary_right(x):
return np.isclose(x[0], l_x)
facets = []
for boundary_locator in [boundary_bottom, boundary_top, boundary_left, boundary_right]:
facet = dolfinx.mesh.locate_entities_boundary(mesh, fdim, boundary_locator)
facets.append(facet)
facets_indexes_dict = {
"bottom": 1,
"top": 2,
"left": 3,
"right": 4
}
num_facets = mesh.topology.index_map(fdim).size_local
markers = np.zeros(num_facets, dtype=np.int32)
for facet, facet_index in zip(facets, facets_indexes_dict.values()):
markers[facet] = facet_index
facet_tags = dolfinx.mesh.meshtags(mesh, fdim, np.arange(num_facets, dtype=np.int32), markers)
V = fem.functionspace(mesh, ("Lagrange", 1, (gdim,)))
zero = fem.Function(V)
zero.x.array[:] = 0.0
gamma = fem.Constant(mesh, gamma_max)
u_ext = fem.Function(V)
u_ext_expression = fem.Expression(ufl.as_vector([0.0, gamma]), V.element.interpolation_points()) # u_ext on the top boundary
dofs_bottom_x = fem.locate_dofs_topological(V.sub(0), fdim, facet_tags.find(facets_indexes_dict["bottom"]))
dofs_bottom_y = fem.locate_dofs_topological(V.sub(1), fdim, facet_tags.find(facets_indexes_dict["bottom"]))
dofs_top_x = fem.locate_dofs_topological(V.sub(0), fdim, facet_tags.find(facets_indexes_dict["top"]))
dofs_top_y = fem.locate_dofs_topological(V.sub(1), fdim, facet_tags.find(facets_indexes_dict["top"]))
bc_bottom_x = fem.dirichletbc(zero.sub(0), dofs_bottom_x)
bc_bottom_y = fem.dirichletbc(zero.sub(1), dofs_bottom_y)
bc_top_x = fem.dirichletbc(zero.sub(0), dofs_top_x)
bc_top_y = fem.dirichletbc(u_ext.sub(1), dofs_top_y)
bc_all_disp_ctrl_list = [bc_bottom_y, bc_top_y, bc_bottom_x, bc_top_x]
bc_all_load_ctrl_list = [bc_bottom_y, bc_bottom_x, bc_top_x]
dofsets_known_reaction_list = [dofs_top_y]
# MATERIAL MODEL
U_disp_ctrl = fem.Function(V, name="U_disp_ctrl")
U_load_ctrl = fem.Function(V, name="U_load_ctrl")
E = dolfinx.default_scalar_type(1.0e4)
nu = dolfinx.default_scalar_type(0.3)
mu = fem.Constant(mesh, E / (2 * (1 + nu)))
lmbda = fem.Constant(mesh, E * nu / ((1 + nu) * (1 - 2 * nu)))
I = ufl.Identity(gdim)
sigma_disp_ctrl = 2.0 * mu * ufl.sym(ufl.grad(U_disp_ctrl)) + lmbda * ufl.tr(ufl.sym(ufl.grad(U_disp_ctrl))) * I
sigma_load_ctrl = 2.0 * mu * ufl.sym(ufl.grad(U_load_ctrl)) + lmbda * ufl.tr(ufl.sym(ufl.grad(U_load_ctrl))) * I
# RESIDUALS OF THE FINITE ELEMENT PROBLEMS
dU = ufl.TrialFunction(V)
deltaU = ufl.TestFunction(V)
dS = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tags)
dX = ufl.Measure("dx", domain=mesh)
F_int_disp_ctrl = ufl.inner(sigma_disp_ctrl, ufl.grad(deltaU)) * dX
F_int_load_ctrl = ufl.inner(sigma_load_ctrl, ufl.grad(deltaU)) * dX
resFE_disp_ctrl = F_int_disp_ctrl
N = ufl.FacetNormal(mesh)
t_N = ufl.dot(sigma_disp_ctrl, N)
F_ext_load_ctrl = ufl.inner(t_N, deltaU) * dS(facets_indexes_dict["top"])
# F_ext_load_ctrl = resFE_disp_ctrl
resFE_load_ctrl = F_int_load_ctrl - F_ext_load_ctrl
# JACOBIANS OF THE FINITE ELEMENT RESIDUALS
DresFE_disp_ctrl_DU_disp_ctrl = ufl.derivative(resFE_disp_ctrl, U_disp_ctrl, dU)
DresFE_load_ctrl_DU = ufl.derivative(resFE_load_ctrl, U_load_ctrl, dU)
# DISPLACEMENT-CONTROLLED PROBLEM
problem_disp_ctrl = dolfinx.fem.petsc.NonlinearProblem(F=resFE_disp_ctrl, u=U_disp_ctrl, bcs=bc_all_disp_ctrl_list, J=DresFE_disp_ctrl_DU_disp_ctrl)
solver_disp_ctrl = dolfinx.nls.petsc.NewtonSolver(comm=mesh.comm, problem=problem_disp_ctrl)
# LOAD-CONTROLLED PROBLEM
problem_load_ctrl = dolfinx.fem.petsc.NonlinearProblem(F=resFE_load_ctrl, u=U_load_ctrl, bcs=bc_all_load_ctrl_list, J=DresFE_load_ctrl_DU)
solver_load_ctrl = dolfinx.nls.petsc.NewtonSolver(comm=mesh.comm, problem=problem_load_ctrl)
# SOLVE DISPLACEMENT-CONTROLLED PROBLEM
u_ext.interpolate(u_ext_expression)
solver_disp_ctrl.solve(U_disp_ctrl)
# SOLVE LOAD-CONTROLLED PROBLEM
solver_load_ctrl.solve(U_load_ctrl)
# L-2 NORM OF DIFFERENCE BETWEEN DISPLACEMENT FIELDS
U_diff = U_load_ctrl - U_disp_ctrl
U_diff_L2_norm_form = fem.form(ufl.inner(U_diff, U_diff) * dX)
l2_norm_U_diff = np.sqrt(mesh.comm.allreduce(fem.assemble_scalar(U_diff_L2_norm_form), op=MPI.SUM))
if rank == 0:
print(f"\nL-2 norm of difference between displacement fields: {l2_norm_U_diff:.6e}")
Unfortunately, using the traction vector t_N I am not able to reach consistency:
L-2 norm of difference between displacement fields: 9.492333e-01
However, if I uncomment the line
# F_ext_load_ctrl = resFE_disp_ctrl
and directly impose the residual of the displacement-controlled problem as the external force in the load-controlled FEA, this inconsistency is resolved:
L-2 norm of difference between displacement fields: 8.206510e-13
Any ideas and solutions to solve this problem are very much appreciated.
P.S.: I am aware of using a non-linear solver for a linear elasticity problem; it is just because of the original problem that I need to use the non-linear solver.