Hello,
I have query regarding the renumbering of element mesh data and original cell index obtained giving wrong information. Can you relate the reorder data and help me if I can avoid this renumbering.
I am facing issues when the new order is not captured by mesh.topology.original_cell_index
Please find the MWE here:
the new connectivity and coordinates are:
from dolfinx.io import gmshio
from mpi4py import MPI
mesh, subdomains, boundaries = gmshio.read_from_msh("mh104.msh", MPI.COMM_WORLD,0, gdim=3)
o_cell_idx= mesh.topology.original_cell_index
c_to_v = mesh.topology.connectivity(mesh.topology.dim,0)
coord=mesh.geometry.x
print(c_to_v)
Info : Reading 'mh104.msh'...
Info : 3434 nodes
Info : 6125 elements
Info : Done reading 'mh104.msh'
<AdjacencyList> with 6125 nodes
0: [0 1 2 ]
1: [3 4 5 ]
2: [6 7 8 ]
3: [8 7 9 ]
4: [8 9 1 ]
The element connectivity for second element (after reordering) to compare with original mesh:
Please provide more than partial snippets of code.
This does not show how you got coord.
Please also provide the script for how you generated the msh file, instead of giving the msh file itself.
No, you cannot. Dolfinx is made to work in parallel, and then it is natural that the nodes gets a reordering depending on the result of the mesh partitioner.
@dokken , I apologize for the confusion. I have provided the coord in MWE as
From the obtain reordered mesh connectivity (c_to_v), I took second element connectivity as printed in main MWE given above. My concern is that I cannot find this element reordered element index/number, in original mesh using mesh.topology.original_cell_index
At present, I have the mesh file generated from a different mesh generator using node and element connectivity data. I can generate a separate mesh in dolfinx using the same data. But, if the mesh is working and extracting all different subdomains (in disarranged order), how should the mesh generation process affects mesh.topology.original_cell_index to not work correctly.
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.