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
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)
