gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
meshGFacePack.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 <stdio.h>
7 #include <stdlib.h>
8 #include <queue>
9 #include <map>
10 #include <vector>
11 #include "gmsh.h"
12 #include "SPoint2.h"
13 #include "SVector3.h"
14 #include "meshTriangulation.h"
15 #include "robustPredicates.h"
16 
17 static void analyze2dMetric(std::vector<double> &val, double &C, double &S,
18  double &h1, double &h2)
19 {
20  double a = val[0];
21  double c = val[1];
22  double b = val[4];
23  double l1 = 0.5 * (a + b + sqrt((a - b) * (a - b) + 4 * c * c));
24  double l2 = 0.5 * (a + b - sqrt((a - b) * (a - b) + 4 * c * c));
25  // printf("%12.5E %12.5E %12.5E %12.5E %12.5E\n",l1,l2,a,b,c);
26  h1 = 1. / sqrt(l1);
27  h2 = 1. / sqrt(l2);
28  double theta = 0.5 * atan2(2 * c, a - b);
29  C = cos(theta);
30  S = sin(theta);
31 }
32 
33 double bestParabola(double x0, double y0, double x1, double y1, double &xmid,
34  double &ymid, int VIEW_TAG);
35 double objectiveFunction_L(double x0, double y0, double xt, double yt,
36  double x1, double y1, int VIEW_TAG);
37 
38 static double maxDir(double &c, double &s, double &h, const double &C,
39  const double &S, const double &H1, const double &H2)
40 {
41  double PV1 = c * C + s * S;
42  double PV2 = s * C - c * S;
43  double PV3 = -c * C - s * S;
44  double PV4 = -s * C + c * S;
45  if(PV1 > PV2 && PV1 > PV3 && PV1 > PV4) {
46  c = C;
47  s = S;
48  h = H1;
49  return PV1;
50  }
51  else if(PV2 > PV3 && PV2 > PV4) {
52  c = -S;
53  s = C;
54  h = H2;
55  return PV2;
56  }
57  else if(PV3 > PV4) {
58  c = -C;
59  s = -S;
60  h = H1;
61  return PV3;
62  }
63  else {
64  c = S;
65  s = -C;
66  h = H2;
67  return PV4;
68  }
69 }
70 
71 static int check_triangle_validity_2d(double *xa, double *xb, double *xc)
72 {
73  return -robustPredicates::orient2d(xa, xb, xc) > 1.e-16;
74 }
75 
76 static double interpolateTriangleP2(double xi, double eta, double v0, double v1,
77  double v2, double v3, double v4, double v5)
78 {
79  double L0 = (1 - xi - eta) * (1 - 2 * xi - 2 * eta);
80  double L1 = xi * (2 * xi - 1);
81  double L2 = eta * (2 * eta - 1);
82  double L3 = 4 * xi * (1 - xi - eta);
83  double L4 = 4 * xi * eta;
84  double L5 = 4 * eta * (1 - xi - eta);
85 
86  return L0 * v0 + L1 * v1 + L2 * v2 + L3 * v3 + L4 * v4 + L5 * v5;
87 }
88 
89 static void derivativeP2(double u, double v, double dfdu[6], double dfdv[6])
90 {
91  /* const double sf[6][6] = { { 1, -3, -3, 4, 2, 2},
92  { 0, -1, 0, 0, 2, 0},
93  { 0, 0, -1, 0, 0, 2},
94  { 0, 4, 0, -4, -4, 0},
95  { 0, 0, 0, 4, 0, 0},
96  { 0, 0, 4, -4, 0, -4} };
97  */
98  dfdu[0] = -3 + 4 * v + 4 * u;
99  dfdu[1] = -1 + 4 * u;
100  dfdu[2] = 0;
101  dfdu[3] = 4 - 4 * v - 8 * u;
102  dfdu[4] = 4 * v;
103  dfdu[5] = -4 * v;
104 
105  dfdv[0] = -3 + 4 * u + 4 * v;
106  dfdv[1] = 0;
107  dfdv[2] = -1 + 4 * v;
108  dfdv[3] = -4 * u;
109  dfdv[4] = 4 * u;
110  dfdv[5] = 4 - 4 * u - 8 * v;
111 }
112 
113 static void jac_corners_p2(double *xa, double *xb, double *xc, double *xab,
114  double *xbc, double *xca, double J[6])
115 {
116  double *x[6] = {xa, xb, xc, xab, xbc, xca};
117  double nodes[6][2] = {{0, 0}, {1, 0}, {0, 1}, {0.5, 0}, {0.5, 0.5}, {0, 0.5}};
118  for(int i = 0; i < 6; i++) {
119  double u = nodes[i][0];
120  double v = nodes[i][1];
121  double dfdu[6], dfdv[6];
122  derivativeP2(u, v, dfdu, dfdv);
123  double dxdu = 0, dxdv = 0, dydu = 0, dydv = 0;
124  for(int j = 0; j < 6; j++) {
125  dxdu += x[j][0] * dfdu[j];
126  dxdv += x[j][0] * dfdv[j];
127  dydu += x[j][1] * dfdu[j];
128  dydv += x[j][1] * dfdv[j];
129  }
130  J[i] = -(dxdu * dydv - dydu * dxdv);
131  }
132 }
133 
134 static void jac_points_p2(double *xa, double *xb, double *xc, double *xab,
135  double *xbc, double *xca, double *xis, double *etas,
136  double detJ[6], double J[6][4])
137 {
138  double *x[6] = {xa, xb, xc, xab, xbc, xca};
139  for(int i = 0; i < 6; i++) {
140  double u = xis[i];
141  double v = etas[i];
142  double dfdu[6], dfdv[6];
143  derivativeP2(u, v, dfdu, dfdv);
144  double dxdu = 0, dxdv = 0, dydu = 0, dydv = 0;
145  for(int j = 0; j < 6; j++) {
146  dxdu += x[j][0] * dfdu[j];
147  dxdv += x[j][0] * dfdv[j];
148  dydu += x[j][1] * dfdu[j];
149  dydv += x[j][1] * dfdv[j];
150  }
151  J[i][0] = dxdu;
152  J[i][1] = dxdv;
153  J[i][2] = dydu;
154  J[i][3] = dydv;
155  detJ[i] = -(dxdu * dydv - dydu * dxdv);
156  }
157 }
158 
159 static double validity_p2triangle_formula(double *xa, double *xb, double *xc,
160  double *xab, double *xbc, double *xca)
161 {
162  double J[6];
163  jac_corners_p2(xa, xb, xc, xab, xbc, xca, J);
164  double bez[6] = {J[0],
165  J[1],
166  J[2],
167  2 * J[3] - 0.5 * (J[0] + J[1]),
168  2 * J[4] - 0.5 * (J[1] + J[2]),
169  2 * J[5] - 0.5 * (J[0] + J[2])};
170  double _MIN = 1.e22;
171  for(int i = 0; i < 6; i++) _MIN = bez[i] < _MIN ? bez[i] : _MIN;
172  return _MIN;
173 }
174 
175 FILE *FF = NULL;
176 
177 static double closest(double t, double u)
178 {
179  const double D = M_PI / 2;
180  double U0 = fabs(u - 2 * D - t);
181  double U1 = fabs(u - D - t);
182  double U2 = fabs(u - t);
183  double U3 = fabs(u + D - t);
184  double U4 = fabs(u + 2 * D - t);
185  if(U0 < U1 && U0 < U2 && U0 < U3 && U0 < U4) return u - 2 * D;
186  if(U1 < U2 && U1 < U3 && U1 < U4) return u - D;
187  if(U2 < U3 && U2 < U4) return u;
188  if(U3 < U4) return u + D;
189  return u + 2 * D;
190 }
191 
192 static double p2triangle_alignment_quality_measure(double *xa, double *xb,
193  double *xc, double *xab,
194  double *xbc, double *xca,
195  int VIEW_TAG)
196 {
197  double xis[6] = {0, .2, .4, .6, .8, 1};
198  double etas[6] = {0, 0, 0, 0., 0., 0.};
199  double J[6];
200  double Js[6][4];
201  jac_points_p2(xa, xb, xc, xab, xbc, xca, xis, etas, J, Js);
202 
203  double SUM = 0;
204 
205  double THETA[6];
206  double THETAM[6];
207  double x[6][2];
208 
209  for(uint32_t j = 0; j < 6; j++) {
210  double xi = xis[j];
211  double eta = etas[j];
212 
213  x[j][0] = interpolateTriangleP2(xi, eta, xa[0], xb[0], xc[0], xab[0],
214  xbc[0], xca[0]);
215  x[j][1] = interpolateTriangleP2(xi, eta, xa[1], xb[1], xc[1], xab[1],
216  xbc[1], xca[1]);
217 
218  std::vector<double> val(9);
219  double dist;
220  gmsh::view::probe(VIEW_TAG, x[j][0], x[j][1], 0.0, val, dist);
221  double dxdu = Js[j][0];
222  double dydu = Js[j][2];
223 
224  double C, S, h1, h2;
225  analyze2dMetric(val, C, S, h1, h2);
226 
227  if(j == 0)
228  THETAM[j] = atan2(S, C);
229  else
230  THETAM[j] = closest(THETAM[j - 1], atan2(S, C));
231 
232  if(j == 0)
233  THETA[j] = atan2(dydu, dxdu);
234  else
235  THETA[j] = closest(THETA[j - 1], atan2(dydu, dxdu));
236  }
237  for(uint32_t j = 1; j < 5; j++) {
238  double dtheta = THETA[j + 1] - THETA[j - 1];
239  double dthetam = THETAM[j + 1] - THETAM[j - 1];
240  double D2 = (dtheta - dthetam) * (dtheta - dthetam);
241  SUM += exp(-(D2));
242  }
243  return SUM;
244 }
245 
246 static double p1triangle_quality_metric(int VIEW_TAG, PolyMesh::Vertex *v0,
247  PolyMesh::Vertex *v1,
248  PolyMesh::Vertex *v2)
249 {
250  double *xa = v0->position;
251  double *xb = v1->position;
252  double *xc = v2->position;
253 
254  std::vector<double> val(9);
255  double dist;
256  gmsh::view::probe(VIEW_TAG, (xa[0] + xb[0] + xc[0]) / 3,
257  (xa[1] + xb[1] + xc[1]) / 3, 0.0, val, dist);
258  double M[3] = {val[0], val[4], val[1]};
259  val.clear();
260 
261  // printf("%g %g %g\n", M[0],M[1],M[2]);
262 
263  double area = fabs(robustPredicates::orient2d(xa, xb, xc));
264 
265  double det = sqrt(M[0] * M[1] - M[2] * M[2]);
266 
267  double d1[2] = {xb[0] - xa[0], xb[1] - xa[1]};
268  double d2[2] = {xc[0] - xb[0], xc[1] - xb[1]};
269  double d3[2] = {xa[0] - xc[0], xa[1] - xc[1]};
270 
271  double tmp1 = d1[0] * d1[0] + d2[0] * d2[0] + d3[0] * d3[0];
272  double tmp2 = d1[1] * d1[1] + d2[1] * d2[1] + d3[1] * d3[1];
273  double tmp3 = 2 * (d1[0] * d1[1] + d2[0] * d2[1] + d3[0] * d3[1]);
274 
275  double l2 = (tmp1 * M[0] + tmp2 * M[1] + tmp3 * M[2]);
276 
277  // printf("%g\n",det);
278 
279  return det * area / 3.0 / l2;
280 }
281 
282 int lengthPathInMetricField(double lagrangianPoints[3][2],
283  double *lengthInMetricField, int VIEW_TAG)
284 {
285  *lengthInMetricField = 0.0;
286 
287  size_t n = 8;
288 
289  static double GL_pt8[8] = {-9.602898564975363e-01, -7.966664774136268e-01,
290  -5.255324099163290e-01, -1.834346424956498e-01,
291  1.834346424956498e-01, 5.255324099163290e-01,
292  7.966664774136268e-01, 9.602898564975363e-01};
293 
294  static double GL_wt8[8] = {1.012285362903768e-01, 2.223810344533745e-01,
295  3.137066458778874e-01, 3.626837833783620e-01,
296  3.626837833783620e-01, 3.137066458778874e-01,
297  2.223810344533745e-01, 1.012285362903768e-01};
298 
299  // double gaussPoints[2][2] = {{-1/sqrt(3.0),1},{1/sqrt(3.0),1}};
300 
301  for(size_t i = 0; i < n; ++i) {
302  double t = GL_pt8[i];
303  double w = GL_wt8[i];
304 
305  double x[2] = {-0.5 * t * (1 - t) * lagrangianPoints[0][0] +
306  (1 + t) * (1 - t) * lagrangianPoints[1][0] +
307  0.5 * t * (1 + t) * lagrangianPoints[2][0],
308  -0.5 * t * (1 - t) * lagrangianPoints[0][1] +
309  (1 + t) * (1 - t) * lagrangianPoints[1][1] +
310  0.5 * t * (1 + t) * lagrangianPoints[2][1]};
311 
312  double dx = -0.5 * (1 - 2 * t) * lagrangianPoints[0][0] -
313  2 * t * lagrangianPoints[1][0] +
314  0.5 * (1 + 2 * t) * lagrangianPoints[2][0];
315  double dy = -0.5 * (1 - 2 * t) * lagrangianPoints[0][1] -
316  2 * t * lagrangianPoints[1][1] +
317  0.5 * (1 + 2 * t) * lagrangianPoints[2][1];
318 
319  std::vector<double> val(9);
320  double dist;
321  gmsh::view::probe(VIEW_TAG, x[0], x[1], 0.0, val, dist);
322  if(val.empty()) { *lengthInMetricField = 1.e22; }
323  else {
324  double M[3] = {val[0], val[4], val[1]};
325  *lengthInMetricField =
326  *lengthInMetricField +
327  w * sqrt(M[0] * dx * dx + 2 * M[2] * dx * dy + M[1] * dy * dy);
328  }
329  }
330  return 0;
331 }
332 
333 double objectiveFunction_L(double x0, double y0, double xt, double yt,
334  double x1, double y1, int VIEW_TAG)
335 {
336  double lagrangianPoints[3][2] = {{x0, y0}, {xt, yt}, {x1, y1}};
337  double lengthInMetricField;
338 
339  lengthPathInMetricField(lagrangianPoints, &lengthInMetricField, VIEW_TAG);
340 
341  return lengthInMetricField;
342 }
343 
344 /*
345  An straight edge
346 
347 */
348 
349 double objectiveFunction(double x0, double y0, double x1, double y1, double x2,
350  double y2, int VIEW_TAG)
351 {
352  return objectiveFunction_L(x0, y0, x1, y1, x2, y2, VIEW_TAG);
353 }
354 
355 double goldenSection(double x0, double y0, double &xt, double &yt, double x1,
356  double y1, double astart, double bstart, double &gamma1,
357  double direction[2], // rename direction instead of jF?
358  int VIEW_TAG)
359 {
360  double tol = 1.e-7;
361  double fstart, f, f1, f2;
362 
363  fstart = objectiveFunction(x0, y0, xt, yt, x1, y1, VIEW_TAG);
364  double Lstart = objectiveFunction_L(x0, y0, xt, yt, x1, y1, VIEW_TAG);
365  // return Lstart;
366  if(Lstart > 2) return Lstart;
367 
368  double lambda = 1.0 / (0.5 * (sqrt(5) + 1.0));
369 
370  double a = astart;
371  double b = bstart;
372 
373  double xi = 0.0;
374  double xi1 = b - lambda * (b - a);
375  double xi2 = a + lambda * (b - a);
376 
377  f1 = objectiveFunction(x0, y0, xt - xi1 * direction[0],
378  yt - xi1 * direction[1], x1, y1, VIEW_TAG);
379  f2 = objectiveFunction(x0, y0, xt - xi2 * direction[0],
380  yt - xi2 * direction[1], x1, y1, VIEW_TAG);
381 
382  while(fabs(xi2 - xi1) > tol) {
383  if(f1 < f2) {
384  b = xi2;
385  xi2 = xi1;
386  f2 = f1;
387  xi1 = b - lambda * (b - a);
388  f1 = objectiveFunction(x0, y0, xt - xi1 * direction[0],
389  yt - xi1 * direction[1], x1, y1, VIEW_TAG);
390  }
391  else {
392  a = xi1;
393  xi1 = xi2;
394  f1 = f2;
395  xi2 = a + lambda * (b - a);
396  f2 = objectiveFunction(x0, y0, xt - xi2 * direction[0],
397  yt - xi2 * direction[1], x1, y1, VIEW_TAG);
398  }
399  }
400 
401  if(f1 < f2) {
402  xi = xi1;
403  f = f1;
404  }
405  else {
406  xi = xi2;
407  f = f2;
408  }
409 
410  if(fstart <= f) { // we don't want to move if we are on the minimum
411  xi = 0;
412  }
413  gamma1 = xi;
414 
415  // double Lend = objectiveFunction_L( x0, y0, xt-xi*direction[0],
416  // yt-xi*direction[1], x1, y1, VIEW_TAG); if (Lstart < 1.2) printf("xi =
417  // %12.5E fstart %12.5E f %12.5E\n",xi, Lstart, Lend);
418  return f;
419 }
420 
421 int shortestParabola(double *xstart, double *xend, int VIEW_TAG, double *xmid,
422  double *lengthGeodesicsTwoPoints)
423 {
424  double x0 = xstart[0];
425  double x1 = xend[0];
426  double y0 = xstart[1];
427  double y1 = xend[1];
428 
429  // double XM = (x0 + x1)/2;
430  // double YM = (y0 + y1)/2;
431 
432  // double Lstart = objectiveFunction_L( x0, y0, XM, YM, x1, y1, VIEW_TAG );
433  // double KK;
434  // if (Lstart > 3) KK = Lstart ;
435  // else KK = bestParabola( x0, y0, x1, y1, XM, YM, VIEW_TAG);
436 
437  // xmid[0] = XM;
438  // xmid[1] = YM;
439 
440  // *lengthGeodesicsTwoPoints = KK;
441  // return 0;
442 
443  xmid[0] = (x0 + x1) / 2;
444  xmid[1] = (y0 + y1) / 2;
445 
446  double direction[2] = {y1 - y0, x0 - x1};
447  double gamma1 = 0.0;
448  // double astart = -.01;
449  // double bstart = .01;
450 
451  // precompute (x0, y0, x1, y1, VIEW_TAG);
452 
453  // goldenSection(x0, y0, xmid[0], xmid[1], x1, y1, astart, bstart, gamma1,
454  // direction, VIEW_TAG);
455 
456  xmid[0] = xmid[0] - gamma1 * direction[0];
457  xmid[1] = xmid[1] - gamma1 * direction[1];
458 
459  *lengthGeodesicsTwoPoints =
460  objectiveFunction_L(x0, y0, xmid[0], xmid[1], x1, y1, VIEW_TAG);
461 
462  return 0;
463 
464  // printf("%12.5E %12.5E %12.5E %12.5E %12.5E\n",gamma1,
465  // direction[0],direction[1],xmid[0],xmid[1]); return 0;
466 }
467 
468 static double smallestLengthParabola(int VIEW_TAG, SPoint2 &p1, SPoint2 &p2,
469  SPoint2 &p3)
470 {
471  double lengthGeodesicsTwoPoints;
472  shortestParabola(p1, p2, VIEW_TAG, p3, &lengthGeodesicsTwoPoints);
473  // printf("LENGTH %12.5E\n",lengthGeodesicsTwoPoints);
474  return lengthGeodesicsTwoPoints;
475 }
476 
477 static double _LIMIT = 0.7;
478 static double closestPoint(int VIEW_TAG, std::vector<SPoint2> &_points,
479  SPoint2 &p, SPoint2 &c)
480 {
481  double closest = 2 * _LIMIT;
482  for(size_t i = 0; i < _points.size(); i++) {
483  SPoint2 pp = _points[i];
484  SPoint2 pmid;
485  double dist = smallestLengthParabola(VIEW_TAG, pp, p, pmid);
486  if(dist < closest) {
487  closest = dist;
488  c = pp;
489  }
490  }
491  return closest;
492 }
493 
494 static bool computeNeighbor(int VIEW_TAG, const SPoint2 &p, int DIR,
495  double adimensionalLength, SPoint2 &n)
496 {
497  std::vector<double> val;
498  double dist;
499  gmsh::view::probe(VIEW_TAG, p.x(), p.y(), 0.0, val, dist);
500  if(val.empty()) return false;
501  double C, S, h1, h2;
502  analyze2dMetric(val, C, S, h1, h2);
503  double A[4] = {C, -C, S, -S};
504  double B[4] = {S, -S, -C, C};
505  double H[4] = {h1, h1, h2, h2};
506 
507  size_t N = 50;
508  double dx = A[DIR];
509  double dy = B[DIR];
510  double h = H[DIR];
511  std::vector<SPoint2> path;
512  SPoint2 PP = p;
513  for(size_t iter = 0; iter < N; iter++) {
514  SPoint2 pp(PP.x() + dx * h * adimensionalLength / N,
515  PP.y() + dy * h * adimensionalLength / N);
516  double Cpp, Spp, h1pp, h2pp;
517  val.clear();
518  gmsh::view::probe(VIEW_TAG, pp.x(), pp.y(), 0.0, val, dist);
519  if(val.empty()) { iter = N + 1; }
520  else {
521  analyze2dMetric(val, Cpp, Spp, h1pp, h2pp);
522  maxDir(dx, dy, h, Cpp, Spp, h1pp, h2pp);
523  path.push_back(pp);
524  PP = pp;
525  }
526  }
527  if(path.size() == N) {
528  n = path[N - 1];
529  return true;
530  }
531  return false;
532 }
533 
534 static int delaunayEdgeCriterionAnIsotropic(PolyMesh::HalfEdge *he,
535  int VIEW_TAG)
536 {
537  if(!he->f) return -1;
538  if(he->opposite == NULL) return -1;
539  if(he->opposite->data != he->data) return -1;
540  if(he->f->data < 0) return -1;
541 
542  PolyMesh::Vertex *v0 = he->v;
543  PolyMesh::Vertex *v1 = he->next->v;
544  PolyMesh::Vertex *v2 = he->next->next->v;
545  PolyMesh::Vertex *v3 = he->opposite->next->next->v;
546 
547  if(!check_triangle_validity_2d(v3->position, v1->position, v2->position) ||
548  !check_triangle_validity_2d(v3->position, v2->position, v0->position))
549  return -1;
550 
551  double q012 = p1triangle_quality_metric(VIEW_TAG, v0, v1, v2);
552  double q103 = p1triangle_quality_metric(VIEW_TAG, v1, v0, v3);
553 
554  double q312 = p1triangle_quality_metric(VIEW_TAG, v3, v1, v2);
555  double q320 = p1triangle_quality_metric(VIEW_TAG, v3, v2, v0);
556 
557  if(std::min(q312, q320) > std::min(q012, q103)) return 1;
558  return 0;
559 }
560 
561 double
562 triangleValidityP2(PolyMesh::HalfEdge *hea,
563  std::map<PolyMesh::HalfEdge *, SPoint2> *midPoints = nullptr)
564 {
565  PolyMesh::HalfEdge *heb = hea->next;
566  PolyMesh::HalfEdge *hec = heb->next;
567  SPoint2 ab(0.5 * (hea->v->position.x() + heb->v->position.x()),
568  0.5 * (hea->v->position.y() + heb->v->position.y()));
569  SPoint2 bc(0.5 * (heb->v->position.x() + hec->v->position.x()),
570  0.5 * (heb->v->position.y() + hec->v->position.y()));
571  SPoint2 ca(0.5 * (hec->v->position.x() + hea->v->position.x()),
572  0.5 * (hec->v->position.y() + hea->v->position.y()));
573  if(midPoints) {
574  ab = (*midPoints)[hea];
575  bc = (*midPoints)[heb];
576  ca = (*midPoints)[hec];
577  }
578  return validity_p2triangle_formula(hea->v->position, heb->v->position,
579  hec->v->position, ab, bc, ca);
580 }
581 
582 double
583 triangleQualityP2(int VIEW_TAG, PolyMesh::HalfEdge *hea,
584  std::map<PolyMesh::HalfEdge *, SPoint2> *midPoints = nullptr)
585 {
586  // hea = hea->f->he->next->next;
587  PolyMesh::HalfEdge *heb = hea->next;
588  PolyMesh::HalfEdge *hec = heb->next;
589  SPoint2 ab(0.5 * (hea->v->position.x() + heb->v->position.x()),
590  0.5 * (hea->v->position.y() + heb->v->position.y()));
591  SPoint2 bc(0.5 * (heb->v->position.x() + hec->v->position.x()),
592  0.5 * (heb->v->position.y() + hec->v->position.y()));
593  SPoint2 ca(0.5 * (hec->v->position.x() + hea->v->position.x()),
594  0.5 * (hec->v->position.y() + hea->v->position.y()));
595  if(midPoints) {
596  ab = (*midPoints)[hea];
597  bc = (*midPoints)[heb];
598  ca = (*midPoints)[hec];
599  }
600 
601  if(FF) {
602  std::vector<double> val0(9), val1(9), val2(9), val3(9), val4(9), val5(9);
603  double dist;
604  gmsh::view::probe(VIEW_TAG, hea->v->position.x(), hea->v->position.y(), 0.0,
605  val0, dist);
606  gmsh::view::probe(VIEW_TAG, heb->v->position.x(), heb->v->position.y(), 0.0,
607  val1, dist);
608  gmsh::view::probe(VIEW_TAG, hec->v->position.x(), hec->v->position.y(), 0.0,
609  val2, dist);
610  gmsh::view::probe(VIEW_TAG, ab.x(), ab.y(), 0.0, val3, dist);
611  gmsh::view::probe(VIEW_TAG, bc.x(), bc.y(), 0.0, val4, dist);
612  gmsh::view::probe(VIEW_TAG, ca.x(), ca.y(), 0.0, val5, dist);
613 
614  fprintf(FF,
615  "ST2(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,0,%g,%g,0,%g,%g,0){%g,%g,%g,%"
616  "g,%g,%g};\n",
617  hea->v->position.x(), hea->v->position.y(), hea->v->position.z(),
618  heb->v->position.x(), heb->v->position.y(), heb->v->position.z(),
619  hec->v->position.x(), hec->v->position.y(), hec->v->position.z(),
620  ab.x(), ab.y(), bc.x(), bc.y(), ca.x(), ca.y(), val0[8], val1[8],
621  val2[8], val3[8], val4[8], val5[8]);
622  }
623 
624  double validity = triangleValidityP2(hea, midPoints);
625  if(validity < 0) return validity;
626  return p2triangle_alignment_quality_measure(
627  hea->v->position, heb->v->position, hec->v->position, ab, bc, ca, VIEW_TAG);
628 }
629 
630 double triangleValidityP1(PolyMesh::HalfEdge *hea)
631 {
632  PolyMesh::HalfEdge *heb = hea->next;
633  PolyMesh::HalfEdge *hec = heb->next;
634  return -robustPredicates::orient2d(hea->v->position, heb->v->position,
635  hec->v->position);
636 }
637 
638 double LENGTH(PolyMesh::HalfEdge *he,
639  std::map<PolyMesh::HalfEdge *, SPoint2> &midPoints, int VIEW_TAG,
640  bool straight = false)
641 {
642  double x0 = he->v->position[0];
643  double x1 = he->next->v->position[0];
644  double y0 = he->v->position[1];
645  double y1 = he->next->v->position[1];
646  SPoint2 mid =
647  straight ? SPoint2(0.5 * (x0 + x1), 0.5 * (y0 + y1)) : midPoints[he];
648  return objectiveFunction_L(x0, y0, mid.x(), mid.y(), x1, y1, VIEW_TAG);
649 }
650 
651 static double bestParabola(PolyMesh::HalfEdge *he, int VIEW_TAG,
652  std::map<PolyMesh::HalfEdge *, SPoint2> &midPoints,
653  double ximax = .2)
654 {
655  double XI_MIN = -ximax;
656  double XI_MAX = ximax;
657 
658  double x0 = he->v->position[0];
659  double x1 = he->next->v->position[0];
660  double y0 = he->v->position[1];
661  double y1 = he->next->v->position[1];
662 
663  double xmid = (x0 + x1) / 2;
664  double ymid = (y0 + y1) / 2;
665 
666  double direction[2] = {y1 - y0, x0 - x1};
667 
668  const int N = 31;
669 
670  SPoint2 CUR = midPoints[he->opposite];
671  double Q1 = triangleQualityP2(VIEW_TAG, he, &midPoints);
672  double Q2 = triangleQualityP2(VIEW_TAG, he->opposite, &midPoints);
673  double QMIN = std::min(Q1, Q2);
674  // printf("----------------- iterating START %4d %4d %12.5E %12.5E
675  // --------------------- \n",he->v->data, he->next->v->data, Q1, Q2); return;
676  for(int i = 0; i < N; i++) {
677  double t = (double)i / (N - 1);
678  double xi = XI_MIN + (XI_MAX - XI_MIN) * t;
679  SPoint2 X(xmid - xi * direction[0], ymid - xi * direction[1]);
680  midPoints[he] = midPoints[he->opposite] = X;
681  Q1 = triangleQualityP2(VIEW_TAG, he, &midPoints);
682  Q2 = triangleQualityP2(VIEW_TAG, he->opposite, &midPoints);
683  if(std::min(Q1, Q2) > QMIN) {
684  QMIN = std::min(Q1, Q2);
685  CUR = X;
686  }
687  }
688  midPoints[he] = midPoints[he->opposite] = CUR;
689  return QMIN;
690  // printf("%12.5E %12.5E %12.5E %12.5E\n",ximin,QMIN,triangleQualityP2
691  // (VIEW_TAG, he, &midPoints), triangleQualityP2 (VIEW_TAG, he->opposite,
692  // &midPoints));
693 }
694 
695 static double bestParabola(PolyMesh::HalfEdge *he,
696  std::map<PolyMesh::HalfEdge *, SPoint2> &midPoints,
697  double er(double *, double *, double *, double *,
698  double *, double *),
699  double ximax = 0.2)
700 {
701  double XI_MIN = -ximax;
702  double XI_MAX = ximax;
703 
704  double x0 = he->v->position[0];
705  double x1 = he->next->v->position[0];
706  double y0 = he->v->position[1];
707  double y1 = he->next->v->position[1];
708 
709  double xmid = (x0 + x1) / 2;
710  double ymid = (y0 + y1) / 2;
711 
712  double direction[2] = {y1 - y0, x0 - x1};
713 
714  const int N = 81;
715 
716  SPoint2 CUR = midPoints[he->opposite];
717  double Q1 =
718  er(he->v->position, he->next->v->position, he->next->next->v->position,
719  midPoints[he], midPoints[he->next], midPoints[he->next->next]);
720  double Q2 =
721  er(he->opposite->v->position, he->opposite->next->v->position,
722  he->opposite->next->next->v->position, midPoints[he->opposite],
723  midPoints[he->opposite->next], midPoints[he->opposite->next->next]);
724  double V1 = triangleValidityP2(he, &midPoints);
725  double V2 = triangleValidityP2(he->opposite, &midPoints);
726  if(V1 <= 0) Q1 = 1.e22;
727  if(V2 <= 0) Q2 = 1.e22;
728 
729  double QSUM = Q1 + Q2;
730  // printf("----------------- iterating START %4d %4d %12.5E %12.5E
731  // --------------------- \n",he->v->data, he->next->v->data, Q1, Q2); return;
732  for(int i = 0; i < N; i++) {
733  double t = (double)i / (N - 1);
734  double xi = XI_MIN + (XI_MAX - XI_MIN) * t;
735  SPoint2 X(xmid - xi * direction[0], ymid - xi * direction[1]);
736  midPoints[he] = midPoints[he->opposite] = X;
737 
738  double Q1 =
739  er(he->v->position, he->next->v->position, he->next->next->v->position,
740  midPoints[he], midPoints[he->next], midPoints[he->next->next]);
741  double Q2 =
742  er(he->opposite->v->position, he->opposite->next->v->position,
743  he->opposite->next->next->v->position, midPoints[he->opposite],
744  midPoints[he->opposite->next], midPoints[he->opposite->next->next]);
745 
746  double V1 = triangleValidityP2(he, &midPoints);
747  double V2 = triangleValidityP2(he->opposite, &midPoints);
748 
749  if(V1 <= 0) Q1 = 1.e22;
750  if(V2 <= 0) Q2 = 1.e22;
751 
752  if(Q1 + Q2 < QSUM) {
753  QSUM = Q1 + Q2;
754  CUR = X;
755  }
756  }
757  midPoints[he] = midPoints[he->opposite] = CUR;
758  return QSUM;
759  // printf("%12.5E %12.5E %12.5E %12.5E\n",ximin,QMIN,triangleQualityP2
760  // (VIEW_TAG, he, &midPoints), triangleQualityP2 (VIEW_TAG, he->opposite,
761  // &midPoints));
762 }
763 
764 double bestParabola(double x0, double y0, double x1, double y1, double &xmid,
765  double &ymid, int VIEW_TAG)
766 {
767  double XI_MIN = -0.2;
768  double XI_MAX = 0.2;
769 
770  SPoint2 P0(x0, y0);
771  SPoint2 P1(x1, y1);
772 
773  xmid = (x0 + x1) / 2;
774  ymid = (y0 + y1) / 2;
775  SPoint2 PM(xmid, ymid);
776  SPoint2 CUR = PM;
777  double direction[2] = {y1 - y0, x0 - x1};
778 
779  const int N = 11;
780 
781  double QMAX = -1.e22;
782  // printf("----------------- iterating START %4d %4d %12.5E %12.5E
783  // --------------------- \n",he->v->data, he->next->v->data, Q1, Q2); return;
784  for(int i = 0; i < N; i++) {
785  double t = (double)i / (N - 1);
786  double xi = XI_MIN + (XI_MAX - XI_MIN) * t;
787  SPoint2 X(xmid - xi * direction[0], ymid - xi * direction[1]);
788  double Q =
789  p2triangle_alignment_quality_measure(P0, P1, P0, PM, PM, PM, VIEW_TAG);
790  if(Q > QMAX) {
791  QMAX = Q;
792  CUR = X;
793  }
794  }
795  xmid = CUR.x();
796  ymid = CUR.y();
797  return objectiveFunction_L(x0, y0, xmid, ymid, x1, y1, VIEW_TAG);
798  // printf("%12.5E %12.5E %12.5E %12.5E\n",ximin,QMIN,triangleQualityP2
799  // (VIEW_TAG, he, &midPoints), triangleQualityP2 (VIEW_TAG, he->opposite,
800  // &midPoints));
801 }
802 
803 int computePointsUsingScaledCrossFieldPlanarP2(
804  const char *modelForMetric, const char *modelForMesh, int VIEW_TAG,
805  int faceTag, std::vector<double> &pts,
806  double er(double *, double *, double *, double *, double *, double *))
807 {
808  PolyMesh *pm;
809 
810  FILE *f = fopen("PTS.pos", "w");
811  fprintf(f, "View \" \"{\n");
812 
813  gmsh::model::setCurrent(modelForMesh);
814  int result = GFace2PolyMesh(faceTag, &pm);
815  if(result) return result;
816 
817  // std::queue<SPoint2> _q;
818  std::stack<SPoint2> _q;
819  std::vector<SPoint2> _points;
820  for(auto v : pm->vertices) {
821  if(pm->degree(v) < 0) {
822  SPoint2 pp(v->position.x(), v->position.y());
823  _q.push(pp);
824  _points.push_back(pp);
825  }
826  }
827 
828  gmsh::model::setCurrent(modelForMetric);
829 
830  std::vector<double> additional;
831 
832  while(!_q.empty()) {
833  // SPoint2 p = _q.front();
834  SPoint2 p = _q.top();
835  _q.pop();
836  // std::vector<SPoint2> n;
837  // computeNeighbors (VIEW_TAG, p , n);
838  for(int DIR = 0; DIR < 4; DIR++) {
839  SPoint2 pp;
840  if(computeNeighbor(VIEW_TAG, p, DIR, sqrt(_LIMIT), pp)) {
841  // printf("%12.5E %12.5E\n",pp.x(), pp.y());
842  SPoint2 close;
843  double dist = closestPoint(VIEW_TAG, _points, pp, close);
844  if(dist > _LIMIT) {
845  additional.push_back(pp.x());
846  additional.push_back(pp.y());
847  _q.push(pp);
848  fprintf(f, "SL(%22.15E,%22.15E,0,%22.15E,%22.15E,0){%g,%g};\n",
849  pp.x(), pp.y(), p.x(), p.y(), dist, dist);
850  _points.push_back(pp);
851  fprintf(f, "SP(%22.15E,%22.15E,0){%lu};\n", pp.x(), pp.y(),
852  _points.size());
853  }
854  else {
855  // fprintf(f,"SL(%22.15E,%22.15E,0,%22.15E,%22.15E,0){%g,%g};\n",close.x(),close.y(),pp.x(),pp.y(),dist,
856  //dist);
857  }
858  }
859  // if (_points.size() > 1900)break;
860  if(_points.size() % 100 == 0)
861  printf("q size %lu p size %lu\n", _q.size(), _points.size());
862  }
863  }
864  fprintf(f, "};\n");
865  fclose(f);
866 
867  gmsh::model::setCurrent(modelForMesh);
868  PolyMesh *newMesh = GFaceInitialMesh(faceTag, 1, &additional);
869  newMesh->print4debug(100);
870  gmsh::model::setCurrent(modelForMetric);
871 
872  while(1) {
873  int count = 0;
874  for(auto he : newMesh->hedges) {
875  if(delaunayEdgeCriterionAnIsotropic(he, VIEW_TAG) == 1) {
876  newMesh->swap_edge(he);
877  count++;
878  }
879  }
880  // printf("count %3d\n",count);
881  if(count == 0) break;
882  }
883  newMesh->print4debug(200);
884 
885  std::map<PolyMesh::HalfEdge *, SPoint2> midPoints;
886  std::map<PolyMesh::HalfEdge *, double> LS;
887  for(auto he : newMesh->hedges) {
888  SPoint2 p_0(he->v->position.x(), he->v->position.y());
889  SPoint2 p_1(he->next->v->position.x(), he->next->v->position.y());
890  SPoint2 p_m = (p_0 + p_1) * .5;
891  midPoints[he] = p_m;
892  }
893 
894  int iter = 0;
895  double ximax = 0.3 / 5;
896  while(iter++ < 6) {
897  for(auto he : newMesh->hedges) {
898  SPoint2 p_0(he->v->position.x(), he->v->position.y());
899  SPoint2 p_1(he->next->v->position.x(), he->next->v->position.y());
900  SPoint2 p_m = (p_0 + p_1) * .5;
901  if(he->opposite && he->f->data >= 0 && he->opposite->f->data >= 0) {
902  if(er && iter > 14)
903  bestParabola(he, midPoints, er, ximax);
904  else
905  bestParabola(he, VIEW_TAG, midPoints, ximax);
906  double L = LENGTH(he, midPoints, VIEW_TAG);
907  LS[he] = L;
908  LS[he->opposite] = L;
909  }
910  }
911  ximax += 0.3 / 5;
912  }
913 
914  iter = 0;
915  while(iter++ < 3) {
916  int count = 0;
917  for(auto he : newMesh->hedges) {
918  if(!he->f || he->opposite == NULL || he->opposite->data != he->data ||
919  he->f->data < 0 || he->opposite->f->data < 0)
920  continue;
921  double LA = LENGTH(he, midPoints, VIEW_TAG);
922  SPoint2 old = midPoints[he];
923  if(LA > 1) {
924  double A = er ? bestParabola(he, midPoints, er) :
925  bestParabola(he, VIEW_TAG, midPoints);
926  newMesh->swap_edge(he);
927  double B = er ? bestParabola(he, midPoints, er) :
928  bestParabola(he, VIEW_TAG, midPoints);
929 
930  double LB = LENGTH(he, midPoints, VIEW_TAG);
931  bool invalid = false;
932  bool better = false;
933  if(er && B > 1.e10) invalid = true;
934  if(er && B > A) better = true;
935  if(!er && B < 0) invalid = true;
936  if(!er && B < A) better = true;
937  if(invalid) {
938  newMesh->swap_edge(he);
939  midPoints[he] = midPoints[he->opposite] = old;
940  }
941  else if(better) {
942  newMesh->swap_edge(he);
943  if(er)
944  bestParabola(he, midPoints, er);
945  else
946  bestParabola(he, VIEW_TAG, midPoints);
947  }
948  else {
949  count++;
950  printf("%12.5E %12.5E %12.5E %12.5E %12.5E \n", A, LA, B, LB,
951  LENGTH(he, midPoints, VIEW_TAG));
952  }
953  }
954  }
955  printf("count %3d\n", count);
956  if(count == 0) break;
957  }
958 
959  {
960  FILE *FFF = fopen("meshMetric.pos", "w");
961  fprintf(FFF, "View \" metric \"{\n");
962 
963  for(auto he : newMesh->hedges) {
964  if(he->opposite && he->f->data >= 0 && he->opposite->f->data >= 0) {
965  FF = FFF;
966  triangleQualityP2(VIEW_TAG, he, &midPoints);
967  FF = NULL;
968  }
969  }
970  fprintf(FFF, "};\n");
971  fclose(FFF);
972  }
973 
974  {
975  char name[256];
976  sprintf(name, "polyMeshCurved%d.pos", 200);
977  FILE *f = fopen(name, "w");
978  sprintf(name, "polyMeshCurved%d.pos", 300);
979  FILE *f2 = fopen(name, "w");
980  fprintf(f, "View \" %s \"{\n", name);
981  fprintf(f2, "View \" %s \"{\n", name);
982  for(auto he : newMesh->hedges) {
983  SPoint2 p = midPoints[he];
984  newMesh->high_order_nodes.push_back(SVector3(p.x(), p.y(), 0.0));
985  if(he->f->data >= 0) {
986  // printf("%g %g %g\n",p[0], p[1],LS[he]);
987  double Lhe = LENGTH(he, midPoints, VIEW_TAG);
988  double Lhe2 = LENGTH(he, midPoints, VIEW_TAG, true);
989  fprintf(f, "SL2(%g,%g,0,%g,%g,0,%g,%g,0){%g,%g,%g};\n",
990  he->v->position.x(), he->v->position.y(),
991  he->next->v->position.x(), he->next->v->position.y(), p.x(),
992  p.y(), Lhe, Lhe, Lhe);
993  fprintf(f2, "SL(%g,%g,0,%g,%g,0){%g,%g};\n", he->v->position.x(),
994  he->v->position.y(), he->next->v->position.x(),
995  he->next->v->position.y(), Lhe2, Lhe2);
996  }
997  }
998  fprintf(f, "};\n");
999  fclose(f);
1000  fprintf(f2, "};\n");
1001  fclose(f2);
1002  }
1003 
1004  gmsh::model::setCurrent(modelForMesh);
1005  PolyMesh2GFace(newMesh, faceTag);
1006 
1007  return 0;
1008 }
D
#define D
Definition: DefaultOptions.h:24
objectiveFunction
double objectiveFunction(double x0, double y0, double x1, double y1, double x2, double y2, int VIEW_TAG)
Definition: meshGFacePack.cpp:349
analyze2dMetric
static void analyze2dMetric(std::vector< double > &val, double &C, double &S, double &h1, double &h2)
Definition: meshGFacePack.cpp:17
check_triangle_validity_2d
static int check_triangle_validity_2d(double *xa, double *xb, double *xc)
Definition: meshGFacePack.cpp:71
PolyMesh::Face::data
int data
Definition: meshPolyMesh.h:56
PolyMesh::swap_edge
int swap_edge(HalfEdge *he0)
Definition: meshPolyMesh.h:197
computePointsUsingScaledCrossFieldPlanarP2
int computePointsUsingScaledCrossFieldPlanarP2(const char *modelForMetric, const char *modelForMesh, int VIEW_TAG, int faceTag, std::vector< double > &pts, double er(double *, double *, double *, double *, double *, double *))
Definition: meshGFacePack.cpp:803
jac_points_p2
static void jac_points_p2(double *xa, double *xb, double *xc, double *xab, double *xbc, double *xca, double *xis, double *etas, double detJ[6], double J[6][4])
Definition: meshGFacePack.cpp:134
LENGTH
double LENGTH(PolyMesh::HalfEdge *he, std::map< PolyMesh::HalfEdge *, SPoint2 > &midPoints, int VIEW_TAG, bool straight=false)
Definition: meshGFacePack.cpp:638
SPoint2
Definition: SPoint2.h:12
_LIMIT
static double _LIMIT
Definition: meshGFacePack.cpp:477
PolyMesh::HalfEdge::v
Vertex * v
Definition: meshPolyMesh.h:38
maxDir
static double maxDir(double &c, double &s, double &h, const double &C, const double &S, const double &H1, const double &H2)
Definition: meshGFacePack.cpp:38
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
LegendrePolynomials::f
void f(int n, double u, double *val)
Definition: orthogonalBasis.cpp:77
goldenSection
double goldenSection(double x0, double y0, double &xt, double &yt, double x1, double y1, double astart, double bstart, double &gamma1, double direction[2], int VIEW_TAG)
Definition: meshGFacePack.cpp:355
SVector3
Definition: SVector3.h:16
shortestParabola
int shortestParabola(double *xstart, double *xend, int VIEW_TAG, double *xmid, double *lengthGeodesicsTwoPoints)
Definition: meshGFacePack.cpp:421
SVector3.h
lengthPathInMetricField
int lengthPathInMetricField(double lagrangianPoints[3][2], double *lengthInMetricField, int VIEW_TAG)
Definition: meshGFacePack.cpp:282
PolyMesh::HalfEdge::next
HalfEdge * next
Definition: meshPolyMesh.h:41
closest
static double closest(double t, double u)
Definition: meshGFacePack.cpp:177
closestPoint
static double closestPoint(int VIEW_TAG, std::vector< SPoint2 > &_points, SPoint2 &p, SPoint2 &c)
Definition: meshGFacePack.cpp:478
smallestLengthParabola
static double smallestLengthParabola(int VIEW_TAG, SPoint2 &p1, SPoint2 &p2, SPoint2 &p3)
Definition: meshGFacePack.cpp:468
validity_p2triangle_formula
static double validity_p2triangle_formula(double *xa, double *xb, double *xc, double *xab, double *xbc, double *xca)
Definition: meshGFacePack.cpp:159
FF
FILE * FF
Definition: meshGFacePack.cpp:175
GL_pt8
static double GL_pt8[8]
Definition: GaussLegendre1D.h:66
PolyMesh::Vertex
Definition: meshPolyMesh.h:21
GFace2PolyMesh
int GFace2PolyMesh(int faceTag, PolyMesh **pm)
Definition: meshTriangulation.cpp:195
GL_wt8
static double GL_wt8[8]
Definition: GaussLegendre1D.h:72
derivativeP2
static void derivativeP2(double u, double v, double dfdu[6], double dfdv[6])
Definition: meshGFacePack.cpp:89
SVector3::x
double x(void) const
Definition: SVector3.h:30
triangleValidityP1
double triangleValidityP1(PolyMesh::HalfEdge *hea)
Definition: meshGFacePack.cpp:630
S
#define S
Definition: DefaultOptions.h:21
bestParabola
double bestParabola(double x0, double y0, double x1, double y1, double &xmid, double &ymid, int VIEW_TAG)
Definition: meshGFacePack.cpp:764
GFaceInitialMesh
PolyMesh * GFaceInitialMesh(int faceTag, int recover, std::vector< double > *additional)
Definition: meshTriangulation.cpp:753
SPoint2::x
double x(void) const
Definition: SPoint2.h:86
jac_corners_p2
static void jac_corners_p2(double *xa, double *xb, double *xc, double *xab, double *xbc, double *xca, double J[6])
Definition: meshGFacePack.cpp:113
PolyMesh::hedges
std::vector< HalfEdge * > hedges
Definition: meshPolyMesh.h:60
p1triangle_quality_metric
static double p1triangle_quality_metric(int VIEW_TAG, PolyMesh::Vertex *v0, PolyMesh::Vertex *v1, PolyMesh::Vertex *v2)
Definition: meshGFacePack.cpp:246
meshTriangulation.h
SVector3::y
double y(void) const
Definition: SVector3.h:31
triangleValidityP2
double triangleValidityP2(PolyMesh::HalfEdge *hea, std::map< PolyMesh::HalfEdge *, SPoint2 > *midPoints=nullptr)
Definition: meshGFacePack.cpp:562
PolyMesh::print4debug
void print4debug(const int debugTag)
Definition: meshPolyMesh.h:73
SVector3::z
double z(void) const
Definition: SVector3.h:32
PolyMesh::HalfEdge::data
int data
Definition: meshPolyMesh.h:43
PolyMesh::HalfEdge::opposite
HalfEdge * opposite
Definition: meshPolyMesh.h:42
direction
static void direction(Vertex *v1, Vertex *v2, double d[3])
Definition: Geo.cpp:218
computeNeighbor
static bool computeNeighbor(int VIEW_TAG, const SPoint2 &p, int DIR, double adimensionalLength, SPoint2 &n)
Definition: meshGFacePack.cpp:494
PolyMesh::high_order_nodes
std::vector< SVector3 > high_order_nodes
Definition: meshPolyMesh.h:62
delaunayEdgeCriterionAnIsotropic
static int delaunayEdgeCriterionAnIsotropic(PolyMesh::HalfEdge *he, int VIEW_TAG)
Definition: meshGFacePack.cpp:534
interpolateTriangleP2
static double interpolateTriangleP2(double xi, double eta, double v0, double v1, double v2, double v3, double v4, double v5)
Definition: meshGFacePack.cpp:76
PolyMesh
Definition: meshPolyMesh.h:15
triangleQualityP2
double triangleQualityP2(int VIEW_TAG, PolyMesh::HalfEdge *hea, std::map< PolyMesh::HalfEdge *, SPoint2 > *midPoints=nullptr)
Definition: meshGFacePack.cpp:583
PolyMesh::degree
int degree(const Vertex *v) const
Definition: meshPolyMesh.h:102
PolyMesh::Vertex::position
SVector3 position
Definition: meshPolyMesh.h:27
robustPredicates::orient2d
REAL orient2d(REAL *pa, REAL *pb, REAL *pc)
Definition: robustPredicates.cpp:1633
PolyMesh::vertices
std::vector< Vertex * > vertices
Definition: meshPolyMesh.h:59
objectiveFunction_L
double objectiveFunction_L(double x0, double y0, double xt, double yt, double x1, double y1, int VIEW_TAG)
Definition: meshGFacePack.cpp:333
PolyMesh::HalfEdge
Definition: meshPolyMesh.h:32
PolyMesh2GFace
int PolyMesh2GFace(PolyMesh *pm, int faceTag)
Definition: meshTriangulation.cpp:88
PolyMesh::HalfEdge::f
Face * f
Definition: meshPolyMesh.h:39
SPoint2.h
p2triangle_alignment_quality_measure
static double p2triangle_alignment_quality_measure(double *xa, double *xb, double *xc, double *xab, double *xbc, double *xca, int VIEW_TAG)
Definition: meshGFacePack.cpp:192
SPoint2::y
double y(void) const
Definition: SPoint2.h:88