XDMF -> GMSH Physical Groups -> Attributes

[jorgensd]

Uses DOLFINx to read in a mesh in xdmf format, converts the physical groups (cell markers) to attributes that can be visualized in gmsh msh files.

"""Convert DOLFINx mesh to msh format using gmsh

Author: Jørgen S. Dokken
SPDX-License-Identifier: MIT
"""

from mpi4py import MPI
import dolfinx
import numpy as np


assert MPI.COMM_WORLD.size == 1
with dolfinx.io.XDMFFile(MPI.COMM_SELF, "../../wild-web/marked_brain.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh()
    ct = xdmf.read_meshtags(mesh, name="mesh_tags")

num_cells_local = mesh.topology.index_map(mesh.topology.dim).size_local
g_indices = dolfinx.mesh.entities_to_geometry(mesh, mesh.topology.dim, np.arange(num_cells_local, dtype=np.int32))

import gmsh
gmsh.initialize()

element_tags = np.arange(num_cells_local, dtype=np.int32) + 1

unique_vals = np.unique(ct.values)
node_tags = np.arange(mesh.geometry.x.shape[0], dtype=np.int32) + 1
gmsh.model.addDiscreteEntity(0, 1)
gmsh.model.mesh.addNodes(0, 1, node_tags, mesh.geometry.x.flatten())

t1 = gmsh.view.add("TestElementData")
unique_vals = np.unique(ct.values)
data = []
num_elements = 0
for val in unique_vals:
    pos = ct.find(val)
    x_dofs = g_indices[pos].flatten()+1
    data.append(np.full(len(pos), val))
    # node_coordinates = mesh.geometry.x[g_indices[pos]]
    # node_coordinates = np.swapaxes(node_coordinates, 1, 2) 
    # node_coordinates = np.append(node_coordinates.reshape(node_coordinates.shape[0], -1), np.full((len(pos),4), val), axis=1).flatten()

    # data.append(node_coordinates)
    # num_elements += len(pos)
    volume = gmsh.model.addDiscreteEntity(3, tag=val)
    gmsh.model.mesh.addElementsByType(val, 4, element_tags[pos], x_dofs)
    gmsh.model.addPhysicalGroup(3, [val], val, f"{val}")

gmsh.view.addModelData(t1, 1, gmsh.model.getCurrent(), "ElementData", ct.indices+1, ct.values.reshape(-1,1),
                       numComponents=1)
#gmsh.view.add_list_data(t1, "SS", num_elements, np.hstack(data))
gmsh.view.option.setNumber(t1, "TimeStep", 1)
gmsh.view.write(t1, "test_view.msh")
gmsh.model.mesh.remove_duplicate_nodes()
gmsh.model.mesh.remove_duplicate_elements()

breakpoint()

gmsh.write("brain.msh")
gmsh.finalize()
gmsh.initialize()
mesh_data = dolfinx.io.gmshio.read_from_msh("brain.msh", MPI.COMM_WORLD, 0)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "new_brain.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh_data.mesh)
    xdmf.write_meshtags(mesh_data.cell_tags, mesh_data.mesh.geometry)