Get boundary points from custom cell tags

mnotazio · September 13, 2023, 7:59am

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.addCurveLoop(rect1_wire, 1)
gmsh.model.occ.addCurveLoop(rect2_wire, 2)
gmsh.model.occ.synchronize()
gmsh.model.occ.addPlaneSurface([1], 1)
gmsh.model.occ.addPlaneSurface([1, 2], 2)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(2, [1], 1)
gmsh.model.addPhysicalGroup(2, [2], 2)
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)):
        triangle = Polygon(mesh.geometry.x[cells.links(ii)])
        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.add_mesh(grid, show_edges=True, show_scalar_bar=False)
    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)

dokken · September 13, 2023, 1:14pm

When you say boundary points, do you mean

the points on the interface between subdomain 0 and subdomain 1?
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.

mnotazio · September 13, 2023, 1:27pm

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)

dokken · September 13, 2023, 2:08pm

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.addCurveLoop(rect1_wire, 1)
gmsh.model.occ.addCurveLoop(rect2_wire, 2)
gmsh.model.occ.synchronize()
gmsh.model.occ.addPlaneSurface([1], 1)
gmsh.model.occ.addPlaneSurface([1, 2], 2)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(2, [1], 1)
gmsh.model.addPhysicalGroup(2, [2], 2)
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)):
        triangle = Polygon(mesh.geometry.x[cells.links(ii)])
        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.add_mesh(grid, show_edges=True, show_scalar_bar=False)
    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):
    cells = f_to_c.links(facet)
    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:
        vertices = f_to_v.links(facet)
        vertex_marker[vertices] = 1
    # If on interface between the two domains
    if 1 in cell_markers and 0 in cell_markers:
        vertices = f_to_v.links(facet)
        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.

mnotazio · September 13, 2023, 2:13pm

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

mnotazio · September 18, 2023, 1:55pm

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?

dokken · September 18, 2023, 2:43pm

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.

mnotazio · September 18, 2023, 2:50pm

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?

Topic		Replies	Views
Find coordinates of cells in a subdomain General	3	135	January 11, 2024
Tag entities between two subdomains mesh	3	353	March 21, 2023
What is the more succint way to get the "cell" and "facet" tags? mesh	5	232	October 27, 2023
Using facet tags to define boundary conditions - GMSH - FEniCSx - dolfinx	4	1381	October 3, 2022
Defining subdomains from gmsh on CG elements dolfinx	1	718	June 27, 2022

Get boundary points from custom cell tags

Related Topics