I am working to create a mesh with subdomains (based on here), and I am visualizing it to make sure the subdomain assignments are correct.

The above example uses the `model_to_mesh`

function from `dolfinx.io.gmshio`

but this appears to have been removed (I get the error `dolfinx.io is not a package`

and there isn’t reference to `gmshio`

in the docs). Because of this, I switched to adapting some other examples where the gmsh geometry is directly used in FEniCSx.

When I specify the cell tags from gmsh in the pyvista code, the colors of the different subdomains are wrong. Looking at the cell tags, the number of cells associated with the different subdomains seems appropriate, so it looks like it might just be an ordering issue. **Are there ways to ensure that the cell ordering is the same between the gmsh mesh and the ufl/FEniCSx mesh?** I have the gmsh MWE below and what it generates, and what the coloring should look like (hand-colored).

```
import gmsh
import numpy as np
from mpi4py import MPI
import matplotlib.pyplot as plt
rank = MPI.COMM_WORLD.rank
gmsh.initialize()
gdim = 2 # Geometric dimension of the mesh
model_rank = 0
mesh_comm = MPI.COMM_WORLD
spacing = 0.01
body_dx = 1
body_dy = body_dx
floor_dx = 1.6
floor_dy = 0.1
indenter_dx = 0.5
indenter_dy = 0.5
if mesh_comm.rank == model_rank:
floor = gmsh.model.occ.addRectangle(-floor_dx/2,0,0,floor_dx,floor_dy)
body1 = gmsh.model.occ.addRectangle(-body_dx/2,floor_dy+spacing,0,body_dx,body_dy)
indenter = gmsh.model.occ.addRectangle(-indenter_dx/2,body_dy+floor_dy+2*spacing,0,indenter_dx,indenter_dy)
background = gmsh.model.occ.addRectangle(-floor_dx/2,0,0,floor_dx,indenter_dy+body_dy+floor_dy+2*spacing)
gmsh.model.occ.synchronize()
bodies = [(2,floor),(2,body1),(2,indenter)]
gmsh.model.occ.synchronize()
whole_domain = gmsh.model.occ.fragment([(2,background)],bodies)
gmsh.model.occ.synchronize()
background_surfaces = []
other_surfaces = []
for i,domain in enumerate(whole_domain[0]):
com = gmsh.model.occ.getCenterOfMass(domain[0],domain[1])
gmsh.model.addPhysicalGroup(domain[0], [domain[1]], tag=i)
if np.isclose(com[1],1.6002/2):
background_surfaces.append(domain)
else:
other_surfaces.append(domain)
gmsh.model.mesh.generate(2)
from dolfinx import fem, mesh, plot, io, cpp, nls, la
import pyvista
if MPI.COMM_WORLD.rank == 0:
# Get mesh geometry
geometry_data = io.extract_gmsh_geometry(gmsh.model)
# Get mesh topology for each element
topology_data = io.extract_gmsh_topology_and_markers(gmsh.model)
if MPI.COMM_WORLD.rank == 0:
# Extract the cell type and number of nodes per cell and broadcast
# it to the other processors
gmsh_cell_type = list(topology_data.keys())[0]
properties = gmsh.model.mesh.getElementProperties(gmsh_cell_type)
name, dim, order, num_nodes, local_coords, _ = properties
cells = topology_data[gmsh_cell_type]["topology"]
cell_id, num_nodes = MPI.COMM_WORLD.bcast([gmsh_cell_type, num_nodes], root=0)
else:
cell_id, num_nodes = MPI.COMM_WORLD.bcast([None, None], root=0)
cells, geometry_data = np.empty([0, num_nodes]), np.empty([0, gdim])
ufl_domain = io.ufl_mesh_from_gmsh(cell_id, 2)
gmsh_cell_perm = io.cell_perm_gmsh(cpp.mesh.to_type(str(ufl_domain.ufl_cell())), num_nodes)
cells = cells[:, gmsh_cell_perm]
# Create distributed mesh
domain = mesh.create_mesh(MPI.COMM_WORLD, cells, geometry_data[:, :gdim], ufl_domain)
import pyvista
pyvista.set_jupyter_backend("pythreejs")
from dolfinx.plot import create_vtk_mesh
plotter = pyvista.Plotter()
topology, cell_types, geometry = create_vtk_mesh(domain, domain.topology.dim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry )
num_local_cells = domain.topology.index_map(domain.topology.dim).size_local
grid.cell_data["Marker"] = topology_data[2]['cell_data']
grid.set_active_scalars("Marker")
actor = plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
pyvista.start_xvfb()
cell_tag_fig = plotter.screenshot("cell_tags.png")
```

**What the Code Generates:**

**What I am Looking to Generate:**