Hi. I want to create a mesh like the one below and define the domains of the unit squares separately while the top and middle rectangles are another domain. However, during this process, I encounter the info message “Cannot bind existing OpenCASCADE surface 6050 to second tag 7040,” which results in the entire surface being defined as one unit. Consequently, I am unable to utilize the frag boolean operator. I would greatly appreciate an explanation for this issue.
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")
gmsh.initialize()
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)
fuse1 = gmsh.model.occ.fuse([(2, SurfaceRb)], [(2, SurfaceRt)])
print("fuse1", fuse1)
fuse11 = gmsh.model.occ.fuse([(2, SurfaceB)], [(2, SurfaceD)])
print("fuse11", fuse11)
fuse111 = gmsh.model.occ.fuse([(2, SurfaceC)], [(2, SurfaceA)])
print("fuse111", fuse111)
fuse1111 = gmsh.model.occ.fuse(fuse11[0],
fuse111[0])
print("fuse1111", fuse1111)
fuse11111 = gmsh.model.occ.fuse(fuse1[0],
fuse1111[0])
print("fuse11111", fuse11111)
gmsh.model.occ.synchronize()
main_entities = fuse1[0].copy()
gmsh.model.occ.synchronize()
# Create physical groups
pg1 = gmsh.model.addPhysicalGroup(2, [entity[1] for entity in main_entities])
gmsh.model.setPhysicalName(2, pg1, "main_geometry")
pg1 = gmsh.model.addPhysicalGroup(2, [SurfaceRt, SurfaceRb])
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)
