Using mpc with ufl.MixedFunctionSpace

General
1

I am enforcing continuity across non-conformal interface using MPC and need to enforce further constraints with real elements. Following the advice here related to the main branch, I constructed a model but it seems to be running into problems with form compilation. I will appreciate if someone familiar could point out my mistake. The MWE below (dolfinx 0.10, conda build) reproduces the error:

import dolfinx, gmsh, ufl, mpi4py, petsc4py, pyvista, scifem
import numpy as np
import matplotlib.pyplot as plt
import pyvista, dolfinx_mpc

gmsh.finalize()
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
gmsh.clear()
# gmsh.model.add("generator_2d")

occ = gmsh.model.occ
gdim = 2

b1 = occ.addRectangle(0, 0, 0, 0.5, 1)
b2 = occ.addRectangle(0.5, 0, 0, 0.5, 1)
occ.synchronize()

all_doms = gmsh.model.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

# number all boundaries
all_edges = gmsh.model.getEntities(gdim - 1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

gmsh.model.mesh.generate(gdim)
# for _ in range(refine):
#     gmsh.model.mesh.refine()
# gmsh.model.mesh.setOrder(fem_order)
gmsh.model.mesh.optimize()

mpi_rank = 0
dolfinx_model = dolfinx.io.gmsh.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, mpi_rank, gdim)
mesh, ct, ft = dolfinx_model.mesh, dolfinx_model.cell_tags, dolfinx_model.facet_tags

interface_bnd_idx = [2, 8] # by visual inspection of the mesh

Omega_A = dolfinx.fem.functionspace(mesh, ("CG", 2))
Omega_r = scifem.create_real_functionspace(mesh)

# mpc
mpc = dolfinx_mpc.MultiPointConstraint(Omega_A)
tol = float(1e-8)
mpc.create_contact_inelastic_condition(
    ft, interface_bnd_idx[0], interface_bnd_idx[1], eps2=tol)
mpc.finalize()

W = ufl.MixedFunctionSpace(mpc.function_space, Omega_r)
A, Vc = ufl.TrialFunctions(W)
At, Vct = ufl.TestFunctions(W)

dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct)

zero = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0))
F = ufl.inner(A, At)*dx + ufl.inner(Vc, Vct)*dx + ufl.inner(A, Vct)*dx + ufl.inner(Vc, At)*dx
F += ufl.inner(zero, Vct)*dx + ufl.inner(zero, At)*dx

left_dofs = dolfinx.fem.locate_dofs_geometrical(Omega_A, lambda x: np.isclose(x[0], 0))
dbc = dolfinx.fem.dirichletbc(1.0, left_dofs, Omega_A)

a_block, L_block = ufl.extract_blocks(F)
a_form = dolfinx.fem.form(a_block)
L_form = dolfinx.fem.form(L_block)

# solve the linear model
petsc_options = {"pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
problem = dolfinx_mpc.problem.LinearProblem(a_form, L_form, [mpc, None], bcs=[dbc], petsc_options=petsc_options)
A_sol = problem.solve()

The form compiler leads to a long error where I suspect the following is relevant:

---------------------------------------------------------------------------
ArityMismatch                             Traceback (most recent call last)
Cell In[17], line 64
     61 dbc = dolfinx.fem.dirichletbc(1.0, left_dofs, Omega_A)
     63 a_block, L_block = ufl.extract_blocks(F)
---> 64 a_form = dolfinx.fem.form(a_block)
     65 L_form = dolfinx.fem.form(L_block)
     67 # solve the linear model

File ~/install/miniforge3/envs/fenics-0.10/lib/python3.11/site-packages/dolfinx/fem/forms.py:449, in form(form, dtype, form_compiler_options, jit_options, entity_maps)
    446     else:
    447         return form
--> 449 return _create_form(form)

File ~/install/miniforge3/envs/fenics-0.10/lib/python3.11/site-packages/dolfinx/fem/forms.py:445, in form.<locals>._create_form(form)
    443     return _zero_form(form)
    444 elif isinstance(form, collections.abc.Iterable):
--> 445     return list(map(lambda sub_form: _create_form(sub_form), form))
    446 else:
    447     return form

File ~/install/miniforge3/envs/fenics-0.10/lib/python3.11/site-packages/dolfinx/fem/forms.py:445, in form.<locals>._create_form.<locals>.<lambda>(sub_form)
    443     return _zero_form(form)
    444 elif isinstance(form, collections.abc.Iterable):
--> 445     return list(map(lambda sub_form: _create_form(sub_form), form))
...
     60         f"{tuple(map(_afmt, a))} vs {tuple(map(_afmt, b))}."
     61     )
     62 return a

ArityMismatch: Adding expressions with non-matching form arguments ('v_0^0',) vs ('v_0^0', 'v_1^0').
2

There are several mistakes in your code.

This is not correct, should use ufl.system and then ufl.extract_blocks.

Need to pass in an empty constraint for second variable.

Here is something that runs:

import dolfinx, gmsh, ufl, mpi4py, scifem
import numpy as np
import  dolfinx_mpc

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
gmsh.clear()
# gmsh.model.add("generator_2d")

occ = gmsh.model.occ
gdim = 2

b1 = occ.addRectangle(0, 0, 0, 0.5, 1)
b2 = occ.addRectangle(0.5, 0, 0, 0.5, 1)
occ.synchronize()

all_doms = gmsh.model.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

# number all boundaries
all_edges = gmsh.model.getEntities(gdim - 1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

gmsh.model.mesh.generate(gdim)
# for _ in range(refine):
#     gmsh.model.mesh.refine()
# gmsh.model.mesh.setOrder(fem_order)
gmsh.model.mesh.optimize()

mpi_rank = 0
dolfinx_model = dolfinx.io.gmsh.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, mpi_rank, gdim)
mesh, ct, ft = dolfinx_model.mesh, dolfinx_model.cell_tags, dolfinx_model.facet_tags

interface_bnd_idx = [2, 8] # by visual inspection of the mesh

Omega_A = dolfinx.fem.functionspace(mesh, ("CG", 2))
Omega_r = scifem.create_real_functionspace(mesh)

# mpc
mpc = dolfinx_mpc.MultiPointConstraint(Omega_A)
tol = float(1e-8)
mpc.create_contact_inelastic_condition(
    ft, interface_bnd_idx[0], interface_bnd_idx[1], eps2=tol)
mpc.finalize()
mpc_r  = dolfinx_mpc.MultiPointConstraint(Omega_r)
mpc_r.finalize()

W = ufl.MixedFunctionSpace(mpc.function_space, Omega_r)
A, Vc = ufl.TrialFunctions(W)
At, Vct = ufl.TestFunctions(W)

dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct)

zero = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0))
F = ufl.inner(A, At)*dx + ufl.inner(Vc, Vct)*dx + ufl.inner(A, Vct)*dx + ufl.inner(Vc, At)*dx
one = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(1))
F += ufl.inner(one, Vct)*dx + ufl.inner(zero, At)*dx

left_dofs = dolfinx.fem.locate_dofs_geometrical(Omega_A, lambda x: np.isclose(x[0], 0))
dbc = dolfinx.fem.dirichletbc(1.0, left_dofs, Omega_A)

a_block, L_block = ufl.system(F)
a_form =  ufl.extract_blocks(a_block)
L_form = ufl.extract_blocks(L_block)

# solve the linear model
petsc_options = {"pc_type": "lu", "pc_factor_mat_solver_type": "mumps", "ksp_error_if_not_converged": True,"ksp_monitor": None,
                 "ksp_type": "preonly"}
problem = dolfinx_mpc.problem.LinearProblem(a_form, L_form, [mpc,mpc_r], bcs=[dbc], petsc_options=petsc_options)
A_sol = problem.solve()
print(A_sol[0].x.array, A_sol[1].x.array)

Note that I change the volume force from 0 to 1.
Still not sure if this is the expected solution of the problem.
I would need the full variational form to verify what you are doing (in mathematical syntax).

3

Thanks a lot for the quick fix! The weak form in the shared MWE is just a filler to reproduce the error. The actual problem is that of a rotating electric generator which is working with this correction.

While we are on this, can I quickly ask if there are alternatives for enforcing field (dis)continuity across non-conformal interfaces?

4

There probably is, like nitsche/mortar methods. But they are quite hard to implement. What specific problem are you considering?

5

I am eyeing moving (rotating, sliding) electric machines in 3D. Generally the sliding interface is in air (or non-conducting) domain which is modeled with scalar magnetic potential (CG, Poisson equation). I believe MPC is useful here to enforce continuity. But I would like to push the boundaries by considering sliding contact between two metals, a problem beyond the commercial tools I know. Due to current flow, the sliding domains have to be considered with Ampere’s law and edge elements. I would appreciate if I could get referral to existing fenics examples involving such interfaces irrespective of the context. Once I understand how to work with fenics I can search for appropriate mathematical formulations on my own.