Manually apply Dirichlet BC after assembly for mixed problem

[fem_fatale]

Hello all,

I would like to know, if it is possible to apply Dirichlet boundary conditions to a mixed problem after assembly. I’ve seen this tutorial Application of Dirichlet boundary conditions — FEniCS Tutorial @ Sorbonne, where the identity row approach is described. I would like to do something like the non-blocked direct solver in this demo Stokes equations with Taylor-Hood elements — DOLFINx 0.6.0 documentation (I’m using v0.6.0), but also manually apply Dirichlet BC after assembly.

The reason is I’m using the Customquad library, where I cannot pass boundary conditions for assembly. A level set function defines two subdomains, for each I want to define a test and trial function (linear scalar Lagrange). On the interface I want to set conditions to couple both. Since I’m using the Cut FEM method, the functions should ideally be restricted to their subdomain, so the current approach would be to manually set zero Dirichlet bc on all respective outside dofs.

Would this be possible? Would it work in parallel?

I appreciate any help!

[dokken]

@fem_fatale

I would like to know, if it is possible to apply Dirichlet boundary conditions to a mixed problem after assembly. I’ve seen this tutorial Application of Dirichlet boundary conditions — FEniCS Tutorial @ Sorbonne, where the identity row approach is described. I would like to do something like the non-blocked direct solver in this demo Stokes equations with Taylor-Hood elements — DOLFINx 0.6.0 documentation (I’m using v0.6.0), but also manually apply Dirichlet BC after assembly.

Yes, this is is possible and should work in parallel.

[fem_fatale]

Dear dokken,

thank you very much for your answer. Can you tell me how the matrix is shaped in principle, so that I can apply zeroRowsLocal to the correct dofs? I have a FunctionSpace like this:

el = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
El = ufl.MixedElement([el,el])
V = dolfinx.fem.FunctionSpace(mesh, El)
V0 ,_ = V.sub(0).collapse()
V1,_ = V.sub(1).collapse()

Am I right, that I need the dofs of the respective collapsed subspace or is it the non collapsed one?
Assuming a boundary condition like this:

# inactive_dofs_V0,  inactive_dofs_V1 are np.arrays with the respective dofs in the collapsed space
zero = dolfinx.fem.Constant(mesh,0.0)
bcs = [dolfinx.fem.dirichletbc(zero,inactive_dofs_V0,V.sub(0)) , dolfinx.fem.dirichletbc(zero,inactive_dofs_V1,V.sub(1))]

what would the identity row approach look like? I guess since the matrix has a different shape, I can’t just do this exactly like in the tutorial:

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

Do I need an offset for the dofs for each subspace?

And since the subspaces are already in the DirichletBC object, would this work as is?

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)

Again thank you very much for your help, I appreciate it a lot!

Best,
fem_fatale

[dokken]

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)