You can start from:
The most tricky bit (when in parallel) is to remove the duplicate vertices at the shared boundary.
Even in serial, for your example, as you are actually reflecting the mesh, it is not the same nodes that should match eachother, which would mean that one would have to do a node by node matching, that is expensive.
Will the assumption always be that you have one component that can be translated, rotated etc, to match the next one? If so, I would rather send in a single grid, and instructions on how each block should be merged into the next.
A general stitching command (which would have to check for node-to-node matching) is quite involved. If the nodes do not match, you should use dolfinx_mpc: GitHub - jorgensd/dolfinx_mpc: Extension for dolfinx to handle multi-point constraints.
if you assume that you have a single mesh, where the same interface always should merge with the same interface of the same mesh, it becomes easier (but still a bit tedious):
# # Create multicomponent mesh by copying, translating, rotating and merging a single component
# Author: Jørgen S. Dokken
# SPDX License Identitfier: MIT
from mpi4py import MPI
from dolfinx import mesh, io as _io, graph
import numpy as np
import ufl
M = 5
assert MPI.COMM_WORLD.size == 1, "This script is for serial execution only."
domain = mesh.create_unit_cube(MPI.COMM_WORLD, M, M, M)
# Tag facets prior to shifting the mesh
def interface(x):
return np.isclose(x[0], 1.0)
def other_boundary(x):
return np.isclose(x[1], 1.0)
def some_boundary(x):
return np.isclose(x[2], 0.0)
def yet_another_boundary(x):
return np.isclose(x[0], 0.0)
f_right = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, interface)
f_top = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, other_boundary)
f_some = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, some_boundary)
f_yet_another = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, yet_another_boundary)
domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim)
all_values = np.full(domain.topology.index_map(domain.topology.dim - 1).size_local, 0, dtype=np.int32)
all_values[f_top] = 2
all_values[f_right] = 1
all_values[f_some] = 3
all_values[f_yet_another] = 4
indices = np.flatnonzero(all_values)
ft = mesh.meshtags(domain, domain.topology.dim - 1, indices, all_values[indices])
# Copy geometry prior to initial translation
domain.geometry.x[:, 0] -= 1.0
domain.geometry.x[:, 1] -= 0.5
domain.geometry.x[:, 1] *= - 0.5 * domain.geometry.x[:, 0] + 0.5
new_geometry = domain.geometry.x.copy()
# Move middle of geometry to origin to flip
new_geometry[:, 0] += 0.5
# Flip x-coordinate
new_geometry[:, 0] *= -1
# Translate to common boundary
new_geometry[:, 0] += 0.5
new_node_start = domain.geometry.x.shape[0]
cells = mesh.entities_to_geometry(domain, domain.topology.dim, np.arange(domain.topology.index_map(domain.topology.dim).size_local, dtype=np.int32))
new_cells = cells.copy() + new_node_start
replacement_nodes = np.unique(mesh.entities_to_geometry(domain, domain.topology.dim - 1, ft.find(1)).flatten())
keep_indices = np.full(new_geometry.shape[0], True, dtype=bool)
keep_indices[replacement_nodes] = False
new_geometry = new_geometry[keep_indices]
replacement_map = np.arange(new_node_start + domain.geometry.x.shape[0], dtype=np.int32)
# Figure out how many nodes have been removed prior to each node
all_nodes_subtractor = np.zeros(domain.geometry.x.shape[0], dtype=np.int32)
all_nodes_subtractor[replacement_nodes] = 1
all_nodes_subtractor = np.cumsum(all_nodes_subtractor)
replacement_map[new_node_start:] -= all_nodes_subtractor
# Replace those that should actually be replaced
replacement_map[replacement_nodes+new_node_start] = replacement_nodes
mapped_new_cells = replacement_map[new_cells]
all_cells = np.vstack([cells, mapped_new_cells]).astype(np.int64)
all_geometry = np.vstack([domain.geometry.x, new_geometry])
new_mesh = mesh.create_mesh(MPI.COMM_WORLD, cells=all_cells, x=all_geometry[:, :domain.geometry.dim], e=ufl.Mesh(domain.ufl_domain().ufl_coordinate_element()))
# Map facet tags to new mesh
facets_as_vertices = mesh.entities_to_geometry(domain, domain.topology.dim - 1, ft.indices)
facet_values = ft.values
new_facets = facets_as_vertices.copy() + new_node_start
new_facets = replacement_map[new_facets]
all_facets = np.vstack([facets_as_vertices, new_facets]).astype(np.int64)
all_facet_values = np.hstack([facet_values, facet_values])
local_facet, local_val = _io.distribute_entity_data(new_mesh, new_mesh.topology.dim-1, all_facets, all_facet_values)
new_mesh.topology.create_connectivity(new_mesh.topology.dim - 1, new_mesh.topology.dim)
new_ft = mesh.meshtags_from_entities(new_mesh, new_mesh.topology.dim - 1, graph.adjacencylist(local_facet), local_val)
with _io.XDMFFile(MPI.COMM_WORLD, "mwe_merge.xdmf", "w") as xdmf:
xdmf.write_mesh(new_mesh)
xdmf.write_meshtags(new_ft, new_mesh.geometry)
yields:
