Create Mesh through combination of two xdmf files

mesh
1

Hello everyone,

I have created two parts, saved them as a xdmf file and would now like to create a single mesh out of them, so that I can define my FunctionSpaces on the total domain.

I know how to do it for a mesh stored in a single xdmf.file:

with io.XDMFFile(MPI.COMM_WORLD, “Meshfile.xdmf”, “r”) as xdmf:
domain = xdmf.read_mesh(name=“Grid”)
cell_tag = xdmf.read_meshtags(domain, name=“Grid”)

But I do not know, how to combine two dolfinx.mesh.Mesh types together.

I tried combining the meshes using meshio, but there the difficulty lies in keeping the definition of the subdomains.

Thank you in advance for any assistance you can provide. I’m open to any suggestions or solutions that might help resolve this issue. Your help is greatly appreciated!

Martin

2

Hi Martin,

What do you mean by

the difficulty lies in keeping the definition of the subdomains

when combining the meshes using meshio. Couldn’t you use a cell_data field that would be, say 1 on the cells of your first mesh and 2 on the cells of the second, so that you could then extract meshtags in FEniCS ?

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)
4

Dear Dokken,

Thank you very much for your fast reply, I will try this!