Thank you for your answer. Unfortunately, this didn't help.
I don't think it is related to the Lagrange term because problem persists even when I use only "m" and "phi".
I've created another MWE that includes two domains: in one, linear elasticity is coupled with magnetostatics, and in the other, it's just linear elasticity. I created the function space for displacement on the full mesh, and a separate function space for the magnetic potential on the submesh. When I use the coupling across the whole domain without a mixed function space, everything works. But when I try to couple them across multiple domains, I get the same error.
MWE
import ufl
import numpy as np
from dolfinx import mesh, fem, default_scalar_type
from mpi4py import MPI
from dolfinx import io
from scifem import BlockedNewtonSolver
from customNewtonSolver import *
import numpy.typing as npt
nx, ny, nz = 20, 20, 20 # Number of elements in x and y directions
points = [[0., 0., 0.], [1., 1., 1.]] # Domain corners
nElement = [nx, ny, nz]
msh = mesh.create_box(MPI.COMM_WORLD, points, nElement, mesh.CellType.tetrahedron)
def up(x, tol=1e-14):
return x[2] >= 0.5 - tol
tdim = msh.topology.dim
cell_map = msh.topology.index_map(tdim)
num_cells_local = cell_map.size_local + cell_map.num_ghosts
marker = np.ones(num_cells_local, dtype=np.int32)
marker[mesh.locate_entities(msh, tdim, up)] = 2
cell_tags = mesh.meshtags(msh, msh.topology.dim, np.arange(num_cells_local, dtype=np.int32), marker)
sm, sm_cell_map, sm_vertex_map, _ = mesh.create_submesh(msh, cell_tags.dim, cell_tags.find(2))
sm.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)
mesh_to_heat_entity = np.full(num_cells_local, -1, dtype=np.int32)
mesh_to_heat_entity[sm_cell_map] = np.arange(len(sm_cell_map), dtype=np.int32)
entity_maps = {sm: mesh_to_heat_entity}
# Material properties
E = fem.Constant(msh, default_scalar_type(21000.0))
nu = fem.Constant(msh, default_scalar_type(3.0e-1))
mu = E / (2.0 * (1.0 + nu))
lmbda = E * nu / ((1 - nu) * (1 + nu))
fs_u = fem.functionspace(msh, ("Lagrange", 1, (msh.geometry.dim, )))
fs_phi = fem.functionspace(sm, ("Lagrange", 1))
V = ufl.MixedFunctionSpace(fs_u, fs_phi)
u = fem.Function(fs_u, name="Magnetization")
phi = fem.Function(fs_phi, name="Magnetic Potential")
du, dphi = ufl.TrialFunctions(V)
eps = ufl.sym(ufl.grad(u))
eps_var = ufl.variable(ufl.sym(ufl.grad(u)))
deps = ufl.sym(ufl.grad(du))
H = -ufl.grad(phi)
H_var = ufl.variable(-ufl.grad(phi))
dH = -ufl.grad(dphi)
a=ufl.as_vector([1, 0, 0])
psi_phi = - 0.5 * ufl.dot(H, H)
psi_u = (0.5 * lmbda * ufl.tr(eps) * ufl.tr(eps) + mu * ufl.inner(eps, eps))
psi_mix = ufl.tr(eps)*ufl.dot(H, H) + ufl.tr(ufl.dot(eps, ufl.outer(a,a)))*ufl.dot(H, H) + ufl.tr(ufl.dot(eps, ufl.outer(H, a)))
psi = psi_phi + psi_u + psi_mix
sig = ufl.diff(psi, eps_var)
B = -ufl.diff(psi, H_var)
fdim = msh.topology.dim - 1
facets_left = mesh.locate_entities_boundary(msh, fdim, lambda x: np.isclose(x[0], 0.0))
facets_right = mesh.locate_entities_boundary(msh, fdim, lambda x: np.isclose(x[0], 1.0))
dofs_left_u = fem.locate_dofs_topological(fs_u.sub(0), fdim, facets_left)
dofs_right_u = fem.locate_dofs_topological(fs_u.sub(0), fdim, facets_right)
bc_u_left = fem.dirichletbc(0.0, dofs_left_u, fs_u.sub(0))
bc_u_right = fem.dirichletbc(1.0, dofs_right_u, fs_u.sub(0))
bcs_u = [bc_u_left, bc_u_right]
dx = ufl.Measure("dx", domain=msh, subdomain_data=cell_tags, subdomain_id=2)
F = ufl.dot(dH, B) * dx(2) + ufl.inner(deps, sig) * dx
residual = fem.form(ufl.extract_blocks(F), entity_maps=entity_maps)
a, b = ufl.TestFunctions(V)
jac = ufl.derivative(F, u, a) + ufl.derivative(F, phi, b)
J = fem.form(ufl.extract_blocks(jac), entity_maps=entity_maps)
solver = customNewtonSolver(residual, J, [u, phi], bcs_u)
solver.solve()