Setting scalar function for BDM boundary condition with controlled flux [mixed Poisson]

I’m interested in systems with controlled flux at the boundary. As I understand the divergence-conforming elements such as BDM (or Raviart-Thomas) are a natural choice for this. Following the example of the mixed Poisson demo, the boundary condition for vector-valued function \sigma is

\sigma \cdot n = g \; \mathrm{on\ } \Gamma_N

so the flux (normal component) at the boundary is determined by scalar function g.

What I want is a symmetric solution with the same scalar function g applied at top and bottom.
But in the mixed Poisson demo the boundary condition is implemented as a vector-valued Dirichlet BC, rather than the scalar function g directly. To try to understand how it is applied, I modified the demo to set a constant value 1 at the top and bottom, via

    values[1, :] = 1

for both f1 and f2, equivalent to applying the same f1 to both the top and bottom boundaries. This does not obtain a symmetric solution. A cross section of u in the y-direction is


The vector-valued Dirichlet BC is applied literally, giving a vector value +1 in the y direction both at top and bottom.

What I want is the same flux “into” the domain, so the vector value at the bottom should be opposite (negative) to the top. That is, I want to control the boundary with scalar function g only, with direction determined by the normal vector.

In the case of the mixed Poisson demo, I can obtain the desired symmetric solution by manually setting the value to -1 in f2 (modified code below).

But that does not extend well to the case of applying a simple scalar boundary condition (for the normal component) over a complex boundary?

How should a scalar function g be applied as a BDM boundary condition?

Modified mixed Poisson demo:

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np

from basix.ufl import element, mixed_element
from dolfinx import fem, io, mesh
from dolfinx.fem import Function
from dolfinx.fem.petsc import LinearProblem
from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions,
                 div, exp, inner)

domain = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32, mesh.CellType.quadrilateral)

k = 1
Q_el = element("BDMCF", domain.basix_cell(), k)
P_el = element("DG", domain.basix_cell(), k - 1)
V_el = mixed_element([Q_el, P_el])
V = fem.functionspace(domain, V_el)

(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)

x = SpatialCoordinate(domain)
f = 10.0 * exp(-((x[0] - 0.5) * (x[0] - 0.5) + (x[1] - 0.5) * (x[1] - 0.5)) / 0.02)

dx = Measure("dx", domain)
a = inner(sigma, tau) * dx + inner(u, div(tau)) * dx + inner(div(sigma), v) * dx
L = -inner(f, v) * dx


fdim = domain.topology.dim - 1
facets_top = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[1], 1.0))
Q, _ = V.sub(0).collapse()
dofs_top = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_top)


def f1(x):
    values = np.zeros((2, x.shape[1]))
    values[1, :] = 1
    return values


f_h1 = fem.Function(Q)
f_h1.interpolate(f1)
bc_top = fem.dirichletbc(f_h1, dofs_top, V.sub(0))


facets_bottom = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[1], 0.0))
dofs_bottom = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_bottom)


def f2(x):
    values = np.zeros((2, x.shape[1]))

    # setting the same value as f2
    # gives an asymmetric solution, not desired
    #values[1, :] = 1

    # gives the desired symmetric solution
    # but can this be done using a simple scalar value
    # of the normal component in any direction?
    values[1, :] = -1
    return values


f_h2 = fem.Function(Q)
f_h2.interpolate(f2)
bc_bottom = fem.dirichletbc(f_h2, dofs_bottom, V.sub(0))


bcs = [bc_top, bc_bottom]

problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu",
                                                      "pc_factor_mat_solver_type": "mumps"})
try:
    w_h = problem.solve()
except PETSc.Error as e:  # type: ignore
    if e.ierr == 92:
        print("The required PETSc solver/preconditioner is not available. Exiting.")
        print(e)
        exit(0)
    else:
        raise e

sigma_h, u_h = w_h.split()

with io.XDMFFile(domain.comm, "out_mixed_poisson/u.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_function(u_h)

Vl = fem.functionspace(domain, ("Lagrange", k, (domain.geometry.dim,)))
sigmal_h = Function(Vl)
sigmal_h.interpolate(sigma_h)
with io.VTXWriter(domain.comm, "out_mixed_poisson/sigmal.bp", sigmal_h) as file:
    file.write(0)

Since the dofs are associated with an edge, I would probably project the Facet normal into the space, using code similar to: dolfinx_mpc/python/dolfinx_mpc/utils/mpc_utils.py at afa9fd51eb7deb8098469fcb3b9185c4ce3099c0 · jorgensd/dolfinx_mpc · GitHub

Here is a version (i bashed it together quickly, so there could be typos):

from mpi4py import MPI
from petsc4py import PETSc
from basix.ufl import element, mixed_element
from dolfinx import fem, io, mesh, fem, cpp
from dolfinx.fem import Function
from dolfinx.fem.petsc import LinearProblem
from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions,
                 div, exp, inner)
import ufl
import numpy as np


def facet_normal_approximation(V):
    """
    Approximate the facet normal by projecting it into the function space for a set of facets

    Args:
        V: The function space to project into
    """
    n = ufl.FacetNormal(V.mesh)
    nh = fem.Function(V)
    u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
    ds = ufl.ds(domain=V.mesh)
    a = (ufl.inner(u, v) * ds)
    L = ufl.inner(n, v) * ds

    # Find all dofs that are not boundary dofs
    imap = V.dofmap.index_map
    all_blocks = np.arange(imap.size_local, dtype=np.int32)
    boundary_facets = mesh.exterior_facet_indices(V.mesh.topology)
    top_blocks = fem.locate_dofs_topological(
        V, V.mesh.topology.dim - 1, boundary_facets)
    deac_blocks = all_blocks[np.isin(all_blocks, top_blocks, invert=True)]

    # Note there should be a better way to do this
    # Create sparsity pattern only for constraint + bc
    bilinear_form = fem.form(a)
    pattern = fem.create_sparsity_pattern(bilinear_form)
    pattern.insert_diagonal(deac_blocks)
    pattern.finalize()
    u_0 = fem.Function(V)
    u_0.vector.set(0)

    bc_deac = fem.dirichletbc(u_0, deac_blocks)
    A = cpp.la.petsc.create_matrix(V.mesh.comm, pattern)
    A.zeroEntries()

    # Assemble the matrix with all entries
    form_coeffs = cpp.fem.pack_coefficients(bilinear_form._cpp_object)
    form_consts = cpp.fem.pack_constants(bilinear_form._cpp_object)
    cpp.fem.petsc.assemble_matrix(A, bilinear_form._cpp_object,
                                  form_consts, form_coeffs, [bc_deac._cpp_object])
    if bilinear_form.function_spaces[0] is bilinear_form.function_spaces[1]:
        A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH)  # type: ignore
        A.assemblyEnd(PETSc.Mat.AssemblyType.FLUSH)  # type: ignore
        cpp.fem.petsc.insert_diagonal(A, bilinear_form.function_spaces[0], [
            bc_deac._cpp_object], 1.0)
    A.assemble()
    linear_form = fem.form(L)
    b = fem.petsc.assemble_vector(linear_form)

    fem.petsc.apply_lifting(b, [bilinear_form], [[bc_deac]])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES,
                  mode=PETSc.ScatterMode.REVERSE)  # type: ignore
    fem.petsc.set_bc(b, [bc_deac])

    # Solve Linear problem
    solver = PETSc.KSP().create(V.mesh.comm)  # type: ignore
    solver.setType("cg")
    solver.rtol = 1e-8
    solver.setOperators(A)
    solver.solve(b, nh.vector)
    nh.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT,
                          mode=PETSc.ScatterMode.FORWARD)  # type: ignore
    return nh


domain = mesh.create_unit_square(
    MPI.COMM_WORLD, 32, 32, mesh.CellType.quadrilateral)

k = 1
Q_el = element("BDMCF", domain.basix_cell(), k)
P_el = element("DG", domain.basix_cell(), k - 1)
V_el = mixed_element([Q_el, P_el])
V = fem.functionspace(domain, V_el)

(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)

x = SpatialCoordinate(domain)
f = 10.0 * exp(-((x[0] - 0.5) * (x[0] - 0.5) +
               (x[1] - 0.5) * (x[1] - 0.5)) / 0.02)

dx = Measure("dx", domain)
a = inner(sigma, tau) * dx + inner(u, div(tau)) * \
    dx + inner(div(sigma), v) * dx
L = -inner(f, v) * dx


fdim = domain.topology.dim - 1
facets_top = mesh.locate_entities_boundary(
    domain, fdim, lambda x: np.isclose(x[1], 1.0))
Q, _ = V.sub(0).collapse()
dofs_top = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_top)

nh = facet_normal_approximation(Q)
bc_top = fem.dirichletbc(nh, dofs_top, V.sub(0))


facets_bottom = mesh.locate_entities_boundary(
    domain, fdim, lambda x: np.isclose(x[1], 0.0))
dofs_bottom = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_bottom)
bc_bottom = fem.dirichletbc(nh, dofs_bottom, V.sub(0))


bcs = [bc_top, bc_bottom]

problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu",
                                                      "pc_factor_mat_solver_type": "mumps"})
try:
    w_h = problem.solve()
except PETSc.Error as e:  # type: ignore
    if e.ierr == 92:
        print("The required PETSc solver/preconditioner is not available. Exiting.")
        print(e)
        exit(0)
    else:
        raise e

sigma_h, u_h = w_h.split()

with io.XDMFFile(domain.comm, "out_mixed_poisson/u.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_function(u_h)

Vl = fem.functionspace(domain, ("DG", k+1, (domain.geometry.dim,)))
sigmal_h = Function(Vl)
sigmal_h.interpolate(sigma_h)
with io.VTXWriter(domain.comm, "out_mixed_poisson/sigmal.bp", sigmal_h, engine="BP4") as file:
    file.write(0)

giving


Please let me know if this helps.

The issue atm is the facet normal approximated at the corner. that would need some better treatment.

@dparsons I improved the facet representations a bit, but it doesn’t seem to be sufficient:

from mpi4py import MPI
from petsc4py import PETSc
from basix.ufl import element, mixed_element
from dolfinx import fem, io, mesh, fem, cpp
from dolfinx.fem import Function
from dolfinx.fem.petsc import LinearProblem
from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions,
                 div, exp, inner)
import ufl
import numpy as np


def facet_normal_approximation(V, mt: mesh.MeshTags, mt_id: int):
    """
    Approximate the facet normal by projecting it into the function space for a set of facets

    Args:
        V: The function space to project into
    """
    n = ufl.FacetNormal(V.mesh)
    nh = fem.Function(V)
    u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
    ds = ufl.ds(domain=V.mesh, subdomain_data=mt, subdomain_id=mt_id)
    a = ufl.inner(u, v) * ds
    L = ufl.inner(n, v) * ds

    # Find all dofs that are not boundary dofs
    imap = V.dofmap.index_map
    all_blocks = np.arange(imap.size_local, dtype=np.int32)
    top_blocks = fem.locate_dofs_topological(
        V, V.mesh.topology.dim - 1, mt.find(mt_id))
    deac_blocks = all_blocks[np.isin(all_blocks, top_blocks, invert=True)]

    # Note there should be a better way to do this
    # Create sparsity pattern only for constraint + bc
    bilinear_form = fem.form(a)
    pattern = fem.create_sparsity_pattern(bilinear_form)
    pattern.insert_diagonal(deac_blocks)
    pattern.finalize()
    u_0 = fem.Function(V)
    u_0.vector.set(0)

    bc_deac = fem.dirichletbc(u_0, deac_blocks)
    A = cpp.la.petsc.create_matrix(V.mesh.comm, pattern)
    A.zeroEntries()

    # Assemble the matrix with all entries
    form_coeffs = cpp.fem.pack_coefficients(bilinear_form._cpp_object)
    form_consts = cpp.fem.pack_constants(bilinear_form._cpp_object)
    cpp.fem.petsc.assemble_matrix(A, bilinear_form._cpp_object,
                                  form_consts, form_coeffs, [bc_deac._cpp_object])
    if bilinear_form.function_spaces[0] is bilinear_form.function_spaces[1]:
        A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH)  # type: ignore
        A.assemblyEnd(PETSc.Mat.AssemblyType.FLUSH)  # type: ignore
        cpp.fem.petsc.insert_diagonal(A, bilinear_form.function_spaces[0], [
            bc_deac._cpp_object], 1.0)
    A.assemble()

    linear_form = fem.form(L)
    b = fem.petsc.assemble_vector(linear_form)
    fem.petsc.apply_lifting(b, [bilinear_form], [[bc_deac]])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES,
                  mode=PETSc.ScatterMode.REVERSE)  # type: ignore
    fem.petsc.set_bc(b, [bc_deac])
    b.ghostUpdate(addv=PETSc.InsertMode.INSERT,
                  mode=PETSc.ScatterMode.FORWARD)  # type: ignore
    # Solve Linear problem
    solver = PETSc.KSP().create(V.mesh.comm)  # type: ignore
    solver.setType("cg")
    solver.rtol = 1e-8
    solver.setOperators(A)
    solver.solve(b, nh.vector)
    nh.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT,
                          mode=PETSc.ScatterMode.FORWARD)  # type: ignore
    return nh


domain = mesh.create_unit_square(
    MPI.COMM_WORLD, 32, 32, mesh.CellType.quadrilateral)

k = 1
Q_el = element("BDMCF", domain.basix_cell(), k)
P_el = element("DG", domain.basix_cell(), k - 1)
V_el = mixed_element([Q_el, P_el])
V = fem.functionspace(domain, V_el)

(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)

x = SpatialCoordinate(domain)
f = 10.0 * exp(-((x[0] - 0.5) * (x[0] - 0.5) +
               (x[1] - 0.5) * (x[1] - 0.5)) / 0.02)

dx = Measure("dx", domain)
a = inner(sigma, tau) * dx + inner(u, div(tau)) * \
    dx + inner(div(sigma), v) * dx
L = -inner(f, v) * dx


fdim = domain.topology.dim - 1
facets_top = mesh.locate_entities_boundary(
    domain, fdim, lambda x: np.isclose(x[1], 1.0))
facets_bottom = mesh.locate_entities_boundary(
    domain, fdim, lambda x: np.isclose(x[1], 0.0))
vals_top = np.full_like(facets_top, 3, dtype=np.int32)
vals_bottom = np.full_like(facets_bottom, 7, dtype=np.int32)
all_facets = np.hstack([facets_top, facets_bottom])
sort = np.argsort(all_facets)
all_values = np.hstack([vals_top, vals_bottom])[sort]
mt = mesh.meshtags(domain, fdim, all_facets[sort], all_values)

Q, _ = V.sub(0).collapse()
dofs_top = fem.locate_dofs_topological((V.sub(0), Q), fdim, mt.find(3))
nh_top = facet_normal_approximation(Q, mt, 3)
bc_top = fem.dirichletbc(nh_top, dofs_top, V.sub(0))


dofs_bottom = fem.locate_dofs_topological((V.sub(0), Q), fdim, mt.find(7))
nh_bot = facet_normal_approximation(Q, mt, 7)
bc_bottom = fem.dirichletbc(nh_bot, dofs_bottom, V.sub(0))


bcs = [bc_top, bc_bottom]

problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu",
                                                      "pc_factor_mat_solver_type": "mumps"})
try:
    w_h = problem.solve()
except PETSc.Error as e:  # type: ignore
    if e.ierr == 92:
        print("The required PETSc solver/preconditioner is not available. Exiting.")
        print(e)
        exit(0)
    else:
        raise e

sigma_h, u_h = w_h.split()

with io.XDMFFile(domain.comm, "out_mixed_poisson/u.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_function(u_h)

Vl = fem.functionspace(domain, ("DG", k+1, (domain.geometry.dim,)))
sigmal_h = Function(Vl)
sigmal_h.interpolate(sigma_h)

n_out = Function(Vl)
n_out.name = "nh_top"
n_out.interpolate(nh_top)
n_b = Function(Vl)
n_b.name = "nh_bottom"
n_b.interpolate(nh_bot)

with io.VTXWriter(domain.comm, "out_mixed_poisson/sigmal.bp", [sigmal_h, n_out, n_b], engine="BP4") as file:
    file.write(0)

I guess a more robust approach would be to push the facet normal forward to the physical domain with Piola map, as shown in: Setting boundary conditions on H(div) spaces defined on general meshes - #2 by dokken

Let’s see if I have further time to look at this soon.

Thanks for that. The idea makes sense, generating the complex boundary vector function which is to be used as the Direchlet BC for the BDM vector function.

I gather an alternative approach could be to impose the condition in weak form via Nitsche’s method. Would be interesting to see how the two approaches compare.

I apologise for being a bit obtuse, but I have to ask, where is the actual scalar function here? There is the function u_0, but it is set to 0. There is the constant 1.0 in

cpp.fem.petsc.insert_diagonal(A, bilinear_form.function_spaces[0], [
            bc_deac._cpp_object], 1.0)

but changing it does not change the value of the vector solution.

What the code above does is that it sets the boundary condition sigma = n where n is the facet normal (in physical space).

The Dirichlet bc in that function just inserts diagonal entries for all dofs not on a given boundary.

The justification behind my thoughts here is by looking at the basis functions: DefElement

I see it now. Multiplying that L by the value I need gives the control I’m looking for.