Hi,
I am quite new to mesh generation and dolfinx, so I apologize if this is a somewhat trivial matter. If anyone has any suggestions or tips, that would be greatly appreciated.
Goal: Convert 2D triangulization meshes created by Triangle into mesh and facet XDMF files which can then be used with dolfinx for a variety of BVPs.
Current Approach: I am currently parsing the .poly, .node, and .ele files generated by Triangle into a python dictionary. Then, using pygmsh and meshio, I am creating a .msh file from the information extracted into the python dictionary. Finally, I am converting the .msh file into two XDMF files (one for facets in order to define boundary conditions later) using the create_mesh() function defined here.
Current Issue: Currently, a .msh file and two XDMF files are being generated from my code. Using paraview and gmsh, I am able to check that the meshes do at least seem to be created correctly (i.e. they look the same as the original triangle mesh). However, I am running into several errors (bad file descriptor or value error: vector) when trying to read the XDMF files as described in the " Defining subdomains for different materials" tutorial by Jørgen S. Dokken. Specifically, I am having trouble with the read_meshtags function.
Code:
import meshio
import pygmsh
import gmsh
import pyvista
import numpy as np
import collections
def read_triangle(poly_file, ele_file, node_file = None):
"""
simple poly-file reader, that creates a python dictionary
with information about vertices, edges and holes.
it assumes that vertices have no attributes
no regional attributes or area constraints are parsed.
if verticies are not given in the poly-file, they must be given in a node-file.
if node-file is not given, and vertices are not given in the poly-file, the function will search for a node-file with the same name as the poly-file, but with the extension .node
"""
output = {'vertices': None, 'holes': None, 'segments': None, 'triangles': None}
# open file and store lines in a list
file = open(poly_file, 'r')
lines = file.readlines()
file.close()
lines = [x.strip('\n').split() for x in lines]
# clean up lines (remove empty lines and comments)
lines = [x for x in lines if x != []]
lines = [x for x in lines if x[0] != '#']
# Divide lines into vertices, segments and holes
vertex_lines, segment_lines, hole_lines = [], [], []
# handle vertices
n_vertices, dimension, attr, bdry_markers = [int(x) for x in lines[0]]
if n_vertices == 0:
if node_file is None:
# no node-file given, so search for a node-file with the same name as the poly-file, but with the extension .node
node_file = poly_file[:-4] + 'node'
# open file and store lines in a list
try:
file = open(node_file, 'r')
vertex_lines = file.readlines()
file.close()
vertex_lines = [x.strip('\n').split() for x in vertex_lines]
vertex_lines = [x for x in vertex_lines if x != []]
vertex_lines = [x for x in vertex_lines if x[0] != '#']
except:
raise Exception(f"no vertices given in poly-file and no node-file found at {node_file}")
# append vertex lines to lines, so that the rest of the code can be used as is.
lines = vertex_lines + lines[1:]
# vertex stats
n_vertices, dimension, attr, bdry_markers = [int(x) for x in lines[0]]
vertex_lines = lines[1:n_vertices+1]
# segment stats
n_segments, bdry_markers = [int(x) for x in lines[n_vertices+1]]
segment_lines = lines[n_vertices+2:n_vertices+n_segments+2]
# hole stats
n_holes = int(lines[n_segments+n_vertices+2][0])
hole_lines = lines[n_segments + n_vertices + 3:n_segments + n_vertices + 3 + n_holes]
# store vertices
vertices = []
for line in vertex_lines:
values = [float(x) for x in line] # read as tuple in case of boundary markers
vertices.append(values[1:]) # ignore label
output['vertices'] = np.array(vertices)
# store segments
segments = []
for line in segment_lines:
values = [int(x) for x in line] # read as tuple in case of boundary markers
values[1] -= 1 # subtract 1 to get 0-based indexing
values[2] -= 1 # subtract 1 to get 0-based indexing
segments.append(values[1:]) # ignore label
output['segments'] = np.array(segments)
# store holes
holes = []
for line in hole_lines:
values = [float(x) for x in line]
holes.append(values[1:]) # ignore label
output['holes'] = np.array(holes)
# Read triangles
try:
# open file and store lines in a list
file = open(ele_file, 'r')
triangle_lines = file.readlines()
file.close()
triangle_lines = [x.strip('\n').split() for x in triangle_lines]
triangle_lines = [x for x in triangle_lines if x != []]
triangle_lines = [x for x in triangle_lines if x[0] != '#']
# store triangles
triangles = []
for line in triangle_lines[1:]:
values = [int(x) for x in line] # read as tuple in case of boundary markers
values[1] -= 1 # subtract 1 to get 0-based indexing
values[2] -= 1 # subtract 1 to get 0-based indexing
values[3] -= 1 # subtract 1 to get 0-based indexing
triangles.append(values[1:]) # ignore label
output['triangles'] = np.array(triangles)
except:
raise Exception(f"no triangles given in ele-file")
return output
def create_xdmf(cell_type, mesh_name):
"""
creates an xdmf file from a mesh file
"""
# read mesh
mesh = meshio.read(mesh_name)
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data('gmsh:physical', cell_type=cell_type)
points = mesh.points[:, :2]
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data]})
return out_mesh
def create_mesh(output_name, poly_dict):
"""
creates a mesh using pygmsh and meshio from the verticies, segments and holes in the poly_dict
"""
# initialize empty geometry
geometry = pygmsh.geo.Geometry()
model = geometry.__enter__()
# add vertices
points = [model.add_point([vertex[0], vertex[1], 0]) for vertex in poly_dict['vertices']]
# add segments
lines = collections.defaultdict(list)
for segment in poly_dict['segments']:
start = points[segment[0]]
end = points[segment[1]]
boundary = None if len(segment) == 2 else segment[2]
line = model.add_line(start, end)
if boundary:
lines[boundary].append(line)
# add triangles
triangles = []
for triangle in poly_dict['triangles']:
p1 = points[triangle[0]].x
p2 = points[triangle[1]].x
p3 = points[triangle[2]].x
#model.add_polygon([p1, p2, p3])
triangles.append(model.add_polygon([p1, p2, p3]))
# call gmsh kernel before adding physical groups
model.synchronize()
# add physical groups
for boundary, line_list in lines.items():
model.add_physical(line_list, label=f"boundary_{boundary}")
# add compound surface
geometry.add_physical(triangles, label="compound_surface")
# write mesh to file
geometry.generate_mesh(dim=2)
gmsh.write(f"{output_name}.msh")
gmsh.clear()
geometry.__exit__()
# create xdmf file
meshio.write(f"{output_name}.xdmf", create_xdmf("triangle", f"{output_name}.msh"))
meshio.write(f"{output_name}_facet.xdmf", create_xdmf("line", f"{output_name}.msh"))
def plot_mesh(mesh):
"""
plots a mesh using pyvista
"""
pyvista.set_plot_theme("document")
p = pyvista.Plotter(window_size=(800, 800))
p.add_mesh(
mesh=pyvista.from_meshio(mesh),
scalar_bar_args={"title": "Materials"},
show_scalar_bar=False,
show_edges=True,
)
p.view_xy()
p.show()
# run file
if __name__ == '__main__':
# print versions of imports
print("-------------------------------------")
print('meshio version: ', meshio.__version__)
print('pygmsh version: ', pygmsh.__version__)
print('gmsh version: ', gmsh.__version__)
print("-------------------------------------\n")
# read poly-file
poly_file = 'meshes/two_holes/two_holes.1.poly'
node_file = 'meshes/two_holes/two_holes.1.node'
ele_file = 'meshes/two_holes/two_holes.1.ele'
mesh = read_triangle(poly_file, ele_file, node_file)
print()
print(f"vertices: {mesh['vertices'].shape}")
print(f"segments: {mesh['segments'].shape}")
print(f"holes: {mesh['holes'].shape}")
print(f"triangles: {mesh['triangles'].shape}")
print()
# create mesh
create_mesh('test', mesh)
# plot mesh
plot_mesh(meshio.read('test.msh'))
# Code used to read XDMF files
with XDMFFile(MPI.COMM_WORLD, "test.xdmf", "r") as xdmf:
mesh = xdmf.read_mesh(name="Grid")
ct = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)
with XDMFFile(MPI.COMM_WORLD, "test_facet.xdmf", "r") as xdmf:
ft = xdmf.read_mesh(mesh, name="Grid")
If there there is any more information that I can provide to help, please let me know!
Thank you.
