Hi all!
With reference to the post from some time ago (Get boundary points from custom cell tags), I created custom cell tags and identified the facets that are at the interface between the domain with cell tag 0 and domain with cell tag 1. I’m using dolfinx 0.8.0.
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 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, _, _ = 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):
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), dtype=np.int32)
values[id_cells_inpoly] = 1
meshtags = dolfinx.mesh.meshtags_from_entities(
mesh, mesh.topology.dim, cells, values
)
return meshtags
# create custom cell tags
ct = create_meshtags(mesh, Polygon(rect1))
def find_domain_boundary(mesh, ct, domains):
tdim = mesh.topology.dim
fdim = tdim - 1
mesh.topology.create_connectivity(fdim, tdim)
f_to_c = mesh.topology.connectivity(fdim, tdim)
c_to_f = mesh.topology.connectivity(tdim, fdim)
f_to_v = mesh.topology.connectivity(fdim, fdim - 1)
mesh.topology.create_connectivity(0, fdim)
v_to_f = mesh.topology.connectivity(0, fdim)
# create facet tags for the exterior facet of the plasma subdomain
num_facets = (
mesh.topology.index_map(fdim).size_local
+ mesh.topology.index_map(fdim).num_ghosts
)
exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
all_boundary_pts = []
all_facets = []
for value in domains:
boundary_ordered_vertexes = []
facets = np.array([], dtype=np.int32)
if value in np.unique(ct.values):
# all cells related to value
c_i = np.where(ct.values == value)
# all facets of c_i
f_i = np.unique(c_to_f.array.reshape((-1, 3))[c_i].flatten())
# add exterior facets
facets = np.append(facets, np.intersect1d(exterior_facets, f_i))
temp = np.array([len(np.unique(ct.values[f_to_c.links(i)])) for i in f_i])
facets = np.append(facets, f_i[np.where(temp > 1)[0]])
vertexes = f_to_v.links(facets[0])
boundary_ordered_vertexes.append(vertexes[0])
ind = 0
for i in range(len(facets) - 1):
ind_v = np.where(~np.isin(vertexes, boundary_ordered_vertexes))[0]
assert len(ind_v) == 1
boundary_ordered_vertexes.append(vertexes[ind_v[0]])
linked_facets = v_to_f.links(vertexes[ind_v[0]])
new_ind = np.isin(facets, linked_facets)
new_ind[ind] = False
new_ind = np.where(new_ind)[0]
if len(new_ind) == 1:
ind = new_ind[0]
vertexes = f_to_v.links(facets[ind])
boundary_pts = mesh.geometry.x[boundary_ordered_vertexes]
else:
boundary_pts = []
all_boundary_pts.append(boundary_pts)
all_facets.append(np.array(facets, dtype=np.int32))
return all_facets, all_boundary_pts
boundary_ft, boundary_pts = find_domain_boundary(mesh, ct, [1])
num_ft = mesh.topology.index_map(mesh.topology.dim - 1).size_local
ft_values = np.zeros((num_ft,), dtype=np.int32)
ft_values[boundary_ft[0]] = 1
ft = dolfinx.mesh.meshtags(mesh, 1, np.arange(num_ft, dtype=np.int32), ft_values)
ft.name = "FacetTags"
ct.name = "CellTags"
with dolfinx.io.XDMFFile(mesh.comm, "mesh.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
mesh.topology.create_connectivity(0, mesh.topology.dim)
xdmf.write_meshtags(ft, mesh.geometry)
xdmf.write_meshtags(ct, mesh.geometry)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
mesh = xdmf.read_mesh(name="mesh")
ft_orig = xdmf.read_meshtags(mesh, name="FacetTags")
ct_orig = xdmf.read_meshtags(mesh, name="CellTags")
I exported the mesh to xdmf format and checked on paraview that the result is as expected.
Having the need to re-import the mesh and tags, I tried to re-import everything but I get as error:
ft_orig = xdmf.read_meshtags(mesh, name="FacetTags")
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
line 280, in read_meshtags
mt = super().read_meshtags(mesh._cpp_object, name, xpath)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: Mesh entities of dimension 1have not been created.
Does anyone know how to fix this?
I leave attached the reference script and the image with the correct facet tags created.
Thanks in advance.