I’m been playing around with revising constraining across all axis (ref: Convergence difficulties with 3D periodic BCs - #2 by dokken)
which I’ve applied to your example:
import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import fem, io, mesh
import petsc4py.PETSc as PETSc
from ufl import grad, ln, tr, det, variable, derivative, TestFunction, TrialFunction
from dolfinx_mpc import NonlinearProblem, MultiPointConstraint
# ------------------------------------------------------------------------------
# Mesh Generation
# ------------------------------------------------------------------------------
comm = MPI.COMM_WORLD
# Generate simple rectangular mesh
length_plate = 1.0
n_elem = 30
domain = mesh.create_rectangle(
comm=comm,
points=((0.0, 0.0), (length_plate, length_plate)),
n=(n_elem, n_elem),
cell_type=mesh.CellType.triangle,
)
dim = domain.topology.dim
# ------------------------------------------------------------------------------
# Boundary Conditions and Function Spaces
# ------------------------------------------------------------------------------
def bot_left_corner(x):
return np.logical_and(np.isclose(x[0], 0.0), np.isclose(x[1], 0.0))
# Space for the periodic fluctuation field 'u'
V = fem.functionspace(domain, ("P", 1, (dim,)))
# Fix displacement at bottom-left corner (to remove rigid body motions)
pinned_dofs = fem.locate_dofs_geometrical(V, bot_left_corner)
u_zero = fem.Constant(domain, np.array([0.0, 0.0], dtype=PETSc.ScalarType))
bcs = [] # [fem.dirichletbc(u_zero, pinned_dofs, V)]
# Initialize MultiPointConstraint
tol = 1e-10
def PeriodicBoundary(axis):
"""
Full surface minus dofs constrained by BCs
"""
def complicated_mpc(x):
condition = np.isclose(x[axis], length_plate, atol=tol)
if axis == 0:
return condition
elif axis == 1:
return np.logical_and(
condition, np.logical_and(x[0] > tol, x[0] < length_plate - tol)
)
elif axis == 2:
return np.logical_and(
condition,
np.logical_and(
np.logical_and(x[0] > tol, x[0] < length_plate - tol),
np.logical_and(x[1] > tol, x[1] < length_plate - tol),
),
)
else:
raise RuntimeError("Axis should be 0, 1, or 2")
return complicated_mpc
def create_periodic_map(axis):
def periodic_relation(x):
out_x = x.copy()
out_x[axis] = 0
return out_x
return periodic_relation
mpc = MultiPointConstraint(V)
for axis in range(2):
mpc.create_periodic_constraint_geometrical(
V, PeriodicBoundary(axis), create_periodic_map(axis), bcs, tol=1e-10
)
mpc.finalize()
# ------------------------------------------------------------------------------
# Variational Formulation (Fluctuation Method)
# ------------------------------------------------------------------------------
V_mpc = mpc.function_space
u = fem.Function(V_mpc, name="Fluctuation")
# Macroscopic Deformation Gradient
F_macro_data = np.array([[1.0, 0.0], [0.0, 1.0]], dtype=PETSc.ScalarType)
F_macro = fem.Constant(domain, F_macro_data)
# Kinematics: Total F = F_macro + grad(u_fluctuation)
F = variable(F_macro + grad(u))
# Material Properties
E_1 = 20.0
nu_1 = 0.3
mu_1_val = E_1 / (2 * (1 + nu_1))
lmbda_1_val = E_1 * nu_1 / ((1 + nu_1) * (1 - 2 * nu_1))
# Strain Energy Density Function (Neo-Hookean)
C = F.T * F
I1 = tr(C)
J = det(F)
psi = (mu_1_val / 2) * (I1 - 3 - 2 * ln(J)) + (lmbda_1_val / 2) * (J - 1) ** 2
# Total Potential Energy
dx = ufl.Measure("dx", domain=domain)
E_pot = psi * dx
# Derivatives
v = TestFunction(V)
du = TrialFunction(V)
Residual = derivative(E_pot, u, v)
Jacobian = derivative(Residual, u, du)
# Problem Setup
problem = NonlinearProblem(
Residual,
u,
mpc,
bcs=bcs,
J=Jacobian,
petsc_options={
"snes_monitor": None,
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"mat_mumps_icntl_24": 1,
"mat_mumps_icntl_25": 0,
"snes_type": "newtonls",
"snes_linesearch_type": "none",
"snes_monitor": None,
"snes_rtol": 1e-8,
"snes_atol": 1e-10,
"snes_max_it": 50,
},
)
solver = problem.solver
# specific solvers omitted from MWE - they appear to give the same issues as default
# ------------------------------------------------------------------------------
# Simulation Loop and Output
# ------------------------------------------------------------------------------
# Aux space for visualization of total displacement
W = fem.functionspace(domain, ("P", 1, (dim,)))
u_total_vis = fem.Function(W, name="Total_Displacement")
x_coords = fem.Function(W)
x_coords.interpolate(lambda x: x[:2]) # Capture reference coords
# Output file
out_file = "MWE_output.xdmf"
with io.XDMFFile(domain.comm, out_file, "w") as xdmf:
xdmf.write_mesh(domain)
for n, eps in enumerate(np.linspace(0, 0.4, 21)):
if n == 0:
continue
# Update Macroscopic Gradient (Uniaxial Compression in Y)
current_F_macro = np.array(
[[1.0, 0.0], [0.0, 1.0 - eps]], dtype=PETSc.ScalarType
)
F_macro.value[:] = current_F_macro
print(f"Step {n}: eps = {eps:.4f}")
# Solve the nonlinear problem
u_mpc, converged, iters = problem.solve()
if converged > 0:
print(f"Converged reason: {converged} in {iters} iterations.")
else:
print(f"Failed to converge. Reason: {converged}, after {iters} iterations.")
# Update ghost values
u.x.scatter_forward()
# Reconstruct Total Displacement: u_total = u_fluc + (F_macro - I) * X
F_minus_I = current_F_macro - np.eye(2)
x_array = x_coords.x.array.reshape((-1, 2))
affine_disp = x_array @ F_minus_I.T
u_total_vis.x.array[:] = u.x.array + affine_disp.flatten()
# Save to file
xdmf.write_function(u_total_vis, n)
which yields a solution for every step and the final solution
which looks better, but something weird (although periodic) happens at the corners. I don’t have more time to look at this atm.