Manually apply Dirichlet BC after assembly for mixed problem

dolfinx
4

For future posts, I would suggest you test such modifications before asking for advice, as the code presented in the Sorbonne tutorial works for both mixed and non-mixed problems, as illustrated with the following MWE:

from mpi4py import MPI
from petsc4py import PETSc
import dolfinx.fem.petsc
import ufl
import basix.ufl
from importlib.metadata import version

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

if version("fenics-dolfinx") < "0.8.0":
    el = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
    El = ufl.MixedElement([el, el])
    V = dolfinx.fem.FunctionSpace(mesh, El)
else:
    el = basix.ufl.element("Lagrange", mesh.basix_cell(), 1)
    El = basix.ufl.mixed_element([el, el])
    V = dolfinx.fem.functionspace(mesh, El)

V0, _ = V.sub(0).collapse()
V1, _ = V.sub(1).collapse()

fdim = mesh.topology.dim - 1
mesh.topology.create_connectivity(fdim, fdim + 1)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)

u0 = dolfinx.fem.Function(V0)
u0.interpolate(lambda x: x[0])
u1 = dolfinx.fem.Function(V1)
u1.interpolate(lambda x: x[1])
dofs0 = dolfinx.fem.locate_dofs_topological((V.sub(0), V0), fdim, boundary_facets)
dofs1 = dolfinx.fem.locate_dofs_topological((V.sub(1), V1), fdim, boundary_facets)

bc0 = dolfinx.fem.dirichletbc(u0, dofs0, V.sub(0))
bc1 = dolfinx.fem.dirichletbc(u1, dofs1, V.sub(1))
bcs = [bc0, bc1]

u, v = ufl.TrialFunctions(V)
p, q = ufl.TestFunctions(V)

a = ufl.inner(u, p) * ufl.dx + ufl.inner(v, q) * ufl.dx
L = ufl.inner(0.1, p) * ufl.dx + ufl.inner(0.3, q) * ufl.dx

a_compiled = dolfinx.fem.form(a)
L_compiled = dolfinx.fem.form(L)

A = dolfinx.fem.petsc.assemble_matrix(a_compiled)
A.assemble()

for bc in bcs:
    dofs, _ = bc._cpp_object.dof_indices()
    A.zeroRowsLocal(dofs, diag=1)

b = dolfinx.fem.petsc.assemble_vector(L_compiled)
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(b, bcs)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)

uh = dolfinx.fem.Function(V)
solver = PETSc.KSP().create(mesh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
solver.getPC().setFactorSolverType("mumps")
solver.solve(b, uh.vector)
uh.vector.ghostUpdate(
    addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD
)

u0_out = uh.sub(0).collapse()
u0_out.name = "u0"
u1_out = uh.sub(1).collapse()
u1_out.name = "u1"
with dolfinx.io.VTXWriter(mesh.comm, "u0.bp", [u0_out], engine="BP5") as bp:
    bp.write(0.0)

with dolfinx.io.VTXWriter(mesh.comm, "u1.bp", [u1_out], engine="BP5") as bp:
    bp.write(0.0)

image
image1107×768 66.7 KB

1 Like