gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
meshGRegionDelaunayInsertion.h
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 #ifndef MESH_GREGION_DELAUNAY_INSERTION_H
7 #define MESH_GREGION_DELAUNAY_INSERTION_H
8 
9 #include <list>
10 #include <set>
11 #include <map>
12 #include <stack>
13 #include "MTetrahedron.h"
14 #include "Numeric.h"
15 #include "BackgroundMeshTools.h"
16 #include "qualityMeasures.h"
17 #include "robustPredicates.h"
18 
19 //#define GMSH_PRE_ALLOCATE_STRATEGY 1
20 class GRegion;
21 class GFace;
22 class GModel;
23 class splitQuadRecovery;
24 
25 double tetcircumcenter(double a[3], double b[3], double c[3], double d[3],
26  double circumcenter[3], double *xi, double *eta,
27  double *zeta);
28 
29 class MTet4Factory;
30 
31 // Memory usage for 1 million tets:
32 //
33 // * sizeof(MTet4) = 36 Bytes and sizeof(MTetrahedron) = 28 Bytes
34 // -> 64 MB
35 // * rb tree containing all pointers sorted with respect to tet
36 // radius: each bucket of the tree contains 4 pointers (16 Bytes)
37 // plus the data -> 20 MB
38 // * sizeof(MVertex) = 44 Bytes and there are about 200000 verts per
39 // million tet -> 9MB
40 // * vector of char lengths per vertex -> 1.6Mb
41 // * vectors in GEntities to store the element and vertex pointers
42 // -> 5Mb
43 //
44 // Grand total should thus be about 100 MB.
45 
46 class MTet4 {
47  friend class MTet4Factory;
48 
49 private:
50  bool deleted;
51  double circum_radius;
52  MTetrahedron *base;
53  MTet4 *neigh[4];
54  GRegion *gr;
55 
56 public:
57  static int radiusNorm; // 2 is euclidian norm, -1 is infinite norm
58  ~MTet4() {}
59  MTet4() : deleted(false), circum_radius(0.0), base(nullptr), gr(nullptr)
60  {
61  neigh[0] = neigh[1] = neigh[2] = neigh[3] = nullptr;
62  }
63  MTet4(MTetrahedron *t, double qual)
64  : deleted(false), circum_radius(qual), base(t), gr(nullptr)
65  {
66  neigh[0] = neigh[1] = neigh[2] = neigh[3] = nullptr;
67  }
68  MTet4(MTetrahedron *t, const qmTetrahedron::Measures &qm)
69  : deleted(false), base(t), gr(nullptr)
70  {
71  neigh[0] = neigh[1] = neigh[2] = neigh[3] = nullptr;
72  double vol;
73  circum_radius = qmTetrahedron::qm(t, qm, &vol);
74  }
75  void circumcenter(double *res)
76  {
77  MVertex *v0 = base->getVertex(0);
78  MVertex *v1 = base->getVertex(1);
79  MVertex *v2 = base->getVertex(2);
80  MVertex *v3 = base->getVertex(3);
81  double A[4] = {v0->x(), v0->y(), v0->z()};
82  double B[4] = {v1->x(), v1->y(), v1->z()};
83  double C[4] = {v2->x(), v2->y(), v2->z()};
84  double D[4] = {v3->x(), v3->y(), v3->z()};
85  tetcircumcenter(A, B, C, D, res, nullptr, nullptr, nullptr);
86  }
87 
88  void setup(MTetrahedron *t, std::vector<double> &sizes,
89  std::vector<double> &sizesBGM)
90  {
91  base = t;
92  neigh[0] = neigh[1] = neigh[2] = neigh[3] = nullptr;
93  double center[3];
94  circumcenter(center);
95  const double dx = base->getVertex(0)->x() - center[0];
96  const double dy = base->getVertex(0)->y() - center[1];
97  const double dz = base->getVertex(0)->z() - center[2];
98  circum_radius = std::sqrt(dx * dx + dy * dy + dz * dz);
99  /*
100  if (base->getVertex(0)->getIndex() >= sizes.size() ||
101  base->getVertex(1)->getIndex() >= sizes.size() ||
102  base->getVertex(2)->getIndex() >= sizes.size() ||
103  base->getVertex(3)->getIndex() >= sizes.size()){
104  printf("ERROR %d vs %d %d %d %d\n",sizes.size() ,
105  base->getVertex(0)->getIndex(),
106  base->getVertex(1)->getIndex(),
107  base->getVertex(2)->getIndex(),
108  base->getVertex(3)->getIndex());
109 
110  }
111  */
112  double lc1 = 0.25 * (sizes[base->getVertex(0)->getIndex()] +
113  sizes[base->getVertex(1)->getIndex()] +
114  sizes[base->getVertex(2)->getIndex()] +
115  sizes[base->getVertex(3)->getIndex()]);
116  double lcBGM = 0.25 * (sizesBGM[base->getVertex(0)->getIndex()] +
117  sizesBGM[base->getVertex(1)->getIndex()] +
118  sizesBGM[base->getVertex(2)->getIndex()] +
119  sizesBGM[base->getVertex(3)->getIndex()]);
120  double lc = Extend2dMeshIn3dVolumes() ? std::min(lc1, lcBGM) : lcBGM;
121  circum_radius /= lc;
122  deleted = false;
123  }
124 
125  void setup(MTetrahedron *t, std::vector<double> &sizes,
126  std::vector<double> &sizesBGM, double lcA, double lcB)
127  {
128  base = t;
129  neigh[0] = neigh[1] = neigh[2] = neigh[3] = nullptr;
130  double center[3];
131  circumcenter(center);
132  const double dx = base->getVertex(0)->x() - center[0];
133  const double dy = base->getVertex(0)->y() - center[1];
134  const double dz = base->getVertex(0)->z() - center[2];
135  circum_radius = std::sqrt(dx * dx + dy * dy + dz * dz);
136 
137  /*
138  if (base->getVertex(0)->getIndex() >= sizes.size() ||
139  base->getVertex(1)->getIndex() >= sizes.size() ||
140  base->getVertex(2)->getIndex() >= sizes.size()){
141  printf("ERROR %d vs %d %d %d %d\n",sizes.size() ,
142  base->getVertex(0)->getIndex(),
143  base->getVertex(1)->getIndex(),
144  base->getVertex(2)->getIndex(),
145  base->getVertex(3)->getIndex());
146 
147  }
148  */
149  double lc1 = 0.25 * (sizes[base->getVertex(0)->getIndex()] +
150  sizes[base->getVertex(1)->getIndex()] +
151  sizes[base->getVertex(2)->getIndex()] + lcA);
152  double lcBGM = 0.25 * (sizesBGM[base->getVertex(0)->getIndex()] +
153  sizesBGM[base->getVertex(1)->getIndex()] +
154  sizesBGM[base->getVertex(2)->getIndex()] + lcB);
155  double lc = Extend2dMeshIn3dVolumes() ? std::min(lc1, lcBGM) : lcBGM;
156  circum_radius /= lc;
157  deleted = false;
158  }
159 
160  GRegion *onWhat() const { return gr; }
161  void setOnWhat(GRegion *g) { gr = g; }
162  bool isDeleted() const { return deleted; }
163  void forceRadius(double r) { circum_radius = r; }
164  double getRadius() const { return circum_radius; }
165  double getQuality() const { return circum_radius; }
166  void setQuality(const double &q) { circum_radius = q; }
167  MTetrahedron *tet() const { return base; }
168  MTetrahedron *&tet() { return base; }
169  void setTet(MTetrahedron *t) { base = t; }
170  void setNeigh(int iN, MTet4 *n) { neigh[iN] = n; }
171  MTet4 *getNeigh(int iN) const { return neigh[iN]; }
172  int inCircumSphere(const double *p) const;
173  int inCircumSphere(double x, double y, double z) const
174  {
175  const double p[3] = {x, y, z};
176  return inCircumSphere(p);
177  }
178  int inCircumSphere(const MVertex *v) const
179  {
180  return inCircumSphere(v->x(), v->y(), v->z());
181  }
182  double getVolume() const
183  {
184  double pa[3] = {base->getVertex(0)->x(), base->getVertex(0)->y(),
185  base->getVertex(0)->z()};
186  double pb[3] = {base->getVertex(1)->x(), base->getVertex(1)->y(),
187  base->getVertex(1)->z()};
188  double pc[3] = {base->getVertex(2)->x(), base->getVertex(2)->y(),
189  base->getVertex(2)->z()};
190  double pd[3] = {base->getVertex(3)->x(), base->getVertex(3)->y(),
191  base->getVertex(3)->z()};
192  return std::abs(robustPredicates::orient3d(pa, pb, pc, pd)) / 6.0;
193  }
194  void setDeleted(bool const d) { deleted = d; }
195  bool assertNeigh() const
196  {
197  if(deleted) return true;
198  for(int i = 0; i < 4; i++)
199  if(neigh[i] && (neigh[i]->isNeigh(this) == false)) return false;
200  return true;
201  }
202  inline bool isNeigh(const MTet4 *t) const
203  {
204  for(int i = 0; i < 4; i++)
205  if(neigh[i] == t) return true;
206  return false;
207  }
208 };
209 
210 void connectTets(std::list<MTet4 *> &,
211  const std::set<MFace, MFaceLessThan> * = nullptr);
212 void connectTets(std::vector<MTet4 *> &,
213  const std::set<MFace, MFaceLessThan> * = nullptr);
214 void delaunayMeshIn3D(std::vector<MVertex *> &, std::vector<MTetrahedron *> &,
215  bool removeBox = false);
216 void insertVerticesInRegion(GRegion *gr, int maxIter,
217  double worstTetRadiusTarget, bool _classify = true,
218  splitQuadRecovery *sqr = nullptr);
219 void bowyerWatsonFrontalLayers(GRegion *gr, bool hex);
220 
221 struct compareTet4Ptr {
222  bool operator()(MTet4 const *const a, MTet4 const *const b) const
223  {
224  if(a->getRadius() > b->getRadius()) return true;
225  if(a->getRadius() < b->getRadius()) return false;
226  return a->tet()->getNum() < b->tet()->getNum();
227  }
228 };
229 
230 class MTet4Factory {
231 public:
232  typedef std::set<MTet4 *, compareTet4Ptr> container;
233  typedef container::iterator iterator;
234 
235 private:
236  container allTets;
237 #ifdef GMSH_PRE_ALLOCATE_STRATEGY
238  MTet4 *allSlots;
239  int s_last, s_alloc;
240  std::stack<MTet4 *> emptySlots;
241  inline MTet4 *getANewSlot()
242  {
243  if(s_last >= s_alloc) return 0;
244  MTet4 *t = &(allSlots[s_last]);
245  s_last++;
246  return t;
247  }
248  inline MTet4 *getAnEmptySlot()
249  {
250  if(!emptySlots.empty()) {
251  MTet4 *t = emptySlots.top();
252  emptySlots.pop();
253  return t;
254  }
255  return getANewSlot();
256  };
257 #endif
258 public:
259  MTet4Factory(int _size = 1000000)
260  {
261 #ifdef GMSH_PRE_ALLOCATE_STRATEGY
262  s_last = 0;
263  s_alloc = _size;
264  allSlots = new MTet4[s_alloc];
265 #endif
266  }
267  ~MTet4Factory()
268  {
269 #ifdef GMSH_PRE_ALLOCATE_STRATEGY
270  delete[] allSlots;
271 #endif
272  }
273  MTet4 *Create(MTetrahedron *t, std::vector<double> &sizes,
274  std::vector<double> &sizesBGM)
275  {
276 #ifdef GMSH_PRE_ALLOCATE_STRATEGY
277  MTet4 *t4 = getAnEmptySlot();
278 #else
279  MTet4 *t4 = new MTet4;
280 #endif
281  t4->setup(t, sizes, sizesBGM);
282  return t4;
283  }
284  MTet4 *Create(MTetrahedron *t, std::vector<double> &sizes,
285  std::vector<double> &sizesBGM, double lc1, double lc2)
286  {
287 #ifdef GMSH_PRE_ALLOCATE_STRATEGY
288  MTet4 *t4 = getAnEmptySlot();
289 #else
290  MTet4 *t4 = new MTet4;
291 #endif
292  t4->setup(t, sizes, sizesBGM, lc1, lc2);
293  return t4;
294  }
295 
296  void Free(MTet4 *t)
297  {
298  if(t->tet()) delete t->tet();
299  t->tet() = nullptr;
300 #ifdef GMSH_PRE_ALLOCATE_STRATEGY
301  emptySlots.push(t);
302  t->setDeleted(true);
303 #else
304  delete t;
305 #endif
306  }
307  void changeTetRadius(iterator it, double r)
308  {
309  MTet4 *t = *it;
310  allTets.erase(it);
311  t->forceRadius(r);
312  allTets.insert(t);
313  }
314  container &getAllTets() { return allTets; }
315 };
316 
317 void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm);
318 
319 #endif
MElement::getNum
virtual std::size_t getNum() const
Definition: MElement.h:68
D
#define D
Definition: DefaultOptions.h:24
qualityMeasures.h
sqr
static int sqr(int x)
Definition: gl2gif.cpp:266
MTet4::assertNeigh
bool assertNeigh() const
Definition: meshGRegionDelaunayInsertion.h:195
MTet4::inCircumSphere
int inCircumSphere(const MVertex *v) const
Definition: meshGRegionDelaunayInsertion.h:178
MTet4::base
MTetrahedron * base
Definition: meshGRegionDelaunayInsertion.h:52
compareTet4Ptr::operator()
bool operator()(MTet4 const *const a, MTet4 const *const b) const
Definition: meshGRegionDelaunayInsertion.h:222
MTetrahedron
Definition: MTetrahedron.h:34
GFace
Definition: GFace.h:33
MTet4Factory::Create
MTet4 * Create(MTetrahedron *t, std::vector< double > &sizes, std::vector< double > &sizesBGM, double lc1, double lc2)
Definition: meshGRegionDelaunayInsertion.h:284
insertVerticesInRegion
void insertVerticesInRegion(GRegion *gr, int maxIter, double worstTetRadiusTarget, bool _classify=true, splitQuadRecovery *sqr=nullptr)
Definition: meshGRegionDelaunayInsertion.cpp:1217
robustPredicates.h
c
static double c(int i, int j, fullMatrix< double > &CA, const std::vector< SPoint3 > &P, const std::vector< SPoint3 > &Q)
Definition: discreteFrechetDistance.cpp:15
MTet4::tet
MTetrahedron *& tet()
Definition: meshGRegionDelaunayInsertion.h:168
MVertex
Definition: MVertex.h:24
MVertex::z
double z() const
Definition: MVertex.h:62
MTet4::setup
void setup(MTetrahedron *t, std::vector< double > &sizes, std::vector< double > &sizesBGM, double lcA, double lcB)
Definition: meshGRegionDelaunayInsertion.h:125
tetcircumcenter
double tetcircumcenter(double a[3], double b[3], double c[3], double d[3], double circumcenter[3], double *xi, double *eta, double *zeta)
Definition: meshGRegionDelaunayInsertion.cpp:1058
MTet4::setNeigh
void setNeigh(int iN, MTet4 *n)
Definition: meshGRegionDelaunayInsertion.h:170
MTet4::setQuality
void setQuality(const double &q)
Definition: meshGRegionDelaunayInsertion.h:166
MTet4Factory::container
std::set< MTet4 *, compareTet4Ptr > container
Definition: meshGRegionDelaunayInsertion.h:232
MTet4Factory::iterator
container::iterator iterator
Definition: meshGRegionDelaunayInsertion.h:233
MTet4::setup
void setup(MTetrahedron *t, std::vector< double > &sizes, std::vector< double > &sizesBGM)
Definition: meshGRegionDelaunayInsertion.h:88
MTet4Factory::Free
void Free(MTet4 *t)
Definition: meshGRegionDelaunayInsertion.h:296
MTet4
Definition: meshGRegionDelaunayInsertion.h:46
MTet4::getVolume
double getVolume() const
Definition: meshGRegionDelaunayInsertion.h:182
delaunayMeshIn3D
void delaunayMeshIn3D(std::vector< MVertex * > &, std::vector< MTetrahedron * > &, bool removeBox=false)
Definition: meshGRegionDelaunayInsertion.cpp:1559
MTet4::tet
MTetrahedron * tet() const
Definition: meshGRegionDelaunayInsertion.h:167
MTet4::inCircumSphere
int inCircumSphere(double x, double y, double z) const
Definition: meshGRegionDelaunayInsertion.h:173
MTetrahedron::getVertex
virtual MVertex * getVertex(int num)
Definition: MTetrahedron.h:67
MTet4::onWhat
GRegion * onWhat() const
Definition: meshGRegionDelaunayInsertion.h:160
Numeric.h
GModel
Definition: GModel.h:44
qmTetrahedron::qm
static double qm(MTetrahedron *t, const Measures &cr, double *volume=nullptr)
Definition: qualityMeasures.cpp:616
MTet4::forceRadius
void forceRadius(double r)
Definition: meshGRegionDelaunayInsertion.h:163
Extend2dMeshIn3dVolumes
bool Extend2dMeshIn3dVolumes()
Definition: BackgroundMeshTools.cpp:345
MTet4::MTet4
MTet4(MTetrahedron *t, double qual)
Definition: meshGRegionDelaunayInsertion.h:63
MTet4Factory::getAllTets
container & getAllTets()
Definition: meshGRegionDelaunayInsertion.h:314
splitQuadRecovery
Definition: meshGRegion.h:72
GRegion
Definition: GRegion.h:28
MTet4Factory::~MTet4Factory
~MTet4Factory()
Definition: meshGRegionDelaunayInsertion.h:267
MTet4::MTet4
MTet4()
Definition: meshGRegionDelaunayInsertion.h:59
MTet4::deleted
bool deleted
Definition: meshGRegionDelaunayInsertion.h:50
MVertex::getIndex
long int getIndex() const
Definition: MVertex.h:93
MTet4::getQuality
double getQuality() const
Definition: meshGRegionDelaunayInsertion.h:165
MTetrahedron.h
MTet4Factory
Definition: meshGRegionDelaunayInsertion.h:230
MTet4::MTet4
MTet4(MTetrahedron *t, const qmTetrahedron::Measures &qm)
Definition: meshGRegionDelaunayInsertion.h:68
z
const double z
Definition: GaussQuadratureQuad.cpp:56
MTet4Factory::MTet4Factory
MTet4Factory(int _size=1000000)
Definition: meshGRegionDelaunayInsertion.h:259
MTet4::setOnWhat
void setOnWhat(GRegion *g)
Definition: meshGRegionDelaunayInsertion.h:161
MTet4::circum_radius
double circum_radius
Definition: meshGRegionDelaunayInsertion.h:51
MTet4::setTet
void setTet(MTetrahedron *t)
Definition: meshGRegionDelaunayInsertion.h:169
optimizeMesh
void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
Definition: meshGRegionDelaunayInsertion.cpp:876
MTet4::circumcenter
void circumcenter(double *res)
Definition: meshGRegionDelaunayInsertion.h:75
MTet4::~MTet4
~MTet4()
Definition: meshGRegionDelaunayInsertion.h:58
BackgroundMeshTools.h
robustPredicates::orient3d
REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
Definition: robustPredicates.cpp:2321
connectTets
void connectTets(std::list< MTet4 * > &, const std::set< MFace, MFaceLessThan > *=nullptr)
Definition: meshGRegionDelaunayInsertion.cpp:277
MTet4::inCircumSphere
int inCircumSphere(const double *p) const
Definition: meshGRegionDelaunayInsertion.cpp:150
MVertex::y
double y() const
Definition: MVertex.h:61
MTet4Factory::allTets
container allTets
Definition: meshGRegionDelaunayInsertion.h:236
MTet4Factory::Create
MTet4 * Create(MTetrahedron *t, std::vector< double > &sizes, std::vector< double > &sizesBGM)
Definition: meshGRegionDelaunayInsertion.h:273
qmTetrahedron::Measures
Measures
Definition: qualityMeasures.h:68
MTet4::getRadius
double getRadius() const
Definition: meshGRegionDelaunayInsertion.h:164
MTet4::gr
GRegion * gr
Definition: meshGRegionDelaunayInsertion.h:54
bowyerWatsonFrontalLayers
void bowyerWatsonFrontalLayers(GRegion *gr, bool hex)
MTet4::isNeigh
bool isNeigh(const MTet4 *t) const
Definition: meshGRegionDelaunayInsertion.h:202
MTet4Factory::changeTetRadius
void changeTetRadius(iterator it, double r)
Definition: meshGRegionDelaunayInsertion.h:307
MTet4::neigh
MTet4 * neigh[4]
Definition: meshGRegionDelaunayInsertion.h:53
MTet4::getNeigh
MTet4 * getNeigh(int iN) const
Definition: meshGRegionDelaunayInsertion.h:171
MTet4::radiusNorm
static int radiusNorm
Definition: meshGRegionDelaunayInsertion.h:57
MTet4::isDeleted
bool isDeleted() const
Definition: meshGRegionDelaunayInsertion.h:162
MTet4::setDeleted
void setDeleted(bool const d)
Definition: meshGRegionDelaunayInsertion.h:194
MVertex::x
double x() const
Definition: MVertex.h:60
compareTet4Ptr
Definition: meshGRegionDelaunayInsertion.h:221