gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
meshGRegionNetgen.cpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 
6 #include <stdlib.h>
7 #include <vector>
8 #include "GmshConfig.h"
9 #include "GmshMessage.h"
10 #include "meshGRegion.h"
11 #include "GModel.h"
12 #include "GRegion.h"
13 #include "GFace.h"
14 #include "MTriangle.h"
15 #include "MTetrahedron.h"
16 #include "ExtrudeParams.h"
17 #include "BackgroundMeshTools.h"
18 #include "Context.h"
19 
20 #if defined(HAVE_NETGEN)
21 
22 namespace nglib {
23 #include "nglib_gmsh.h"
24 }
25 using namespace nglib;
26 
27 static void getAllBoundingVertices(
28  GRegion *gr, std::set<MVertex *, MVertexPtrLessThan> &allBoundingVertices)
29 {
30  std::vector<GFace *> faces = gr->faces();
31  auto it = faces.begin();
32 
33  while(it != faces.end()) {
34  GFace *gf = (*it);
35  for(std::size_t i = 0; i < gf->triangles.size(); i++) {
36  MTriangle *t = gf->triangles[i];
37  for(int k = 0; k < 3; k++)
38  if(allBoundingVertices.find(t->getVertex(k)) ==
39  allBoundingVertices.end())
40  allBoundingVertices.insert(t->getVertex(k));
41  }
42  ++it;
43  }
44 }
45 
46 static Ng_Mesh *buildNetgenStructure(GRegion *gr, bool importVolumeMesh,
47  std::vector<MVertex *> &numberedV)
48 {
49  Ng_Init();
50  Ng_Mesh *ngmesh = Ng_NewMesh();
51 
52  std::set<MVertex *, MVertexPtrLessThan> allBoundingVertices;
53  getAllBoundingVertices(gr, allBoundingVertices);
54 
55  auto itv = allBoundingVertices.begin();
56  int I = 1;
57  while(itv != allBoundingVertices.end()) {
58  double tmp[3];
59  tmp[0] = (*itv)->x();
60  tmp[1] = (*itv)->y();
61  tmp[2] = (*itv)->z();
62  (*itv)->setIndex(I++);
63  numberedV.push_back(*itv);
64  Ng_AddPoint(ngmesh, tmp);
65  ++itv;
66  }
67 
68  if(importVolumeMesh) {
69  for(std::size_t i = 0; i < gr->mesh_vertices.size(); i++) {
70  double tmp[3];
71  tmp[0] = gr->mesh_vertices[i]->x();
72  tmp[1] = gr->mesh_vertices[i]->y();
73  tmp[2] = gr->mesh_vertices[i]->z();
74  gr->mesh_vertices[i]->setIndex(I++);
75  Ng_AddPoint(ngmesh, tmp);
76  }
77  }
78  std::vector<GFace *> faces = gr->faces();
79  auto it = faces.begin();
80  while(it != faces.end()) {
81  GFace *gf = (*it);
82  for(std::size_t i = 0; i < gf->triangles.size(); i++) {
83  MTriangle *t = gf->triangles[i];
84  int tmp[3];
85  tmp[0] = t->getVertex(0)->getIndex();
86  tmp[1] = t->getVertex(1)->getIndex();
87  tmp[2] = t->getVertex(2)->getIndex();
88  Ng_AddSurfaceElement(ngmesh, NG_TRIG, tmp);
89  }
90  ++it;
91  }
92 
93  if(importVolumeMesh) {
94  for(std::size_t i = 0; i < gr->tetrahedra.size(); i++) {
95  MTetrahedron *t = gr->tetrahedra[i];
96  // Netgen expects tet with negative volume
97  if(t->getVolumeSign() > 0) t->reverse();
98  int tmp[4];
99  tmp[0] = t->getVertex(0)->getIndex();
100  tmp[1] = t->getVertex(1)->getIndex();
101  tmp[2] = t->getVertex(2)->getIndex();
102  tmp[3] = t->getVertex(3)->getIndex();
103  Ng_AddVolumeElement(ngmesh, NG_TET, tmp);
104  }
105  }
106 
107  return ngmesh;
108 }
109 
110 static void TransferVolumeMesh(GRegion *gr, Ng_Mesh *ngmesh,
111  std::vector<MVertex *> &numberedV)
112 {
113  // get total number of vertices and elements of Netgen's mesh
114  int nbv = Ng_GetNP(ngmesh);
115  int nbe = Ng_GetNE(ngmesh);
116  if(!nbv || !nbe) return;
117 
118  int nbpts = numberedV.size();
119 
120  // record which vertices are actually used by tetrahedra (Netgen can return
121  // additional "isolated" vertices)
122  std::set<int> used;
123  for(int i = 0; i < nbe; i++) {
124  int tmp[4];
125  Ng_GetVolumeElement(ngmesh, i + 1, tmp);
126  for(int j = 0; j < 4; j++) {
127  if(tmp[j] > nbpts) used.insert(tmp[j]);
128  }
129  }
130 
131  // create new vertices
132  for(int i = nbpts; i < nbv; i++) {
133  double tmp[3];
134  Ng_GetPoint(ngmesh, i + 1, tmp);
135  if(used.count(i + 1)) {
136  MVertex *v = new MVertex(tmp[0], tmp[1], tmp[2], gr);
137  numberedV.push_back(v);
138  gr->mesh_vertices.push_back(v);
139  }
140  else{
141  numberedV.push_back(nullptr);
142  }
143  }
144 
145  // create new tets
146  for(int i = 0; i < nbe; i++) {
147  int tmp[4];
148  Ng_GetVolumeElement(ngmesh, i + 1, tmp);
149  MVertex *v[4];
150  for(int j = 0; j < 4; j++)
151  v[j] = numberedV[tmp[j] - 1];
152  if(v[0] && v[1] && v[2] && v[3])
153  gr->tetrahedra.push_back(new MTetrahedron(v[0], v[1], v[2], v[3]));
154  else
155  Msg::Error("Tetrahedron with unknown node - should never be here!");
156  }
157 }
158 
159 // X_1 (1-u-v) + X_2 u + X_3 v = P_x + t N_x
160 // Y_1 (1-u-v) + Y_2 u + Y_3 v = P_y + t N_y
161 // Z_1 (1-u-v) + Z_2 u + Z_3 v = P_z + t N_z
162 
163 static int intersectLineTriangle(double X[3], double Y[3], double Z[3],
164  double P[3], double N[3],
165  const double eps_prec)
166 {
167  double mat[3][3], det;
168  double b[3], res[3];
169 
170  mat[0][0] = X[1] - X[0];
171  mat[0][1] = X[2] - X[0];
172  mat[0][2] = N[0];
173 
174  mat[1][0] = Y[1] - Y[0];
175  mat[1][1] = Y[2] - Y[0];
176  mat[1][2] = N[1];
177 
178  mat[2][0] = Z[1] - Z[0];
179  mat[2][1] = Z[2] - Z[0];
180  mat[2][2] = N[2];
181 
182  b[0] = P[0] - X[0];
183  b[1] = P[1] - Y[0];
184  b[2] = P[2] - Z[0];
185 
186  if(!sys3x3_with_tol(mat, b, res, &det)) { return 0; }
187  // printf("coucou %g %g %g\n",res[0],res[1],res[2]);
188  if(res[0] >= eps_prec && res[0] <= 1.0 - eps_prec && res[1] >= eps_prec &&
189  res[1] <= 1.0 - eps_prec && 1 - res[0] - res[1] >= eps_prec &&
190  1 - res[0] - res[1] <= 1.0 - eps_prec) {
191  // the line clearly intersects the triangle
192  return (res[2] > 0) ? 1 : 0;
193  }
194  else if(res[0] < -eps_prec || res[0] > 1.0 + eps_prec || res[1] < -eps_prec ||
195  res[1] > 1.0 + eps_prec || 1 - res[0] - res[1] < -eps_prec ||
196  1 - res[0] - res[1] > 1.0 + eps_prec) {
197  // the line clearly does NOT intersect the triangle
198  return 0;
199  }
200  else {
201  // printf("non robust stuff\n");
202  // the intersection is not robust, try another triangle
203  return -10000;
204  }
205 }
206 
207 static void setRand(double r[6])
208 {
209  for(int i = 0; i < 6; i++)
210  r[i] = 0.0001 * ((double)rand() / (double)RAND_MAX);
211 }
212 
213 static void meshNormalsPointOutOfTheRegion(GRegion *gr)
214 {
215  std::vector<GFace *> faces = gr->faces();
216  auto it = faces.begin();
217 
218  // perform intersection check in normalized coordinates
219  SBoundingBox3d bbox = gr->bounds();
220  double scaling = norm(SVector3(bbox.max(), bbox.min()));
221  if(!scaling) {
222  Msg::Warning("Bad scaling in meshNormalsPointOutOfTheRegion");
223  scaling = 1.;
224  }
225 
226  double rrr[6];
227  setRand(rrr);
228 
229  while(it != faces.end()) {
230  GFace *gf = (*it);
231  int nb_intersect = 0;
232  for(std::size_t i = 0; i < gf->triangles.size(); i++) {
233  MTriangle *t = gf->triangles[i];
234  double X[3] = {t->getVertex(0)->x(), t->getVertex(1)->x(),
235  t->getVertex(2)->x()};
236  double Y[3] = {t->getVertex(0)->y(), t->getVertex(1)->y(),
237  t->getVertex(2)->y()};
238  double Z[3] = {t->getVertex(0)->z(), t->getVertex(1)->z(),
239  t->getVertex(2)->z()};
240  for(int j = 0; j < 3; j++) {
241  X[j] /= scaling;
242  Y[j] /= scaling;
243  Z[j] /= scaling;
244  }
245 
246  double P[3] = {(X[0] + X[1] + X[2]) / 3., (Y[0] + Y[1] + Y[2]) / 3.,
247  (Z[0] + Z[1] + Z[2]) / 3.};
248  double v1[3] = {X[0] - X[1], Y[0] - Y[1], Z[0] - Z[1]};
249  double v2[3] = {X[2] - X[1], Y[2] - Y[1], Z[2] - Z[1]};
250  double N[3];
251  prodve(v1, v2, N);
252  norme(v1);
253  norme(v2);
254  norme(N);
255  N[0] += rrr[0] * v1[0] + rrr[1] * v2[0];
256  N[1] += rrr[2] * v1[1] + rrr[3] * v2[1];
257  N[2] += rrr[4] * v1[2] + rrr[5] * v2[2];
258  norme(N);
259  auto it_b = faces.begin();
260  while(it_b != faces.end()) {
261  GFace *gf_b = (*it_b);
262  for(std::size_t i_b = 0; i_b < gf_b->triangles.size(); i_b++) {
263  MTriangle *t_b = gf_b->triangles[i_b];
264  if(t_b != t) {
265  double X_b[3] = {t_b->getVertex(0)->x(), t_b->getVertex(1)->x(),
266  t_b->getVertex(2)->x()};
267  double Y_b[3] = {t_b->getVertex(0)->y(), t_b->getVertex(1)->y(),
268  t_b->getVertex(2)->y()};
269  double Z_b[3] = {t_b->getVertex(0)->z(), t_b->getVertex(1)->z(),
270  t_b->getVertex(2)->z()};
271  for(int j = 0; j < 3; j++) {
272  X_b[j] /= scaling;
273  Y_b[j] /= scaling;
274  Z_b[j] /= scaling;
275  }
276  int inters = intersectLineTriangle(X_b, Y_b, Z_b, P, N, 1.e-9);
277  nb_intersect += inters;
278  }
279  }
280  ++it_b;
281  }
282  Msg::Info("Region %d Face %d, %d intersect", gr->tag(), gf->tag(),
283  nb_intersect);
284  if(nb_intersect >= 0)
285  break; // negative value means intersection is not "robust"
286  }
287 
288  if(nb_intersect < 0) { setRand(rrr); }
289  else {
290  if(nb_intersect % 2 == 1) {
291  // odd nb of intersections: the normal points inside the region
292  for(std::size_t i = 0; i < gf->triangles.size(); i++) {
293  gf->triangles[i]->reverse();
294  }
295  }
296  ++it;
297  }
298  }
299 
300  // FILE *fp = Fopen("debug.pos", "w");
301  // if(fp){
302  // fprintf(fp, "View \"debug\" {\n");
303  // for(auto it = faces.begin(); it != faces.end(); it++)
304  // for(std::size_t i = 0; i < (*it)->triangles.size(); i++)
305  // (*it)->triangles[i]->writePOS(fp, 1., (*it)->tag());
306  // fprintf(fp, "};\n");
307  // fclose(fp);
308  // }
309 }
310 
311 #endif
312 
313 void meshGRegionNetgen(GRegion *gr)
314 {
315 #if !defined(HAVE_NETGEN)
316  Msg::Error("Frontal algorithm requires Netgen");
317 #else
318  // sanity check for frontal algo
319  std::vector<GFace *> faces = gr->faces();
320  for(auto it = faces.begin(); it != faces.end(); it++) {
321  if((*it)->quadrangles.size()) {
322  Msg::Error(
323  "Cannot use frontal 3D algorithm with quadrangles on boundary");
324  return;
325  }
326  }
327 
328  Msg::Info("Meshing volume %d (Frontal)", gr->tag());
329  // orient the triangles of with respect to this region
330  meshNormalsPointOutOfTheRegion(gr);
331  std::vector<MVertex *> numberedV;
332  Ng_Mesh *ngmesh = buildNetgenStructure(gr, false, numberedV);
333  SPoint3 pt = gr->bounds().center();
334  double lc = BGM_MeshSize(gr, 0, 0, pt.x(), pt.y(), pt.z());
335  Ng_GenerateVolumeMesh(ngmesh, lc);
336  TransferVolumeMesh(gr, ngmesh, numberedV);
337  Ng_DeleteMesh(ngmesh);
338  Ng_Exit();
339 #endif
340 }
341 
342 void optimizeMeshGRegionNetgen::operator()(GRegion *gr, bool always)
343 {
344  gr->model()->setCurrentMeshEntity(gr);
345 
346  if(!always && gr->geomType() == GEntity::DiscreteVolume) return;
347 
348  // don't optimize transfinite or extruded meshes
349  if(gr->meshAttributes.method == MESH_TRANSFINITE) return;
350  ExtrudeParams *ep = gr->meshAttributes.extrude;
351  if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == EXTRUDED_ENTITY) return;
352 
353  if(gr->prisms.size() || gr->hexahedra.size() || gr->pyramids.size()) {
354  Msg::Info("Skipping Netgen optimizer for hybrid mesh");
355  return;
356  }
357 
358 #if !defined(HAVE_NETGEN)
359  Msg::Error("Netgen optimizer is not compiled in this version of Gmsh");
360 #else
361  Msg::Info("Optimizing volume %d", gr->tag());
362  // import mesh into netgen, including volume tets
363  std::vector<MVertex *> numberedV;
364  Ng_Mesh *ngmesh = buildNetgenStructure(gr, true, numberedV);
365  // delete volume vertices and tets
366  deMeshGRegion dem;
367  dem(gr);
368  // optimize mesh
369  SPoint3 pt = gr->bounds().center();
370  double lc = BGM_MeshSize(gr, 0, 0, pt.x(), pt.y(), pt.z());
371  Ng_OptimizeVolumeMesh(ngmesh, lc);
372  TransferVolumeMesh(gr, ngmesh, numberedV);
373  Ng_DeleteMesh(ngmesh);
374  Ng_Exit();
375 #endif
376 }
MTriangle.h
GRegion::method
char method
Definition: GRegion.h:146
GFace.h
meshGRegionNetgen
void meshGRegionNetgen(GRegion *gr)
Definition: meshGRegionNetgen.cpp:313
MTetrahedron
Definition: MTetrahedron.h:34
GFace
Definition: GFace.h:33
GFace::triangles
std::vector< MTriangle * > triangles
Definition: GFace.h:428
Msg::Info
static void Info(const char *fmt,...)
Definition: GmshMessage.cpp:587
GEntity::model
GModel * model() const
Definition: GEntity.h:277
MESH_TRANSFINITE
#define MESH_TRANSFINITE
Definition: GmshDefines.h:259
MVertex
Definition: MVertex.h:24
Msg::Warning
static void Warning(const char *fmt,...)
Definition: GmshMessage.cpp:543
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
MVertex::z
double z() const
Definition: MVertex.h:62
SPoint3
Definition: SPoint3.h:14
SBoundingBox3d::min
SPoint3 min() const
Definition: SBoundingBox3d.h:90
SVector3
Definition: SVector3.h:16
MTetrahedron::reverse
virtual void reverse()
Definition: MTetrahedron.h:118
MTriangle::getVertex
virtual MVertex * getVertex(int num)
Definition: MTriangle.h:62
GmshMessage.h
ExtrudeParams::geo
struct ExtrudeParams::@15 geo
deMeshGRegion
Definition: meshGRegion.h:48
SBoundingBox3d::center
SPoint3 center() const
Definition: SBoundingBox3d.h:92
GRegion::faces
virtual std::vector< GFace * > faces() const
Definition: GRegion.h:64
SPoint3::x
double x(void) const
Definition: SPoint3.h:125
GRegion::hexahedra
std::vector< MHexahedron * > hexahedra
Definition: GRegion.h:164
GRegion::extrude
ExtrudeParams * extrude
Definition: GRegion.h:148
GRegion.h
ExtrudeParams::mesh
struct ExtrudeParams::@14 mesh
GEntity::mesh_vertices
std::vector< MVertex * > mesh_vertices
Definition: GEntity.h:56
EXTRUDED_ENTITY
#define EXTRUDED_ENTITY
Definition: ExtrudeParams.h:17
optimizeMeshGRegionNetgen::operator()
void operator()(GRegion *, bool always=false)
Definition: meshGRegionNetgen.cpp:342
MTetrahedron::getVertex
virtual MVertex * getVertex(int num)
Definition: MTetrahedron.h:67
GEntity::DiscreteVolume
@ DiscreteVolume
Definition: GEntity.h:120
norm
void norm(const double *vec, double *norm)
Definition: gmshLevelset.cpp:202
BGM_MeshSize
double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y, double Z)
Definition: BackgroundMeshTools.cpp:255
ExtrudeParams.h
GEntity::geomType
virtual GeomType geomType() const
Definition: GEntity.h:238
SPoint3::y
double y(void) const
Definition: SPoint3.h:127
ExtrudeParams::Mode
int Mode
Definition: ExtrudeParams.h:49
GEntity::tag
int tag() const
Definition: GEntity.h:280
MTetrahedron::getVolumeSign
virtual int getVolumeSign()
Definition: MTetrahedron.h:137
prodve
void prodve(double a[3], double b[3], double c[3])
Definition: Numeric.h:105
intersectLineTriangle
static int intersectLineTriangle(double X[3], double Y[3], double Z[3], double P[3], double N[3], const double eps_prec)
Definition: GRegion.cpp:797
GRegion
Definition: GRegion.h:28
ExtrudeParams::ExtrudeMesh
bool ExtrudeMesh
Definition: ExtrudeParams.h:36
GRegion::prisms
std::vector< MPrism * > prisms
Definition: GRegion.h:165
MTriangle
Definition: MTriangle.h:26
GRegion::pyramids
std::vector< MPyramid * > pyramids
Definition: GRegion.h:166
MVertex::getIndex
long int getIndex() const
Definition: MVertex.h:93
Context.h
MTetrahedron.h
GRegion::meshAttributes
struct GRegion::@20 meshAttributes
setRand
static void setRand(double r[6])
Definition: GRegion.cpp:787
SPoint3::z
double z(void) const
Definition: SPoint3.h:129
meshGRegion.h
BackgroundMeshTools.h
GModel.h
MVertex::y
double y() const
Definition: MVertex.h:61
norme
double norme(double a[3])
Definition: Numeric.h:123
SBoundingBox3d::max
SPoint3 max() const
Definition: SBoundingBox3d.h:91
ExtrudeParams
Definition: ExtrudeParams.h:26
SBoundingBox3d
Definition: SBoundingBox3d.h:21
GModel::setCurrentMeshEntity
void setCurrentMeshEntity(GEntity *e)
Definition: GModel.cpp:2202
GRegion::bounds
virtual SBoundingBox3d bounds(bool fast=false)
Definition: GRegion.cpp:179
sys3x3_with_tol
int sys3x3_with_tol(double mat[3][3], double b[3], double res[3], double *det)
Definition: Numeric.cpp:177
GRegion::tetrahedra
std::vector< MTetrahedron * > tetrahedra
Definition: GRegion.h:163
MVertex::x
double x() const
Definition: MVertex.h:60
faces
static int faces[4][3]
Definition: meshGRegionDelaunayInsertion.cpp:165