# How to assign boundary conditions based on values of dolfinx.Function

Hello everybody,

I would like to set boundary conditions based on values of dolfinx.Function. For instance, let’s define the Function in the following manner (as suggested by Dokken):

``````import dolfinx
from mpi4py import MPI
from dolfinx.io import XDMFFile

mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 9, 9)
V = dolfinx.FunctionSpace(mesh, ("DG", 0))
v = dolfinx.Function(V)
x = V.tabulate_dof_coordinates()
for i in range(x.shape):
midpoint = x[i,:]
if midpoint> 0.5 and midpoint>0.25:
v.vector.setValueLocal(i, 2)
else:
v.vector.setValueLocal(i, 1)

xdmf = XDMFFile(MPI.COMM_WORLD, "test.xdmf", "w")
xdmf.write_mesh(mesh)
xdmf.write_function(v)
``````

so the dolfinx.Function v has different values at different locations. Now I want to set e.g. Neumann bc on facets where v equals 2 and Robin bc where v equals 1 (visualized in figure below).

I could use marker functions as suggested throughout the Dokken’s tutorial, however, in my problem the distribution of values in dolfinx.Function is more complex and changes with every time step, so it is not feasible to define locations for Robin and Neumann geometrically (numpy.isclose and so on).

I have adapted the code suggested by Dokken so that the last defined boundary conditions (in “boundaries”, see the code below) rewrite previously defined ones. Thus I can e.g. first set Robin bc on the whole boundary and then substitute Robin with Neumann where for instance x ==1 (as shown in the code below), but I want to substitute where dolfinx.Function has values of two. Unfortunately, I don’t know how I can use values of dolfinx.Function instead of marker function.

``````import numpy as np

boundaries = [(1, lambda x: np.full(x.shape, True, dtype=bool)),  # Set robin on the whole boundary
(2, lambda x: np.isclose(x, 1))]   # Now Robin will be rewritten by Neumann  for x == 1
# I want the same but for dolfinx.Function values == 2

facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in reversed(boundaries):
facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full(len(facets), marker))
facet_indices, pos = np.unique(np.array(np.hstack(facet_indices), dtype=np.int32), return_index=True)
facet_markers = np.array(np.hstack(facet_markers)[pos], dtype=np.int32)
facet_tag = dolfinx.MeshTags(mesh, fdim, facet_indices, facet_markers)

import dolfinx.io
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "facet_tags.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_meshtags(facet_tag)
``````

Any help would be really appreciated. I think there could be a better way to achieve what I want, so it would nice if someone could point out a correct direction.

Thanks

How about using a UFL conditional (Form language — Unified Form Language (UFL) 2021.1.0 documentation)
i.e.

``````v = Function(V)
# Assign values 1 and 2 to v
# ...
# Add boundary condition terms to variational form
F += ufl.conditional(ufl.eq(v, 1), 1, 0) *  robin integrand * ds
F += ufl.conditional(ufl.eq(v, 2), 1, 0) * neumann_integrand * ds
``````
1 Like

Hi Jorgen,

yes, it works. Thank you.

Accidentally I notice one weird thing. In the following code, the solution should be the same independent on the value of “value_in_question”. And it is the same except for value_in_question = 3. In this case the right boundary is not recognized (shown in figure below). Why “3” is so special in this case?

``````import dolfinx
import ufl
from mpi4py import MPI
from dolfinx.io import XDMFFile
import numpy as np
import math

value_in_question = 3

mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 10, 10)
V = dolfinx.FunctionSpace(mesh, ("CG", 1))

bc_func = dolfinx.Function(V)
x = V.tabulate_dof_coordinates()
# Define boundary control function (bc_func)
for i in range(x.shape):
if x[i,0] == 0:
bc_func.vector.setValueLocal(i, 2)
elif math.isclose(x[i,0], 1):
bc_func.vector.setValueLocal(i, value_in_question)
else:
bc_func.vector.setValueLocal(i, 1)
# print(bc_func.vector.array)

u_n = dolfinx.Function(V)
uh = dolfinx.Function(V)
v = ufl.TestFunction(V)
f = dolfinx.Constant(mesh, 0)

xdmf = XDMFFile(MPI.COMM_WORLD, "function.xdmf", "w")
xdmf.write_mesh(mesh)
xdmf.write_function(bc_func)
xdmf.write_function(uh, 0)

F = 0.54 * 0.0078 * uh * v * ufl.dx + 0.042 * ufl.dot(ufl.grad(uh), ufl.grad(v)) * ufl.dx - (0.54 * 0.0078 * u_n + f) * v * ufl.dx

F += ufl.conditional(ufl.eq(bc_func, 1), 1, 0) *  0.4 * ufl.inner(uh-0, v) * ufl.ds
F += ufl.conditional(ufl.eq(bc_func, 2), 1, 0) * ufl.inner(-5, v) * ufl.ds
F += ufl.conditional(ufl.eq(bc_func, value_in_question), 1, 0) * ufl.inner(-10, v) * ufl.ds

problem = dolfinx.fem.NonlinearProblem(F, uh)
solver = dolfinx.NewtonSolver(MPI.COMM_WORLD, problem)
n, converged = solver.solve(uh)
assert(converged)
print(f"Number of interations: {n:d}")

xdmf.write_function(uh, 1)

``````

This has to do with machine precision in assembly. I would suggest replacing `ufl.eq` with the following:

``````def cond(value, tol=1e-13):
return ufl.conditional(ufl.lt(abs(bc_func-value), tol), 1, 0)
``````

This is a more stable way of checking if it is equal or not.
The full modified code can be found below.

``````from dolfinx.fem import assemble_scalar
from dolfinx.mesh import create_unit_square
from dolfinx.fem import FunctionSpace, Function, DirichletBC, Constant, NonlinearProblem
from dolfinx.nls import NewtonSolver
import ufl
from mpi4py import MPI
from dolfinx.io import XDMFFile
import numpy as np
import math

value_in_question = 3

mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
V = FunctionSpace(mesh, ("CG", 1))

bc_func = Function(V)
x = V.tabulate_dof_coordinates()
# Define boundary control function (bc_func)
for i in range(x.shape):
if x[i, 0] == 0:
bc_func.vector.setValueLocal(i, 2)
elif math.isclose(x[i, 0], 1):
bc_func.vector.setValueLocal(i, value_in_question)
else:
bc_func.vector.setValueLocal(i, 1)

def cond(value, tol=1e-13):
return ufl.conditional(ufl.lt(abs(bc_func-value), tol), 1, 0)

u_n = Function(V)
uh = Function(V)
v = ufl.TestFunction(V)
f = Constant(mesh, 0)

xdmf = XDMFFile(MPI.COMM_WORLD, "function.xdmf", "w")
xdmf.write_mesh(mesh)
xdmf.write_function(bc_func)
xdmf.write_function(uh, 0)

F = 0.54 * 0.0078 * uh * v * ufl.dx + 0.042 * \
(0.54 * 0.0078 * u_n + f) * v * ufl.dx

F += cond(1) * \
0.4 * ufl.inner(uh-0, v) * ufl.ds
F += cond(2) * ufl.inner(-5, v) * ufl.ds
F += cond(value_in_question) * ufl.inner(-10, v) * ufl.ds
print()
problem = NonlinearProblem(F, uh)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
n, converged = solver.solve(uh)
assert(converged)
print(f"Number of interations: {n:d}")

xdmf.write_function(uh, 1)
``````
1 Like