# Fail to read tagged facets from mesh generated with GMSH **URL:** \ **Category:** mesh\ **Tags:** io, gmsh\ **Created:** [February 25, 2025, 7:58pm UTC](https://fenicsproject.discourse.group/t/fail-to-read-tagged-facets-from-mesh-generated-with-gmsh/17082 "2025-02-25T19:58:57Z")\ **Posts on this page:** 1\ **Showing post:** 2
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [February 25, 2025, 8:36pm UTC](https://fenicsproject.discourse.group/t/fail-to-read-tagged-facets-from-mesh-generated-with-gmsh/17082/2 "2025-02-25T20:36:05Z")
You are here trying to integrate over the vertices of your mesh? Then you are not integrating over a facet (“ds”) but a vertex (“dp”). However, lets take 10 steps backwards, and consider: ```python import meshio import gmsh gmsh.initialize() gmsh.model.add("my model") c_r = [0, 0, 0] c_R = [0, 0, 0] r = 1 R = 2 resolution = 0.1 print(f"Mesh resolution = {resolution}") disk_r = gmsh.model.occ.addDisk(c_r[0], c_r[1], c_r[2], r, r) disk_R = gmsh.model.occ.addDisk(c_R[0], c_R[1], c_R[2], R, R) ring = gmsh.model.occ.cut([(2, disk_R)], [(2, disk_r)]) p_1 = gmsh.model.occ.addPoint((r + R) / 2.0, 0, 0) p_2 = gmsh.model.occ.addPoint((r + R) / 2.0, r / 10.0, 0) line_p_1_p_2 = gmsh.model.occ.addLine(p_1, p_2) gmsh.model.occ.synchronize() # add 2-dimensional objects surfaces = gmsh.model.occ.getEntities(dim=2) assert surfaces == ring[0] disk_subdomain_id = 6 gmsh.model.addPhysicalGroup(surfaces[0][0], [surfaces[0][1]], disk_subdomain_id) gmsh.model.setPhysicalName(surfaces[0][0], disk_subdomain_id, "disk") # add 1-dimensional objects lines = gmsh.model.occ.getEntities(dim=1) circle_r_subdomain_id = 1 circle_R_subdomain_id = 2 line_p_1_p_2_subdomain_id = 3 gmsh.model.addPhysicalGroup(lines[0][0], [lines[0][1]], circle_r_subdomain_id) gmsh.model.setPhysicalName(lines[0][0], circle_r_subdomain_id, "circle_r") gmsh.model.addPhysicalGroup(lines[1][0], [lines[1][1]], circle_R_subdomain_id) gmsh.model.setPhysicalName(lines[1][0], circle_R_subdomain_id, "circle_R") gmsh.model.addPhysicalGroup(lines[2][0], [lines[2][1]], line_p_1_p_2_subdomain_id) gmsh.model.setPhysicalName(lines[2][0], line_p_1_p_2_subdomain_id, "line_p_1_p_2") # add 0-dimensional objects vertices = gmsh.model.occ.getEntities(dim=0) p_1_subdomain_id = 4 p_2_subdomain_id = 5 gmsh.model.addPhysicalGroup(vertices[0][0], [vertices[0][1]], p_1_subdomain_id) gmsh.model.setPhysicalName(vertices[0][0], p_1_subdomain_id, "p_1") gmsh.model.addPhysicalGroup(vertices[1][0], [vertices[1][1]], p_2_subdomain_id) gmsh.model.setPhysicalName(vertices[1][0], p_2_subdomain_id, "p_2") gmsh.model.occ.synchronize() gmsh.model.mesh.generate(2) gmsh.write("mesh.msh") def create_mesh(mesh, cell_type, prune_z=False): cells = mesh.get_cells_type(cell_type) cell_data = mesh.get_cell_data("gmsh:physical", cell_type) out_mesh = meshio.Mesh( points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data]}, ) return out_mesh mesh_from_file = meshio.read("mesh.msh") # create a triangle mesh in which the surfaces will be stored triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True) meshio.write("triangle_mesh.xdmf", triangle_mesh) # create a line mesh line_mesh = create_mesh(mesh_from_file, "line", True) meshio.write("line_mesh.xdmf", line_mesh) # create a vertex mesh import numpy as np line_vertices = np.unique(line_mesh.cells_dict["line"].flatten()) triangle_vertices = np.unique(triangle_mesh.cells_dict["triangle"].flatten()) print(line_vertices) print(triangle_vertices) print(np.isin(line_vertices, triangle_vertices)) ``` This yields a list of what vertices are in your cells, and what vertices are in your marked facets. This is: ```python [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36] [ 0 1 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60] [ True True False False True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True] ``` Note that the last line indicates that there are some vertices in your tagged facets that are not in the cells. (2, 3). This is because you have not correctly embedded then in your mesh with gmsh. By simple trial and failure, we can identify the lines: ```python gmsh.model.addPhysicalGroup(lines[2][0], [lines[2][1]], line_p_1_p_2_subdomain_id) gmsh.model.setPhysicalName(lines[2][0], line_p_1_p_2_subdomain_id, "line_p_1_p_2") ``` as the problematic ones, which goes back to: ```python p_1 = gmsh.model.occ.addPoint((r + R) / 2.0, 0, 0) p_2 = gmsh.model.occ.addPoint((r + R) / 2.0, r / 10.0, 0) line_p_1_p_2 = gmsh.model.occ.addLine(p_1, p_2) ``` which is not part of your mesh. This can thus be fixed by adding: ```python gmsh.model.occ.synchronize() ring = gmsh.model.occ.cut([(2, disk_R)], [(2, disk_r)]) gmsh.model.occ.synchronize() p_1 = gmsh.model.occ.addPoint((r + R) / 2.0, 0, 0) p_2 = gmsh.model.occ.addPoint((r + R) / 2.0, r / 10.0, 0) line_p_1_p_2 = gmsh.model.occ.addLine(p_1, p_2) gmsh.model.occ.synchronize() gmsh.model.mesh.embed(1, [line_p_1_p_2], 2, ring[0][0][0]) gmsh.model.occ.synchronize() ``` Complete mesh generating script is: ```python import meshio import gmsh gmsh.initialize() gmsh.model.add("my model") c_r = [0, 0, 0] c_R = [0, 0, 0] r = 1 R = 2 resolution = 0.1 print(f"Mesh resolution = {resolution}") disk_r = gmsh.model.occ.addDisk(c_r[0], c_r[1], c_r[2], r, r) disk_R = gmsh.model.occ.addDisk(c_R[0], c_R[1], c_R[2], R, R) gmsh.model.occ.synchronize() ring = gmsh.model.occ.cut([(2, disk_R)], [(2, disk_r)]) gmsh.model.occ.synchronize() p_1 = gmsh.model.occ.addPoint((r + R) / 2.0, 0, 0) p_2 = gmsh.model.occ.addPoint((r + R) / 2.0, r / 10.0, 0) line_p_1_p_2 = gmsh.model.occ.addLine(p_1, p_2) gmsh.model.occ.synchronize() gmsh.model.mesh.embed(1, [line_p_1_p_2], 2, ring[0][0][0]) gmsh.model.occ.synchronize() # add 2-dimensional objects surfaces = gmsh.model.occ.getEntities(dim=2) assert surfaces == ring[0] disk_subdomain_id = 6 gmsh.model.addPhysicalGroup(surfaces[0][0], [surfaces[0][1]], disk_subdomain_id) gmsh.model.setPhysicalName(surfaces[0][0], disk_subdomain_id, "disk") # add 1-dimensional objects lines = gmsh.model.occ.getEntities(dim=1) circle_r_subdomain_id = 1 circle_R_subdomain_id = 2 line_p_1_p_2_subdomain_id = 3 gmsh.model.addPhysicalGroup(lines[0][0], [lines[0][1]], circle_r_subdomain_id) gmsh.model.setPhysicalName(lines[0][0], circle_r_subdomain_id, "circle_r") gmsh.model.addPhysicalGroup(lines[1][0], [lines[1][1]], circle_R_subdomain_id) gmsh.model.setPhysicalName(lines[1][0], circle_R_subdomain_id, "circle_R") gmsh.model.addPhysicalGroup(lines[2][0], [lines[2][1]], line_p_1_p_2_subdomain_id) gmsh.model.setPhysicalName(lines[2][0], line_p_1_p_2_subdomain_id, "line_p_1_p_2") # add 0-dimensional objects vertices = gmsh.model.occ.getEntities(dim=0) p_1_subdomain_id = 4 p_2_subdomain_id = 5 gmsh.model.addPhysicalGroup(vertices[0][0], [vertices[0][1]], p_1_subdomain_id) gmsh.model.setPhysicalName(vertices[0][0], p_1_subdomain_id, "p_1") gmsh.model.addPhysicalGroup(vertices[1][0], [vertices[1][1]], p_2_subdomain_id) gmsh.model.setPhysicalName(vertices[1][0], p_2_subdomain_id, "p_2") gmsh.model.occ.synchronize() gmsh.model.mesh.generate(2) gmsh.write("mesh.msh") def create_mesh(mesh, cell_type, prune_z=False): cells = mesh.get_cells_type(cell_type) cell_data = mesh.get_cell_data("gmsh:physical", cell_type) out_mesh = meshio.Mesh( points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data]}, ) return out_mesh mesh_from_file = meshio.read("mesh.msh") # create a triangle mesh in which the surfaces will be stored triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True) meshio.write("triangle_mesh.xdmf", triangle_mesh) # create a line mesh line_mesh = create_mesh(mesh_from_file, "line", True) meshio.write("line_mesh.xdmf", line_mesh) # create a vertex mesh import numpy as np line_vertices = np.unique(line_mesh.cells_dict["line"].flatten()) triangle_vertices = np.unique(triangle_mesh.cells_dict["triangle"].flatten()) ``` For future posts, please note the following: 1. Your error happens way before the end of your code. This indicates that you shouldn’t post the rest of the code. 2. You should visually inspect your markers in paraview to see if they make sense. For the problem at hand, it was quite clear that the facet that you had marked is not in the mesh (use coarse resolution). --- _[View the full topic](https://fenicsproject.discourse.group/t/fail-to-read-tagged-facets-from-mesh-generated-with-gmsh/17082)._