After the fuse operation, there is only one surface in the model. This is exactly what a fuse operation is supposed to do.
So of course it is impossible to preserve the tags after the fuse operation. However, based on your other post, I assume you want something like this?
from dolfin import *
import gmsh
import numpy as np
import numpy as np
import matplotlib.pyplot as plt
import meshio
a = 0.25
b = 1
x_dir = 1.5
y_dir = 1 + a
scale = 1
commit = "rectangular"
lc = 1e-1
def createGeometryAndMesh(a, b, lc):
# Clear all models and create a new one
gmsh.initialize()
gmsh.clear()
gmsh.model.add("t3")
def geopoints(dx, dy, surf_point):
P1 = gmsh.model.occ.addPoint(scale * 0 + dx, dy, 0, lc)
P2 = gmsh.model.occ.addPoint(scale * b + dx, scale * 0 + dy, 0, lc)
P3 = gmsh.model.occ.addPoint(0 + dx, scale * b + dy, 0, lc)
P4 = gmsh.model.occ.addPoint(b + dx, scale * b + dy, 0, lc)
L1 = gmsh.model.occ.addLine(P1, P2)
L2 = gmsh.model.occ.addLine(P2, P4)
L3 = gmsh.model.occ.addLine(P4, P3)
L4 = gmsh.model.occ.addLine(P3, P1)
Curve1 = gmsh.model.occ.addCurveLoop([L1, L2, L3, L4])
Surface1 = gmsh.model.occ.addPlaneSurface([Curve1], 11 + surf_point)
return Surface1
def geopoints2(dx, dy, surf_point):
P11 = gmsh.model.occ.addPoint(scale * 0 + dx, dy, 0, lc)
P22 = gmsh.model.occ.addPoint(scale * b + dx, scale * 0 + dy, 0, lc)
P33 = gmsh.model.occ.addPoint(0 + dx, scale * b + dy, 0, lc)
P44 = gmsh.model.occ.addPoint(b + dx, scale * b + dy, 0, lc)
L11 = gmsh.model.occ.addLine(P11, P22)
L22 = gmsh.model.occ.addLine(P22, P44)
L33 = gmsh.model.occ.addLine(P44, P33)
L44 = gmsh.model.occ.addLine(P33, P11)
Curve11 = gmsh.model.occ.addCurveLoop([L11, L22, L33, L44])
Surface2 = gmsh.model.occ.addPlaneSurface([Curve11], 12 + surf_point)
return Surface2
def geopoints3(dx, dy, surf_point):
P111 = gmsh.model.occ.addPoint(scale * 0 + dx, dy, 0, lc)
P222 = gmsh.model.occ.addPoint(scale * b + dx, scale * 0 + dy, 0, lc)
P333 = gmsh.model.occ.addPoint(0 + dx, scale * b + dy, 0, lc)
P444 = gmsh.model.occ.addPoint(b + dx, scale * b + dy, 0, lc)
L111 = gmsh.model.occ.addLine(P111, P222)
L222 = gmsh.model.occ.addLine(P222, P444)
L333 = gmsh.model.occ.addLine(P444, P333)
L444 = gmsh.model.occ.addLine(P333, P111)
Curve111 = gmsh.model.occ.addCurveLoop([L111, L222, L333, L444])
Surface3 = gmsh.model.occ.addPlaneSurface([Curve111], 13 + surf_point)
return Surface3
def geopoints4(dx, dy, surf_point):
P1111 = gmsh.model.occ.addPoint(scale * 0 + dx, dy, 0, lc)
P2222 = gmsh.model.occ.addPoint(scale * b + dx, scale * 0 + dy, 0, lc)
P3333 = gmsh.model.occ.addPoint(0 + dx, scale * b + dy, 0, lc)
P4444 = gmsh.model.occ.addPoint(b + dx, scale * b + dy, 0, lc)
L1111 = gmsh.model.occ.addLine(P1111, P2222)
L2222 = gmsh.model.occ.addLine(P2222, P4444)
L3333 = gmsh.model.occ.addLine(P4444, P3333)
L4444 = gmsh.model.occ.addLine(P3333, P1111)
Curve1111 = gmsh.model.occ.addCurveLoop([L1111, L2222, L3333, L4444])
Surface4 = gmsh.model.occ.addPlaneSurface([Curve1111], 14 + surf_point)
return Surface4
def rectangulartp(dx, dy, surf_point):
P5 = gmsh.model.occ.addPoint(-a, scale * 0 + b + dy, 0, lc)
P6 = gmsh.model.occ.addPoint(scale * -a, scale * 0 + b + a + dy, 0, lc)
P7 = gmsh.model.occ.addPoint(
scale * (a + dx + b), scale * 0 + b + a + dy, 0, lc
)
P8 = gmsh.model.occ.addPoint(scale * (a + dx + b), scale * 0 + b + dy, 0, lc)
L5 = gmsh.model.occ.addLine(P5, P8)
L6 = gmsh.model.occ.addLine(P8, P7)
L7 = gmsh.model.occ.addLine(P7, P6)
L8 = gmsh.model.occ.addLine(P6, P5)
Curve5 = gmsh.model.occ.addCurveLoop([L5, L6, L7, L8])
Surface5 = gmsh.model.occ.addPlaneSurface([Curve5], 40 + surf_point)
return Surface5
def rectangularbt(dx, dy, surf_point):
P50 = gmsh.model.occ.addPoint(-a, scale * 0 + b, 0, lc)
P60 = gmsh.model.occ.addPoint(scale * -a, scale * 0 + b + a, 0, lc)
P70 = gmsh.model.occ.addPoint(scale * (a + dx + b), scale * 0 + b + a, 0, lc)
P80 = gmsh.model.occ.addPoint(scale * (a + dx + b), scale * 0 + b, 0, lc)
L50 = gmsh.model.occ.addLine(P50, P80)
L60 = gmsh.model.occ.addLine(P80, P70)
L70 = gmsh.model.occ.addLine(P70, P60)
L80 = gmsh.model.occ.addLine(P60, P50)
Curve6 = gmsh.model.occ.addCurveLoop([L50, L60, L70, L80])
Surface6 = gmsh.model.occ.addPlaneSurface([Curve6], 50 + surf_point)
return Surface6
# Entity creation operations
SurfaceA = geopoints(0, 0, 2000)
SurfaceB = geopoints2(0, y_dir, 3000)
SurfaceC = geopoints3(x_dir, 0, 4000)
SurfaceD = geopoints4(x_dir, y_dir, 5000)
SurfaceRb = rectangularbt(x_dir, y_dir, 6000)
SurfaceRt = rectangulartp(x_dir, y_dir, 7000)
# SurfaceB = geopoints(x_dir+10, y_dir+10, 400)
gmsh.model.occ.synchronize()
# gmsh.fltk.run()
frag, frag_map = gmsh.model.occ.fragment(
[(2, SurfaceRb)],
[(2, SurfaceRt), (2, SurfaceB), (2, SurfaceD), (2, SurfaceC), (2, SurfaceA)],
)
print("frag, frag_map:", frag, frag_map)
gmsh.model.occ.synchronize()
# Create physical groups
# One group for each of the unit squares
for i in range(2, len(frag_map)):
pg1 = gmsh.model.addPhysicalGroup(2, [frag_map[i][0][1]])
gmsh.model.setPhysicalName(2, pg1, f"main_geometry{i}")
# One group for the top and middle rectangles (frag_map[0] and frag_map[1])
pg1 = gmsh.model.addPhysicalGroup(2, [frag_map[0][0][1], frag_map[1][0][1]])
gmsh.model.setPhysicalName(2, pg1, "filled_geometry")
# Set mesh size at all points
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc)
gmsh.model.occ.synchronize()
gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
gmsh.model.mesh.generate()
gmsh.write("./Result--" + str(commit) + "geoforstl" + ".geo_unrolled")
gmsh.write("./Result--" + str(commit) + "geo.msh")
# colors:
gmsh.option.setNumber("Geometry.PointNumbers", 1)
gmsh.option.setColor("Geometry.Points", 255, 165, 0)
gmsh.option.setColor("General.Text", 255, 255, 255)
gmsh.option.setColor("Mesh.Points", 255, 0, 0)
rr, gg, bb, aaa = gmsh.option.getColor("Geometry.Points")
gmsh.option.setColor("Geometry.Surfaces", rr, gg, bb, aaa)
gmsh.onelab.set(
"""[
{"type":"number",
"name":"Parameters/Twisting angle",
"values":[90],
"min":0,
"max":120,
"step":1}]"""
)
def checkForEvent():
action = gmsh.onelab.getString("ONELAB/Action")
if len(action) and action[0] == "check":
gmsh.onelab.setString("ONELAB/Action", [""])
createGeometryAndMesh()
gmsh.graphics.draw()
return True
# Refine the mesh
gmsh.model.mesh.refine()
gmsh.write("t3.opt")
gmsh.finalize()
createGeometryAndMesh(a, b, lc)
# Importing mesh from gmsh and defining surface and boundary markers
msh = meshio.read("./Result--" + str(commit) + "geo.msh")
for key in msh.cell_data_dict["gmsh:physical"].keys():
if key == "triangle":
triangle_data = msh.cell_data_dict["gmsh:physical"][key]
for cell in msh.cells:
if cell.type == "tetra":
tetra_cells = cell.data
elif cell.type == "triangle":
triangle_cells = cell.data
triangle_mesh = meshio.Mesh(
points=msh.points,
cells=[("triangle", triangle_cells)],
cell_data={"name_to_read": [triangle_data]},
)
meshio.write("./Result--" + str(commit) + "mesh.xdmf", triangle_mesh)
meshio.write("mesh.xdmf", triangle_mesh)
