Gmsh is used to create very complex meshes in my applications. The whole process of mesh generation is very time consuming as there are many boolean operations, such as holes and cuts. But luckily, the needed mesh has a lot of repeating units. So I am wondering if I can only create a single repeating unit with matching interface nodes, and translate it along the repeating direction. Can the multiple repeating units be then stitched to a single mesh.Mesh object via DOLFINx?
The following figure demonstrates the process: (a) two connected mesh blocks are created with matching interfaces. (They are deliberately set apart) (b) the two mesh blocks are imported to FEniCSx and stitched to a single mesh object. The stitch_meshesin the MWE is the needed function.
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)
FEniCSx is used as the final computational tool and also an intermediate tool to help build stitched mesh. And yes @Ekrem_Ekici, it is great if you can handle all the stitching processing inside Gmsh in a convenient and efficient manner. See the following problem description. Please do help if you have some good ideas on my problem.
The following figure shows the needed mesh in my use case: there are a lot of repeating units in the structure and they are connected by some specific parts. It involves lot of holes inside each repeating unit, where those Boolean difference operations make the mesh generation processing extremely time consuming. Hence I am seeking a mesh stitching method here.
