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:
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:
[ 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:
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:
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:
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:
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:
- Your error happens way before the end of your code. This indicates that you shouldn’t post the rest of the code.
- 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).