Thank you very much! @conpierce8
It solves the error but my mesh still looks like this:
Also, I want my whole mesh to be continuous (main and filled part). Do I need to fuse all them together?
A few points:
- It seemed to me that creating surfaces with holes, e.g. using
gmsh.model.occ.addPlaneSurface([Curve1, Curve2, Curve3], ...), was creating problems for the boolean operations. Creating solid surfaces, i.e.gmsh.model.occ.addPlaneSurface([Curve1], ...)and using boolean operations to make the holes worked better for me. - Use
gmsh.model.occ.fragmentto create conforming internal interfaces. - You should synchronize the OpenCascade kernel before adding physical groups.
- You should apply the mesh size to points after performing the boolean operations, since the boolean operations add new points. This is the reason your mesh size was so distorted.
Does the following do what you want?
import gmsh
import math
import sys
gmsh.initialize()
gmsh.clear()
gmsh.model.add("t3")
P1 = gmsh.model.occ.addPoint(-0.5, -1, 0, 0.01)
P2 = gmsh.model.occ.addPoint(0.5, -1, 0, 0.01)
P3 = gmsh.model.occ.addPoint(-0.5, -0.5, 0, 0.01)
P4 = gmsh.model.occ.addPoint(0.5, -0.5, 0, 0.01)
P5= gmsh.model.occ.addPoint(0.5, -0.1, 0, 0.01)
P6 = gmsh.model.occ.addPoint(-0.5, -0.1, 0, 0.01)
P7 = gmsh.model.occ.addPoint(-0.4, -0.2, 0, 0.01)
P8 = gmsh.model.occ.addPoint(-0.4, -0.4, 0, 0.01)
P9 = gmsh.model.occ.addPoint(0.3, -0.4, 0, 0.01)
P10 = gmsh.model.occ.addPoint(0.3, -0.2, 0, 0.01)
P11 = gmsh.model.occ.addPoint(-0.4, -0.6, 0, 0.01)
P12 = gmsh.model.occ.addPoint(-0.4, -0.8, 0, 0.01)
P13= gmsh.model.occ.addPoint(0.3, -0.8, 0,0.01)
P14 = gmsh.model.occ.addPoint(0.3, -0.6, 0, 0.01)
P15 = gmsh.model.occ.addPoint(0.05, 0.2, 0, 0.01)
P16 = gmsh.model.occ.addPoint(-0.05, 0.2, 0, 0.01)
P17 = gmsh.model.occ.addPoint(-0.0, 0.2, 0, 0.01)
P18 = gmsh.model.occ.addPoint(-0.0, 0.25, 0, 0.01)
P19 = gmsh.model.occ.addPoint(-0.05, 0.3, 0, 0.01)
P20 = gmsh.model.occ.addPoint(0.05, 0.3, 0, 0.01)
L1 = gmsh.model.occ.addLine(P1, P2)
L2 =gmsh.model.occ.addLine(P2, P4)
L4 = gmsh.model.occ.addLine(P3, P1)
L9 = gmsh.model.occ.addLine(P3, P6)
L10 = gmsh.model.occ.addLine(P4, P5)
L11 = gmsh.model.occ.addLine(P8, P9)
L12 = gmsh.model.occ.addLine(P9, P10)
L13= gmsh.model.occ.addLine(P10, P7)
L14= gmsh.model.occ.addLine(P7, P8)
L15 =gmsh.model.occ.addLine(P16, P19)
L16= gmsh.model.occ.addLine(P19, P20)
L17 =gmsh.model.occ.addLine(P20, P15)
L18 = gmsh.model.occ.addLine(P15, P16)
P21 = gmsh.model.occ.addPoint(-0.4, -0.1, 0, 0.01)
P22 = gmsh.model.occ.addPoint(0.4, -0.1, 0, 0.01)
L19= gmsh.model.occ.addSpline([P6, P18, P5])
L20 = gmsh.model.occ.addSpline([P22, P17, P21])
L21 = gmsh.model.occ.addLine(P21, P22)
Curve1 = gmsh.model.occ.addCurveLoop([L19, -L10, -L2, -L1, -L4, L9])
Curve2 = gmsh.model.occ.addCurveLoop([L20, L21])
Curve3 = gmsh.model.occ.addCurveLoop([L13, L14, L11, L12])
Curve4 = gmsh.model.occ.addCurveLoop([L16, L17, L18, L15])
Surface1 = gmsh.model.occ.addPlaneSurface([Curve1], 1)
Surface2 = gmsh.model.occ.addPlaneSurface([Curve2], 2)
Surface3 = gmsh.model.occ.addPlaneSurface([Curve3], 3)
Surface4 = gmsh.model.occ.addPlaneSurface([Curve4], 4)
fuse1 = gmsh.model.occ.fuse([(2, Surface1)], [(2, Surface4)])
print(fuse1, "\n")
frag1 = gmsh.model.occ.fragment(fuse1[0], [(2, Surface2)])
print(frag1)
main_entities = frag1[1][0].copy()
for entity in frag1[1][1]:
main_entities.remove(entity)
print(" ", main_entities, "\n")
diff1 = gmsh.model.occ.cut(frag1[1][0], [(2, Surface3)])
print(diff1, "\n")
# Synchronize the CAD model
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")
pg2 = gmsh.model.addPhysicalGroup(2, [entity[1] for entity in frag1[1][1]])
gmsh.model.setPhysicalName(2, pg2, "filled_geometry")
# Set mesh size at all points
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), 0.01)
# elementary volumes:
#gmsh.model.addPhysicalGroup(2, [1, 2], 0)
gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
gmsh.model.mesh.generate()
gmsh.fltk.run()
gmsh.finalize()
2 Likes
Yes, thanks a lot! It is working very well!
Sorry I have one other question, when I try to implement the same logic to the code below, I am not having a good result. When I print(frag1) it is creating a new surface. @conpierce8
a = 1.3
b = 0.75
c = 2
aa = 8.65
d2 = 0.0
e = 6.1
f = 0.4
g = 2
h = 5.3
j = 0.1
k = 7.0
m = 0.05
gmsh.clear()
gmsh.model.add("t3")
lc = 1e-1 # mesh
P1 = gmsh.model.occ.addPoint(b, 0, 0, lc)
P2 = gmsh.model.occ.addPoint(b, g, 0, lc)
P3 = gmsh.model.occ.addPoint(0, 3 + (h - 5.3), 0, lc)
P4 = gmsh.model.occ.addPoint(0, e + (h - 5.3), 0, lc)
P5 = gmsh.model.occ.addPoint(0, e + f + (h - 5.3), 0, lc)
P6 = gmsh.model.occ.addPoint(a, h, 0, lc)
P7 = gmsh.model.occ.addPoint(a + 1.5 * f / 3, h + 2.5 * f / 3, 0, lc)
P8 = gmsh.model.occ.addPoint(a + 0.95, k, 0, lc)
P9 = gmsh.model.occ.addPoint(a + 1.05, k, 0, lc)
P10 = gmsh.model.occ.addPoint(b + c, k, 0, lc)
P11 = gmsh.model.occ.addPoint(0, aa + d2, 0, lc)
P12 = gmsh.model.occ.addPoint(0.5, aa + d2 + 2 * j, 0, lc)
P13 = gmsh.model.occ.addPoint(0.5, 10, 0, lc)
P14 = gmsh.model.occ.addPoint(-b, 0, 0, lc)
P15 = gmsh.model.occ.addPoint(-b, g, 0, lc)
P17 = gmsh.model.occ.addPoint(-a, h, 0, lc)
P18 = gmsh.model.occ.addPoint(-a - 1.5 * f / 3, h + 2.5 * f / 3, 0, lc)
P19 = gmsh.model.occ.addPoint(-a - 0.95, k, 0, lc)
P20 = gmsh.model.occ.addPoint(-a - 1.05, k, 0, lc)
P21 = gmsh.model.occ.addPoint(-b - c, k, 0, lc)
P23 = gmsh.model.occ.addPoint(-0.5, aa + d2 + 2 * j, 0, lc)
P24 = gmsh.model.occ.addPoint(-0.5, 10, 0, lc)
P25 = gmsh.model.occ.addPoint(1.55 + j / 2 + d2 / 2, 8.1 + j + d2 + (k - 7) / 2, 0, lc) # depends on k??
P26 = gmsh.model.occ.addPoint(-1.55 - j / 2 - d2 / 2, 8.1 + j + d2 + (k - 7) / 2, 0, lc)
# spline points
P27 = gmsh.model.occ.addPoint(a + 0.95, k + j, 0, lc)
P28 = gmsh.model.occ.addPoint(-a - 0.95, k + j, 0, lc)
P29 = gmsh.model.occ.addPoint(a + 1.05, k + 2 * j, 0, lc)
P30 = gmsh.model.occ.addPoint(-a - 1.05, k + 2 * j, 0, lc)
# spline points and rectangular
P31 = gmsh.model.occ.addPoint(0.2, aa + d2, 0, lc)
P32 = gmsh.model.occ.addPoint(-0.2, aa + d2, 0, lc)
P33 = gmsh.model.occ.addPoint(-0.2, 10, 0, lc)
P34 = gmsh.model.occ.addPoint(0.2, 10, 0, lc)
P35 = gmsh.model.occ.addPoint(0, 8.8, 0, lc)
L1 = gmsh.model.occ.addLine(P14, P1)
L2 = gmsh.model.occ.addLine(P1, P2)
L3 = gmsh.model.occ.addLine(P14, P15)
L4 = gmsh.model.occ.addLine(P2, P10)
L5 = gmsh.model.occ.addLine(P15, P21)
L6 = gmsh.model.occ.addLine(P3, P6)
L7 = gmsh.model.occ.addLine(P3, P17)
L8 = gmsh.model.occ.addLine(P7, P8)
L9 = gmsh.model.occ.addLine(P18, P19)
L10 = gmsh.model.occ.addLine(P10, P9)
L11 = gmsh.model.occ.addLine(P21, P20)
L20 = gmsh.model.occ.addLine(P19, P28)
L21 = gmsh.model.occ.addLine(P20, P30)
L22 = gmsh.model.occ.addLine(P8, P27)
L23 = gmsh.model.occ.addLine(P9, P29)
L24 = gmsh.model.occ.addLine(P32, P31)
L25 = gmsh.model.occ.addLine(P31, P34)
L26 = gmsh.model.occ.addLine(P34, P33)
L27 = gmsh.model.occ.addLine(P33, P32)
L15 = gmsh.model.occ.addSpline([P6, P4, P17])
L16 = gmsh.model.occ.addSpline([P7, P5, P18])
L17 = gmsh.model.occ.addSpline([P27, P11, P28])
L28 = gmsh.model.occ.addSpline([P30, P35, P29])
Curve1 = gmsh.model.occ.addCurveLoop([L28, -L23, -L10, -L4, -L2, -L1, L3, L5, L11, L21])
Curve2 = gmsh.model.occ.addCurveLoop([L17, -L20, -L9, -L16, L8, L22])
Curve3 = gmsh.model.occ.addCurveLoop([L7, -L15, -L6])
Curve4 = gmsh.model.occ.addCurveLoop([L27, L24, L25, L26])
Surface1 = gmsh.model.occ.addPlaneSurface([Curve1], 1)
Surface2 = gmsh.model.occ.addPlaneSurface([Curve2], 2)
Surface3 = gmsh.model.occ.addPlaneSurface([Curve3], 3)
Surface4 = gmsh.model.occ.addPlaneSurface([Curve4], 4)
fuse1 = gmsh.model.occ.fuse([(2, Surface1)], [(2, Surface4)])
print(fuse1, "\n")
frag1 = gmsh.model.occ.fragment(fuse1[1][0], [(2, Surface2)])
print(frag1)
main_entities = frag1[1][0].copy()
for entity in frag1[1][1]:
main_entities.remove(entity)
print(" ", main_entities, "\n")
print(entity)
diff1 = gmsh.model.occ.cut(frag1[1][0], [(2, Surface3)])
print(diff1, "\n")
# Synchronize the CAD model
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")
pg2 = gmsh.model.addPhysicalGroup(2, [entity[1] for entity in frag1[1][1]])
gmsh.model.setPhysicalName(2, pg2, "filled_geometry")
# Set mesh size at all points
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), 0.01)
# elementary volumes:
# gmsh.model.addPhysicalGroup(2, [1, 2], 0)
gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
gmsh.model.mesh.generate()
It looks like the boolean operations are not very robust, and the results of the fuse operation are not always the same (i.e. sometimes you get one surface, sometimes two). With that being said, it is possible to use the outputs of the fragment function to identify which surfaces should be part of main_geometry and which should be part of filled_geometry.
The documentation is not very clear, but it looks like the outDimTagsMap output of fragment identifies the output entity tags that were contained within each of the input entities. In the code below, frag1[0] lists all the surfaces created by the fragment command, and frag1[1][-1] lists the subset of surfaces that were contained within Surface2. Thus frag1[1][-1] lists all the surfaces that should be inside filled_geometry, and the remaining surfaces make up main_geometry.
import gmsh
import math
import sys
a = 1.3
b = 0.75
c = 2
aa = 8.65
d2 = 0.0
e = 6.1
f = 0.4
g = 2
h = 5.3
j = 0.1
k = 7.0
m = 0.05
gmsh.initialize()
gmsh.clear()
gmsh.model.add("t3")
lc = 1e-1 # mesh
P1 = gmsh.model.occ.addPoint(b, 0, 0, lc)
P2 = gmsh.model.occ.addPoint(b, g, 0, lc)
P3 = gmsh.model.occ.addPoint(0, 3 + (h - 5.3), 0, lc)
P4 = gmsh.model.occ.addPoint(0, e + (h - 5.3), 0, lc)
P5 = gmsh.model.occ.addPoint(0, e + f + (h - 5.3), 0, lc)
P6 = gmsh.model.occ.addPoint(a, h, 0, lc)
P7 = gmsh.model.occ.addPoint(a + 1.5 * f / 3, h + 2.5 * f / 3, 0, lc)
P8 = gmsh.model.occ.addPoint(a + 0.95, k, 0, lc)
P9 = gmsh.model.occ.addPoint(a + 1.05, k, 0, lc)
P10 = gmsh.model.occ.addPoint(b + c, k, 0, lc)
P11 = gmsh.model.occ.addPoint(0, aa + d2, 0, lc)
P12 = gmsh.model.occ.addPoint(0.5, aa + d2 + 2 * j, 0, lc)
P13 = gmsh.model.occ.addPoint(0.5, 10, 0, lc)
P14 = gmsh.model.occ.addPoint(-b, 0, 0, lc)
P15 = gmsh.model.occ.addPoint(-b, g, 0, lc)
P17 = gmsh.model.occ.addPoint(-a, h, 0, lc)
P18 = gmsh.model.occ.addPoint(-a - 1.5 * f / 3, h + 2.5 * f / 3, 0, lc)
P19 = gmsh.model.occ.addPoint(-a - 0.95, k, 0, lc)
P20 = gmsh.model.occ.addPoint(-a - 1.05, k, 0, lc)
P21 = gmsh.model.occ.addPoint(-b - c, k, 0, lc)
P23 = gmsh.model.occ.addPoint(-0.5, aa + d2 + 2 * j, 0, lc)
P24 = gmsh.model.occ.addPoint(-0.5, 10, 0, lc)
P25 = gmsh.model.occ.addPoint(1.55 + j / 2 + d2 / 2, 8.1 + j + d2 + (k - 7) / 2, 0, lc) # depends on k??
P26 = gmsh.model.occ.addPoint(-1.55 - j / 2 - d2 / 2, 8.1 + j + d2 + (k - 7) / 2, 0, lc)
# spline points
P27 = gmsh.model.occ.addPoint(a + 0.95, k + j, 0, lc)
P28 = gmsh.model.occ.addPoint(-a - 0.95, k + j, 0, lc)
P29 = gmsh.model.occ.addPoint(a + 1.05, k + 2 * j, 0, lc)
P30 = gmsh.model.occ.addPoint(-a - 1.05, k + 2 * j, 0, lc)
# spline points and rectangular
P31 = gmsh.model.occ.addPoint(0.2, aa + d2, 0, lc)
P32 = gmsh.model.occ.addPoint(-0.2, aa + d2, 0, lc)
P33 = gmsh.model.occ.addPoint(-0.2, 10, 0, lc)
P34 = gmsh.model.occ.addPoint(0.2, 10, 0, lc)
P35 = gmsh.model.occ.addPoint(0, 8.8, 0, lc)
P100 = gmsh.model.occ.addPoint(-a - 1.05, 8.8, 0, lc)
P101 = gmsh.model.occ.addPoint(a + 1.05, 8.8, 0, lc)
L1 = gmsh.model.occ.addLine(P14, P1)
L2 = gmsh.model.occ.addLine(P1, P2)
L3 = gmsh.model.occ.addLine(P14, P15)
L4 = gmsh.model.occ.addLine(P2, P10)
L5 = gmsh.model.occ.addLine(P15, P21)
L6 = gmsh.model.occ.addLine(P3, P6)
L7 = gmsh.model.occ.addLine(P3, P17)
L8 = gmsh.model.occ.addLine(P7, P8)
L9 = gmsh.model.occ.addLine(P18, P19)
L10 = gmsh.model.occ.addLine(P10, P9)
L11 = gmsh.model.occ.addLine(P21, P20)
L20 = gmsh.model.occ.addLine(P19, P28)
L21 = gmsh.model.occ.addLine(P20, P30)
L22 = gmsh.model.occ.addLine(P8, P27)
L23 = gmsh.model.occ.addLine(P9, P29)
L24 = gmsh.model.occ.addLine(P32, P31)
L25 = gmsh.model.occ.addLine(P31, P34)
L26 = gmsh.model.occ.addLine(P34, P33)
L27 = gmsh.model.occ.addLine(P33, P32)
L15 = gmsh.model.occ.addSpline([P6, P4, P17])
L16 = gmsh.model.occ.addSpline([P7, P5, P18])
L17 = gmsh.model.occ.addSpline([P27, P11, P28])
L28 = gmsh.model.occ.addSpline([P30, P35, P29])
Curve1 = gmsh.model.occ.addCurveLoop([L28, -L23, -L10, -L4, -L2, -L1, L3, L5, L11, L21])
Curve2 = gmsh.model.occ.addCurveLoop([L17, -L20, -L9, -L16, L8, L22])
Curve3 = gmsh.model.occ.addCurveLoop([L7, -L15, -L6])
Curve4 = gmsh.model.occ.addCurveLoop([L27, L24, L25, L26])
Surface1 = gmsh.model.occ.addPlaneSurface([Curve1], 1)
Surface2 = gmsh.model.occ.addPlaneSurface([Curve2], 2)
Surface3 = gmsh.model.occ.addPlaneSurface([Curve3], 3)
Surface4 = gmsh.model.occ.addPlaneSurface([Curve4], 4)
fuse1 = gmsh.model.occ.fuse([(2, Surface1)], [(2, Surface4)])
print(fuse1, "\n")
diff1 = gmsh.model.occ.cut(fuse1[0], [(2, Surface3)])
print(diff1, "\n")
frag1 = gmsh.model.occ.fragment(diff1[0], [(2, Surface2)])
print(frag1)
filled_entities = frag1[1][-1]
main_entities = frag1[0].copy()
for entity in filled_entities:
if entity in main_entities:
main_entities.remove(entity)
print(" ", main_entities, "\n")
# Synchronize the CAD model
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")
pg2 = gmsh.model.addPhysicalGroup(2, [entity[1] for entity in filled_entities])
gmsh.model.setPhysicalName(2, pg2, "filled_geometry")
# Set mesh size at all points
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc)
# elementary volumes:
# gmsh.model.addPhysicalGroup(2, [1, 2], 0)
gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
gmsh.model.mesh.generate()
gmsh.fltk.run()
gmsh.finalize()
1 Like
Thanks! @conpierce8
With created mesh I am trying to obtain different domains, main and filled domains.
Here is the code;
from dolfin import *
import meshio
# Importing mesh from gmsh and defining surface and boundary markers
msh = meshio.read("geo.msh")
print(msh.cell_data_dict)
for key in msh.cell_data_dict["gmsh:geometrical"].keys():
if key == "triangle":
triangle_data = msh.cell_data_dict["gmsh:geometrical"][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("mesh.xdmf", triangle_mesh)
# Create mesh and define function space
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
print(mesh.topology().dim() - 1)
File("MSH.pvd") << mesh
File("MSH2.pvd") << mf
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1, 0)
class Boundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0.0)
b_c = Boundary()
b_c.mark(boundary_markers, 3)
File("MSH3.pvd") << boundary_markers
# Compiling subdomains
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx_filled = Measure("dx", domain=mesh, subdomain_data=mf, subdomain_id=2)
dx_main = Measure("dx", domain=mesh, subdomain_data=mf, subdomain_id=1)
V = VectorFunctionSpace(mesh, "CG", 2)
W_tensor = TensorFunctionSpace(mesh, "CG", 2)
...
It defines my physicals 0 and 1 I am not sure if that is a problem. However, I want to have my elementaries only 1 and 2 (for main and filled) but it is going up to 5, I could not change it.
I want my msh file to be like
…
4199 2 2 0 1 2301 2336 763
4200 2 2 0 1 2248 2250 17
4201 2 2 0 2 409 3171 2380
4202 2 2 0 2 2389 3389 16
4203 2 2 0 2 2914 3234 2948
…
Do I need to fix the msh file while reading or creating it?
I get physical entities numbered 1 and 2, and elementary entities numbered 2, 4, and 5. I don’t know of any way in this case to force Gmsh to have only two elementary entities. I guess if this is the way you would like to proceed, then you would have to edit the msh file.
However, is there a reason why you cannot convert your mesh like the following?
for key in msh.cell_data_dict["gmsh:physical"].keys():
if key == "triangle":
triangle_data = msh.cell_data_dict["gmsh:physical"][key]
When I use the above, I obtain:
2 Likes
Thank you so much, here is the code
import gmsh
import math
import sys
a = 1.3
b = 0.75
c = 2
aa = 8.65
d2 = 0.0
e = 6.1
f = 0.4
g = 2
h = 5.3
j = 0.1
k = 7.0
m = 0.05
gmsh.initialize()
gmsh.clear()
gmsh.model.add("t3")
lc = 1e-1 # mesh
P1 = gmsh.model.occ.addPoint(b, 0, 0, lc)
P2 = gmsh.model.occ.addPoint(b, g, 0, lc)
P3 = gmsh.model.occ.addPoint(0, 3 + (h - 5.3), 0, lc)
P4 = gmsh.model.occ.addPoint(0, e + (h - 5.3), 0, lc)
P5 = gmsh.model.occ.addPoint(0, e + f + (h - 5.3), 0, lc)
P6 = gmsh.model.occ.addPoint(a, h, 0, lc)
P7 = gmsh.model.occ.addPoint(a + 1.5 * f / 3, h + 2.5 * f / 3, 0, lc)
P8 = gmsh.model.occ.addPoint(a + 0.95, k, 0, lc)
P9 = gmsh.model.occ.addPoint(a + 1.05, k, 0, lc)
P10 = gmsh.model.occ.addPoint(b + c, k, 0, lc)
P11 = gmsh.model.occ.addPoint(0, aa + d2, 0, lc)
P12 = gmsh.model.occ.addPoint(0.5, aa + d2 + 2 * j, 0, lc)
P13 = gmsh.model.occ.addPoint(0.5, 10, 0, lc)
P14 = gmsh.model.occ.addPoint(-b, 0, 0, lc)
P15 = gmsh.model.occ.addPoint(-b, g, 0, lc)
P17 = gmsh.model.occ.addPoint(-a, h, 0, lc)
P18 = gmsh.model.occ.addPoint(-a - 1.5 * f / 3, h + 2.5 * f / 3, 0, lc)
P19 = gmsh.model.occ.addPoint(-a - 0.95, k, 0, lc)
P20 = gmsh.model.occ.addPoint(-a - 1.05, k, 0, lc)
P21 = gmsh.model.occ.addPoint(-b - c, k, 0, lc)
P23 = gmsh.model.occ.addPoint(-0.5, aa + d2 + 2 * j, 0, lc)
P24 = gmsh.model.occ.addPoint(-0.5, 10, 0, lc)
P25 = gmsh.model.occ.addPoint(1.55 + j / 2 + d2 / 2, 8.1 + j + d2 + (k - 7) / 2, 0, lc) # depends on k??
P26 = gmsh.model.occ.addPoint(-1.55 - j / 2 - d2 / 2, 8.1 + j + d2 + (k - 7) / 2, 0, lc)
# spline points
P27 = gmsh.model.occ.addPoint(a + 0.95, k + j, 0, lc)
P28 = gmsh.model.occ.addPoint(-a - 0.95, k + j, 0, lc)
P29 = gmsh.model.occ.addPoint(a + 1.05, k + 2 * j, 0, lc)
P30 = gmsh.model.occ.addPoint(-a - 1.05, k + 2 * j, 0, lc)
# spline points and rectangular
P31 = gmsh.model.occ.addPoint(0.2, aa + d2, 0, lc)
P32 = gmsh.model.occ.addPoint(-0.2, aa + d2, 0, lc)
P33 = gmsh.model.occ.addPoint(-0.2, 10, 0, lc)
P34 = gmsh.model.occ.addPoint(0.2, 10, 0, lc)
P35 = gmsh.model.occ.addPoint(0, 8.8, 0, lc)
P100 = gmsh.model.occ.addPoint(-a - 1.05, 8.8, 0, lc)
P101 = gmsh.model.occ.addPoint(a + 1.05, 8.8, 0, lc)
L1 = gmsh.model.occ.addLine(P14, P1)
L2 = gmsh.model.occ.addLine(P1, P2)
L3 = gmsh.model.occ.addLine(P14, P15)
L4 = gmsh.model.occ.addLine(P2, P10)
L5 = gmsh.model.occ.addLine(P15, P21)
L6 = gmsh.model.occ.addLine(P3, P6)
L7 = gmsh.model.occ.addLine(P3, P17)
L8 = gmsh.model.occ.addLine(P7, P8)
L9 = gmsh.model.occ.addLine(P18, P19)
L10 = gmsh.model.occ.addLine(P10, P9)
L11 = gmsh.model.occ.addLine(P21, P20)
L20 = gmsh.model.occ.addLine(P19, P28)
L21 = gmsh.model.occ.addLine(P20, P30)
L22 = gmsh.model.occ.addLine(P8, P27)
L23 = gmsh.model.occ.addLine(P9, P29)
L24 = gmsh.model.occ.addLine(P32, P31)
L25 = gmsh.model.occ.addLine(P31, P34)
L26 = gmsh.model.occ.addLine(P34, P33)
L27 = gmsh.model.occ.addLine(P33, P32)
L15 = gmsh.model.occ.addSpline([P6, P4, P17])
L16 = gmsh.model.occ.addSpline([P7, P5, P18])
L17 = gmsh.model.occ.addSpline([P27, P11, P28])
L28 = gmsh.model.occ.addSpline([P30, P35, P29])
Curve1 = gmsh.model.occ.addCurveLoop([L28, -L23, -L10, -L4, -L2, -L1, L3, L5, L11, L21])
Curve2 = gmsh.model.occ.addCurveLoop([L17, -L20, -L9, -L16, L8, L22])
Curve3 = gmsh.model.occ.addCurveLoop([L7, -L15, -L6])
Curve4 = gmsh.model.occ.addCurveLoop([L27, L24, L25, L26])
Surface1 = gmsh.model.occ.addPlaneSurface([Curve1], 1)
Surface2 = gmsh.model.occ.addPlaneSurface([Curve2], 2)
Surface3 = gmsh.model.occ.addPlaneSurface([Curve3], 3)
Surface4 = gmsh.model.occ.addPlaneSurface([Curve4], 4)
fuse1 = gmsh.model.occ.fuse([(2, Surface1)], [(2, Surface4)])
print(fuse1, "\n")
diff1 = gmsh.model.occ.cut(fuse1[0], [(2, Surface3)])
print(diff1, "\n")
frag1 = gmsh.model.occ.fragment(diff1[0], [(2, Surface2)])
print(frag1)
filled_entities = frag1[1][-1]
main_entities = frag1[0].copy()
for entity in filled_entities:
if entity in main_entities:
main_entities.remove(entity)
print(" ", main_entities, "\n")
# Synchronize the CAD model
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")
pg2 = gmsh.model.addPhysicalGroup(2, [entity[1] for entity in filled_entities])
gmsh.model.setPhysicalName(2, pg2, "filled_geometry")
# Set mesh size at all points
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc)
# elementary volumes:
# gmsh.model.addPhysicalGroup(2, [1, 2], 0)
gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
gmsh.model.mesh.generate()
gmsh.fltk.run()
gmsh.finalize()
from dolfin import *
# 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]
elif key == "quad":
quad_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
elif cell.type == "quad":
quad_cells = cell.data
quad_mesh = meshio.Mesh(points=msh.points,
cells=[("quad", quad_cells)],
cell_data={"name_to_read": [quad_data]})
meshio.write("./Result--" + str(commit)+ "mesh.xdmf",
quad_mesh)
meshio.write("mesh.xdmf", quad_mesh)
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True,
"eliminate_zeros": True,
"precompute_basis_const": True,
"precompute_ip_const": True}
# Create mesh and define function space
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
print(mesh.topology().dim() - 1)
File("./Result--" + str(commit)+ "MSH.pvd") << mesh
File("./Result--"+ str(commit) + "/MSH.pvd") << mesh
File("./Result--"+ str(commit) + "MSH2.pvd") << mf
File("./Result--" + str(commit)+ "/MSH2" ".pvd") << mf
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1, 0)
class Boundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0.0)
b_c = Boundary()
b_c.mark(boundary_markers, 3)
File("./Result--" + str(commit)+ "MSH3.pvd") << boundary_markers
# Compiling subdomains
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx_filled = Measure("dx", domain=mesh, subdomain_data=mf, subdomain_id=2)
dx_main = Measure("dx", domain=mesh, subdomain_data=mf, subdomain_id=1)
and I am having
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** fenics-support@googlegroups.com
*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error: Unable to order quadrilateral cell.
*** Reason: Cell is not orderable.
*** Where: This error was encountered inside QuadrilateralCell.cpp.
*** Process: 0
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: ubuntu
*** -------------------------------------------------------------------------
Can you share your complete code and the error you receive?
Thank you, I edited the previous post above.
So, as the error says, the mesh of quadrilateral elements that you have created cannot be ordered in a consistent way (this issue is related to: Unable to order quadrilateral cell)
which lead to the paper by @mscroggs : https://dl.acm.org/doi/full/10.1145/3524456 (arxiv version: [2102.11901] Construction of arbitrary order finite element degree-of-freedom maps on polygonal and polyhedral cell meshes) which has been implemented in DOLFINx.
2 Likes
Thank you so much @dokken. I have been actually looking for a solution for my contacting problem. Above, the red area is my air mesh, so I am using very soft material for that area, however when the young modulus of the void material is lower than a certain amount (1e-3 MPa for my case) when the top curve reaches the bottom one the solver stops converging. If my void is stiffer, than my results are not accurate. Do you have any idea why I am having this problem?
Hello @conpierce8,
I have a question about fusing the surfaces again,
Here I translate my mesh and I want to fuse it with the previous one and name it main_geometry again,
my trial is not working. I would like to hear if you have any idea. Thanks!
Here is my code:
def createGeometryAndMesh(a, b, c, aa, d2, e, f, g, h, j, k, lc):
# Clear all models and create a new one
gmsh.initialize()
gmsh.clear()
gmsh.model.add("t3")
# in units mm
# cube points:
def geopoints(dx, dy, surf_point):
P1 = gmsh.model.occ.addPoint(scale *b + dx, 0+dy, 0, lc)
P2 = gmsh.model.occ.addPoint(scale *b + dx, scale *g+dy, 0, lc)
P3 = gmsh.model.occ.addPoint(0+ dx,scale * m +scale * (h - 5.3) +dy, 0, lc)
P4 = gmsh.model.occ.addPoint(0+ dx, scale *e + scale *(h - 5.3)+dy, 0, lc)
P5 = gmsh.model.occ.addPoint(0+ dx,scale * e +scale * f + scale *(h - 5.3)+dy, 0, lc)
P6 = gmsh.model.occ.addPoint(scale * a+ dx, scale *h+dy, 0, lc)
P7 = gmsh.model.occ.addPoint(scale *a+ dx + scale *1.5 * f / 3, scale *h +scale * 2.5 * f / 3 +dy, 0, lc)
P8 = gmsh.model.occ.addPoint(scale *a+ dx + scale *0.95, scale *k+dy, 0, lc)
P9 = gmsh.model.occ.addPoint(scale *a+ dx + scale *1.05, scale *k+dy, 0, lc)
P10 = gmsh.model.occ.addPoint(scale *b+ dx + scale *c, scale * k+dy, 0, lc)
P11 = gmsh.model.occ.addPoint(scale *0+ dx, scale *aa +scale * d2+dy, 0, lc)
# P12 = gmsh.model.occ.addPoint(0.5, aa + d2 + 2 * j, 0, lc)
# P13 = gmsh.model.occ.addPoint(0.5, 10, 0, lc)
P14 = gmsh.model.occ.addPoint(-scale *b+ dx, 0+dy, 0, lc)
P15 = gmsh.model.occ.addPoint(-scale *b+ dx, scale *g+dy, 0, lc)
P17 = gmsh.model.occ.addPoint(-scale *a+ dx, scale *h+dy, 0, lc)
P18 = gmsh.model.occ.addPoint(-scale *a -scale * 1.5 * f / 3+ dx, scale *h +scale * 2.5 * f / 3+dy, 0, lc)
P19 = gmsh.model.occ.addPoint(-scale *a -scale * 0.95 + dx, scale *k+dy, 0, lc)
P20 = gmsh.model.occ.addPoint(-scale *a - scale *1.05 + dx, scale *k+dy, 0, lc)
P21 = gmsh.model.occ.addPoint(-scale *b - scale *c + dx, scale *k+dy, 0, lc)
# P23 = gmsh.model.occ.addPoint(-0.5, aa + d2 + 2 * j, 0, lc)
# P24 = gmsh.model.occ.addPoint(-0.5, 10, 0, lc)
# P25 = gmsh.model.occ.addPoint(1.55 + j / 2 + d2 / 2, 8.1 + j + d2 + (k - 7) / 2, 0, lc) # depends on k??
# P26 = gmsh.model.occ.addPoint(-1.55 - j / 2 - d2 / 2, 8.1 + j + d2 + (k - 7) / 2, 0, lc)
# spline points
P27 = gmsh.model.occ.addPoint(scale *a + scale *0.95 + dx, scale *k + scale *j+dy, 0, lc)
P28 = gmsh.model.occ.addPoint(-scale *a - scale *0.95 + dx, scale *k + scale *j+dy, 0, lc)
P29 = gmsh.model.occ.addPoint(scale *a + scale *1.05+ dx, scale *k +scale * 2 * j+dy,0, lc)
P30 = gmsh.model.occ.addPoint(-scale *a - scale *1.05+ dx, scale *k +scale * 2 * j+dy, 0, lc)
# spline points and rectangular
P31 = gmsh.model.occ.addPoint(scale *0.3+ dx, scale *aa +scale * d2+dy, 0, lc)
P32 = gmsh.model.occ.addPoint(-scale *0.3+ dx, scale *aa +scale * d2+dy, 0, lc)
P33 = gmsh.model.occ.addPoint(-scale *0.3+ dx,scale * 10+dy, 0, lc)
P34 = gmsh.model.occ.addPoint(scale *0.3+ dx,scale * 10+dy, 0, lc)
P35 = gmsh.model.occ.addPoint(scale *0+ dx, scale *aa +scale * d2 + scale *j+dy, 0, lc)
L1 = gmsh.model.occ.addLine(P14, P1)
L2 = gmsh.model.occ.addLine(P1, P2)
L3 = gmsh.model.occ.addLine(P14, P15)
L4 = gmsh.model.occ.addLine(P2, P10)
L5 = gmsh.model.occ.addLine(P15, P21)
L6 = gmsh.model.occ.addLine(P3, P6)
L7 = gmsh.model.occ.addLine(P3, P17)
L8 = gmsh.model.occ.addLine(P7, P8)
L9 = gmsh.model.occ.addLine(P18, P19)
L10 = gmsh.model.occ.addLine(P10, P9)
L11 = gmsh.model.occ.addLine(P21, P20)
# L12 = gmsh.model.occ.addLine(P24, P13)
# L13 = gmsh.model.occ.addLine(P12, P13)
# L14 = gmsh.model.occ.addLine(P23, P24)
L20 = gmsh.model.occ.addLine(P19, P28)
L21 = gmsh.model.occ.addLine(P20, P30)
L22 = gmsh.model.occ.addLine(P8, P27)
L23 = gmsh.model.occ.addLine(P9, P29)
L24 = gmsh.model.occ.addLine(P32, P31)
L25 = gmsh.model.occ.addLine(P31, P34)
L26 = gmsh.model.occ.addLine(P34, P33)
L27 = gmsh.model.occ.addLine(P33, P32)
L15 = gmsh.model.occ.addSpline([P6, P4, P17])
L16 = gmsh.model.occ.addSpline([P7, P5, P18])
L17 = gmsh.model.occ.addSpline([P27, P11, P28])
# L18 = gmsh.model.occ.addSpline([P29, P25, P12])
# L19 = gmsh.model.occ.addSpline([P30, P26, P23])
L28 = gmsh.model.occ.addSpline([P30, P35, P29])
Curve1 = gmsh.model.occ.addCurveLoop([L28, -L23, -L10, -L4, -L2, -L1, L3, L5, L11, L21])
Curve2 = gmsh.model.occ.addCurveLoop([L17, -L20, -L9, -L16, L8, L22])
Curve3 = gmsh.model.occ.addCurveLoop([L7, -L15, -L6])
Curve4 = gmsh.model.occ.addCurveLoop([L27, L24, L25, L26])
Surface1 = gmsh.model.occ.addPlaneSurface([Curve1], 1 + surf_point)
Surface2 = gmsh.model.occ.addPlaneSurface([Curve2], 2 + surf_point)
Surface3 = gmsh.model.occ.addPlaneSurface([Curve3], 3 + surf_point)
Surface4 = gmsh.model.occ.addPlaneSurface([Curve4], 4 + surf_point)
return Surface1, Surface2, Surface3, Surface4
Surface1, Surface2, Surface3, Surface4 = geopoints(0, 0, 100)
###############################
fuse1 = gmsh.model.occ.fuse([(2, Surface1)], [(2, Surface4)])
print("fuse1")
print(fuse1, "\n")
diff1 = gmsh.model.occ.cut(fuse1[0], [(2, Surface3)])
print(diff1, "\n")
frag1 = gmsh.model.occ.fragment(diff1[0], [(2, Surface2)])
print(frag1[1][-1])
print("fraggg1111")
print(frag1)
filled_entities = frag1[1][-1]
main_entities = frag1[0].copy()
print("qqqqqqqqqqqqmainentities1")
print(main_entities)
print(filled_entities)
for entity in filled_entities:
if entity in main_entities:
main_entities.remove(entity)
print(" ", main_entities, "\n")
print(entity[1])
# Synchronize the CAD model
gmsh.model.occ.synchronize()
print("entity111111")
print(entity[1])
# Create physical groups
pg1 = gmsh.model.addPhysicalGroup(2, [entity[1] for entity in main_entities])
gmsh.model.setPhysicalName(2, pg1, "main_geometry")
pg2 = gmsh.model.addPhysicalGroup(2, [entity[1] for entity in filled_entities])
gmsh.model.setPhysicalName(2, pg2, "filled_geometry")
pg10, pg20 = pg1, pg2
# Set mesh size at all points
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc)
gmsh.model.occ.synchronize()
###############################
Surface1, Surface2, Surface3, Surface4 = geopoints(0, 10, 200)
fuse1 = gmsh.model.occ.fuse([(2, Surface1)], [(2, Surface4)])
print(fuse1, "\n")
diff1 = gmsh.model.occ.cut(fuse1[0], [(2, Surface3)])
print(diff1, "\n")
frag1 = gmsh.model.occ.fragment(diff1[0], [(2, Surface2)])
print(frag1[1][-1])
print(frag1)
filled_entities = frag1[1][-1]
main_entities = frag1[0].copy()
print(main_entities)
print(filled_entities)
for entity in filled_entities:
if entity in main_entities:
main_entities.remove(entity)
print(" ", main_entities, "\n")
print(entity[1])
# Synchronize the CAD model
gmsh.model.occ.synchronize()
print(entity[1])
# Create physical groups
pg1 = gmsh.model.addPhysicalGroup(2, [entity[1] for entity in main_entities])
gmsh.model.setPhysicalName(2, pg1, "main_geometry2")
pg2 = gmsh.model.addPhysicalGroup(2, [entity[1] for entity in filled_entities])
gmsh.model.setPhysicalName(2, pg2, "filled_geometry2")
# Set mesh size at all points
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc)
gmsh.model.occ.synchronize()
#fuse2 = gmsh.model.occ.fuse([(2, pg1)], [(2, pg10)])
# Set mesh size at all points
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc)
gmsh.model.occ.synchronize()
###############################
Surface1, Surface2, Surface3, Surface4 = geopoints(10, 10, 300)
fuse1 = gmsh.model.occ.fuse([(2, Surface1)], [(2, Surface4)])
print(fuse1, "\n")
diff1 = gmsh.model.occ.cut(fuse1[0], [(2, Surface3)])
print(diff1, "\n")
frag1 = gmsh.model.occ.fragment(diff1[0], [(2, Surface2)])
print(frag1[1][-1])
print(frag1)
filled_entities = frag1[1][-1]
main_entities = frag1[0].copy()
print(main_entities)
print(filled_entities)
for entity in filled_entities:
if entity in main_entities:
main_entities.remove(entity)
print(" ", main_entities, "\n")
print(entity[1])
# Synchronize the CAD model
gmsh.model.occ.synchronize()
print(entity[1])
# Create physical groups
pg1 = gmsh.model.addPhysicalGroup(2, [entity[1] for entity in main_entities])
gmsh.model.setPhysicalName(2, pg1, "main_geometry2")
pg2 = gmsh.model.addPhysicalGroup(2, [entity[1] for entity in filled_entities])
gmsh.model.setPhysicalName(2, pg2, "filled_geometry2")
# Set mesh size at all points
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc)
gmsh.model.occ.synchronize()
#fuse2 = gmsh.model.occ.fuse([(2, pg1)], [(2, pg10)])
###############################
Surface1, Surface2, Surface3, Surface4 = geopoints(10, 0, 400)
fuse1 = gmsh.model.occ.fuse([(2, Surface1)], [(2, Surface4)])
print(fuse1, "\n")
diff1 = gmsh.model.occ.cut(fuse1[0], [(2, Surface3)])
print(diff1, "\n")
frag1 = gmsh.model.occ.fragment(diff1[0], [(2, Surface2)])
print(frag1[1][-1])
print(frag1)
filled_entities = frag1[1][-1]
main_entities = frag1[0].copy()
print(main_entities)
print(filled_entities)
for entity in filled_entities:
if entity in main_entities:
main_entities.remove(entity)
print(" ", main_entities, "\n")
print(entity[1])
# Synchronize the CAD model
gmsh.model.occ.synchronize()
print(entity[1])
# Create physical groups
pg1 = gmsh.model.addPhysicalGroup(2, [entity[1] for entity in main_entities])
gmsh.model.setPhysicalName(2, pg1, "main_geometry2")
pg2 = gmsh.model.addPhysicalGroup(2, [entity[1] for entity in filled_entities])
gmsh.model.setPhysicalName(2, pg2, "filled_geometry2")
#fuse2 = gmsh.model.occ.fuse([(2, pg1)], [(2, pg10)])
# 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)
# https://gitlab.onelab.info/gmsh/gmsh/blob/master/examples/api/bspline_filling.py
# Note that color options are special: setting a color option of "X.Y"
# actually sets the option "X.Color.Y".
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}]""")
# Create the geometry and mesh it:
# createGeometryAndMesh(a, b, c, aa, d2, e, f, g, h, j, k, m)
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()
Maybe you should perform all entity creation operations first, followed by all boolean operations?
E.g.
# Entity creation operations
SurfaceA = geopoints(0, 0, 100)
SurfaceB = geopoints(0, 10, 100)
# Boolean operations
fuse1 = gmsh.model.occ.fuse(
[(2, SurfaceA[0]), (2, SurfaceA[3]), (2, SurfaceB[0])],
[(2, SurfaceB[3])]
)
diff1 = gmsh.model.occ.cut(
fuse1[0],
[(2, SurfaceA[2]), (2, SurfaceB[2])]
)
frag1 = gmsh.model.occ.fragment(
diff1[0],
[(2, SurfaceA[1]), (2, SurfaceB[1])]
)
main_entities = frag1[0].copy()
for filled_entity in frag1[1][-2:]:
for entity in filled_entity:
if entity in main_entities:
main_entities.remove(entity)
# Synchronize...
# Create physical groups...
This is just a guess since I’m unable to test your code.
1 Like
Thank you, yes this works!
Hello @conpierce8 ,
I have a small issue with this method. When I use boolean operator for 2 surfaces, it does not preserve the surface tags after fusing. It says Cannot bind existing OpenCASCADE surface, so they replace the tag and that is why I cannot cut the surfaces later. Here I shared the code and I could not find a solution, I would be very happy if you have an idea.
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)
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)
1 Like
Yes! Thank you so much!
1 Like