Boolean operator gmsh

mesh
1

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)

WhatsApp Image 2022-11-25 at 11.03.53 PM
WhatsApp Image 2022-11-25 at 11.03.53 PM656Ă—809 252 KB

WhatsApp Image 2022-11-25 at 11.03.52 PM
WhatsApp Image 2022-11-25 at 11.03.52 PM685Ă—780 276 KB 

// Gmsh project created on Fri Nov 25 16:59:48 2022
SetFactory("OpenCASCADE");
//+
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()
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])

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

# Synchronize the CAD model
gmsh.model.occ.synchronize()
gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)


# elementary volumes:
gmsh.model.addPhysicalGroup(2, [1, 2], 0)
gmsh.model.mesh.generate()
gmsh.fltk.run()
gmsh.finalize()

When I just have the fuse part, I am having a result like,
WhatsApp Image 2022-11-25 at 11.03.52 PM(1)

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

2

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

Thank you very much! @conpierce8
It solves the error but my mesh still looks like this:
Screenshot 2022-11-26 111217

4

Also, I want my whole mesh to be continuous (main and filled part). Do I need to fuse all them together?

5

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.
  3. You should synchronize the OpenCascade kernel before adding physical groups.
  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()
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
6

Yes, thanks a lot! It is working very well!

7

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

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
9

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?

10

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:

screenshot
screenshot1023Ă—789 17.9 KB

2 Likes
11

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

1 Like
12

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

13

Can you share your complete code and the error you receive?

14

Thank you, I edited the previous post above.

15

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
16

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?

17

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

image
image686Ă—789 51.1 KB

18

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
19

Thank you, yes this works!

20

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)