Create Mesh through combination of two xdmf files

mesh
3

Here is an example of merging two meshes which has both facet and cell markers.
In the following example I use gmsh to make two separate xdmf files with marker information in them. Next I read them into dolfinx and combine the meshes there.

from mpi4py import MPI

import gmsh
import numpy as np
import dolfinx.io.gmshio

# Generate two separate meshes with gmsh
gmsh.initialize()
gmsh.model.add("Mesh1")
gmsh.model.setCurrent("Mesh1")

box_1 = gmsh.model.occ.addBox(0, 0, 0, 1, 1, 1)
gmsh.model.occ.synchronize()
volumes = gmsh.model.getEntities(dim=3)
assert len(volumes) == 1
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], 1)
boundaries = gmsh.model.getBoundary(volumes, oriented=False)
c = 3
for boundary in boundaries:
    gmsh.model.addPhysicalGroup(boundary[0], [boundary[1]], c)
    c += 1
gmsh.model.mesh.generate(3)

mesh1, ct1, ft1 = dolfinx.io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh1.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh1)
    xdmf.write_meshtags(ct1, mesh1.geometry)
    xdmf.write_meshtags(ft1, mesh1.geometry)


gmsh.model.add("Mesh2")
gmsh.model.setCurrent("Mesh2")
box_2 = gmsh.model.occ.addBox(2, 2, 2, 3, 3, 3)
gmsh.model.occ.synchronize()
volumes = gmsh.model.getEntities(dim=3)
assert len(volumes) == 1
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], 2)
boundaries = gmsh.model.getBoundary(volumes, oriented=False)
for boundary in boundaries:
    gmsh.model.addPhysicalGroup(boundary[0], [boundary[1]], c)
    c += 1
gmsh.model.mesh.generate(3)
mesh2, ct2, ft2 = dolfinx.io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh2.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh2)
    xdmf.write_meshtags(ct2, mesh2.geometry)
    xdmf.write_meshtags(ft2, mesh2.geometry)


gmsh.finalize()


# Read meshes data into dolfinx
with dolfinx.io.XDMFFile(MPI.COMM_SELF, "mesh1.xdmf", "r") as xdmf:
    mesh_1 = xdmf.read_mesh()
    mesh_1.topology.create_connectivity(mesh_1.topology.dim - 1, mesh_1.topology.dim)
    ct_1 = xdmf.read_meshtags(mesh_1, name="Cell tags")
    ft_1 = xdmf.read_meshtags(mesh_1, name="Facet tags")


with dolfinx.io.XDMFFile(MPI.COMM_SELF, "mesh2.xdmf", "r") as xdmf:
    mesh_2 = xdmf.read_mesh()
    ct_2 = xdmf.read_meshtags(mesh_2, name="Cell tags")
    mesh_2.topology.create_connectivity(mesh_2.topology.dim - 1, mesh_2.topology.dim)
    ft_2 = xdmf.read_meshtags(mesh_2, name="Facet tags")


def extract_local_information(mesh, ct, ft):
    tdim = mesh.topology.dim
    num_cells = mesh.topology.index_map(tdim).size_local
    mesh.topology.create_connectivity(tdim, tdim)
    x_indices = dolfinx.mesh.entities_to_geometry(
        mesh, tdim, np.arange(num_cells, dtype=np.int32), False
    )
    num_nodes = mesh.geometry.index_map().size_local
    x = mesh.geometry.x[:num_nodes]

    ct_indices = dolfinx.mesh.entities_to_geometry(mesh, ct.dim, ct.indices, False)
    ct_values = ct.values
    mesh.topology.create_connectivity(ft.dim, tdim)

    ft_indices = dolfinx.mesh.entities_to_geometry(mesh, ft.dim, ft.indices, False)
    ft_values = ft.values
    return (
        x_indices,
        x,
        mesh.ufl_domain().ufl_coordinate_element(),
        ct_indices,
        ct_values,
        ft_indices,
        ft_values,
    )


# Extract information from each mesh
data_1 = extract_local_information(mesh_1, ct_1, ft_1)
data_2 = extract_local_information(mesh_2, ct_2, ft_2)

# Combine data
cells = np.vstack([data_1[0], data_2[0] + data_1[1].shape[0]])
points = np.vstack([data_1[1], data_2[1]])
ct_indices = np.vstack([data_1[3], data_2[3] + data_1[1].shape[0]])
ct_values = np.hstack([data_1[4], data_2[4]])
ft_indices = np.vstack([data_1[5], data_2[5] + data_1[1].shape[0]])
ft_values = np.hstack([data_1[6], data_2[6]])
assert data_1[2] == data_2[2]
# Create combined mesh and tags
mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells, points, data_1[2])

c, c_values = dolfinx.io.distribute_entity_data(
    mesh, mesh.topology.dim, ct_indices, ct_values
)
f, f_values = dolfinx.io.distribute_entity_data(
    mesh, mesh.topology.dim - 1, ft_indices, ft_values
)

ct = dolfinx.mesh.meshtags_from_entities(
    mesh,
    mesh.topology.dim,
    dolfinx.graph.adjacencylist(c.astype(np.int32)),
    c_values.astype(np.int32),
)
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
ft = dolfinx.mesh.meshtags_from_entities(
    mesh,
    mesh.topology.dim - 1,
    dolfinx.graph.adjacencylist(f.astype(np.int32)),
    f_values.astype(np.int32),
)
ct.name = "Cell tags"
ft.name = "Facet tags"

# Write to file
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "combined_mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(ct, mesh.geometry)
    xdmf.write_meshtags(ft, mesh.geometry)
2 Likes