There are a few things to consider here.
Python has a zero index as a starting point, GMSH starts with index 1. That means that you need to add 1 to any result that you get from DOLFINx (as the input index relates to the index in a python data structure.
Secondly, we do not associate the original_cell_index with the index it had in GMSH, as in the msh format, there might be cells that are not part of the actual mesh.
Following is a code that highlights that the vertices have the correct global numbering, but it then also shows you that for DOLFINx to read a msh file, we extract cells by their physical group, and order them in an array (which the input global index is relative to):
from mpi4py import MPI
import dolfinx
import numpy as np
mesh, subdomains, boundaries = dolfinx.io.gmshio.read_from_msh("mh104.msh", MPI.COMM_WORLD,0, gdim=3)
o_cell_idx= mesh.topology.original_cell_index
c_to_node = dolfinx.mesh.entities_to_geometry(mesh, mesh.topology.dim, np.arange(len(o_cell_idx), dtype=np.int32))
# Sort entities to make them comparable
c_to_node = np.sort(c_to_node, axis=1)
coord=mesh.geometry.x
o_node_idx = mesh.geometry.input_global_indices
# Inspecting gmsh:
# 1116 2 2 2 276 2140 2145 2141
# 2140 0 5.740535104259405e-05 0.1428603951628746
# 2145 0 8.756520587749979e-05 0.1428477232139085
# 2141 0 0.0001658069490290491 0.1428511737801536
# Get vertices of given cell (subtract 1 from gmsh to get correct python index)
v0 = np.flatnonzero(o_node_idx == 2139)
v1 = np.flatnonzero(o_node_idx == 2144)
v2 = np.flatnonzero(o_node_idx == 2140)
input_vertices = np.sort(np.hstack([v0,v1,v2]))
# Find cell with matching indices
cell_index_python = np.flatnonzero(np.all(input_vertices == c_to_node, axis=1))
print(o_cell_idx[cell_index_python], coord[input_vertices])
print(coord[v0], coord[v1], coord[v2])
# Inspect pure gmsh
import gmsh
gmsh.initialize()
gmsh.model.add("Mesh from file")
gmsh.merge("mh104.msh")
breakpoint()
# Extract all cells as done by DOLFINx when reading in the mesh (but do not subtract the 1 from the gmsh topology)
# for easy comparison
phys_grps = gmsh.model.getPhysicalGroups()
topologies = []
for dim, tag in phys_grps:
if dim == 2:
entities = gmsh.model.getEntitiesForPhysicalGroup(dim, tag)
for entity in entities:
(entity_types, entity_tags, entity_topologies) = gmsh.model.mesh.getElements(dim, tag=entity)
assert len(entity_types) == 1
properties = gmsh.model.mesh.getElementProperties(entity_types[0])
name, dim, _, num_nodes, _, _ = properties
topology = entity_topologies[0].reshape(-1, num_nodes)
topologies.append(topology)
print(np.vstack(topologies)[o_cell_idx[cell_index_python]])
You can create your own gmsh reader that respect the ordering in the msh file,
import gmsh
gmsh.initialize()
gmsh.model.add("Mesh from file")
gmsh.merge("mh104.msh")
# Extract cells in original order
topologies_2 = []
entities = gmsh.model.getEntities()
entities_sorted = np.sort([e[1] for e in entities])
for obj in entities_sorted:
(entity_types, entity_tags, entity_topologies) = gmsh.model.mesh.getElements(2, tag=obj)
assert len(entity_types) == 1
properties = gmsh.model.mesh.getElementProperties(entity_types[0])
name, dim, _, num_nodes, _, _ = properties
topology = entity_topologies[0].reshape(-1, num_nodes)
topologies_2.append(topology)
topologies_2 = np.vstack(topologies_2)
print(topologies_2[1115])
and then replicate everything in
(if you want facet tags you need to call get entities(1) as well.