Dolfinx_mpc on a curved interface

dolfinx
1

I am trying to rework the disconnected elasticiy demo of dolfinx_mpc (0.7) on a curved boundary. The mesh across the boundary is non-conformal on purpose. Is the mpc expected to work for this geometry? The result I am getting with my attempt is clearly not working. But then I might have done it wrong. Below is the code

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

def create_geom(ro=50, fov=None, disp=False):

    if fov is None:
        fov = (400, 300)
    fov_x, fov_y = fov
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0) # disable output messages
    gmsh.clear()
    gdim = 2
    model_rank = 0
    
    t_gap = 5 # thickness of airgap
    
    occ = gmsh.model.occ

    # first we create the outer domain
    outer_box = occ.add_rectangle(-fov_x/2, -fov_y/2, 0, fov_x, fov_y)
    occ.synchronize()
    air_shell = occ.addDisk(0, 0, 0, ro+t_gap, ro+t_gap)
    dom1 = occ.cut([(gdim, outer_box)], [(gdim, air_shell)]) # one domain

    # now the inner domain which is the wheel
    dom2 = occ.addDisk(0, 0, 0, ro+t_gap, ro+t_gap)

    occ.synchronize()
    # number all domains
    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]], j + 1)  # create the main group/node

    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.refine()
    gmsh.model.mesh.refine()
    # gmsh.write('geom/mwe_mixed_domains.msh')
    model_rank = 0
    mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)
    if disp:
        gmsh.fltk.run() # visible inspection, for debugging
    gmsh.finalize()
    return mesh, ct, ft

mesh, ct, ft = create_geom(disp=False)

dom_idx_rot = 2 # rotating/circular domain
dom_idx_stat = 1 # stationary domain
bnd_idx_cont = (1, 6) # boundary where continuity is imposed, stationary and rotating sides respectively
bnd_idx_l = (3, ) # left boundary
bnd_idx_r = (4, ) # right boundary

fem_degree = 2
Omega = dolfinx.fem.functionspace(mesh, ("Lagrange", fem_degree))

u = ufl.TrialFunction(Omega)
v = ufl.TestFunction(Omega)

F = ufl.dot(ufl.grad(u), ufl.grad(v))*ufl.dx - ufl.dot(dolfinx.fem.Constant(mesh, 0.0), v)*ufl.dx

a = ufl.lhs(F)
L = ufl.rhs(F)

# dirichlet boundary conditions at the left and right edges

dofs_l = dolfinx.fem.locate_dofs_topological(Omega, ft.dim, ft.find(bnd_idx_l[0]))
dbc_l = dolfinx.fem.dirichletbc(1.0, dofs_l, Omega)

dofs_r = dolfinx.fem.locate_dofs_topological(Omega, ft.dim, ft.find(bnd_idx_r[0]))
dbc_r = dolfinx.fem.dirichletbc(0., dofs_r, Omega)

def gather_dof_coordinates(V, dofs: np.ndarray):
    """
    Distributes the dof coordinates of this subset of dofs to all processors
    """
    x = V.tabulate_dof_coordinates()
    local_dofs = dofs[dofs < V.dofmap.index_map.size_local * V.dofmap.index_map_bs]
    coords = x[local_dofs]
    num_nodes = len(coords)
    glob_num_nodes = mpi4py.MPI.COMM_WORLD.allreduce(num_nodes, op=mpi4py.MPI.SUM)
    recvbuf = None
    if mpi4py.MPI.COMM_WORLD.rank == 0:
        recvbuf = np.zeros(3 * glob_num_nodes, dtype=V.mesh.geometry.x.dtype)
    sendbuf = coords.reshape(-1)
    sendcounts = np.array(mpi4py.MPI.COMM_WORLD.gather(len(sendbuf), 0))
    mpi4py.MPI.COMM_WORLD.Gatherv(sendbuf, (recvbuf, sendcounts), root=0)
    glob_coords = mpi4py.MPI.COMM_WORLD.bcast(recvbuf, root=0).reshape((-1, 3))
    return glob_coords

dofs_stat = dolfinx.fem.locate_dofs_topological(Omega, ft.dim, ft.find(bnd_idx_cont[0]))
dofs_rot = dolfinx.fem.locate_dofs_topological(Omega, ft.dim, ft.find(bnd_idx_cont[1]))

# Given the local coordinates of the dofs, distribute them on all processors
nodes = [gather_dof_coordinates(Omega, dofs_stat), gather_dof_coordinates(Omega, dofs_rot)]
pairs = []
# points are on the circumference across the interface
for x in nodes[0]:
    theta1 = np.arctan2(x[1], x[0])
    for y in nodes[1]:
        theta2 = np.arctan2(y[1], y[0])
        if np.isclose(theta1, theta2):
            pairs.append([x, y])
            break

from dolfinx_mpc.utils import create_point_to_point_constraint
mpc = dolfinx_mpc.MultiPointConstraint(Omega)
for i, pair in enumerate(pairs):
    sl, ms, co, ow, off = create_point_to_point_constraint(
        Omega, pair[0], pair[1])
    mpc.add_constraint(Omega, sl, ms, co, ow, off)
mpc.finalize()

bcs = [dbc_l, dbc_r]
# petsc_options = {"ksp_type": "preonly", "pc_type": "lu"} # does not work without dolfinx_mpc
petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}

from dolfinx_mpc import LinearProblem
problem = LinearProblem(a, L, mpc, bcs=bcs, petsc_options=petsc_options)
uh = problem.solve()

import matplotlib.tri as tri
import matplotlib.animation as animation

fig1, ax1 = plt.subplots(constrained_layout=True, figsize=(4.5, 3.5))
ax1.set_aspect('equal')
triang = tri.Triangulation(uh.function_space.tabulate_dof_coordinates()[:,0], 
                           uh.function_space.tabulate_dof_coordinates()[:,1])
tcf = ax1.tricontourf(triang, uh.x.array, 100, cmap='hot')
plt.colorbar(tcf)
2

First of all, this is wrong, as you are mapping the top to the bottom, not the interfaces to eachother

i.e. you are using 1 and 6 and not 5 and 6 which is your mesh-tags.

Secondly, using create_point_to_point_constraint with identical coordinates is not going to work out.
See for instance: dolfinx_mpc/python/demos/demo_contact_3D.py at main · jorgensd/dolfinx_mpc · GitHub
which creates a 1-1 constraint between to decoupled meshes if you need vector elements, or:
dolfinx_mpc/python/demos/demo_periodic3d_topological.py at main · jorgensd/dolfinx_mpc · GitHub
for mapping facet to facet.

Edit
Testing the constraints for your code, I need to make some minor adjustments to your code to make it work. The easiest fix would be to fix: `FiniteElement::is_mixed` is wrong for 1-dimensional vectors · Issue #2944 · FEniCS/dolfinx · GitHub
and once this is fixed, use something along the lines of:

from mpi4py import MPI
import basix.ufl
from dolfinx.io import XDMFFile
from dolfinx.io import VTXWriter
from dolfinx_mpc import LinearProblem
import dolfinx
import mpi4py
import ufl
import gmsh
import dolfinx_mpc
from dolfinx.io import gmshio
import numpy as np
import matplotlib.pyplot as plt


def create_geom(ro=50, fov=None, disp=False):

    if fov is None:
        fov = (400, 300)
    fov_x, fov_y = fov
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)  # disable output messages
    gmsh.clear()
    gdim = 2
    model_rank = 0

    t_gap = 5  # thickness of airgap

    occ = gmsh.model.occ

    # first we create the outer domain
    outer_box = occ.add_rectangle(-fov_x/2, -fov_y/2, 0, fov_x, fov_y)
    occ.synchronize()
    air_shell = occ.addDisk(0, 0, 0, ro+t_gap, ro+t_gap)
    dom1 = occ.cut([(gdim, outer_box)], [(gdim, air_shell)])  # one domain

    # now the inner domain which is the wheel
    dom2 = occ.addDisk(0, 0, 0, ro+t_gap, ro+t_gap)

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

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

    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.refine()
    gmsh.model.mesh.refine()
    # gmsh.write('geom/mwe_mixed_domains.msh')
    model_rank = 0
    mesh, ct, ft = gmshio.model_to_mesh(
        gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)
    if disp:
        gmsh.fltk.run()  # visible inspection, for debugging
    gmsh.finalize()
    return mesh, ct, ft


mesh, ct, ft = create_geom(disp=False)

dom_idx_rot = 2  # rotating/circular domain
dom_idx_stat = 1  # stationary domain
# boundary where continuity is imposed, stationary and rotating sides respectively
bnd_idx_cont = (5, 6)
bnd_idx_l = (3, )  # left boundary
bnd_idx_r = (4, )  # right boundary

fem_degree = 2
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
el = basix.ufl.element(
    "Lagrange", mesh.topology.cell_name(), fem_degree, shape=(1, ))
Omega = dolfinx.fem.functionspace(mesh, el)
u = ufl.TrialFunction(Omega)
v = ufl.TestFunction(Omega)

F = ufl.inner(ufl.grad(u), ufl.grad(v))*ufl.dx - \
    ufl.inner(dolfinx.fem.Constant(mesh, (0.0,)), v)*ufl.dx

a = ufl.lhs(F)
L = ufl.rhs(F)

# dirichlet boundary conditions at the left and right edges

dofs_l = dolfinx.fem.locate_dofs_topological(
    Omega, ft.dim, ft.find(bnd_idx_l[0]))
cl = dolfinx.fem.Constant(mesh, (1.0,))
dbc_l = dolfinx.fem.dirichletbc(cl, dofs_l, Omega)
cr = dolfinx.fem.Constant(mesh, (0.0,))
dofs_r = dolfinx.fem.locate_dofs_topological(
    Omega, ft.dim, ft.find(bnd_idx_r[0]))
dbc_r = dolfinx.fem.dirichletbc(cr, dofs_r, Omega)

mpc = dolfinx_mpc.MultiPointConstraint(Omega)
tol = float(1e-11)
bcs = [dbc_l, dbc_r]
mpc.create_contact_inelastic_condition(
    ft, bnd_idx_cont[0], bnd_idx_cont[1], eps2=tol)
mpc.finalize()

# petsc_options = {"ksp_type": "preonly", "pc_type": "lu"} # does not work without dolfinx_mpc
petsc_options = {"ksp_type": "preonly", "pc_type": "lu",
                 "pc_factor_mat_solver_type": "mumps"}

problem = LinearProblem(a, L, mpc, bcs=bcs, petsc_options=petsc_options)
uh = problem.solve()


with XDMFFile(mesh.comm, "mf.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(ft, mesh.geometry)


with VTXWriter(mesh.comm, "uh.bp", [uh], engine="BP4") as vtx:
    vtx.write(0.0)
3

I’ve also filed an issue in my github repo, and will have a look at it early in January

4

@sbhasan I’ve made a pull request with fixes for your case:

With that branch you can run the supplied code and get the expected results.

I will merge it into main once the CI is passing.

1 Like
5

@sbhasan it is now part of dolfinx_mpc main branch, which is compatible with Dolfinx main branch.

6

Hi,

I honestly cannot thank you enough!

Meanwhile, I am struggling to acquire your fix. I operate under a conda environment in which “pip install .” on the cloned git repo runs into issues with petsc. Is it possible to acquire the main branch through conda?

7

No, conda doesn’t work on live branches, only on tags.

I’ll see if I can backport the fix to 0.7.x and make a new release. However, this will not happen tonight, as we are closing in on midnight, and probably not until sometime next week.

1 Like
8

Whenever it’s convenient, this is way beyond what I could hope for. Please possibly drop a line here as soon as it is retrievable from conda-forge

9

Hi,

Any chance this could be introduced to conda in the coming days? Otherwise, can we acquire it through docker? I searched for an image but the dolfinx_mpc seems to be accessible up to 0.5.0 only.

Thanks a lot in advance!

10

Not sure, dealing with some other conda issues atm.

Soon, yes, you can use the image ghcr.io/fenics/dolfinx/dolfinx:nightly (Remember to pull it first) and then install the main branch of dolfinx_mpc. I need to deal with some api changes that were introduced friday before this will work.

Dolfinx_mpc image are available up to 0.7.1 at Package dolfinx_mpc · GitHub

1 Like
11

@sbhasan I’ve now backported the fix to DOLFINx v0.7.x.

There will be a docker image with the tag v0.7.2 at: Package dolfinx_mpc · GitHub in about an hour.

1 Like
12

Thanks a lot. The scalar version finally works with the docker image!

13

Just for readers in future, here is the working code

from mpi4py import MPI
import basix.ufl
from dolfinx.io import XDMFFile
from dolfinx.io import VTXWriter
from dolfinx_mpc import LinearProblem
import dolfinx
import mpi4py
import ufl
import gmsh
import dolfinx_mpc
from dolfinx.io import gmshio
import numpy as np
import matplotlib.pyplot as plt


def create_geom(ro=50, fov=None, disp=False):

    if fov is None:
        fov = (400, 300)
    fov_x, fov_y = fov
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)  # disable output messages
    gmsh.clear()
    gdim = 2
    model_rank = 0

    t_gap = 5  # thickness of airgap

    occ = gmsh.model.occ

    # first we create the outer domain
    outer_box = occ.add_rectangle(-fov_x/2, -fov_y/2, 0, fov_x, fov_y)
    occ.synchronize()
    air_shell = occ.addDisk(0, 0, 0, ro+t_gap, ro+t_gap)
    dom1 = occ.cut([(gdim, outer_box)], [(gdim, air_shell)])  # one domain

    # now the inner domain which is the wheel
    dom2 = occ.addDisk(0, 0, 0, ro+t_gap, ro+t_gap)

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

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

    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.refine()
    gmsh.model.mesh.refine()
    # gmsh.write('geom/mwe_mixed_domains.msh')
    model_rank = 0
    mesh, ct, ft = gmshio.model_to_mesh(
        gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)
    if disp:
        gmsh.fltk.run()  # visible inspection, for debugging
    gmsh.finalize()
    return mesh, ct, ft


mesh, ct, ft = create_geom(disp=False)

dom_idx_rot = 2  # rotating/circular domain
dom_idx_stat = 1  # stationary domain
# boundary where continuity is imposed, stationary and rotating sides respectively
bnd_idx_cont = (5, 6)
bnd_idx_l = (2, )  # left boundary
bnd_idx_r = (3, )  # right boundary

fem_degree = 2
# mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
el = basix.ufl.element(
    "Lagrange", mesh.topology.cell_name(), fem_degree, shape=(1, ))
Omega = dolfinx.fem.functionspace(mesh, el)
Omega = dolfinx.fem.functionspace(mesh, ("Lagrange",fem_degree))
u = ufl.TrialFunction(Omega)
v = ufl.TestFunction(Omega)

F = ufl.inner(ufl.grad(u), ufl.grad(v))*ufl.dx - \
    ufl.inner(dolfinx.fem.Constant(mesh, 0.0), v)*ufl.dx

a = ufl.lhs(F)
L = ufl.rhs(F)

# dirichlet boundary conditions at the left and right edges

dofs_l = dolfinx.fem.locate_dofs_topological(
    Omega, ft.dim, ft.find(bnd_idx_l[0]))
cl = dolfinx.fem.Constant(mesh, 1.0)

dbc_l = dolfinx.fem.dirichletbc(cl, dofs_l, Omega)
cr = dolfinx.fem.Constant(mesh, 0.0)
dofs_r = dolfinx.fem.locate_dofs_topological(
    Omega, ft.dim, ft.find(bnd_idx_r[0]))
dbc_r = dolfinx.fem.dirichletbc(cr, dofs_r, Omega)

mpc = dolfinx_mpc.MultiPointConstraint(Omega)
tol = float(1e-11)
bcs = [dbc_l, dbc_r]
mpc.create_contact_inelastic_condition(
    ft, bnd_idx_cont[0], bnd_idx_cont[1], eps2=tol)
mpc.finalize()

# petsc_options = {"ksp_type": "preonly", "pc_type": "lu"} # does not work without dolfinx_mpc
petsc_options = {"ksp_type": "preonly", "pc_type": "lu",
                 "pc_factor_mat_solver_type": "mumps"}

problem = LinearProblem(a, L, mpc, bcs=bcs, petsc_options=petsc_options)
uh = problem.solve()

with XDMFFile(mesh.comm, "mf.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(ft, mesh.geometry)

with VTXWriter(mesh.comm, "uh.bp", [uh], engine="BP4") as vtx:
    vtx.write(0.0)
14

A conda-version for v0.7.2 has been released (Dolfinx Mpc :: Anaconda.org) thanks to @minrk.

1 Like
15

Hi @dokken

I need to redo the problem of this thread for Nedelec elements and in 3D if it matters. Does mpc.create_contact_inelastic_condition still work? If the reminder helps, I need to impose field continuity at overlapping boundaries between two non-conformal meshes.

Thanks

16

Nedelec elements are not currently supported by most of the MPC constructors, as the Finite element dual function is not defined through point evaluations. You could construct your own map that restricts one dof to another, but i am not sure it would make sense for non-conforming interfaces.

2 Likes
17

Is there a possibility of doing it using legacy fenics with whatever tools it has?

18

If the meshes are conforming it might work in legacy fenics.

19

Hi, I tried to run your code with dolfinx-0.8.0, dolfinx_mpc-0.8.1 and gmsh-4.13.0. For me it results in the error:

No masters found on contact surface (when executed in serial). Please make sure that the surfaces are in contact, or increase the tolerance eps2.

I have not been able to fix it, could you confirm if your example is still running with the latest version of dolfinx and dolfinx_mpc?

20

The issue is that the gmsh tags have not been consistently marked. Thus changing the GMSH version has changed what facets are marked with which tag. Here is a working mwe for v0.8.0:

from mpi4py import MPI
import basix.ufl
from dolfinx.io import XDMFFile
from dolfinx.io import VTXWriter
from dolfinx_mpc import LinearProblem
import dolfinx
import mpi4py
import ufl
import gmsh
import dolfinx_mpc
from dolfinx.io import gmshio
import numpy as np
import matplotlib.pyplot as plt


def create_geom(ro=50, fov=None, disp=False):
    if fov is None:
        fov = (400, 300)
    fov_x, fov_y = fov
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)  # disable output messages
    gmsh.clear()
    gdim = 2
    model_rank = 0

    t_gap = 5  # thickness of airgap

    occ = gmsh.model.occ

    # first we create the outer domain
    outer_box = occ.add_rectangle(-fov_x / 2, -fov_y / 2, 0, fov_x, fov_y)
    occ.synchronize()
    air_shell = occ.addDisk(0, 0, 0, ro + t_gap, ro + t_gap)
    dom1 = occ.cut([(gdim, outer_box)], [(gdim, air_shell)])  # one domain

    # now the inner domain which is the wheel
    dom2 = occ.addDisk(0, 0, 0, ro + t_gap, ro + t_gap)

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

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

    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.refine()
    gmsh.model.mesh.refine()
    # gmsh.write('geom/mwe_mixed_domains.msh')
    model_rank = 0
    mesh, ct, ft = gmshio.model_to_mesh(
        gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim
    )
    if disp:
        gmsh.fltk.run()  # visible inspection, for debugging
    gmsh.finalize()
    return mesh, ct, ft


mesh, ct, ft = create_geom(disp=False)
# with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "w") as xdmf:
#     xdmf.write_mesh(mesh)
#     xdmf.write_meshtags(ft, mesh.geometry)
# exit()
dom_idx_rot = 2  # rotating/circular domain
dom_idx_stat = 1  # stationary domain
# boundary where continuity is imposed, stationary and rotating sides respectively
bnd_idx_cont = (1, 6)
bnd_idx_l = (3,)  # left boundary
bnd_idx_r = (4,)  # right boundary

fem_degree = 2
# mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
el = basix.ufl.element("Lagrange", mesh.topology.cell_name(), fem_degree, shape=(1,))
Omega = dolfinx.fem.functionspace(mesh, el)
Omega = dolfinx.fem.functionspace(mesh, ("Lagrange", fem_degree))
u = ufl.TrialFunction(Omega)
v = ufl.TestFunction(Omega)

F = (
    ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
    - ufl.inner(dolfinx.fem.Constant(mesh, 0.0), v) * ufl.dx
)

a = ufl.lhs(F)
L = ufl.rhs(F)

# dirichlet boundary conditions at the left and right edges

dofs_l = dolfinx.fem.locate_dofs_topological(Omega, ft.dim, ft.find(bnd_idx_l[0]))
cl = dolfinx.fem.Constant(mesh, 1.0)

dbc_l = dolfinx.fem.dirichletbc(cl, dofs_l, Omega)
cr = dolfinx.fem.Constant(mesh, 0.0)
dofs_r = dolfinx.fem.locate_dofs_topological(Omega, ft.dim, ft.find(bnd_idx_r[0]))
dbc_r = dolfinx.fem.dirichletbc(cr, dofs_r, Omega)

mpc = dolfinx_mpc.MultiPointConstraint(Omega)
tol = float(1e-11)
bcs = [dbc_l, dbc_r]
mpc.create_contact_inelastic_condition(ft, bnd_idx_cont[0], bnd_idx_cont[1], eps2=tol)
mpc.finalize()

# petsc_options = {"ksp_type": "preonly", "pc_type": "lu"} # does not work without dolfinx_mpc
petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
}

problem = LinearProblem(a, L, mpc, bcs=bcs, petsc_options=petsc_options)
uh = problem.solve()

with XDMFFile(mesh.comm, "mf.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(ft, mesh.geometry)

with VTXWriter(mesh.comm, "uh.bp", [uh], engine="BP4") as vtx:
    vtx.write(0.0)
1 Like