gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
meshGRegionLocalMeshMod.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 "meshGRegionLocalMeshMod.h"
7 #include "GEntity.h"
8 #include "GRegion.h"
9 #include "GmshMessage.h"
10 #include "Numeric.h"
11 #include "MHexahedron.h"
12 #include "MPrism.h"
13 #include "MPyramid.h"
14 
15 typedef struct {
16  int nbr_triangles; // number of different triangles
17  int (*triangles)[3]; // triangles array
18  int nbr_trianguls; // number of different triangulations
19  int nbr_triangles_2; // number of triangles / triangulation
20  int (*trianguls)[5]; // retriangulations array
21 } SwapPattern;
22 
23 static int edges[6][2] = {{0, 1}, {0, 2}, {0, 3}, {1, 2}, {1, 3}, {2, 3}};
24 static int efaces[6][2] = {{0, 2}, {0, 1}, {1, 2}, {0, 3}, {2, 3}, {1, 3}};
25 // static int enofaces[6][2] = {{1,3},{2,3},{0,3},{1,2},{0,1},{0,2}};
26 // static int facesXedges[4][3] = {{0,1,3},{1,2,5},{0,2,4},{3,4,5}};
27 static int faces[4][3] = {{0, 1, 2}, {0, 2, 3}, {0, 1, 3}, {1, 2, 3}};
28 static int vnofaces[4] = {3, 1, 2, 0};
29 static int vFac[4][3] = {{0, 1, 2}, {0, 2, 3}, {0, 1, 3}, {1, 2, 3}};
30 
31 // as input, we give a tet and an edge, as return, we get all tets that share
32 // this edge and all vertices that are forming the outer ring of the cavity; we
33 // return true if the cavity is closed and false if it is open
34 
35 void computeNeighboringTetsOfACavity(const std::vector<MTet4 *> &cavity,
36  std::vector<MTet4 *> &outside)
37 {
38  outside.clear();
39  for(std::size_t i = 0; i < cavity.size(); i++) {
40  for(int j = 0; j < 4; j++) {
41  MTet4 *neigh = cavity[i]->getNeigh(j);
42  if(neigh) {
43  bool found = false;
44  for(std::size_t k = 0; k < outside.size(); k++) {
45  if(outside[k] == neigh) {
46  found = true;
47  break;
48  }
49  }
50  if(!found) {
51  for(std::size_t k = 0; k < cavity.size(); k++) {
52  if(cavity[k] == neigh) { found = true; }
53  }
54  }
55  if(!found) outside.push_back(neigh);
56  }
57  }
58  }
59 }
60 
61 bool buildEdgeCavity(MTet4 *t, int iLocalEdge, MVertex **v1, MVertex **v2,
62  std::vector<MTet4 *> &cavity,
63  std::vector<MTet4 *> &outside,
64  std::vector<MVertex *> &ring)
65 {
66  cavity.clear();
67  ring.clear();
68 
69  *v1 = t->tet()->getVertex(edges[iLocalEdge][0]);
70  *v2 = t->tet()->getVertex(edges[iLocalEdge][1]);
71 
72  // the 5 - i th edge contains the other 2 points of the tet
73  MVertex *lastinring = t->tet()->getVertex(edges[5 - iLocalEdge][0]);
74  ring.push_back(lastinring);
75  cavity.push_back(t);
76 
77  while(1) {
78  MVertex *ov1 = t->tet()->getVertex(edges[5 - iLocalEdge][0]);
79  MVertex *ov2 = t->tet()->getVertex(edges[5 - iLocalEdge][1]);
80  int K = ov1 == lastinring ? 1 : 0;
81  lastinring = ov1 == lastinring ? ov2 : ov1;
82  // look in the 2 faces sharing this edge the one that has vertex
83  // ov2 i.e. edges[5-iLocalEdge][1]
84  int iFace;
85  int iFace1 = efaces[iLocalEdge][0];
86  int iFace2 = efaces[iLocalEdge][1];
87  if(faces[iFace1][0] == edges[5 - iLocalEdge][K] ||
88  faces[iFace1][1] == edges[5 - iLocalEdge][K] ||
89  faces[iFace1][2] == edges[5 - iLocalEdge][K])
90  iFace = iFace1;
91  else if(faces[iFace2][0] == edges[5 - iLocalEdge][K] ||
92  faces[iFace2][1] == edges[5 - iLocalEdge][K] ||
93  faces[iFace2][2] == edges[5 - iLocalEdge][K])
94  iFace = iFace2;
95  else {
96  Msg::Error("Error of connexion");
97  return false;
98  }
99  t = t->getNeigh(iFace);
100  if(!t) return false;
101  if(t->isDeleted()) {
102  Msg::Warning("Strange edge cavity (tet is deleted)");
103  return false;
104  }
105  if(t == cavity[0]) break;
106  ring.push_back(lastinring);
107  cavity.push_back(t);
108  iLocalEdge = -1;
109  for(int i = 0; i < 6; i++) {
110  MVertex *a = t->tet()->getVertex(edges[i][0]);
111  MVertex *b = t->tet()->getVertex(edges[i][1]);
112  if((a == *v1 && b == *v2) || (a == *v2 && b == *v1)) {
113  iLocalEdge = i;
114  break;
115  }
116  }
117  if(iLocalEdge == -1) {
118  Msg::Warning("Strange edge cavity (local edge not found)");
119  return false;
120  }
121  // FIXME when hybrid mesh, this loops for ever
122  if(cavity.size() > 1000) {
123  // printf("cavity size gets laaaarge\n");
124  return false;
125  }
126  // printf("%d %d\n",ITER++, cavity.size());
127  }
128  computeNeighboringTetsOfACavity(cavity, outside);
129  return true;
130 }
131 
132 void BuildSwapPattern3(SwapPattern *sc)
133 {
134  static int trgl[][3] = {{0, 1, 2}};
135  static int trgul[][5] = {{0, -1, -1, -1, -1}};
136 
137  sc->nbr_triangles = 1;
138  sc->nbr_triangles_2 = 1;
139  sc->nbr_trianguls = 1;
140  sc->triangles = trgl;
141  sc->trianguls = trgul;
142 }
143 
144 /*
145  0 1
146  +------------+
147  | |
148  | |
149  | |
150  +------------+
151  3 2
152 
153 */
154 
155 void BuildSwapPattern4(SwapPattern *sc)
156 {
157  static int trgl[][3] = {{0, 1, 2}, {2, 3, 0}, {1, 2, 3}, {3, 0, 1}};
158  static int trgul[][5] = {{0, 1, -1, -1, -1}, {2, 3, -1, -1, -1}};
159 
160  sc->nbr_triangles = 4;
161  sc->nbr_triangles_2 = 2;
162  sc->nbr_trianguls = 2;
163  sc->triangles = trgl;
164  sc->trianguls = trgul;
165 }
166 
167 void BuildSwapPattern5(SwapPattern *sc)
168 {
169  static int trgl[][3] = {{0, 1, 2}, {0, 2, 3}, {0, 3, 4}, {0, 1, 4},
170  {1, 3, 4}, {1, 2, 3}, {2, 3, 4}, {0, 2, 4},
171  {0, 1, 3}, {1, 2, 4}};
172  static int trgul[][5] = {{0, 1, 2, -1, -1},
173  {3, 4, 5, -1, -1},
174  {0, 6, 7, -1, -1},
175  {2, 5, 8, -1, -1},
176  {3, 6, 9, -1, -1}};
177 
178  sc->nbr_triangles = 10;
179  sc->nbr_triangles_2 = 3;
180  sc->nbr_trianguls = 5;
181  sc->triangles = trgl;
182  sc->trianguls = trgul;
183 }
184 
185 void BuildSwapPattern6(SwapPattern *sc)
186 {
187  static int trgl[][3] = {
188  {0, 1, 2}, {0, 2, 3}, {0, 3, 4}, {0, 4, 5}, {0, 2, 5}, {2, 4, 5}, {2, 3, 4},
189  {0, 3, 5}, {3, 4, 5}, {0, 2, 4}, {2, 3, 5}, {1, 2, 3}, {0, 1, 3}, {0, 1, 5},
190  {1, 4, 5}, {1, 3, 4}, {0, 1, 4}, {1, 3, 5}, {1, 2, 4}, {1, 2, 5}};
191  static int trgul[][5] = {
192  {0, 1, 2, 3, -1}, {0, 4, 5, 6, -1}, {0, 1, 7, 8, -1},
193  {0, 3, 6, 9, -1}, {0, 4, 8, 10, -1}, {2, 3, 11, 12, -1},
194  {11, 13, 14, 15, -1}, {7, 8, 11, 12, -1}, {3, 11, 15, 16, -1},
195  {8, 11, 13, 17, -1}, {6, 13, 14, 18, -1}, {3, 6, 16, 18, -1},
196  {5, 6, 13, 19, -1}, {8, 10, 13, 19, -1}};
197 
198  sc->nbr_triangles = 20;
199  sc->nbr_triangles_2 = 4;
200  sc->nbr_trianguls = 14;
201  sc->triangles = trgl;
202  sc->trianguls = trgul;
203 }
204 
205 void BuildSwapPattern7(SwapPattern *sc)
206 {
207  static int trgl[][3] = {
208  {0, 1, 2}, {0, 2, 3}, {0, 3, 4}, {0, 4, 5}, {0, 5, 6}, {0, 3, 6},
209  {3, 5, 6}, {3, 4, 5}, {0, 4, 6}, {4, 5, 6}, {0, 3, 5}, {3, 4, 6},
210  {0, 2, 4}, {2, 3, 4}, {0, 2, 6}, {2, 5, 6}, {2, 4, 5}, {0, 2, 5},
211  {2, 4, 6}, {2, 3, 5}, {2, 3, 6}, {0, 1, 3}, {1, 2, 3}, {0, 1, 4},
212  {1, 3, 4}, {0, 1, 6}, {1, 5, 6}, {1, 4, 5}, {0, 1, 5}, {1, 4, 6},
213  {1, 3, 5}, {1, 3, 6}, {1, 2, 4}, {1, 2, 5}, {1, 2, 6}};
214  static int trgul[][5] = {
215  {0, 1, 2, 3, 4}, {0, 1, 5, 6, 7}, {0, 1, 2, 8, 9},
216  {0, 1, 4, 7, 10}, {0, 1, 5, 9, 11}, {0, 3, 4, 12, 13},
217  {0, 13, 14, 15, 16}, {0, 8, 9, 12, 13}, {0, 4, 13, 16, 17},
218  {0, 9, 13, 14, 18}, {0, 7, 14, 15, 19}, {0, 4, 7, 17, 19},
219  {0, 6, 7, 14, 20}, {0, 9, 11, 14, 20}, {2, 3, 4, 21, 22},
220  {5, 6, 7, 21, 22}, {2, 8, 9, 21, 22}, {4, 7, 10, 21, 22},
221  {5, 9, 11, 21, 22}, {3, 4, 22, 23, 24}, {22, 24, 25, 26, 27},
222  {8, 9, 22, 23, 24}, {4, 22, 24, 27, 28}, {9, 22, 24, 25, 29},
223  {7, 22, 25, 26, 30}, {4, 7, 22, 28, 30}, {6, 7, 22, 25, 31},
224  {9, 11, 22, 25, 31}, {3, 4, 13, 23, 32}, {13, 25, 26, 27, 32},
225  {8, 9, 13, 23, 32}, {4, 13, 27, 28, 32}, {9, 13, 25, 29, 32},
226  {13, 16, 25, 26, 33}, {4, 13, 16, 28, 33}, {13, 15, 16, 25, 34},
227  {9, 13, 18, 25, 34}, {7, 19, 25, 26, 33}, {4, 7, 19, 28, 33},
228  {7, 15, 19, 25, 34}, {6, 7, 20, 25, 34}, {9, 11, 20, 25, 34}};
229 
230  sc->nbr_triangles = 35;
231  sc->nbr_triangles_2 = 5;
232  sc->nbr_trianguls = 42;
233  sc->triangles = trgl;
234  sc->trianguls = trgul;
235 }
236 
237 bool edgeSwap(std::vector<MTet4 *> &newTets, MTet4 *tet, int iLocalEdge,
238  const qmTetrahedron::Measures &cr,
239  const std::set<MFace, MFaceLessThan> &embeddedFaces)
240 {
241  // static int edges[6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};
242  int permut[6] = {0, 3, 1, 2, 5, 4};
243  iLocalEdge = permut[iLocalEdge];
244 
245  std::vector<MTet4 *> cavity;
246  std::vector<MTet4 *> outside;
247  std::vector<MVertex *> ring;
248  MVertex *v1, *v2;
249 
250  // printf("a\n");
251  bool closed =
252  buildEdgeCavity(tet, iLocalEdge, &v1, &v2, cavity, outside, ring);
253  // printf("b\n");
254 
255  if(!closed) return false;
256 
257  double volumeRef = 0.0;
258  double tetQualityRef = 1;
259  for(std::size_t i = 0; i < cavity.size(); i++) {
260  double vol = fabs(cavity[i]->tet()->getVolume());
261  tetQualityRef = std::min(tetQualityRef, cavity[i]->getQuality());
262  volumeRef += vol;
263  }
264 
265  // build swap patterns
266 
267  SwapPattern sp;
268  switch(ring.size()) {
269  case 3: BuildSwapPattern3(&sp); break;
270  case 4: BuildSwapPattern4(&sp); break;
271  case 5: BuildSwapPattern5(&sp); break;
272  case 6: BuildSwapPattern6(&sp); break;
273  case 7: BuildSwapPattern7(&sp); break;
274  default: return false;
275  }
276 
277  // compute qualities of all tets that appear in the patterns
278  double tetQuality1[100], tetQuality2[100];
279  double volume1[100], volume2[100];
280  for(int i = 0; i < sp.nbr_triangles; i++) {
281  int p1 = sp.triangles[i][0];
282  int p2 = sp.triangles[i][1];
283  int p3 = sp.triangles[i][2];
284  tetQuality1[i] =
285  qmTetrahedron::qm(ring[p1], ring[p2], ring[p3], v1, cr, &(volume1[i]));
286  tetQuality2[i] =
287  qmTetrahedron::qm(ring[p1], ring[p2], ring[p3], v2, cr, &(volume2[i]));
288  }
289 
290  // look for the best triangulation, i.e. the one that maximize the minimum
291  // element quality
292  double minQuality[100];
293  // for all triangulations
294  for(int i = 0; i < sp.nbr_trianguls; i++) {
295  // for all triangles in a triangulation
296  minQuality[i] = 1;
297  double volume = 0;
298  for(int j = 0; j < sp.nbr_triangles_2; j++) {
299  int iT = sp.trianguls[i][j];
300  minQuality[i] = std::min(minQuality[i], tetQuality1[iT]);
301  minQuality[i] = std::min(minQuality[i], tetQuality2[iT]);
302  volume += (volume1[iT] + volume2[iT]);
303  }
304  // printf("config %3d qual %12.5E volume %12.5E\n",i,minQuality[i],volume);
305  if(fabs(volume - volumeRef) > 1.e-10 * (volume + volumeRef))
306  minQuality[i] = -1;
307  }
308 
309  int iBest = 0;
310  double best = -1.0;
311  for(int i = 0; i < sp.nbr_trianguls; i++) {
312  if(minQuality[i] > best) {
313  best = minQuality[i];
314  iBest = i;
315  }
316  }
317 
318  // if there exist no swap that enhance the quality
319  if(best <= tetQualityRef + 1e-20) return false;
320  // does random swaps if (best < .01) return false;
321 
322  // we have the best configuration, so we swap
323  // printf("outside size = %d\n",outside.size());
324 
325  // printf("a swap with %d tets reconnect %d tets cavity %d tets\n",
326  // ring.size(),outside.size(),cavity.size());
327 
328  for(int j = 0; j < sp.nbr_triangles_2; j++) {
329  int iT = sp.trianguls[iBest][j];
330  int p1 = sp.triangles[iT][0];
331  int p2 = sp.triangles[iT][1];
332  int p3 = sp.triangles[iT][2];
333  MVertex *pv1 = ring[p1];
334  MVertex *pv2 = ring[p2];
335  MVertex *pv3 = ring[p3];
336  MTetrahedron *tr1 = new MTetrahedron(pv1, pv2, pv3, v1);
337  MTetrahedron *tr2 = new MTetrahedron(pv3, pv2, pv1, v2);
338  MTet4 *t41 = new MTet4(tr1, tetQuality1[iT]);
339  MTet4 *t42 = new MTet4(tr2, tetQuality2[iT]);
340  t41->setOnWhat(cavity[0]->onWhat());
341  t42->setOnWhat(cavity[0]->onWhat());
342  outside.push_back(t41);
343  outside.push_back(t42);
344  newTets.push_back(t41);
345  newTets.push_back(t42);
346  }
347 
348  for(std::size_t i = 0; i < cavity.size(); i++) cavity[i]->setDeleted(true);
349 
350  connectTets(outside, &embeddedFaces);
351 
352  return true;
353 }
354 
355 bool edgeSplit(std::vector<MTet4 *> &newTets, MTet4 *tet, MVertex *newVertex,
356  int iLocalEdge, const qmTetrahedron::Measures &cr)
357 {
358  std::vector<MTet4 *> cavity;
359  std::vector<MTet4 *> outside;
360  std::vector<MVertex *> ring;
361  MVertex *v1, *v2;
362 
363  bool closed =
364  buildEdgeCavity(tet, iLocalEdge, &v1, &v2, cavity, outside, ring);
365  if(!closed) return false;
366 
367  for(std::size_t j = 0; j < ring.size(); j++) {
368  MVertex *pv1 = ring[j];
369  MVertex *pv2 = ring[(j + 1) % ring.size()];
370  MTetrahedron *tr1 = new MTetrahedron(pv1, pv2, newVertex, v1);
371  MTetrahedron *tr2 = new MTetrahedron(newVertex, pv2, pv1, v2);
372  MTet4 *t41 = new MTet4(tr1, cr);
373  MTet4 *t42 = new MTet4(tr2, cr);
374  t41->setOnWhat(cavity[0]->onWhat());
375  t42->setOnWhat(cavity[0]->onWhat());
376  outside.push_back(t41);
377  outside.push_back(t42);
378  newTets.push_back(t41);
379  newTets.push_back(t42);
380  }
381 
382  for(std::size_t i = 0; i < cavity.size(); i++) cavity[i]->setDeleted(true);
383 
384  connectTets(outside);
385 
386  return true;
387 }
388 
389 // swap a face i.e. remove a face shared by 2 tets
390 bool faceSwap(std::vector<MTet4 *> &newTets, MTet4 *t1, int iLocalFace,
391  const qmTetrahedron::Measures &cr,
392  const std::set<MFace, MFaceLessThan> &embeddedFaces)
393 {
394  MTet4 *t2 = t1->getNeigh(iLocalFace);
395  if(!t2) return false;
396  if(t1->onWhat() != t2->onWhat()) return false;
397 
398  MVertex *v1 = t1->tet()->getVertex(vnofaces[iLocalFace]);
399  MVertex *f1 = t1->tet()->getVertex(faces[iLocalFace][0]);
400  MVertex *f2 = t1->tet()->getVertex(faces[iLocalFace][1]);
401  MVertex *f3 = t1->tet()->getVertex(faces[iLocalFace][2]);
402  MVertex *v2 = nullptr;
403  for(int i = 0; i < 4; i++) {
404  MVertex *v = t2->tet()->getVertex(i);
405  if(v != f1 && v != f2 && v != f3) {
406  v2 = v;
407  break;
408  }
409  }
410  if(!v2) {
411  Msg::Warning("Impossible to swap face");
412  return false;
413  }
414 
415  // printf("%p %p -- %p %p %p\n",v1,v2,f1,f2,f3);
416 
417  double vol1 = fabs(t1->tet()->getVolume());
418  double q1 = t1->getQuality();
419  double vol2 = fabs(t2->tet()->getVolume());
420  double q2 = t2->getQuality();
421 
422  double vol3;
423  double q3 = qmTetrahedron::qm(f1, f2, v1, v2, cr, &vol3);
424  double vol4;
425  double q4 = qmTetrahedron::qm(f2, f3, v1, v2, cr, &vol4);
426  double vol5;
427  double q5 = qmTetrahedron::qm(f3, f1, v1, v2, cr, &vol5);
428 
429  if(fabs(vol1 + vol2 - vol3 - vol4 - vol5) > 1.e-10 * (vol1 + vol2))
430  return false;
431  if(std::min(q1, q2) > std::min(std::min(q3, q4), q5)) return false;
432  // printf("%g %g %g\n",vol1 + vol2, vol3 + vol4 + vol5,vol1 + vol2 - vol3 -
433  // vol4 - vol5); printf("qs = %g %g vs %g %g %g\n",q1,q2,q3,q4,q5);
434 
435  std::vector<MTet4 *> outside;
436  for(int i = 0; i < 4; i++) {
437  if(t1->getNeigh(i) && t1->getNeigh(i) != t2) {
438  bool found = false;
439  for(std::size_t j = 0; j < outside.size(); j++) {
440  if(outside[j] == t1->getNeigh(i)) {
441  found = true;
442  break;
443  }
444  }
445  if(!found) outside.push_back(t1->getNeigh(i));
446  }
447  }
448  for(int i = 0; i < 4; i++) {
449  if(t2->getNeigh(i) && t2->getNeigh(i) != t1) {
450  bool found = false;
451  for(std::size_t j = 0; j < outside.size(); j++) {
452  if(outside[j] == t2->getNeigh(i)) {
453  found = true;
454  break;
455  }
456  }
457  if(!found) outside.push_back(t2->getNeigh(i));
458  }
459  }
460 
461  // printf("we have a face swap %d\n",outside.size());
462 
463  t1->setDeleted(true);
464  t2->setDeleted(true);
465 
466  MTetrahedron *tr1 = new MTetrahedron(f1, f2, v1, v2);
467  MTetrahedron *tr2 = new MTetrahedron(f2, f3, v1, v2);
468  MTetrahedron *tr3 = new MTetrahedron(f3, f1, v1, v2);
469  MTet4 *t41 = new MTet4(tr1, q3);
470  MTet4 *t42 = new MTet4(tr2, q4);
471  MTet4 *t43 = new MTet4(tr3, q5);
472  t41->setOnWhat(t1->onWhat());
473  t42->setOnWhat(t1->onWhat());
474  t43->setOnWhat(t1->onWhat());
475  outside.push_back(t41);
476  outside.push_back(t42);
477  outside.push_back(t43);
478  newTets.push_back(t41);
479  newTets.push_back(t42);
480  newTets.push_back(t43);
481  connectTets(outside, &embeddedFaces);
482  return true;
483 }
484 
485 bool buildVertexCavity_recur(MTet4 *t, MVertex *v, std::vector<MTet4 *> &cavity)
486 {
487  if(t->isDeleted()) {
488  Msg::Warning("A deleted tet is a neighbor of a non deleted tet"
489  " - skipping cavity");
490  return false;
491  }
492  int iV = -1;
493  for(int i = 0; i < 4; i++) {
494  if(t->tet()->getVertex(i) == v) {
495  iV = i;
496  break;
497  }
498  }
499  if(iV == -1) {
500  Msg::Warning("Trying to build a cavity of tets for a node that does not "
501  "belong to this tet - skipping cavity");
502  return false;
503  }
504  for(int i = 0; i < 3; i++) {
505  MTet4 *neigh = t->getNeigh(vFac[iV][i]);
506  if(neigh) {
507  bool found = false;
508  for(std::size_t j = 0; j < cavity.size(); j++) {
509  if(cavity[j] == neigh) {
510  found = true;
511  j = cavity.size();
512  }
513  }
514  if(!found) {
515  cavity.push_back(neigh);
516  if(!buildVertexCavity_recur(neigh, v, cavity))
517  return false;
518  }
519  }
520  }
521  return true;
522 }
523 
524 // sliver removal by compound mesh modif postulate : the edge cannot be swopped
525 // so we split it, and then collapse the new vertex on another one (of course,
526 // not the other one on the unswappable edge) after that crap, the sliver is
527 // trashed
528 
529 bool collapseVertex(std::vector<MTet4 *> &newTets, MTet4 *t, int iVertex,
530  int iTarget, const qmTetrahedron::Measures &cr,
531  const localMeshModAction action, double *minQual)
532 {
533  if(t->isDeleted()) {
534  Msg::Warning("Impossible to collapse node");
535  return false;
536  }
537 
538  MVertex *v = t->tet()->getVertex(iVertex);
539  MVertex *tg = t->tet()->getVertex(iTarget);
540 
541  if(v->onWhat()->dim() < 3) return false;
542  if(tg->onWhat()->dim() < 3) return false;
543 
544  std::vector<MTet4 *> cavity_v;
545  std::vector<MTet4 *> outside;
546  cavity_v.push_back(t);
547  if(!buildVertexCavity_recur(t, v, cavity_v)) return false;
548 
549  std::vector<MTet4 *> toDelete;
550  std::vector<MTet4 *> toUpdate;
551  double volume = 0;
552  double worst = 1.0;
553  for(std::size_t i = 0; i < cavity_v.size(); i++) {
554  bool found = false;
555  volume += fabs(cavity_v[i]->tet()->getVolume());
556  double q = cavity_v[i]->getQuality();
557  worst = std::min(worst, q);
558  for(int j = 0; j < 4; j++) {
559  if(cavity_v[i]->tet()->getVertex(j) == tg) found = true;
560  }
561  if(found)
562  toDelete.push_back(cavity_v[i]);
563  else
564  toUpdate.push_back(cavity_v[i]);
565  }
566 
567  double x = v->x();
568  double y = v->y();
569  double z = v->z();
570  v->x() = tg->x();
571  v->y() = tg->y();
572  v->z() = tg->z();
573 
574  double volume_update = 0;
575 
576  double worstAfter = 1.0;
577  std::vector<double> newQuals(toUpdate.size());
578  for(std::size_t i = 0; i < toUpdate.size(); i++) {
579  double vv;
580  newQuals[i] = qmTetrahedron::qm(toUpdate[i]->tet(), cr, &vv);
581  worstAfter = std::min(worstAfter, newQuals[i]);
582  volume_update += vv;
583  }
584 
585  // printf("%12.5E %12.5E %12.5E %12.5E %d\n",
586  // volume,volume_update,worstAfter,worst,toUpdate.size());
587 
588  if(fabs(volume - volume_update) > 1.e-10 * volume || worstAfter < worst) {
589  v->x() = x;
590  v->y() = y;
591  v->z() = z;
592  return false;
593  }
594  if(action == GMSH_EVALONLY) {
595  *minQual = worstAfter;
596  return true;
597  }
598  // ok we collapse
599  computeNeighboringTetsOfACavity(cavity_v, outside);
600  for(std::size_t i = 0; i < toUpdate.size(); i++) {
601  MTetrahedron *tr1 = new MTetrahedron(
602  toUpdate[i]->tet()->getVertex(0) == v ? tg :
603  toUpdate[i]->tet()->getVertex(0),
604  toUpdate[i]->tet()->getVertex(1) == v ? tg :
605  toUpdate[i]->tet()->getVertex(1),
606  toUpdate[i]->tet()->getVertex(2) == v ? tg :
607  toUpdate[i]->tet()->getVertex(2),
608  toUpdate[i]->tet()->getVertex(3) == v ? tg :
609  toUpdate[i]->tet()->getVertex(3));
610  MTet4 *t41 = new MTet4(tr1, cr);
611  t41->setOnWhat(cavity_v[0]->onWhat());
612  t41->setQuality(newQuals[i]);
613  outside.push_back(t41);
614  newTets.push_back(t41);
615  }
616  for(std::size_t i = 0; i < cavity_v.size(); i++)
617  cavity_v[i]->setDeleted(true);
618 
619  connectTets(outside);
620 
621  return true;
622 }
623 
624 bool smoothVertex(MTet4 *t, int iVertex, const qmTetrahedron::Measures &cr)
625 {
626  if(t->isDeleted()) {
627  Msg::Warning("Impossible to collapse node");
628  return false;
629  }
630  if(t->tet()->getVertex(iVertex)->onWhat()->dim() < 3) return false;
631 
632  std::vector<MTet4 *> cavity;
633  cavity.push_back(t);
634  if(!buildVertexCavity_recur(t, t->tet()->getVertex(iVertex), cavity))
635  return false;
636 
637  double xcg = 0, ycg = 0, zcg = 0;
638  double vTot = 0;
639  double worst = 1.0;
640 
641  for(std::size_t i = 0; i < cavity.size(); i++) {
642  double volume = fabs(cavity[i]->tet()->getVolume());
643  double q = cavity[i]->getQuality();
644  worst = std::min(worst, q);
645  xcg += 0.25 *
646  (cavity[i]->tet()->getVertex(0)->x() +
647  cavity[i]->tet()->getVertex(1)->x() +
648  cavity[i]->tet()->getVertex(2)->x() +
649  cavity[i]->tet()->getVertex(3)->x()) *
650  volume;
651  ycg += 0.25 *
652  (cavity[i]->tet()->getVertex(0)->y() +
653  cavity[i]->tet()->getVertex(1)->y() +
654  cavity[i]->tet()->getVertex(2)->y() +
655  cavity[i]->tet()->getVertex(3)->y()) *
656  volume;
657  zcg += 0.25 *
658  (cavity[i]->tet()->getVertex(0)->z() +
659  cavity[i]->tet()->getVertex(1)->z() +
660  cavity[i]->tet()->getVertex(2)->z() +
661  cavity[i]->tet()->getVertex(3)->z()) *
662  volume;
663  vTot += volume;
664  }
665  xcg /= (vTot);
666  ycg /= (vTot);
667  zcg /= (vTot);
668  double volumeAfter = 0.0;
669 
670  double x = t->tet()->getVertex(iVertex)->x();
671  double y = t->tet()->getVertex(iVertex)->y();
672  double z = t->tet()->getVertex(iVertex)->z();
673 
674  t->tet()->getVertex(iVertex)->x() = xcg;
675  t->tet()->getVertex(iVertex)->y() = ycg;
676  t->tet()->getVertex(iVertex)->z() = zcg;
677  double worstAfter = 1.0;
678  std::vector<double> newQuals(cavity.size());
679  for(std::size_t i = 0; i < cavity.size(); i++) {
680  double volume;
681  newQuals[i] = qmTetrahedron::qm(cavity[i]->tet(), cr, &volume);
682  volumeAfter += volume;
683  worstAfter = std::min(worstAfter, newQuals[i]);
684  }
685 
686  if(fabs(volumeAfter - vTot) > 1.e-10 * vTot || worstAfter < worst) {
687  t->tet()->getVertex(iVertex)->x() = x;
688  t->tet()->getVertex(iVertex)->y() = y;
689  t->tet()->getVertex(iVertex)->z() = z;
690  return false; // smoothVertexOptimize(t, iVertex, cr);
691  }
692  else {
693  // restore new quality
694  for(std::size_t i = 0; i < cavity.size(); i++) {
695  cavity[i]->setQuality(newQuals[i]);
696  }
697  return true;
698  }
699 }
700 
701 struct smoothVertexData3D {
702  MVertex *v;
703  std::vector<MTet4 *> ts;
704  double LC;
705 };
706 
707 double smoothing_objective_function_3D(double X, double Y, double Z, MVertex *v,
708  std::vector<MTet4 *> &ts)
709 {
710  const double oldX = v->x();
711  const double oldY = v->y();
712  const double oldZ = v->z();
713  v->x() = X;
714  v->y() = Y;
715  v->z() = Z;
716 
717  auto it = ts.begin();
718  auto ite = ts.end();
719  double qMin = 1, vol;
720  while(it != ite) {
721  qMin = std::min(
722  qmTetrahedron::qm((*it)->tet(), qmTetrahedron::QMTET_GAMMA, &vol), qMin);
723  ++it;
724  }
725  v->x() = oldX;
726  v->y() = oldY;
727  v->z() = oldZ;
728  return -qMin;
729 }
730 
731 void deriv_smoothing_objective_function_3D(double *XYZ, double *dF, double &F,
732  void *data)
733 {
734  smoothVertexData3D *svd = (smoothVertexData3D *)data;
735  MVertex *v = svd->v;
736  const double LARGE = svd->LC * 1.e5;
737  const double SMALL = 1. / LARGE;
738  F = smoothing_objective_function_3D(XYZ[0], XYZ[1], XYZ[2], v, svd->ts);
739  double F_X =
740  smoothing_objective_function_3D(XYZ[0] + SMALL, XYZ[1], XYZ[2], v, svd->ts);
741  double F_Y =
742  smoothing_objective_function_3D(XYZ[0], XYZ[1] + SMALL, XYZ[2], v, svd->ts);
743  double F_Z =
744  smoothing_objective_function_3D(XYZ[0], XYZ[1], XYZ[2] + SMALL, v, svd->ts);
745  dF[0] = (F_X - F) * LARGE;
746  dF[1] = (F_Y - F) * LARGE;
747  dF[2] = (F_Z - F) * LARGE;
748 }
749 
750 double smooth_obj_3D(double *XYZ, void *data)
751 {
752  smoothVertexData3D *svd = (smoothVertexData3D *)data;
753  return smoothing_objective_function_3D(XYZ[0], XYZ[1], XYZ[2], svd->v,
754  svd->ts);
755 }
756 
757 bool smoothVertexOptimize(MTet4 *t, int iVertex,
758  const qmTetrahedron::Measures &cr)
759 {
760  if(t->tet()->getVertex(iVertex)->onWhat()->dim() < 3) return false;
761 
762  smoothVertexData3D vd;
763  vd.ts.push_back(t);
764  vd.v = t->tet()->getVertex(iVertex);
765  vd.LC = 1.0; // WRONG
766  if(!buildVertexCavity_recur(t, t->tet()->getVertex(iVertex), vd.ts))
767  return false;
768 
769  double xyzopti[3] = {vd.v->x(), vd.v->y(), vd.v->z()};
770 
771  // double val = 0.;
772  Msg::Error("Fletcher-Reeves minimizer routine must be reimplemented");
773  // minimize_N(3, smooth_obj_3D, deriv_smoothing_objective_function_3D, &vd, 4,
774  // xyzopti, val);
775 
776  double vTot = 0;
777 
778  for(std::size_t i = 0; i < vd.ts.size(); i++) {
779  double volume = fabs(vd.ts[i]->tet()->getVolume());
780  vTot += volume;
781  }
782 
783  double volumeAfter = 0.0;
784 
785  double x = t->tet()->getVertex(iVertex)->x();
786  double y = t->tet()->getVertex(iVertex)->y();
787  double z = t->tet()->getVertex(iVertex)->z();
788 
789  t->tet()->getVertex(iVertex)->x() = xyzopti[0];
790  t->tet()->getVertex(iVertex)->y() = xyzopti[1];
791  t->tet()->getVertex(iVertex)->z() = xyzopti[2];
792 
793  std::vector<double> newQuals(vd.ts.size());
794  for(std::size_t i = 0; i < vd.ts.size(); i++) {
795  double volume;
796  newQuals[i] = qmTetrahedron::qm(vd.ts[i]->tet(), cr, &volume);
797  volumeAfter += volume;
798  }
799 
800  if(fabs(volumeAfter - vTot) > 1.e-10 * vTot) {
801  t->tet()->getVertex(iVertex)->x() = x;
802  t->tet()->getVertex(iVertex)->y() = y;
803  t->tet()->getVertex(iVertex)->z() = z;
804  return false;
805  }
806  else {
807  // restore new quality
808  for(std::size_t i = 0; i < vd.ts.size(); i++) {
809  vd.ts[i]->setQuality(newQuals[i]);
810  }
811  return true;
812  }
813 }
814 
815 template <class ITERATOR>
816 void fillv_(std::multimap<MVertex *, MElement *> &vertexToElement,
817  ITERATOR it_beg, ITERATOR it_end)
818 {
819  for(ITERATOR IT = it_beg; IT != it_end; ++IT) {
820  MElement *el = *IT;
821  for(std::size_t j = 0; j < el->getNumVertices(); j++) {
822  MVertex *e = el->getVertex(j);
823  vertexToElement.insert(std::make_pair(e, el));
824  }
825  }
826 }
827 
828 int LaplaceSmoothing(GRegion *gr)
829 {
830  std::multimap<MVertex *, MElement *> vertexToElement;
831  fillv_(vertexToElement, (gr)->tetrahedra.begin(), (gr)->tetrahedra.end());
832  fillv_(vertexToElement, (gr)->hexahedra.begin(), (gr)->hexahedra.end());
833  fillv_(vertexToElement, (gr)->prisms.begin(), (gr)->prisms.end());
834  fillv_(vertexToElement, (gr)->pyramids.begin(), (gr)->pyramids.end());
835  int N = 0;
836  for(std::size_t i = 0; i < gr->mesh_vertices.size(); i++) {
837  MVertex *v = gr->mesh_vertices[i];
838  auto it = vertexToElement.lower_bound(v);
839  auto it_low = it;
840  auto it_up = vertexToElement.upper_bound(v);
841  double minQual = 1.e22;
842  double volTot = 0.0;
843  double xold = v->x(), yold = v->y(), zold = v->z();
844  SPoint3 pNew(0, 0, 0);
845  for(; it != it_up; ++it) {
846  minQual = std::min(minQual, it->second->minSICNShapeMeasure());
847  double vol = fabs(it->second->getVolume());
848  SPoint3 cog = it->second->barycenter();
849  pNew += cog * vol;
850  volTot += vol;
851  }
852  pNew *= (1. / volTot);
853  v->setXYZ(pNew.x(), pNew.y(), pNew.z());
854  double minQual2 = 1.e22;
855  for(it = it_low; it != it_up; ++it) {
856  minQual2 = std::min(minQual2, it->second->minSICNShapeMeasure());
857  if(minQual2 < minQual) {
858  v->setXYZ(xold, yold, zold);
859  break;
860  }
861  }
862  if(minQual < minQual2) N++;
863  }
864  return N;
865 }
smoothVertexData3D
Definition: meshGRegionLocalMeshMod.cpp:701
smooth_obj_3D
double smooth_obj_3D(double *XYZ, void *data)
Definition: meshGRegionLocalMeshMod.cpp:750
BuildSwapPattern4
void BuildSwapPattern4(SwapPattern *sc)
Definition: meshGRegionLocalMeshMod.cpp:155
BuildSwapPattern3
void BuildSwapPattern3(SwapPattern *sc)
Definition: meshGRegionLocalMeshMod.cpp:132
smoothing_objective_function_3D
double smoothing_objective_function_3D(double X, double Y, double Z, MVertex *v, std::vector< MTet4 * > &ts)
Definition: meshGRegionLocalMeshMod.cpp:707
BuildSwapPattern5
void BuildSwapPattern5(SwapPattern *sc)
Definition: meshGRegionLocalMeshMod.cpp:167
MTetrahedron
Definition: MTetrahedron.h:34
meshGRegionLocalMeshMod.h
vFac
static int vFac[4][3]
Definition: meshGRegionLocalMeshMod.cpp:29
F
#define F
Definition: DefaultOptions.h:23
edgeSplit
bool edgeSplit(std::vector< MTet4 * > &newTets, MTet4 *tet, MVertex *newVertex, int iLocalEdge, const qmTetrahedron::Measures &cr)
Definition: meshGRegionLocalMeshMod.cpp:355
GEntity.h
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
deriv_smoothing_objective_function_3D
void deriv_smoothing_objective_function_3D(double *XYZ, double *dF, double &F, void *data)
Definition: meshGRegionLocalMeshMod.cpp:731
SwapPattern::nbr_triangles_2
int nbr_triangles_2
Definition: meshGRegionLocalMeshMod.cpp:19
smoothVertexData3D::v
MVertex * v
Definition: meshGRegionLocalMeshMod.cpp:702
MVertex::onWhat
GEntity * onWhat() const
Definition: MVertex.h:82
MTet4::setQuality
void setQuality(const double &q)
Definition: meshGRegionDelaunayInsertion.h:166
GmshMessage.h
edges
static int edges[6][2]
Definition: meshGRegionLocalMeshMod.cpp:23
vnofaces
static int vnofaces[4]
Definition: meshGRegionLocalMeshMod.cpp:28
edgeSwap
bool edgeSwap(std::vector< MTet4 * > &newTets, MTet4 *tet, int iLocalEdge, const qmTetrahedron::Measures &cr, const std::set< MFace, MFaceLessThan > &embeddedFaces)
Definition: meshGRegionLocalMeshMod.cpp:237
MElement::getNumVertices
virtual std::size_t getNumVertices() const =0
MTet4
Definition: meshGRegionDelaunayInsertion.h:46
connectTets
void connectTets(ITER beg, ITER end, const std::set< MFace, MFaceLessThan > *allEmbeddedFaces=nullptr)
Definition: meshGRegionDelaunayInsertion.cpp:250
SwapPattern::triangles
int(* triangles)[3]
Definition: meshGRegionLocalMeshMod.cpp:17
MElement::getVertex
virtual const MVertex * getVertex(int num) const =0
SPoint3::x
double x(void) const
Definition: SPoint3.h:125
GRegion.h
GEntity::mesh_vertices
std::vector< MVertex * > mesh_vertices
Definition: GEntity.h:56
MHexahedron.h
MTet4::tet
MTetrahedron * tet() const
Definition: meshGRegionDelaunayInsertion.h:167
faces
static int faces[4][3]
Definition: meshGRegionLocalMeshMod.cpp:27
MTetrahedron::getVertex
virtual MVertex * getVertex(int num)
Definition: MTetrahedron.h:67
MTet4::onWhat
GRegion * onWhat() const
Definition: meshGRegionDelaunayInsertion.h:160
LaplaceSmoothing
int LaplaceSmoothing(GRegion *gr)
Definition: meshGRegionLocalMeshMod.cpp:828
Numeric.h
SwapPattern::nbr_triangles
int nbr_triangles
Definition: meshGRegionLocalMeshMod.cpp:16
qmTetrahedron::qm
static double qm(MTetrahedron *t, const Measures &cr, double *volume=nullptr)
Definition: qualityMeasures.cpp:616
MVertex::setXYZ
void setXYZ(double x, double y, double z)
Definition: MVertex.h:68
localMeshModAction
localMeshModAction
Definition: meshGRegionLocalMeshMod.h:15
MPyramid.h
GMSH_EVALONLY
@ GMSH_EVALONLY
Definition: meshGRegionLocalMeshMod.h:15
SPoint3::y
double y(void) const
Definition: SPoint3.h:127
MElement
Definition: MElement.h:30
SwapPattern::nbr_trianguls
int nbr_trianguls
Definition: meshGRegionLocalMeshMod.cpp:18
efaces
static int efaces[6][2]
Definition: meshGRegionLocalMeshMod.cpp:24
buildEdgeCavity
bool buildEdgeCavity(MTet4 *t, int iLocalEdge, MVertex **v1, MVertex **v2, std::vector< MTet4 * > &cavity, std::vector< MTet4 * > &outside, std::vector< MVertex * > &ring)
Definition: meshGRegionLocalMeshMod.cpp:61
fillv_
void fillv_(std::multimap< MVertex *, MElement * > &vertexToElement, ITERATOR it_beg, ITERATOR it_end)
Definition: meshGRegionLocalMeshMod.cpp:816
MTetrahedron::getVolume
virtual double getVolume()
Definition: MTetrahedron.cpp:114
smoothVertexOptimize
bool smoothVertexOptimize(MTet4 *t, int iVertex, const qmTetrahedron::Measures &cr)
Definition: meshGRegionLocalMeshMod.cpp:757
GRegion
Definition: GRegion.h:28
SwapPattern
Definition: meshGRegionLocalMeshMod.cpp:15
faceSwap
bool faceSwap(std::vector< MTet4 * > &newTets, MTet4 *t1, int iLocalFace, const qmTetrahedron::Measures &cr, const std::set< MFace, MFaceLessThan > &embeddedFaces)
Definition: meshGRegionLocalMeshMod.cpp:390
MTet4::getQuality
double getQuality() const
Definition: meshGRegionDelaunayInsertion.h:165
z
const double z
Definition: GaussQuadratureQuad.cpp:56
MTet4::setOnWhat
void setOnWhat(GRegion *g)
Definition: meshGRegionDelaunayInsertion.h:161
BuildSwapPattern6
void BuildSwapPattern6(SwapPattern *sc)
Definition: meshGRegionLocalMeshMod.cpp:185
BuildSwapPattern7
void BuildSwapPattern7(SwapPattern *sc)
Definition: meshGRegionLocalMeshMod.cpp:205
MPrism.h
XYZ
Definition: Camera.h:13
smoothVertexData3D::ts
std::vector< MTet4 * > ts
Definition: meshGRegionLocalMeshMod.cpp:703
qmTetrahedron::QMTET_GAMMA
@ QMTET_GAMMA
Definition: qualityMeasures.h:68
buildVertexCavity_recur
bool buildVertexCavity_recur(MTet4 *t, MVertex *v, std::vector< MTet4 * > &cavity)
Definition: meshGRegionLocalMeshMod.cpp:485
SPoint3::z
double z(void) const
Definition: SPoint3.h:129
SwapPattern::trianguls
int(* trianguls)[5]
Definition: meshGRegionLocalMeshMod.cpp:20
MVertex::y
double y() const
Definition: MVertex.h:61
smoothVertex
bool smoothVertex(MTet4 *t, int iVertex, const qmTetrahedron::Measures &cr)
Definition: meshGRegionLocalMeshMod.cpp:624
qmTetrahedron::Measures
Measures
Definition: qualityMeasures.h:68
GEntity::dim
virtual int dim() const
Definition: GEntity.h:196
computeNeighboringTetsOfACavity
void computeNeighboringTetsOfACavity(const std::vector< MTet4 * > &cavity, std::vector< MTet4 * > &outside)
Definition: meshGRegionLocalMeshMod.cpp:35
MTet4::getNeigh
MTet4 * getNeigh(int iN) const
Definition: meshGRegionDelaunayInsertion.h:171
MTet4::isDeleted
bool isDeleted() const
Definition: meshGRegionDelaunayInsertion.h:162
smoothVertexData3D::LC
double LC
Definition: meshGRegionLocalMeshMod.cpp:704
MTet4::setDeleted
void setDeleted(bool const d)
Definition: meshGRegionDelaunayInsertion.h:194
MVertex::x
double x() const
Definition: MVertex.h:60
collapseVertex
bool collapseVertex(std::vector< MTet4 * > &newTets, MTet4 *t, int iVertex, int iTarget, const qmTetrahedron::Measures &cr, const localMeshModAction action, double *minQual)
Definition: meshGRegionLocalMeshMod.cpp:529