gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
qualityMeasures.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 "robustPredicates.h"
7 #include "qualityMeasures.h"
8 #include "BDS.h"
9 #include "MVertex.h"
10 #include "MTriangle.h"
11 #include "MQuadrangle.h"
12 #include "MTetrahedron.h"
13 #include "MPrism.h"
14 #include "MHexahedron.h"
15 #include "Numeric.h"
16 #include "polynomialBasis.h"
17 #include "JacobianBasis.h"
18 #include "GmshMessage.h"
19 #include <limits>
20 #include <string.h>
21 
22 namespace {
23 
24  // Compute unit vector and gradients w.r.t. x0, y0, z0
25  // Remark: gradients w.r.t. x1, y1, z1 not computed as they are just the
26  // opposite
27  inline void unitVecAndGrad(const SPoint3 &p0, const SPoint3 &p1,
28  SVector3 &vec, std::vector<SVector3> &grad)
29  {
30  vec = SVector3(p0, p1);
31  const double n = vec.normalize(), invN = 1 / n;
32  grad[0] = invN * vec[0] * vec;
33  grad[0][0] -= invN; // dv01/dx0
34  grad[1] = invN * vec[1] * vec;
35  grad[1][1] -= invN; // dv01/dy0
36  grad[2] = invN * vec[2] * vec;
37  grad[2][2] -= invN; // dv01/dz0
38  }
39 
40  // Given vectors (0, 1) and (0, 2), their gradients and opposite of gradients,
41  // and the unit normal vector, compute NCJ from area of triangle defined by
42  // both vectors and gradients w.r.t. x0, y0, z0, x1, y1, z1, x2, y2, z2
43  inline void
44  NCJAndGrad2D(const SVector3 &v01, const std::vector<SVector3> &dv01dp0,
45  const std::vector<SVector3> &dv01dp1, const SVector3 &v02,
46  const std::vector<SVector3> &dv02dp0,
47  const std::vector<SVector3> &dv02dp2, const SVector3 &normal,
48  double &NCJ, std::vector<double> &dNCJ)
49  {
50  const SVector3 dvndx0 =
51  crossprod(v01, dv02dp0[0]) +
52  crossprod(dv01dp0[0], v02); // v01 x dv02/dx0 + dv01/dx0 x v02
53  const SVector3 dvndy0 =
54  crossprod(v01, dv02dp0[1]) +
55  crossprod(dv01dp0[1], v02); // v01 x dv02/dy0 + dv01/dy0 x v02
56  const SVector3 dvndz0 =
57  crossprod(v01, dv02dp0[2]) +
58  crossprod(dv01dp0[2], v02); // v01 x dv02/dz0 + dv01/dz0 x v02
59  const SVector3 dvndx1 =
60  crossprod(dv01dp1[0], v02); // dv01/dx1 x v02 (= -dv01/dx0 x v02)
61  const SVector3 dvndy1 =
62  crossprod(dv01dp1[1], v02); // dv01/dy1 x v02 (= -dv01/dy0 x v02)
63  const SVector3 dvndz1 =
64  crossprod(dv01dp1[2], v02); // dv01/dz1 x v02 (= -dv01/dz0 x v02)
65  const SVector3 dvndx2 =
66  crossprod(v01, dv02dp2[0]); // v01 x dv02/dx2 (= v01 x -dv02/dx0)
67  const SVector3 dvndy2 =
68  crossprod(v01, dv02dp2[1]); // v01 x dv02/dy2 (= v01 x -dv02/dy0)
69  const SVector3 dvndz2 =
70  crossprod(v01, dv02dp2[2]); // v01 x dv02/dz2 (= v01 x -dv02/dz0)
71 
72  SVector3 vn = crossprod(v01, v02);
73  NCJ = dot(vn, normal); // NCJ
74  dNCJ[0] = dot(dvndx0, normal); // dNCJ/dx0
75  dNCJ[1] = dot(dvndy0, normal); // dNCJ/dy0
76  dNCJ[2] = dot(dvndz0, normal); // dNCJ/dz0
77  dNCJ[3] = dot(dvndx1, normal); // dNCJ/dx1
78  dNCJ[4] = dot(dvndy1, normal); // dNCJ/dy1
79  dNCJ[5] = dot(dvndz1, normal); // dNCJ/dz1
80  dNCJ[6] = dot(dvndx2, normal); // dNCJ/dx2
81  dNCJ[7] = dot(dvndy2, normal); // dNCJ/dy2
82  dNCJ[8] = dot(dvndz2, normal); // dNCJ/dz2
83  }
84 
86  // inline void revertVG(const fullMatrix<double> &vg, fullMatrix<double> &res)
87  //{
88  // res(0, 0) = -vg(0, 3); res(0, 1) = -vg(0, 4); res(0, 2) = -vg(0, 5);
89  // res(0, 6) = -vg(0, 6); res(1, 0) = -vg(1, 3); res(1, 1) = -vg(1, 4);
90  // res(1, 2) = -vg(1, 5); res(1, 6) = -vg(1, 6); res(2, 0) = -vg(2, 3);
91  // res(2, 1) = -vg(2, 4); res(2, 2) = -vg(2, 5); res(2, 6) = -vg(2, 6);
92  //}
93 
94  // Scatter the NCJ gradients at vertex iV w.r.t vertices i0, i1 and i2
95  // in the vector of gradients for 2D element of nV vertices
96  template <int iV, int nV, int i0, int i1, int i2>
97  inline void scatterNCJGrad(const std::vector<double> &dNCJi,
98  std::vector<double> &dNCJ)
99  {
100  dNCJ[(iV * nV + i0) * 3] = dNCJi[0];
101  dNCJ[(iV * nV + i0) * 3 + 1] = dNCJi[1];
102  dNCJ[(iV * nV + i0) * 3 + 2] = dNCJi[2];
103  dNCJ[(iV * nV + i1) * 3] = dNCJi[3];
104  dNCJ[(iV * nV + i1) * 3 + 1] = dNCJi[4];
105  dNCJ[(iV * nV + i1) * 3 + 2] = dNCJi[5];
106  dNCJ[(iV * nV + i2) * 3] = dNCJi[6];
107  dNCJ[(iV * nV + i2) * 3 + 1] = dNCJi[7];
108  dNCJ[(iV * nV + i2) * 3 + 2] = dNCJi[8];
109  }
110 
111  // Scatter the NCJ gradients at vertex iV w.r.t vertices i0, i1, i2 and i3
112  // in the vector of gradients for 3D element of nV vertices
113  template <int iV, int nV, int i0, int i1, int i2, int i3>
114  inline void scatterNCJGrad(const std::vector<double> &dNCJi,
115  std::vector<double> &dNCJ)
116  {
117  dNCJ[iV * nV + i0 * 3] = dNCJi[0];
118  dNCJ[iV * nV + i0 * 3 + 1] = dNCJi[1];
119  dNCJ[iV * nV + i0 * 3 + 2] = dNCJi[2];
120  dNCJ[iV * nV + i1 * 3] = dNCJi[3];
121  dNCJ[iV * nV + i1 * 3 + 1] = dNCJi[4];
122  dNCJ[iV * nV + i1 * 3 + 2] = dNCJi[5];
123  dNCJ[iV * nV + i2 * 3] = dNCJi[6];
124  dNCJ[iV * nV + i2 * 3 + 1] = dNCJi[7];
125  dNCJ[iV * nV + i2 * 3 + 2] = dNCJi[8];
126  dNCJ[iV * nV + i2 * 3] = dNCJi[9];
127  dNCJ[iV * nV + i2 * 3 + 1] = dNCJi[10];
128  dNCJ[iV * nV + i2 * 3 + 2] = dNCJi[11];
129  }
130 
131 } // namespace
132 
133 double qmTriangle::gamma(const BDS_Point *p1, const BDS_Point *p2,
134  const BDS_Point *p3)
135 {
136  return gamma(p1->X, p1->Y, p1->Z, p2->X, p2->Y, p2->Z, p3->X, p3->Y, p3->Z);
137 }
138 
139 double qmTriangle::gamma(BDS_Face *t)
140 {
141  BDS_Point *n[4];
142  t->getNodes(n);
143  return gamma(n[0], n[1], n[2]);
144 }
145 
146 double qmTriangle::gamma(MTriangle *t)
147 {
148  return gamma(t->getVertex(0), t->getVertex(1), t->getVertex(2));
149 }
150 
151 double qmTriangle::gamma(const MVertex *v1, const MVertex *v2,
152  const MVertex *v3)
153 {
154  return gamma(v1->x(), v1->y(), v1->z(), v2->x(), v2->y(), v2->z(), v3->x(),
155  v3->y(), v3->z());
156 }
157 
158 // Triangle abc
159 // quality is between 0 and 1
160 double qmTriangle::gamma(const double &xa, const double &ya, const double &za,
161  const double &xb, const double &yb, const double &zb,
162  const double &xc, const double &yc, const double &zc)
163 {
164  // quality = rho / R = 2 * inscribed radius / circumradius
165  double a[3] = {xc - xb, yc - yb, zc - zb};
166  double b[3] = {xa - xc, ya - yc, za - zc};
167  double c[3] = {xb - xa, yb - ya, zb - za};
168  norme(a);
169  norme(b);
170  norme(c);
171  double pva[3];
172  prodve(b, c, pva);
173  const double sina = norm3(pva);
174  double pvb[3];
175  prodve(c, a, pvb);
176  const double sinb = norm3(pvb);
177  double pvc[3];
178  prodve(a, b, pvc);
179  const double sinc = norm3(pvc);
180 
181  if(sina == 0.0 && sinb == 0.0 && sinc == 0.0)
182  return 0.0;
183  else
184  return 2 * (2 * sina * sinb * sinc / (sina + sinb + sinc));
185 }
186 
187 double qmTriangle::eta(MTriangle *el)
188 {
189  MVertex *_v[3] = {el->getVertex(0), el->getVertex(1), el->getVertex(2)};
190 
191  double a1 = 180 * angle3Vertices(_v[0], _v[1], _v[2]) / M_PI;
192  double a2 = 180 * angle3Vertices(_v[1], _v[2], _v[0]) / M_PI;
193  double a3 = 180 * angle3Vertices(_v[2], _v[0], _v[1]) / M_PI;
194 
195  double amin = std::min(std::min(a1, a2), a3);
196  double angle = std::abs(60. - amin);
197  return 1. - angle / 60;
198 }
199 
200 double qmTriangle::angles(MTriangle *e)
201 {
202  double a = 500;
203  double worst_quality = std::numeric_limits<double>::max();
204  double mat[3][3];
205  double mat2[3][3];
206  double den = atan(a * (M_PI / 9)) + atan(a * (M_PI / 9));
207 
208  // This matrix is used to "rotate" the triangle to get each vertex
209  // as the "origin" of the mapping in turn
210  double rot[3][3];
211  rot[0][0] = -1;
212  rot[0][1] = 1;
213  rot[0][2] = 0;
214  rot[1][0] = -1;
215  rot[1][1] = 0;
216  rot[1][2] = 0;
217  rot[2][0] = 0;
218  rot[2][1] = 0;
219  rot[2][2] = 1;
220  double tmp[3][3];
221 
222  // double minAngle = 120.0;
223  for(std::size_t i = 0; i < e->getNumPrimaryVertices(); i++) {
224  const double u = i == 1 ? 1 : 0;
225  const double v = i == 2 ? 1 : 0;
226  const double w = 0;
227  e->getJacobian(u, v, w, mat);
228  e->getPrimaryJacobian(u, v, w, mat2);
229  for(std::size_t j = 0; j < i; j++) {
230  matmat(rot, mat, tmp);
231  memcpy(mat, tmp, sizeof(mat));
232  }
233  // get angle
234  double v1[3] = {mat[0][0], mat[0][1], mat[0][2]};
235  double v2[3] = {mat[1][0], mat[1][1], mat[1][2]};
236  double v3[3] = {mat2[0][0], mat2[0][1], mat2[0][2]};
237  double v4[3] = {mat2[1][0], mat2[1][1], mat2[1][2]};
238  norme(v1);
239  norme(v2);
240  norme(v3);
241  norme(v4);
242  double v12[3], v34[3];
243  prodve(v1, v2, v12);
244  prodve(v3, v4, v34);
245  norme(v12);
246  norme(v34);
247  double const orientation = prosca(v12, v34);
248 
249  // If the triangle is "flipped" it's no good
250  if(orientation < 0) return -std::numeric_limits<double>::max();
251 
252  double const c = prosca(v1, v2);
253  double x = std::acos(c) - M_PI / 3;
254  // double angle = (x+M_PI/3)/M_PI*180;
255  double quality =
256  (std::atan(a * (x + M_PI / 9)) + std::atan(a * (M_PI / 9 - x))) / den;
257  worst_quality = std::min(worst_quality, quality);
258 
259  // minAngle = std::min(angle, minAngle);
260  // printf("Angle %g ", angle);
261  // printf("Quality %g\n",quality);
262  }
263  // printf("MinAngle %g \n", minAngle);
264  // return minAngle;
265 
266  return worst_quality;
267 }
268 
269 void qmTriangle::NCJRange(const MTriangle *el, double &valMin, double &valMax)
270 {
271  const JacobianBasis *jac = el->getJacobianFuncSpace();
272  fullMatrix<double> primNodesXYZ(3, 3);
273  for(int i = 0; i < jac->getNumPrimMapNodes(); i++) {
274  const MVertex *v = el->getVertex(i);
275  primNodesXYZ(i, 0) = v->x();
276  primNodesXYZ(i, 1) = v->y();
277  primNodesXYZ(i, 2) = v->z();
278  }
279  fullMatrix<double> nM(1, 3);
280  jac->getPrimNormal2D(primNodesXYZ, nM);
281  SVector3 normal(nM(0, 0), nM(0, 1), nM(0, 2));
282 
283  std::vector<double> ncj(3);
284  NCJ(el->getVertex(0)->point(), el->getVertex(1)->point(),
285  el->getVertex(2)->point(), normal, ncj);
286  valMin = *std::min_element(ncj.begin(), ncj.end());
287  valMax = *std::max_element(ncj.begin(), ncj.end());
288 }
289 
290 void qmTriangle::NCJ(const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2,
291  const SVector3 &normal, std::vector<double> &NCJ)
292 {
293  // Compute unit vectors for each edge
294  SVector3 v01n(p0, p1), v12n(p1, p2), v20n(p2, p0);
295  v01n.normalize();
296  v12n.normalize();
297  v20n.normalize();
298 
299  // Compute NCJ at vertex from unit vectors a and b as
300  // 0.5*||a^b||/A_equilateral Factor = 2./sqrt(3.) = 0.5/A_equilateral
301  NCJ[0] = 2.0 / std::sqrt(3.0) * dot(crossprod(v01n, -v20n), normal);
302  NCJ[1] = 2.0 / std::sqrt(3.0) * dot(crossprod(v12n, -v01n), normal);
303  NCJ[2] = 2.0 / std::sqrt(3.0) * dot(crossprod(v20n, -v12n), normal);
304 }
305 
306 // Compute NCJ and its gradients at corners
307 // Gradients packed in vector: (dNCJ0/dx0, dNCJ0/dy0, dNCJ0/dz0,
308 // dNCJ0/dx1, ... dNCJ0/dz3, dNCJ1/dx0, ...,
309 // dNCJ3/dz3)
310 void qmTriangle::NCJAndGradients(const SPoint3 &p0, const SPoint3 &p1,
311  const SPoint3 &p2, const SVector3 &normal,
312  std::vector<double> &NCJ,
313  std::vector<double> &dNCJ)
314 {
315  // Factor = 2./sqrt(3.) = 0.5/A_equilateral
316  static const double fact = 2. / sqrt(3.);
317 
318  // Compute unit vector, its gradients and opposite grandients for edge (0, 1)
319  SVector3 v01n, v10n;
320  std::vector<SVector3> dv01ndp0(3), dv01ndp1(3);
321  unitVecAndGrad(p0, p1, v01n, dv01ndp0);
322  v10n = -v01n;
323  for(int i = 0; i < 3; i++) dv01ndp1[i] = -dv01ndp0[i];
324  const std::vector<SVector3> &dv10ndp1 = dv01ndp0, &dv10ndp0 = dv01ndp1;
325 
326  // Compute unit vector, its gradients and opposite grandients for edge (1, 2)
327  SVector3 v12n, v21n;
328  std::vector<SVector3> dv12ndp1(3), dv12ndp2(3);
329  unitVecAndGrad(p1, p2, v12n, dv12ndp1);
330  v21n = -v12n;
331  for(int i = 0; i < 3; i++) dv12ndp2[i] = -dv12ndp1[i];
332  const std::vector<SVector3> &dv21ndp2 = dv12ndp1, &dv21ndp1 = dv12ndp2;
333 
334  // Compute unit vector, its gradients and opposite grandients for edge (2, 0)
335  SVector3 v20n, v02n;
336  std::vector<SVector3> dv20ndp2(3), dv20ndp0(3);
337  unitVecAndGrad(p2, p0, v20n, dv20ndp2);
338  v02n = -v20n;
339  for(int i = 0; i < 3; i++) dv20ndp0[i] = -dv20ndp2[i];
340  const std::vector<SVector3> &dv02ndp0 = dv20ndp2, &dv02ndp2 = dv20ndp0;
341 
342  // Compute NCJ at vertex 0 as 0.5*||u01^u02||/A_triEqui
343  // and gradients w.r.t. x0, y0, z0, x1, y1, z1, x2, y2, z2
344  std::vector<double> dNCJ0(9);
345  NCJAndGrad2D(v01n, dv01ndp0, dv01ndp1, v02n, dv02ndp0, dv02ndp2, normal,
346  NCJ[0], dNCJ0);
347  // dNCJ[0] = dNCJ0[0]; dNCJ[1] = dNCJ0[1]; dNCJ[2] = dNCJ0[2];
348  // dNCJ[3] = dNCJ0[3]; dNCJ[4] = dNCJ0[4]; dNCJ[5] = dNCJ0[5];
349  // dNCJ[6] = dNCJ0[6]; dNCJ[7] = dNCJ0[7]; dNCJ[8] = dNCJ0[8];
350  scatterNCJGrad<0, 3, 0, 1, 2>(dNCJ0, dNCJ);
351 
352  // Compute NCJ at vertex 1 as 0.5*||u12^u10||/A_triEqui
353  // and gradients w.r.t. x1, y1, z1, x2, y2, z2, x0, y0, z0
354  std::vector<double> dNCJ1(9);
355  NCJAndGrad2D(v12n, dv12ndp1, dv12ndp2, v10n, dv10ndp1, dv10ndp0, normal,
356  NCJ[1], dNCJ1);
357  // dNCJ[9] = dNCJ1[6]; dNCJ[10] = dNCJ1[7]; dNCJ[11] = dNCJ1[8];
358  // dNCJ[10] = dNCJ1[0]; dNCJ[11] = dNCJ1[1]; dNCJ[12] = dNCJ1[2];
359  // dNCJ[13] = dNCJ1[3]; dNCJ[14] = dNCJ1[4]; dNCJ[15] = dNCJ1[5];
360  scatterNCJGrad<1, 3, 1, 2, 0>(dNCJ1, dNCJ);
361 
362  // Compute NCJ at vertex 2 as 0.5*||u20^u21||/A_triEqui
363  // Compute NCJ at vertex 2 and gradients w.r.t. x2, y2, z2, x0, y0, z0, x1,
364  // y1, z1
365  std::vector<double> dNCJ2(9);
366  NCJAndGrad2D(v20n, dv20ndp2, dv20ndp0, v21n, dv21ndp2, dv21ndp1, normal,
367  NCJ[2], dNCJ2);
368  // dNCJ[16] = dNCJ2[3]; dNCJ[17] = dNCJ2[4]; dNCJ[18] = dNCJ2[5];
369  // dNCJ[19] = dNCJ2[6]; dNCJ[20] = dNCJ2[7]; dNCJ[21] = dNCJ2[8];
370  // dNCJ[22] = dNCJ2[0]; dNCJ[23] = dNCJ2[1]; dNCJ[24] = dNCJ2[2];
371  scatterNCJGrad<2, 3, 2, 0, 1>(dNCJ2, dNCJ);
372 
373  for(int i = 0; i < 3; i++) NCJ[i] *= fact;
374  for(int i = 0; i < 27; i++) dNCJ[i] *= fact;
375 
376  // for (int iV=0; iV<3; iV++) {
377  // std::cout << "DBGTT: Vertex " << iV << ":\n";
378  // std::cout << "DBGTT: -> NCJ = " << NCJ[iV] << "\n";
379  // for (unsigned ig=0; ig<3; ig++) {
380  // int ind = iV*9+ig*3;
381  // std::cout << "DBGTT: -> dNCJ/dp" << ig << " = (" << dNCJ[ind] <<
382  // ", " << dNCJ[ind+1] << ", " << dNCJ[ind+2] << ")\n";
388  // }
389  // }
390 }
391 
392 double qmQuadrangle::eta(MQuadrangle *el)
393 {
394  double AR = 1; // pow(el->minEdge()/el->maxEdge(),.25);
395 
396  MVertex *_v[4] = {el->getVertex(0), el->getVertex(1), el->getVertex(2),
397  el->getVertex(3)};
398 
399  SVector3 v01(_v[1]->x() - _v[0]->x(), _v[1]->y() - _v[0]->y(),
400  _v[1]->z() - _v[0]->z());
401  SVector3 v12(_v[2]->x() - _v[1]->x(), _v[2]->y() - _v[1]->y(),
402  _v[2]->z() - _v[1]->z());
403  SVector3 v23(_v[3]->x() - _v[2]->x(), _v[3]->y() - _v[2]->y(),
404  _v[3]->z() - _v[2]->z());
405  SVector3 v30(_v[0]->x() - _v[3]->x(), _v[0]->y() - _v[3]->y(),
406  _v[0]->z() - _v[3]->z());
407 
408  SVector3 a = crossprod(v01, v12);
409  SVector3 b = crossprod(v12, v23);
410  SVector3 c = crossprod(v23, v30);
411  SVector3 d = crossprod(v30, v01);
412 
413  double sign = 1.0;
414  if(dot(a, b) < 0 || dot(a, c) < 0 || dot(a, d) < 0) sign = -1;
415  // FIXME ...
416  // if (a.z() > 0 || b.z() > 0 || c.z() > 0 || d.z() > 0) sign = -1;
417 
418  double a1 = 180 * angle3Vertices(_v[0], _v[1], _v[2]) / M_PI;
419  double a2 = 180 * angle3Vertices(_v[1], _v[2], _v[3]) / M_PI;
420  double a3 = 180 * angle3Vertices(_v[2], _v[3], _v[0]) / M_PI;
421  double a4 = 180 * angle3Vertices(_v[3], _v[0], _v[1]) / M_PI;
422 
423  a1 = std::min(180., a1);
424  a2 = std::min(180., a2);
425  a3 = std::min(180., a3);
426  a4 = std::min(180., a4);
427  double angle = fabs(90. - a1);
428  angle = std::max(fabs(90. - a2), angle);
429  angle = std::max(fabs(90. - a3), angle);
430  angle = std::max(fabs(90. - a4), angle);
431 
432  return sign * (1. - angle / 90) * AR;
433 }
434 
435 double qmQuadrangle::angles(MQuadrangle *e)
436 {
437  double a = 100;
438  double worst_quality = std::numeric_limits<double>::max();
439  double mat[3][3];
440  double mat2[3][3];
441  double den = atan(a * (M_PI / 4)) + atan(a * (2 * M_PI / 4 - (M_PI / 4)));
442 
443  // This matrix is used to "rotate" the triangle to get each vertex
444  // as the "origin" of the mapping in turn
445  // double rot[3][3];
446  // rot[0][0]=-1; rot[0][1]=1; rot[0][2]=0;
447  // rot[1][0]=-1; rot[1][1]=0; rot[1][2]=0;
448  // rot[2][0]= 0; rot[2][1]=0; rot[2][2]=1;
449  // double tmp[3][3];
450 
451  const double u[9] = {-1, -1, 1, 1, 0, 0, 1, -1, 0};
452  const double v[9] = {-1, 1, 1, -1, -1, 1, 0, 0, 0};
453 
454  for(int i = 0; i < 9; i++) {
455  e->getJacobian(u[i], v[i], 0, mat);
456  e->getPrimaryJacobian(u[i], v[i], 0, mat2);
457  // for (int j = 0; j < i; j++) {
458  // matmat(rot,mat,tmp);
459  // memcpy(mat, tmp, sizeof(mat));
460  //}
461 
462  // get angle
463  double v1[3] = {mat[0][0], mat[0][1], mat[0][2]};
464  double v2[3] = {mat[1][0], mat[1][1], mat[1][2]};
465  double v3[3] = {mat2[0][0], mat2[0][1], mat2[0][2]};
466  double v4[3] = {mat2[1][0], mat2[1][1], mat2[1][2]};
467  norme(v1);
468  norme(v2);
469  norme(v3);
470  norme(v4);
471  double v12[3], v34[3];
472  prodve(v1, v2, v12);
473  prodve(v3, v4, v34);
474  norme(v12);
475  norme(v34);
476 
477  // If the if the triangle is "flipped" it's no good
478  // double const orientation = prosca(v12, v34);
479  // if (orientation < 0)
480  // return -std::numeric_limits<double>::max();
481 
482  double const c = prosca(v1, v2);
483  double const x = std::abs(std::acos(c)) - M_PI / 2;
484  // double angle = std::fabs(std::acos(c))*180/M_PI;
485  double const quality = (std::atan(a * (x + M_PI / 4)) +
486  std::atan(a * (2 * M_PI / 4 - (x + M_PI / 4)))) /
487  den;
488  worst_quality = std::min(worst_quality, quality);
489  }
490  return worst_quality;
491 }
492 
493 void qmQuadrangle::NCJRange(const MQuadrangle *el, double &valMin,
494  double &valMax)
495 {
496  const JacobianBasis *jac = el->getJacobianFuncSpace();
497  fullMatrix<double> primNodesXYZ(4, 3);
498  for(int i = 0; i < jac->getNumPrimMapNodes(); i++) {
499  const MVertex *v = el->getVertex(i);
500  primNodesXYZ(i, 0) = v->x();
501  primNodesXYZ(i, 1) = v->y();
502  primNodesXYZ(i, 2) = v->z();
503  }
504  fullMatrix<double> nM(1, 3);
505  jac->getPrimNormal2D(primNodesXYZ, nM);
506  SVector3 normal(nM(0, 0), nM(0, 1), nM(0, 2));
507 
508  std::vector<double> ncj(4);
509  NCJ(el->getVertex(0)->point(), el->getVertex(1)->point(),
510  el->getVertex(2)->point(), el->getVertex(3)->point(), normal, ncj);
511  valMin = *std::min_element(ncj.begin(), ncj.end());
512  valMax = *std::max_element(ncj.begin(), ncj.end());
513 }
514 
515 void qmQuadrangle::NCJ(const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2,
516  const SPoint3 &p3, const SVector3 &normal,
517  std::vector<double> &ncj)
518 {
519  // Compute unit vectors for each edge
520  SVector3 v01n(p0, p1), v12n(p1, p2), v23n(p2, p3), v30n(p3, p0);
521  v01n.normalize();
522  v12n.normalize();
523  v23n.normalize();
524  v30n.normalize();
525 
526  // Compute NCJ at vertex from unit vectors a and b as
527  // 0.5*||a^b||/A_equilateral
528  ncj[0] = dot(crossprod(v01n, -v30n), normal);
529  ncj[1] = dot(crossprod(v12n, -v01n), normal);
530  ncj[2] = dot(crossprod(v23n, -v12n), normal);
531  ncj[3] = dot(crossprod(v30n, -v23n), normal);
532 }
533 
534 void qmQuadrangle::NCJAndGradients(const SPoint3 &p0, const SPoint3 &p1,
535  const SPoint3 &p2, const SPoint3 &p3,
536  const SVector3 &normal,
537  std::vector<double> &NCJ,
538  std::vector<double> &dNCJ)
539 {
540  // Compute unit vector, its gradients and opposite gradients for edge (0,1)
541  SVector3 v01n, v10n;
542  std::vector<SVector3> dv01ndp0(3), dv01ndp1(3);
543  unitVecAndGrad(p0, p1, v01n, dv01ndp0);
544  v10n = -v01n;
545  for(int i = 0; i < 3; i++) dv01ndp1[i] = -dv01ndp0[i];
546  const std::vector<SVector3> &dv10ndp1 = dv01ndp0, &dv10ndp0 = dv01ndp1;
547 
548  // Compute unit vector, its gradients and opposite gradients for edge (1,2)
549  SVector3 v12n, v21n;
550  std::vector<SVector3> dv12ndp1(3), dv12ndp2(3);
551  unitVecAndGrad(p1, p2, v12n, dv12ndp1);
552  v21n = -v12n;
553  for(int i = 0; i < 3; i++) dv12ndp2[i] = -dv12ndp1[i];
554  const std::vector<SVector3> &dv21ndp2 = dv12ndp1, &dv21ndp1 = dv12ndp2;
555 
556  // Compute unit vector, its gradients and opposite gradients for edge (2,3)
557  SVector3 v23n, v32n;
558  std::vector<SVector3> dv23ndp2(3), dv23ndp3(3);
559  unitVecAndGrad(p2, p3, v23n, dv23ndp2);
560  v32n = -v23n;
561  for(int i = 0; i < 3; i++) dv23ndp3[i] = -dv23ndp2[i];
562  const std::vector<SVector3> &dv32ndp3 = dv23ndp2, &dv32ndp2 = dv23ndp3;
563 
564  // Compute unit vector, its gradients and opposite gradients for edge (3,0)
565  SVector3 v30n, v03n;
566  std::vector<SVector3> dv30ndp3(3), dv30ndp0(3);
567  unitVecAndGrad(p3, p0, v30n, dv30ndp3);
568  v03n = -v30n;
569  for(int i = 0; i < 3; i++) dv30ndp0[i] = -dv30ndp3[i];
570  const std::vector<SVector3> &dv03ndp0 = dv30ndp3, &dv03ndp3 = dv30ndp0;
571 
572  // Compute NCJ at vertex 0 as 0.5*||u01^u03||/A_triRect
573  // and gradients w.r.t. x0, y0, z0, x1, y1, z1, x3, y3, z3
574  std::vector<double> dNCJ0(9);
575  NCJAndGrad2D(v01n, dv01ndp0, dv01ndp1, v03n, dv03ndp0, dv03ndp3, normal,
576  NCJ[0], dNCJ0);
577  scatterNCJGrad<0, 4, 0, 1, 3>(dNCJ0, dNCJ);
578 
579  // Compute NCJ at vertex 1 as 0.5*||u12^u10||/A_triRect
580  // and gradients w.r.t. x1, y1, z1, x2, y2, z2, x0, y0, z0
581  std::vector<double> dNCJ1(9);
582  NCJAndGrad2D(v12n, dv12ndp1, dv12ndp2, v10n, dv10ndp1, dv10ndp0, normal,
583  NCJ[1], dNCJ1);
584  scatterNCJGrad<1, 4, 1, 2, 0>(dNCJ1, dNCJ);
585 
586  // Compute NCJ at vertex 2 as 0.5*||u23^u21||/A_triRect
587  // and gradients w.r.t. x2, y2, z2, x3, y3, z3, x1, y1, z1
588  std::vector<double> dNCJ2(9);
589  NCJAndGrad2D(v23n, dv23ndp2, dv23ndp3, v21n, dv21ndp2, dv21ndp1, normal,
590  NCJ[2], dNCJ2);
591  scatterNCJGrad<2, 4, 2, 3, 1>(dNCJ2, dNCJ);
592 
593  // Compute NCJ at vertex 3 as 0.5*||u30^u32||/A_triRect
594  // and gradients w.r.t. x3, y3, z3, x0, y0, z0, x2, y2, z2
595  std::vector<double> dNCJ3(9);
596  NCJAndGrad2D(v30n, dv30ndp3, dv30ndp0, v32n, dv32ndp3, dv32ndp2, normal,
597  NCJ[3], dNCJ3);
598  scatterNCJGrad<3, 4, 3, 0, 2>(dNCJ3, dNCJ);
599 
600  // for (int iV=0; iV<4; iV++) {
601  // std::cout << "DBGTT: Vertex " << iV << ":\n";
602  // std::cout << "DBGTT: -> NCJ = " << NCJ[iV] << "\n";
603  // for (unsigned ig=0; ig<4; ig++) {
604  // int ind = iV*12+ig*3;
605  // std::cout << "DBGTT: -> dNCJ/dp" << ig << " = (" << dNCJ[ind] <<
606  // ", " << dNCJ[ind+1] << ", " << dNCJ[ind+2] << ")\n";
612  // }
613  // }
614 }
615 
616 double qmTetrahedron::qm(MTetrahedron *t, const Measures &cr, double *volume)
617 {
618  return qm(t->getVertex(0), t->getVertex(1), t->getVertex(2), t->getVertex(3),
619  cr, volume);
620 }
621 
622 double qmTetrahedron::qm(const MVertex *v1, const MVertex *v2,
623  const MVertex *v3, const MVertex *v4,
624  const Measures &cr, double *volume)
625 {
626  return qm(v1->x(), v1->y(), v1->z(), v2->x(), v2->y(), v2->z(), v3->x(),
627  v3->y(), v3->z(), v4->x(), v4->y(), v4->z(), cr, volume);
628 }
629 
630 double qmTetrahedron::qm(const double &x1, const double &y1, const double &z1,
631  const double &x2, const double &y2, const double &z2,
632  const double &x3, const double &y3, const double &z3,
633  const double &x4, const double &y4, const double &z4,
634  const Measures &cr, double *volume)
635 {
636  switch(cr) {
637  case QMTET_ONE: return 1.0;
638  case QMTET_ETA:
639  return eta(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, volume);
640  case QMTET_GAMMA: {
641  double G = gamma(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, volume);
642  *volume = fabs(*volume);
643  return G;
644  }
645  case QMTET_COND:
646  return cond(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, volume);
647  default: Msg::Error("Unknown quality measure"); return 0.;
648  }
649 }
650 
651 double qmTetrahedron::eta(const double &x1, const double &y1, const double &z1,
652  const double &x2, const double &y2, const double &z2,
653  const double &x3, const double &y3, const double &z3,
654  const double &x4, const double &y4, const double &z4,
655  double *volume)
656 {
657  double p0[3] = {x1, y1, z1};
658  double p1[3] = {x2, y2, z2};
659  double p2[3] = {x3, y3, z3};
660  double p3[3] = {x4, y4, z4};
661 
662  *volume = fabs(robustPredicates::orient3d(p0, p1, p2, p3)) / 6.0;
663 
664  double l =
665  ((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1) + (z2 - z1) * (z2 - z1));
666  l += ((x3 - x1) * (x3 - x1) + (y3 - y1) * (y3 - y1) + (z3 - z1) * (z3 - z1));
667  l += ((x4 - x1) * (x4 - x1) + (y4 - y1) * (y4 - y1) + (z4 - z1) * (z4 - z1));
668  l += ((x3 - x2) * (x3 - x2) + (y3 - y2) * (y3 - y2) + (z3 - z2) * (z3 - z2));
669  l += ((x4 - x2) * (x4 - x2) + (y4 - y2) * (y4 - y2) + (z4 - z2) * (z4 - z2));
670  l += ((x3 - x4) * (x3 - x4) + (y3 - y4) * (y3 - y4) + (z3 - z4) * (z3 - z4));
671  return 12. * pow(3 * fabs(*volume), 2. / 3.) / l;
672 }
673 
674 double qmTetrahedron::gamma(const double &x1, const double &y1,
675  const double &z1, const double &x2,
676  const double &y2, const double &z2,
677  const double &x3, const double &y3,
678  const double &z3, const double &x4,
679  const double &y4, const double &z4, double *volume)
680 {
681  // quality = rho / R = 3 * inradius / circumradius
682 
683  double p0[3] = {x1, y1, z1};
684  double p1[3] = {x2, y2, z2};
685  double p2[3] = {x3, y3, z3};
686  double p3[3] = {x4, y4, z4};
687 
688  *volume = (robustPredicates::orient3d(p0, p1, p2, p3)) / 6.0;
689 
690  if(fabs(*volume) == 0) return 0;
691 
692  double la =
693  (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1) + (z2 - z1) * (z2 - z1);
694  double lb =
695  (x3 - x1) * (x3 - x1) + (y3 - y1) * (y3 - y1) + (z3 - z1) * (z3 - z1);
696  double lc =
697  (x4 - x1) * (x4 - x1) + (y4 - y1) * (y4 - y1) + (z4 - z1) * (z4 - z1);
698  double lA =
699  (x4 - x3) * (x4 - x3) + (y4 - y3) * (y4 - y3) + (z4 - z3) * (z4 - z3);
700  double lB =
701  (x4 - x2) * (x4 - x2) + (y4 - y2) * (y4 - y2) + (z4 - z2) * (z4 - z2);
702  double lC =
703  (x3 - x2) * (x3 - x2) + (y3 - y2) * (y3 - y2) + (z3 - z2) * (z3 - z2);
704 
705  double lalA = std::sqrt(la * lA);
706  double lblB = std::sqrt(lb * lB);
707  double lclC = std::sqrt(lc * lC);
708 
709  double insideSqrt = (lalA + lblB + lclC) * (lalA + lblB - lclC) *
710  (lalA - lblB + lclC) * (-lalA + lblB + lclC);
711 
712  // This happens when the 4 points are (nearly) co-planar
713  // => R is actually undetermined but the quality is (close to) zero
714  if(insideSqrt <= 0) return 0;
715 
716  double R = std::sqrt(insideSqrt) / 24 / *volume;
717 
718  double s1 = fabs(triangle_area(p0, p1, p2));
719  double s2 = fabs(triangle_area(p0, p2, p3));
720  double s3 = fabs(triangle_area(p0, p1, p3));
721  double s4 = fabs(triangle_area(p1, p2, p3));
722  double rho = 3 * 3 * *volume / (s1 + s2 + s3 + s4);
723 
724  return rho / R;
725 }
726 
727 double qmTetrahedron::cond(const double &x1, const double &y1, const double &z1,
728  const double &x2, const double &y2, const double &z2,
729  const double &x3, const double &y3, const double &z3,
730  const double &x4, const double &y4, const double &z4,
731  double *volume)
732 {
734  double INVW[3][3] = {{1, -1. / sqrt(3.), -1. / sqrt(6.)},
735  {0, 2 / sqrt(3.), -1. / sqrt(6.)},
736  {0, 0, sqrt(1.5)}};
737  double A[3][3] = {{x2 - x1, y2 - y1, z2 - z1},
738  {x3 - x1, y3 - y1, z3 - z1},
739  {x4 - x1, y4 - y1, z4 - z1}};
740  double S[3][3], INVS[3][3];
741  matmat(A, INVW, S);
742  *volume = inv3x3(S, INVS) * 0.70710678118654762; // 2/sqrt(2);
743  double normS = norm2(S);
744  double normINVS = norm2(INVS);
745  return normS * normINVS;
746 }
747 
748 // TODO: Replace this
749 static double prismNCJ(const MVertex *a, const MVertex *b, const MVertex *c,
750  const MVertex *d)
751 {
752  static const double fact = 2. / sqrt(3.);
753 
754  const SVector3 vec1 =
755  SVector3(b->x() - a->x(), b->y() - a->y(), b->z() - a->z());
756  const SVector3 vec2 =
757  SVector3(c->x() - a->x(), c->y() - a->y(), c->z() - a->z());
758  const SVector3 vec3 =
759  SVector3(d->x() - a->x(), d->y() - a->y(), d->z() - a->z());
760 
761  const double l1 = vec1.norm();
762  const double l2 = vec2.norm();
763  const double l3 = vec3.norm();
764 
765  const double val = dot(vec1, crossprod(vec2, vec3));
766  return fact * fabs(val) / (l1 * l2 * l3);
767 }
768 
769 double qmPrism::minNCJ(const MPrism *el)
770 {
771  const MVertex *a = el->getVertex(0), *b = el->getVertex(1),
772  *c = el->getVertex(2);
773  const MVertex *d = el->getVertex(3), *e = el->getVertex(4),
774  *f = el->getVertex(5);
775  const double j[6] = {prismNCJ(a, b, c, d), prismNCJ(b, a, c, e),
776  prismNCJ(c, a, b, f), prismNCJ(d, a, e, f),
777  prismNCJ(e, b, d, f), prismNCJ(f, c, d, e)};
778  const double result = *std::min_element(j, j + 6);
779  return result;
780 }
781 
782 // void qmPrism::NCJ(const double &x0, const double &y0, const double &z0,
783 // const double &x1, const double &y1, const double &z1,
784 // const double &x2, const double &y2, const double &z2,
785 // const double &x3, const double &y3, const double &z3,
786 // const double &x4, const double &y4, const double &z4,
787 // fullVector<double> &ncj)
788 //{
789 // // Compute unit vectors for each edge
790 // double x01n, y01n, z01n, x12n, y12n, z12n, x23n, y23n, z23n, x30n, y30n,
791 // z30n; unitVec(x0, y0, z0, x1, y1, z1, x01n, y01n, z01n); unitVec(x1, y1, z1,
792 // x2, y2, z2, x12n, y12n, z12n); unitVec(x2, y2, z2, x3, y3, z3, x23n, y23n,
793 // z23n); unitVec(x3, y3, z3, x0, y0, z0, x30n, y30n, z30n);
794 //
795 // // Compute NCJ at vertex from unit vectors a and b as
796 // 0.5*||a^b||/A_equilateral
797 // // Factor = 2./sqrt(3.) = 0.5/A_equilateral
798 // static const double fact = 1.1547005383792517;
799 // ncj(0) = triArea(fact, x01n, y01n, z01n, -x20n, -y20n, -z20n);
800 // ncj(1) = triArea(fact, x12n, y12n, z12n, -x01n, -y01n, -z01n);
801 // ncj(2) = triArea(fact, x20n, y20n, z20n, -x12n, -y12n, -z12n);
802 //}
803 
804 // TODO: Remove this (useless as quality measure)
805 double qmHexahedron::angles(MHexahedron *el)
806 {
807  double angleMax = 0.0;
808  double angleMin = M_PI;
809  double zeta = 0.0;
810  for(int i = 0; i < el->getNumFaces(); i++) {
811  std::vector<MVertex *> vv;
812  vv.push_back(el->getFace(i).getVertex(0));
813  vv.push_back(el->getFace(i).getVertex(1));
814  vv.push_back(el->getFace(i).getVertex(2));
815  vv.push_back(el->getFace(i).getVertex(3));
816  // MVertex *v0 = new MVertex(0, 0, 0); vv.push_back(v0);
817  // MVertex *v1 = new MVertex(1., 0, 0);vv.push_back(v1);
818  // MVertex *v2 = new MVertex(2., 1., 0);vv.push_back(v2);
819  // MVertex *v3 = new MVertex(1, 1., 0);vv.push_back(v3);
820  for(int j = 0; j < 4; j++) {
821  SVector3 a(vv[(j + 2) % 4]->x() - vv[(j + 1) % 4]->x(),
822  vv[(j + 2) % 4]->y() - vv[(j + 1) % 4]->y(),
823  vv[(j + 2) % 4]->z() - vv[(j + 1) % 4]->z());
824  SVector3 b(vv[(j + 1) % 4]->x() - vv[(j) % 4]->x(),
825  vv[(j + 1) % 4]->y() - vv[(j) % 4]->y(),
826  vv[(j + 1) % 4]->z() - vv[(j) % 4]->z());
827  double angle = acos(dot(a, b) / (norm(a) * norm(b))); //*180/M_PI;
828  angleMax = std::max(angleMax, angle);
829  angleMin = std::min(angleMin, angle);
830  }
831  // printf("angle max =%g min =%g \n", angleMax*180/M_PI, angleMin*180/M_PI);
832  }
833  zeta = 1. - std::max((angleMax - 0.5 * M_PI) / (0.5 * M_PI),
834  (0.5 * M_PI - angleMin) / (0.5 * M_PI));
835  return zeta;
836 }
crossprod
SVector3 crossprod(const SVector3 &a, const SVector3 &b)
Definition: SVector3.h:150
MTriangle.h
qualityMeasures.h
JacobianBasis::getNumPrimMapNodes
int getNumPrimMapNodes() const
Definition: JacobianBasis.h:85
inv3x3
double inv3x3(double mat[3][3], double inv[3][3])
Definition: Numeric.cpp:232
MElement::getPrimaryJacobian
virtual double getPrimaryJacobian(double u, double v, double w, double jac[3][3]) const
Definition: MElement.cpp:939
qmQuadrangle::eta
static double eta(MQuadrangle *el)
Definition: qualityMeasures.cpp:392
SurfaceProjectorUtils::vec2
std::array< double, 2 > vec2
Definition: meshOctreeLibOL.cpp:27
dot
double dot(const SVector3 &a, const SMetric3 &m, const SVector3 &b)
Definition: STensor3.h:71
qmTetrahedron::gamma
static double gamma(const double &x1, const double &y1, const double &z1, const double &x2, const double &y2, const double &z2, const double &x3, const double &y3, const double &z3, const double &x4, const double &y4, const double &z4, double *volume=nullptr)
Definition: qualityMeasures.cpp:674
MTetrahedron
Definition: MTetrahedron.h:34
MQuadrangle::getVertex
virtual MVertex * getVertex(int num)
Definition: MQuadrangle.h:64
angle
double angle(const SVector3 &a, const SVector3 &b)
Definition: SVector3.h:157
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
MVertex
Definition: MVertex.h:24
a4
#define a4
Definition: GaussQuadratureTet.cpp:11
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
LegendrePolynomials::f
void f(int n, double u, double *val)
Definition: orthogonalBasis.cpp:77
SVector3
Definition: SVector3.h:16
qmQuadrangle::NCJRange
static void NCJRange(const MQuadrangle *e, double &valMin, double &valMax)
Definition: qualityMeasures.cpp:493
MTriangle::getVertex
virtual MVertex * getVertex(int num)
Definition: MTriangle.h:62
MPrism
Definition: MPrism.h:34
MVertex::point
SPoint3 point() const
Definition: MVertex.h:67
MHexahedron::getFace
virtual MFace getFace(int num) const
Definition: MHexahedron.h:92
GmshMessage.h
qmPrism::minNCJ
static double minNCJ(const MPrism *el)
Definition: qualityMeasures.cpp:769
qmTriangle::NCJ
static void NCJ(const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2, const SVector3 &normal, std::vector< double > &NCJ)
Definition: qualityMeasures.cpp:290
qmTriangle::NCJAndGradients
static void NCJAndGradients(const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2, const SVector3 &normal, std::vector< double > &NCJ, std::vector< double > &dNCJ)
Definition: qualityMeasures.cpp:310
qmQuadrangle::NCJAndGradients
static void NCJAndGradients(const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2, const SPoint3 &p3, const SVector3 &normal, std::vector< double > &NCJ, std::vector< double > &dNCJ)
Definition: qualityMeasures.cpp:534
norm3
double norm3(double a[3])
Definition: Numeric.h:119
qmTetrahedron::eta
static double eta(const double &x1, const double &y1, const double &z1, const double &x2, const double &y2, const double &z2, const double &x3, const double &y3, const double &z3, const double &x4, const double &y4, const double &z4, double *volume=nullptr)
Definition: qualityMeasures.cpp:651
fullMatrix< double >
BDS_Point::X
double X
Definition: BDS.h:56
MPrism::getVertex
virtual MVertex * getVertex(int num)
Definition: MPrism.h:71
MHexahedron.h
MVertex.h
qmTriangle::gamma
static double gamma(MTriangle *f)
Definition: qualityMeasures.cpp:146
matmat
void matmat(double mat1[3][3], double mat2[3][3], double res[3][3])
Definition: Numeric.cpp:54
a1
const double a1
Definition: GaussQuadratureHex.cpp:10
MTetrahedron::getVertex
virtual MVertex * getVertex(int num)
Definition: MTetrahedron.h:67
qmTetrahedron::QMTET_ONE
@ QMTET_ONE
Definition: qualityMeasures.h:68
norm
void norm(const double *vec, double *norm)
Definition: gmshLevelset.cpp:202
Numeric.h
qmQuadrangle::NCJ
static void NCJ(const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2, const SPoint3 &p3, const SVector3 &normal, std::vector< double > &ncj)
Definition: qualityMeasures.cpp:515
qmTriangle::NCJRange
static void NCJRange(const MTriangle *e, double &valMin, double &valMax)
Definition: qualityMeasures.cpp:269
qmTetrahedron::qm
static double qm(MTetrahedron *t, const Measures &cr, double *volume=nullptr)
Definition: qualityMeasures.cpp:616
MFace::getVertex
MVertex * getVertex(std::size_t i) const
Definition: MFace.h:30
SVector3::norm
double norm() const
Definition: SVector3.h:33
qmTetrahedron::QMTET_ETA
@ QMTET_ETA
Definition: qualityMeasures.h:68
MHexahedron
Definition: MHexahedron.h:28
prosca
double prosca(double const a[3], double const b[3])
Definition: Numeric.h:112
JacobianBasis.h
angle3Vertices
double angle3Vertices(const MVertex *p1, const MVertex *p2, const MVertex *p3)
Definition: MVertex.cpp:16
BDS.h
prodve
void prodve(double a[3], double b[3], double c[3])
Definition: Numeric.h:105
a2
const double a2
Definition: GaussQuadratureHex.cpp:12
prismNCJ
static double prismNCJ(const MVertex *a, const MVertex *b, const MVertex *c, const MVertex *d)
Definition: qualityMeasures.cpp:749
norm2
double norm2(double a[3][3])
Definition: Numeric.cpp:38
triangle_area
double triangle_area(double p0[3], double p1[3], double p2[3])
Definition: Numeric.cpp:293
SurfaceProjectorUtils::vec3
std::array< double, 3 > vec3
Definition: meshOctreeLibOL.cpp:28
S
#define S
Definition: DefaultOptions.h:21
MTriangle
Definition: MTriangle.h:26
JacobianBasis::getPrimNormal2D
double getPrimNormal2D(const fullMatrix< double > &nodesXYZ, fullMatrix< double > &result, bool ideal=false) const
Definition: JacobianBasis.cpp:385
MHexahedron::getNumFaces
virtual int getNumFaces()
Definition: MHexahedron.h:91
BDS_Face
Definition: BDS.h:164
BDS_Face::getNodes
bool getNodes(BDS_Point *_n[4]) const
Definition: BDS.h:201
polynomialBasis.h
qmQuadrangle::angles
static double angles(MQuadrangle *e)
Definition: qualityMeasures.cpp:435
MTetrahedron.h
qmHexahedron::angles
static double angles(MHexahedron *el)
Definition: qualityMeasures.cpp:805
z
const double z
Definition: GaussQuadratureQuad.cpp:56
JacobianBasis
Definition: JacobianBasis.h:60
MQuadrangle.h
MElement::getNumPrimaryVertices
std::size_t getNumPrimaryVertices() const
Definition: MElement.h:160
BDS_Point::Y
double Y
Definition: BDS.h:56
MElement::getJacobianFuncSpace
virtual const JacobianBasis * getJacobianFuncSpace(int orderElement=-1) const
Definition: MElement.cpp:679
qmTriangle::angles
static double angles(MTriangle *e)
Definition: qualityMeasures.cpp:200
MPrism.h
qmTetrahedron::cond
static double cond(const double &x1, const double &y1, const double &z1, const double &x2, const double &y2, const double &z2, const double &x3, const double &y3, const double &z3, const double &x4, const double &y4, const double &z4, double *volume=nullptr)
Definition: qualityMeasures.cpp:727
qmTetrahedron::QMTET_GAMMA
@ QMTET_GAMMA
Definition: qualityMeasures.h:68
qmTriangle::eta
static double eta(MTriangle *el)
Definition: qualityMeasures.cpp:187
MElement::getJacobian
virtual double getJacobian(const fullMatrix< double > &gsf, double jac[3][3]) const
Definition: MElement.cpp:868
robustPredicates::orient3d
REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
Definition: robustPredicates.cpp:2321
BDS_Point
Definition: BDS.h:49
MVertex::y
double y() const
Definition: MVertex.h:61
norme
double norme(double a[3])
Definition: Numeric.h:123
qmTetrahedron::Measures
Measures
Definition: qualityMeasures.h:68
MQuadrangle
Definition: MQuadrangle.h:26
qmTetrahedron::QMTET_COND
@ QMTET_COND
Definition: qualityMeasures.h:68
MVertex::x
double x() const
Definition: MVertex.h:60
SVector3::normalize
double normalize()
Definition: SVector3.h:38
BDS_Point::Z
double Z
Definition: BDS.h:56