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
|
# -----------------------------------------------------------------------------
#
# Gmsh Julia extended tutorial 2
#
# Mesh import, discrete entities, hybrid models, terrain meshing
#
# -----------------------------------------------------------------------------
import gmsh
# The API can be used to import a mesh without reading it from a file, by
# creating nodes and elements on the fly and storing them in model
# entities. These model entities can be existing CAD entities, or can be
# discrete entities, entirely defined by the mesh.
#
# Discrete entities can be reparametrized (see `t13.jl') so that they can be
# remeshed later on; and they can also be combined with built-in CAD entities to
# produce hybrid models.
#
# We combine all these features in this tutorial to perform terrain meshing,
# where the terrain is described by a discrete surface (that we then
# reparametrize) combined with a CAD representation of the underground.
gmsh.initialize()
gmsh.model.add("x2")
# We will create the terrain surface mesh from N x N input data points:
N = 100
# Helper function to return a node tag given two indices i and j:
function tag(i, j)
return (N + 1) * i + j + 1
end
# The x, y, z coordinates of all the nodes:
coords = []
# The tags of the corresponding nodes:
nodes = []
# The connectivities of the triangle elements (3 node tags per triangle) on the
# terrain surface:
tris = []
# The connectivities of the line elements on the 4 boundaries (2 node tags
# for each line element):
lin = [[], [], [], []]
# The connectivities of the point elements on the 4 corners (1 node tag for each
# point element):
pnt = [tag(0, 0), tag(N, 0), tag(N, N), tag(0, N)]
for i in 0:N
for j in 0:N
push!(nodes, tag(i, j))
push!(coords, float(i) / N, float(j) / N,
0.05 * sin(10 * float(i + j) / N))
if i > 0 && j > 0
push!(tris, tag(i - 1, j - 1), tag(i, j - 1), tag(i - 1, j))
push!(tris, tag(i, j - 1), tag(i, j), tag(i - 1, j))
end
if (i == 0 || i == N) && j > 0
push!(lin[i == 0 ? 4 : 2], tag(i, j - 1), tag(i, j))
end
if (j == 0 || j == N) && i > 0
push!(lin[j == 0 ? 1 : 3], tag(i - 1, j), tag(i, j))
end
end
end
# Create 4 discrete points for the 4 corners of the terrain surface:
for i in 0:3
gmsh.model.addDiscreteEntity(0, i + 1)
end
gmsh.model.setCoordinates(1, 0, 0, coords[3 * tag(0, 0)])
gmsh.model.setCoordinates(2, 1, 0, coords[3 * tag(N, 0)])
gmsh.model.setCoordinates(3, 1, 1, coords[3 * tag(N, N)])
gmsh.model.setCoordinates(4, 0, 1, coords[3 * tag(0, N)])
# Create 4 discrete bounding curves, with their boundary points:
for i in 0:3
gmsh.model.addDiscreteEntity(1, i + 1, [i + 1, i < 3 ? i + 2 : 1])
end
# Create one discrete surface, with its bounding curves:
gmsh.model.addDiscreteEntity(2, 1, [1, 2, -3, -4])
# Add all the nodes on the surface (for simplicity... see below):
gmsh.model.mesh.addNodes(2, 1, nodes, coords)
# Add point elements on the 4 points, line elements on the 4 curves, and
# triangle elements on the surface:
for i in 0:3
# Type 15 for point elements:
gmsh.model.mesh.addElementsByType(i + 1, 15, [], [pnt[i + 1]])
# Type 1 for 2-node line elements:
gmsh.model.mesh.addElementsByType(i + 1, 1, [], lin[i + 1])
end
# Type 2 for 3-node triangle elements:
gmsh.model.mesh.addElementsByType(1, 2, [], tris)
# Reclassify the nodes on the curves and the points (since we put them all on
# the surface before with `addNodes' for simplicity)
gmsh.model.mesh.reclassifyNodes()
# Create a geometry for the discrete curves and surfaces, so that we can remesh
# them later on:
gmsh.model.mesh.createGeometry()
# Note that for more complicated meshes, e.g. for on input unstructured STL
# mesh, we could use `classifySurfaces()' to automatically create the discrete
# entities and the topology; but we would then have to extract the boundaries
# afterwards.
# Create other build-in CAD entities to form one volume below the terrain
# surface. Beware that only built-in CAD entities can be hybrid, i.e. have
# discrete entities on their boundary: OpenCASCADE does not support this
# feature.
p1 = gmsh.model.geo.addPoint(0, 0, -0.5)
p2 = gmsh.model.geo.addPoint(1, 0, -0.5)
p3 = gmsh.model.geo.addPoint(1, 1, -0.5)
p4 = gmsh.model.geo.addPoint(0, 1, -0.5)
c1 = gmsh.model.geo.addLine(p1, p2)
c2 = gmsh.model.geo.addLine(p2, p3)
c3 = gmsh.model.geo.addLine(p3, p4)
c4 = gmsh.model.geo.addLine(p4, p1)
c10 = gmsh.model.geo.addLine(p1, 1)
c11 = gmsh.model.geo.addLine(p2, 2)
c12 = gmsh.model.geo.addLine(p3, 3)
c13 = gmsh.model.geo.addLine(p4, 4)
ll1 = gmsh.model.geo.addCurveLoop([c1, c2, c3, c4])
s1 = gmsh.model.geo.addPlaneSurface([ll1])
ll3 = gmsh.model.geo.addCurveLoop([c1, c11, -1, -c10])
s3 = gmsh.model.geo.addPlaneSurface([ll3])
ll4 = gmsh.model.geo.addCurveLoop([c2, c12, -2, -c11])
s4 = gmsh.model.geo.addPlaneSurface([ll4])
ll5 = gmsh.model.geo.addCurveLoop([c3, c13, 3, -c12])
s5 = gmsh.model.geo.addPlaneSurface([ll5])
ll6 = gmsh.model.geo.addCurveLoop([c4, c10, 4, -c13])
s6 = gmsh.model.geo.addPlaneSurface([ll6])
sl1 = gmsh.model.geo.addSurfaceLoop([s1, s3, s4, s5, s6, 1])
v1 = gmsh.model.geo.addVolume([sl1])
gmsh.model.geo.synchronize()
# Set this to true to build a fully hex mesh:
#transfinite = true
transfinite = false
transfiniteAuto = false
if transfinite
NN = 30
for c in gmsh.model.getEntities(1)
gmsh.model.mesh.setTransfiniteCurve(c[2], NN)
end
for s in gmsh.model.getEntities(2)
gmsh.model.mesh.setTransfiniteSurface(s[2])
gmsh.model.mesh.setRecombine(s[1], s[2])
gmsh.model.mesh.setSmoothing(s[1], s[2], 100)
end
gmsh.model.mesh.setTransfiniteVolume(v1)
elseif transfiniteAuto
gmsh.option.setNumber("Mesh.MeshSizeMin", 0.5)
gmsh.option.setNumber("Mesh.MeshSizeMax", 0.5)
# setTransfiniteAutomatic() uses the sizing constraints to set the number
# of points
gmsh.model.mesh.setTransfiniteAutomatic()
else
gmsh.option.setNumber("Mesh.MeshSizeMin", 0.05)
gmsh.option.setNumber("Mesh.MeshSizeMax", 0.05)
end
gmsh.model.mesh.generate(3)
gmsh.write("x2.msh")
# Launch the GUI to see the results:
if !("-nopopup" in ARGS)
gmsh.fltk.run()
end
gmsh.finalize()
|
