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
|
// -----------------------------------------------------------------------------
//
// Gmsh C++ extended tutorial 2
//
// Mesh import, discrete entities, hybrid models, terrain meshing
//
// -----------------------------------------------------------------------------
#include <set>
#include <iostream>
#include <gmsh.h>
// 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.py') 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.
int main(int argc, char **argv)
{
gmsh::initialize();
gmsh::model::add("x2");
// We will create the terrain surface mesh from N x N input data points:
int N = 100;
// Helper function to return a node tag given two indices i and j:
auto tag = [N](std::size_t i, std::size_t j) { return (N + 1) * i + j + 1; };
// The x, y, z coordinates of all the points:
std::vector<double> coords;
// The tags of the corresponding nodes:
std::vector<std::size_t> nodes;
// The connectivities of the triangle elements (3 node tags per triangle) on
// the terrain surface:
std::vector<std::size_t> tris;
// The connectivities of the line elements on the 4 boundaries (2 node tags
// for each line element):
std::vector<std::size_t> lin[4];
// The connectivities of the point elements on the 4 corners (1 node tag for
// each point element):
std::size_t pnt[4] = {tag(0, 0), tag(N, 0), tag(N, N), tag(0, N)};
for(std::size_t i = 0; i < N + 1; i++) {
for(std::size_t j = 0; j < N + 1; j++) {
nodes.push_back(tag(i, j));
coords.insert(coords.end(), {(double)i / N, (double)j / N,
0.05 * sin(10 * (double)(i + j) / N)});
if(i > 0 && j > 0) {
tris.insert(tris.end(),
{tag(i - 1, j - 1), tag(i, j - 1), tag(i - 1, j)});
tris.insert(tris.end(), {tag(i, j - 1), tag(i, j), tag(i - 1, j)});
}
if((i == 0 || i == N) && j > 0) {
std::size_t s = (i == 0) ? 3 : 1;
lin[s].insert(lin[s].end(), {tag(i, j - 1), tag(i, j)});
}
if((j == 0 || j == N) && i > 0) {
std::size_t s = (j == 0) ? 0 : 2;
lin[s].insert(lin[s].end(), {tag(i - 1, j), tag(i, j)});
}
}
}
// Create 4 discrete points for the 4 corners of the terrain surface:
for(int i = 0; i < 4; i++) gmsh::model::addDiscreteEntity(0, i + 1);
gmsh::model::setCoordinates(1, 0, 0, coords[3 * tag(0, 0) - 1]);
gmsh::model::setCoordinates(2, 1, 0, coords[3 * tag(N, 0) - 1]);
gmsh::model::setCoordinates(3, 1, 1, coords[3 * tag(N, N) - 1]);
gmsh::model::setCoordinates(4, 0, 1, coords[3 * tag(0, N) - 1]);
// Create 4 discrete bounding curves, with their boundary points:
for(int i = 0; i < 4; i++)
gmsh::model::addDiscreteEntity(1, i + 1, {i + 1, (i < 3) ? (i + 2) : 1});
// 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(int i = 0; i < 4; i++) {
// Type 15 for point elements:
gmsh::model::mesh::addElementsByType(i + 1, 15, {}, {pnt[i]});
// Type 1 for 2-node line elements:
gmsh::model::mesh::addElementsByType(i + 1, 1, {}, lin[i]);
}
// 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 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.
int p1 = gmsh::model::geo::addPoint(0, 0, -0.5);
int p2 = gmsh::model::geo::addPoint(1, 0, -0.5);
int p3 = gmsh::model::geo::addPoint(1, 1, -0.5);
int p4 = gmsh::model::geo::addPoint(0, 1, -0.5);
int c1 = gmsh::model::geo::addLine(p1, p2);
int c2 = gmsh::model::geo::addLine(p2, p3);
int c3 = gmsh::model::geo::addLine(p3, p4);
int c4 = gmsh::model::geo::addLine(p4, p1);
int c10 = gmsh::model::geo::addLine(p1, 1);
int c11 = gmsh::model::geo::addLine(p2, 2);
int c12 = gmsh::model::geo::addLine(p3, 3);
int c13 = gmsh::model::geo::addLine(p4, 4);
int ll1 = gmsh::model::geo::addCurveLoop({c1, c2, c3, c4});
int s1 = gmsh::model::geo::addPlaneSurface({ll1});
int ll3 = gmsh::model::geo::addCurveLoop({c1, c11, -1, -c10});
int s3 = gmsh::model::geo::addPlaneSurface({ll3});
int ll4 = gmsh::model::geo::addCurveLoop({c2, c12, -2, -c11});
int s4 = gmsh::model::geo::addPlaneSurface({ll4});
int ll5 = gmsh::model::geo::addCurveLoop({c3, c13, 3, -c12});
int s5 = gmsh::model::geo::addPlaneSurface({ll5});
int ll6 = gmsh::model::geo::addCurveLoop({c4, c10, 4, -c13});
int s6 = gmsh::model::geo::addPlaneSurface({ll6});
int sl1 = gmsh::model::geo::addSurfaceLoop({s1, s3, s4, s5, s6, 1});
int v1 = gmsh::model::geo::addVolume({sl1});
gmsh::model::geo::synchronize();
// Set this to true to build a fully hex mesh:
bool transfinite = false;
bool transfiniteAuto = false;
if(transfinite) {
int NN = 30;
std::vector<std::pair<int, int> > tmp;
gmsh::model::getEntities(tmp, 1);
for(auto c : tmp) { gmsh::model::mesh::setTransfiniteCurve(c.second, NN); }
gmsh::model::getEntities(tmp, 2);
for(auto s : tmp) {
gmsh::model::mesh::setTransfiniteSurface(s.second);
gmsh::model::mesh::setRecombine(s.first, s.second);
gmsh::model::mesh::setSmoothing(s.first, s.second, 100);
}
gmsh::model::mesh::setTransfiniteVolume(v1);
}
else if(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);
}
gmsh::model::mesh::generate(3);
gmsh::write("x2.msh");
// Launch the GUI to see the results:
std::set<std::string> args(argv, argv + argc);
if(!args.count("-nopopup")) gmsh::fltk::run();
gmsh::finalize();
return 0;
}
|
