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[0]):
midpoint = x[i,:]
if midpoint[0]> 0.5 and midpoint[1]>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[1], True, dtype=bool)), # Set robin on the whole boundary
(2, lambda x: np.isclose(x[0], 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
