Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
# ------------------------------------------------------------------------------
#
# Gmsh Python tutorial for meshing a layered sediment basin stratigraphy
#
# Geometry basics, elementary entities, physical groups
#
# ------------------------------------------------------------------------------
# The Python API is entirely defined in the `gmsh.py' module (which contains the
# full documentation of all the functions in the API):
import numpy
import gmsh
# Before using any functions in the Python API, Gmsh must be initialized:
gmsh.initialize()
# By default Gmsh will not print out any messages: in order to output messages
# on the terminal, just set the "General.Terminal" option to 1:
gmsh.option.setNumber("General.Terminal", 1)
# Next we add a new model named "t1" (if gmsh.model.add() is not called a new
# unnamed model will be created on the fly, if necessary):
gmsh.model.add("basinMesh")
# The Python API provides direct access to each supported geometry kernel. The
# built-in kernel is used in this first tutorial: the corresponding API
# functions have the `gmsh.model.geo' prefix.
# The first type of `elementary entity' in Gmsh is a `Point'. To create a point
# with the built-in geometry kernel, the Python API function is
# gmsh.model.geo.addPoint():
# - the first 3 arguments are the point coordinates (x, y, z)
# - the next (optional) argument is the target mesh size (the "characteristic
# length") close to the point
# - the last (optional) argument is the point tag (a stricly positive integer
# that uniquely identifies the point)
# The distribution of the mesh element sizes will be obtained by interpolation
# of these characteristic lengths throughout the geometry. Another method to
# specify characteristic lengths is to use general mesh size Fields (see
# `t10.py'). A particular case is the use of a background mesh (see `t7.py').
#
# If no target mesh size of provided, a default uniform coarse size will be used
# for the model, based on the overall model size.
# Units: meter
km = 1e3 # km in m
w = 120 * km # width
d = 7.5 * km # depth
lc = 1e3 # m
# We can then define some additional points. All points should have different tags:
gmsh.model.geo.addPoint(w / 2, 0, 0, lc, 1)
gmsh.model.geo.addPoint(-w / 2, 0, 0, lc, 2)
gmsh.model.geo.addPoint(-w / 2, -d, 0, lc, 3)
gmsh.model.geo.addPoint(w / 2, -d, 0, lc, 4)
# If the tag is not provided explicitly, a new tag is automatically created, and
# returned by the function:
pLay1centr = gmsh.model.geo.addPoint(0, -d / 6, 0, lc)
pLay2centr = gmsh.model.geo.addPoint(0, -d / 3, 0, lc)
pLay3centr = gmsh.model.geo.addPoint(0, -d / 2, 0, lc)
# Auxiliary points for the sediment layers
pLay1north = gmsh.model.geo.addPoint(-w / 8, 0, 0, lc)
pLay1south = gmsh.model.geo.addPoint(w / 8, 0, 0, lc)
pLay2north = gmsh.model.geo.addPoint(-w / 6, 0, 0, lc)
pLay2south = gmsh.model.geo.addPoint(w / 6, 0, 0, lc)
pLay3north = gmsh.model.geo.addPoint(-w / 4, 0, 0, lc)
pLay3south = gmsh.model.geo.addPoint(w / 4, 0, 0, lc)
# Curves are Gmsh's second type of elementery entities, and, amongst curves,
# straight lines are the simplest. The API to create straight line segments with
# the built-in kernel follows the same conventions: the first 2 arguments are
# point tags (the start and end points of the line), and the last (optional one)
# is the line tag.
#
# Note that curve tags are separate from point tags - hence we can reuse tag `1'
# for our first curve. And as a general rule, elementary entity tags in Gmsh
# have to be unique per geometrical dimension.
# gmsh.model.geo.addLine(1, 2, 1) # (start, end, tag)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)
gmsh.model.geo.addLine(pLay3south, pLay2south, 5)
gmsh.model.geo.addLine(pLay2south, pLay1south, 6)
gmsh.model.geo.addLine(pLay1south, pLay1north, 7)
gmsh.model.geo.addLine(pLay1north, pLay2north, 8)
gmsh.model.geo.addLine(pLay2north, pLay3north, 9)
gmsh.model.geo.addLine(1, pLay3south, 10)
gmsh.model.geo.addLine(pLay3north, 2, 14)
gmsh.model.geo.addSpline([pLay1south, pLay1centr, pLay1north], 11)
gmsh.model.geo.addSpline([pLay2south, pLay2centr, pLay2north], 12)
gmsh.model.geo.addSpline([pLay3south, pLay3centr, pLay3north], 13)
# The third elementary entity is the surface. In order to define a surface
# from the curves defined above, a curve loop has first to be defined.
# A curve loop is defined by an ordered list of connected curves,
# a sign being associated with each curve (depending on the orientation of the
# curve to form a loop). The API function to create curve loops takes a list
# of integers as first argument, and the curve loop tag (which must be unique
# amongst curve loops) as the second (optional) argument:
# e.g. gmsh.model.geo.addCurveLoop([4, 1, 2, 3], 1)
gmsh.model.geo.addCurveLoop([7, -11], 1)
gmsh.model.geo.addCurveLoop([6, 11, 8, -12], 2)
gmsh.model.geo.addCurveLoop([5, 12, 9, -13], 3)
gmsh.model.geo.addCurveLoop([2, 3, 4, 10, 13, 14], 4)
# Add plane surfaces defined by one or more curve loops. The first curve
# loop defines the exterior contour; additional curve loop define holes.
# (only one here, representing the external contour, since there are no holes
# --see `t4.py' for an example of a surface with a hole):
for l in range(1, 5):
gmsh.model.geo.addPlaneSurface([l], l)
# At this level, Gmsh knows everything to display the surfaces and
# to mesh it. An optional step is needed if we want to group elementary
# geometrical entities into more meaningful groups, e.g. to define some
# mathematical ("domain", "boundary"), functional ("left wing", "fuselage") or
# material ("steel", "carbon") properties.
#
# Such groups are called "Physical Groups" in Gmsh. By default, if physical
# groups are defined, Gmsh will export in output files only mesh elements that
# belong to at least one physical group. (To force Gmsh to save all elements,
# whether they belong to physical groups or not, set the `Mesh.SaveAll' option
# to 1.) Physical groups are also identified by tags, i.e. stricly positive
# integers, that should be unique per dimension (0D, 1D, 2D or 3D). Physical
# groups can also be given names.
#
# Here we define physical curves that groups the left, bottom, top and right
# curves in single groups for the one-dimensional boundary meshes
dim = 1
North = gmsh.model.addPhysicalGroup(dim, [2])
gmsh.model.setPhysicalName(dim, North, "Left")
Base = gmsh.model.addPhysicalGroup(dim, [3])
gmsh.model.setPhysicalName(dim, Base, "Bottom")
South = gmsh.model.addPhysicalGroup(dim, [4])
gmsh.model.setPhysicalName(dim, South, "Right")
Ground = gmsh.model.addPhysicalGroup(dim, [10, 5, 6, 7, 8, 9, 14])
gmsh.model.setPhysicalName(dim, Ground, "Top")
# ... and we define physical surfaces for the different
# material domains containing the corresponding geometrical surfaces:
dim = 2
Lay1 = gmsh.model.addPhysicalGroup(dim, [1])
gmsh.model.setPhysicalName(dim, Lay1, "SedimentLayer1")
Lay2 = gmsh.model.addPhysicalGroup(dim, [2])
gmsh.model.setPhysicalName(dim, Lay2, "SedimentLayer2")
Lay3 = gmsh.model.addPhysicalGroup(dim, [3])
gmsh.model.setPhysicalName(dim, Lay3, "SedimentLayer3")
Bed = gmsh.model.addPhysicalGroup(dim, [4])
gmsh.model.setPhysicalName(dim, Bed, "RockBed")
# Before it can be meshed, the internal CAD representation must be synchronized
# with the Gmsh model, which will create the relevant Gmsh data structures. This
# is achieved by the gmsh.model.geo.synchronize() API call for the built-in
# geometry kernel. Synchronizations can be called at any time, but they involve
# a non trivial amount of processing; so while you could synchronize the
# internal CAD data after every CAD command, it is usually better to minimize
# the number of synchronization points.
gmsh.model.geo.synchronize()
# We can then generate a 2D mesh...
gmsh.model.mesh.generate(2)
# ... and save it to disk: msh=current format, msh22 older format
gmsh.write("mesh_basin.msh")
# Remember that by default, if physical groups are defined, Gmsh will export in
# the output mesh file only those elements that belong to at least one physical
# group. To force Gmsh to save all elements, you can use
#
# gmsh.option.setNumber("Mesh.SaveAll", 1)
# By default, Gmsh saves meshes in the latest version of the Gmsh mesh file
# format (the `MSH' format). You can save meshes in other mesh formats by
# specifying a filename with a different extension. For example
# gmsh.write("basinMesh.pdf") only accessible via GUI
# gmsh.write("basinMesh.pvtu")
# gmsh.write("basinMesh.vtk")
# gmsh.write("basinMesh.unv")
# gmsh.write("basinMesh.stl")
# will save the mesh in the vtk, unv and stl format.
# To visualize the model we can run the graphical user interface with:
gmsh.fltk.run()
# Note that starting with Gmsh 3.0, models can be built using other geometry
# kernels than the default "built-in" kernel. To use the OpenCASCADE geometry
# kernel instead of the built-in kernel, you should use the functions with the
# `gmsh.model.occ' prefix.
#
# Different geometry kernels have different features. With OpenCASCADE, instead
# of defining the surface by successively defining 4 points, 4 curves and 1
# curve loop, one can define the rectangular surface directly with
#
# gmsh.model.occ.addRectangle(.2, 0, 0, .1, .3)
#
# Boolean operation using OpenCascade module (occ, see below)
# gmsh.model.occ.cut([(3, 1)], [(3, 2)], 3)
#
# After synchronization with the Gmsh model with
#
# gmsh.model.occ.synchronize()
#
# the underlying curves and points could be accessed with
# gmsh.model.getBoundary().
#
# See e.g. `t16.py', `t18.py', `t19.py' or `t20.py' for complete examples based
# on OpenCASCADE, and `demos/api' for more.
# This should be called when you are done using the Gmsh Python API:
gmsh.finalize()