# Get boundary points from custom cell tags

Hi everyone!

I’m using dolfinx 0.5.2 and I need to find the boundary points of a subdomain defined by custom cell tags (the example code is reported below). Specifically, the cells of the subdomain of interest are marked with 1 (cells of the inner rectangle for this case), all the other cells with 0. Thanks in advance for the help!

``````from shapely import Polygon, Point
import numpy as np
import gmsh
from mpi4py import MPI
from dolfinx.io.gmshio import model_to_mesh
import dolfinx
import dolfinx.plot as plot
import pyvista

rank = MPI.COMM_WORLD.rank

rect1 = np.array([[1, 1], [4, 1], [4, 4], [1, 4], [1, 1]], dtype=np.float64)
rect2 = np.array([[0, 0], [10, 0], [10, 5], [0, 5], [0, 0]], dtype=np.float64)

gmsh.initialize()

rect1_tags = []
rect2_tags = []
for point1, point2 in zip(rect1, rect2):
p1 = gmsh.model.occ.addPoint(point1[0], point1[1], 0, 1)
p2 = gmsh.model.occ.addPoint(point2[0], point2[1], 0, 1)
rect1_tags.append(p1)
rect2_tags.append(p2)
gmsh.model.occ.synchronize()

rect1_wire = []
rect2_wire = []
for ii in range(len(rect1_tags)):
if ii > 0:
l1 = gmsh.model.occ.addLine(rect1_tags[ii - 1], rect1_tags[ii])
l2 = gmsh.model.occ.addLine(rect2_tags[ii - 1], rect2_tags[ii])
rect1_wire.append(l1)
rect2_wire.append(l2)
gmsh.model.occ.synchronize()

gmsh.model.occ.synchronize()
gmsh.model.occ.synchronize()
gmsh.model.occ.synchronize()

# generate 2D mesh
gmsh.option.setNumber("Mesh.Algorithm", 2)
gmsh.model.mesh.generate(2)
gmsh.model.mesh.optimize("Netgen")
gmsh.model.occ.synchronize()
gmsh.write("mesh.geo_unrolled")

model_rank = 0
mesh_comm = MPI.COMM_WORLD
mesh, ct, ft = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
gmsh.finalize()

def create_meshtags(
mesh: dolfinx.mesh.Mesh, poly: Polygon
) -> dolfinx.cpp.mesh.MeshTags_int32:
cells = mesh.topology.connectivity(mesh.topology.dim, 0)

id_cells_inpoly = []
for ii in range(len(cells)):
centroid = triangle.centroid
if poly.contains(Point(centroid)):
id_cells_inpoly.append(ii)

values = np.zeros([len(cells), 1], dtype=np.int32)
values[id_cells_inpoly] = 1
meshtags = dolfinx.mesh.meshtags_from_entities(
mesh, mesh.topology.dim, cells, values
)
return meshtags

def plot_meshtags(
mesh: dolfinx.mesh.Mesh,
meshtags: dolfinx.cpp.mesh.MeshTags_int32,
filename: str = "meshtags.png",
):
# Create VTK mesh
cells, types, x = plot.create_vtk_mesh(mesh)
grid = pyvista.UnstructuredGrid(cells, types, x)

# Attach the cells tag data to the pyvita grid
grid.cell_data["Marker"] = meshtags.values
grid.set_active_scalars("Marker")

plotter = pyvista.Plotter()
plotter.view_xy()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
plotter.screenshot(filename)

# plot original cell tags
plot_meshtags(mesh, ct)

# create custom cell tags
ct = create_meshtags(mesh, Polygon(rect1))
# plot custom cell tags
plot_meshtags(mesh, ct)
``````

When you say boundary points, do you mean

1. the points on the interface between subdomain 0 and subdomain 1?
2. any vertex/facet on the external boundary of your mesh (but only on cells marked with 1)? An external boundary is defined as a facet which is only connected to a single cell.

In your example all exterior facets are connected to cells with the marker 1, so just getting the exterior facets would do the trick in your example.

Sorry, I’m not expert about it so I’ll try to explain my self better. I need the two vertexes for each cell marked with 1 (in yellow) which are part of the boundary for subdomain 1 (In order to define the boundary, for example, of a more complex domain) → (the red points in the figure)

1 Like

Consider the following:

``````from shapely import Polygon, Point
import numpy as np
import gmsh
from mpi4py import MPI
from dolfinx.io.gmshio import model_to_mesh
import dolfinx
import dolfinx.plot as plot
import pyvista
pyvista.set_jupyter_backend("static")
pyvista.start_xvfb()
rank = MPI.COMM_WORLD.rank

rect1 = np.array([[1, 1], [4, 1], [4, 4], [1, 4], [1, 1]], dtype=np.float64)
rect2 = np.array([[0, 0], [10, 0], [10, 5], [0, 5], [0, 0]], dtype=np.float64)

gmsh.initialize()

rect1_tags = []
rect2_tags = []
for point1, point2 in zip(rect1, rect2):
p1 = gmsh.model.occ.addPoint(point1[0], point1[1], 0, 1)
p2 = gmsh.model.occ.addPoint(point2[0], point2[1], 0, 1)
rect1_tags.append(p1)
rect2_tags.append(p2)
gmsh.model.occ.synchronize()

rect1_wire = []
rect2_wire = []
for ii in range(len(rect1_tags)):
if ii > 0:
l1 = gmsh.model.occ.addLine(rect1_tags[ii - 1], rect1_tags[ii])
l2 = gmsh.model.occ.addLine(rect2_tags[ii - 1], rect2_tags[ii])
rect1_wire.append(l1)
rect2_wire.append(l2)
gmsh.model.occ.synchronize()

gmsh.model.occ.synchronize()
gmsh.model.occ.synchronize()
gmsh.model.occ.synchronize()

# generate 2D mesh
gmsh.option.setNumber("Mesh.Algorithm", 2)
gmsh.model.mesh.generate(2)
gmsh.model.mesh.optimize("Netgen")
gmsh.model.occ.synchronize()
gmsh.write("mesh.geo_unrolled")

model_rank = 0
mesh_comm = MPI.COMM_WORLD
cell_partitioner = dolfinx.mesh.create_cell_partitioner(dolfinx.mesh.GhostMode.shared_facet)
mesh, ct, ft = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2, partitioner=cell_partitioner)
gmsh.finalize()

def create_meshtags(
mesh: dolfinx.mesh.Mesh, poly: Polygon
) -> dolfinx.cpp.mesh.MeshTags_int32:
cells = mesh.topology.connectivity(mesh.topology.dim, 0)

id_cells_inpoly = []
for ii in range(len(cells)):
centroid = triangle.centroid
if poly.contains(Point(centroid)):
id_cells_inpoly.append(ii)

values = np.zeros([len(cells), 1], dtype=np.int32)
values[id_cells_inpoly] = 1
meshtags = dolfinx.mesh.meshtags_from_entities(
mesh, mesh.topology.dim, cells, values
)
return meshtags

def plot_meshtags(
mesh: dolfinx.mesh.Mesh,
meshtags: dolfinx.cpp.mesh.MeshTags_int32,
filename: str = "meshtags.png",
):
# Create VTK mesh
cells, types, x = plot.create_vtk_mesh(mesh)
grid = pyvista.UnstructuredGrid(cells, types, x)

# Attach the cells tag data to the pyvita grid
grid.cell_data["Marker"] = meshtags.values
grid.set_active_scalars("Marker")

plotter = pyvista.Plotter()
plotter.view_xy()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
plotter.screenshot(filename)

# plot original cell tags
plot_meshtags(mesh, ct)

# create custom cell tags
ct = create_meshtags(mesh, Polygon(rect1))
# plot custom cell tags
plot_meshtags(mesh, ct)

cell_map = mesh.topology.index_map(mesh.topology.dim)
vertex_map = mesh.topology.index_map(0)
mesh.topology.create_entities(mesh.topology.dim-1)
facet_map = mesh.topology.index_map(mesh.topology.dim-1)
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
f_to_c =  mesh.topology.connectivity(mesh.topology.dim-1, mesh.topology.dim)
f_to_v =  mesh.topology.connectivity(mesh.topology.dim-1, 0)

assert len(ct.indices) == cell_map.size_local + cell_map.num_ghosts
exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
vertex_marker = np.zeros(vertex_map.size_local + vertex_map.num_ghosts, dtype=np.int32)
for facet in range(facet_map.size_local+ facet_map.num_ghosts):
cell_markers = ct.values[cells]
# If exterior facet and connected cell is marked with 1
if 1 in cell_markers and facet in exterior_facets:
vertex_marker[vertices] = 1
# If on interface between the two domains
if 1 in cell_markers and 0 in cell_markers:
vertex_marker[vertices] = 1

vt = dolfinx.mesh.meshtags(mesh, 0, np.arange(len(vertex_marker), dtype=np.int32), vertex_marker)
with dolfinx.io.XDMFFile(mesh.comm, "vertices.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
mesh.topology.create_connectivity(0, mesh.topology.dim)
xdmf.write_meshtags(vt)
``````

generating

Thus you could gather all the vertices from ` vertices = f_to_v.links(facet)` whenever it is called into an array and do whatever you plan with it.

1 Like

Thank you so much, that’s exactly what I need! You’re great @dokken !

If I may, I would like to link back to this issue to ask a further question. In this case, the vertices obtained are not in order (by plotting, you do not in fact form a closed line but a series of intersecting lines). I would like to reorder the vertices found in such a way that I can plot the polyline that describes the boundary of subdomain 1, is it possible to do this by exploiting the dolfinx functions?

You can do this using `f_to_v` and the equivalent `v_to_f` and checking your now existing marker to generate a closed loop.

Ok thanks a lot for your quick answer. I apologize in advance for the silly request, however, I do not quite understand how to apply it to my algorithm. Would it be possible to have an example?