# Boolean operator gmsh

Hello,

I want to add a rectangle to my main geometry as I shared below and combine it with the main shape using the boolean operator. I achieved that with gmsh but I could not have it while using API.
Here is my gmsh code and figure, (before and after boolean operator)

``````// Gmsh project created on Fri Nov 25 16:59:48 2022
//+
Point(1) = {-0.5, -1, 0, 1.0};
//+
Point(2) = {0.5, -1, 0, 1.0};
//+
Point(3) = {-0.5, -0.5, 0, 1.0};
//+
Point(4) = {0.5, -0.5, 0, 1.0};
//+
Point(5) = {0.5, -0.1, 0, 1.0};
//+
Point(6) = {-0.5, -0.1, 0, 1.0};
//+
Point(7) = {-0.4, -0.2, 0, 1.0};
//+
Point(8) = {-0.4, -0.4, 0, 1.0};
//+
Point(9) = {0.3, -0.4, 0, 1.0};
//+
Point(10) = {0.3, -0.2, 0, 1.0};
//+
Point(11) = {-0.4, -0.6, 0, 1.0};
//+
Point(12) = {-0.4, -0.8, 0, 1.0};
//+
Point(13) = {0.3, -0.8, 0, 1.0};
//+
Point(14) = {0.3, -0.6, 0, 1.0};
//+
Point(15) = {0.05, 0.2, 0, 1.0};
//+
Point(16) = {-0.05, 0.2, 0, 1.0};
//+
Point(17) = {-0.0, 0.2, 0, 1.0};
//+
Point(18) = {-0.0, 0.25, 0, 1.0};
//+
Point(19) = {-0.05, 0.3, 0, 1.0};
//+
Point(20) = {0.05, 0.3, 0, 1.0};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 4};
//+
Line(4) = {3, 1};
//+
Line(9) = {3, 6};
//+
Line(10) = {4, 5};
//+
Line(11) = {8, 9};
//+
Line(12) = {9, 10};
//+
Line(13) = {10, 7};
//+
Line(14) = {7, 8};
//+
Line(15) = {16, 19};
//+
Line(16) = {19, 20};
//+
Line(17) = {20, 15};
//+
Line(18) = {15, 16};
//+
Point(21) = {-0.4, -0.1, 0, 1.0};
//+
Point(22) = {0.4, -0.1, 0, 1.0};
//+
Spline(19) = {6, 18, 5};
//+
Spline(20) = {22, 17, 21};
//+
Line(21) = {21, 22};
//+
Curve Loop(1) = {19, -10, -2, -1, -4, 9};
//+
Curve Loop(2) = {20, 21};
//+
Curve Loop(3) = {13, 14, 11, 12};
//+
Plane Surface(1) = {1, 2, 3};
//+
Curve Loop(4) = {16, 17, 18, 15};
//+
Plane Surface(2) = {2};
//+
Plane Surface(3) = {4};
//+
BooleanUnion{ Surface{3}; Delete; }{ Surface{1}; Delete; }
//+
BooleanIntersection{ Surface{4}; Delete; }{ Surface{1}; }
//+
Physical Surface("main_geometry", 1) = {1, 3};
//+
Physical Surface("filled_geometry", 2) = {2};
``````

However, when I use occ.fuse function as below, I am having an error as â€śsetting an array element with a sequence. The requested array has an inhomogeneous shape after 2 dimensions. The detected shape was (1, 2) + inhomogeneous partâ€ť.
Here is the code,

``````import gmsh
import math
import sys

gmsh.initialize()
gmsh.clear()
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)
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)
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)
P21 = gmsh.model.occ.addPoint(-0.4, -0.1, 0, 0.01)
P22 = gmsh.model.occ.addPoint(0.4, -0.1, 0, 0.01)
Curve1 = gmsh.model.occ.addCurveLoop([L19, -L10, -L2, -L1, -L4, L9])
Curve3 = gmsh.model.occ.addCurveLoop([L13, L14, L11, L12])
Curve4 = gmsh.model.occ.addCurveLoop([L16, L17, L18, L15])

sf = gmsh.model.occ.fuse([(2, 1)], [(2, 3)], removeObject=True, removeTool=True)
gmsh.model.setPhysicalName(2, ps, "main_geometry")
gmsh.model.setPhysicalName(2, ps2, "filled_geometry")

gmsh.model.occ.synchronize()
gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)

# elementary volumes:
gmsh.model.mesh.generate()
gmsh.fltk.run()
gmsh.finalize()

``````

When I just have the fuse part, I am having a result like,

I would be so happy if you tell me my mistake.
Thank you so much!

Have a look at the documentation of `fuse`:

It returns a list of the geometric entities created by the fuse operation, and a map of the entities created by the fuse operation. In Python, this looks like:

``````>>> print(sf)
([(2, 1), (2, 3), (2, 4)], [[(2, 1)], [(2, 3), (2, 4)]])
``````

You are then wrapping that data in a list when you pass it to `addPhysicalGroup`; thus you are passing a data structure of the following form as the second argument to `addPhysicalGroup`:

``````list[
tuple[
list[tuple, tuple, tuple],
list[
list[tuple],
list[tuple, tuple]
]
]
]
``````

This is not the expected form of the second argument, which should be a list of `int`s. You should do the following:

``````sf = gmsh.model.occ.fuse([(2, 1)], [(2, 3)], removeObject=True, removeTool=True)
gmsh.model.occ.synchronize()
ps = gmsh.model.addPhysicalGroup(2, [entity[1] for entity in sf[0]])
``````

Note that you need to call `gmsh.model.occ.synchronize()` before you attempt to create physical groups, or else the geometric entities created by the `occ` module will not be accessible to the `gmsh.model` module.

2 Likes

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:

1. 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.
2. Use `gmsh.model.occ.fragment` to create conforming internal interfaces.
4. 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()
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)
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)
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)
P21 = gmsh.model.occ.addPoint(-0.4, -0.1, 0, 0.01)
P22 = gmsh.model.occ.addPoint(0.4, -0.1, 0, 0.01)

Curve1 = gmsh.model.occ.addCurveLoop([L19, -L10, -L2, -L1, -L4, L9])
Curve3 = gmsh.model.occ.addCurveLoop([L13, L14, L11, L12])
Curve4 = gmsh.model.occ.addCurveLoop([L16, L17, L18, L15])

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")

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.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()

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)
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])
Curve4 = gmsh.model.occ.addCurveLoop([L27, L24, L25, L26])
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")
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.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()

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)
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])
Curve4 = gmsh.model.occ.addCurveLoop([L27, L24, L25, L26])

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")

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.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

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)],

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:
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
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

You are right, thanks for your great help! @conpierce8

1 Like

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()

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)
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])
Curve4 = gmsh.model.occ.addCurveLoop([L27, L24, L25, L26])

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")

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.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]
for cell in msh.cells:
if cell.type == "tetra":
tetra_cells = cell.data
elif cell.type == "triangle":
triangle_cells = cell.data

meshio.write("./Result--" + str(commit)+ "mesh.xdmf",

# 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:
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
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

*** 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
*** -------------------------------------------------------------------------

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()
# 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)

# L18 = gmsh.model.occ.addSpline([P29, P25, P12])
# L19 = gmsh.model.occ.addSpline([P30, P26, P23])

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])
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])

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])

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])

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])

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()

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)

Curve1 = gmsh.model.occ.addCurveLoop([L1, L2, L3, L4])

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)

Curve11 = gmsh.model.occ.addCurveLoop([L11, L22, L33, L44])

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)

Curve111 = gmsh.model.occ.addCurveLoop([L111, L222, L333, L444])

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)

Curve1111 = gmsh.model.occ.addCurveLoop([L1111, L2222, L3333, L4444])

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)

Curve5 = gmsh.model.occ.addCurveLoop([L5, L6, L7, L8])

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)

Curve6 = gmsh.model.occ.addCurveLoop([L50, L60, L70, L80])

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")
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)],