From 1674ef93c99efcd199e49eb15203b0e3b4a1d2e1 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Tue, 22 Jan 2008 07:48:12 +0000
Subject: [PATCH] *** empty log message ***
---
Mesh/qualityMeasures.cpp | 210 ++++++++++++++++++++-------------------
Mesh/qualityMeasures.h | 36 ++++---
2 files changed, 126 insertions(+), 120 deletions(-)
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 5f6e2a9e32..4b94e87053 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -4,141 +4,143 @@
#include "MElement.h"
#include "Numeric.h"
-
-double qmTriangle ( const BDS_Point *p1, const BDS_Point *p2, const BDS_Point *p3, const gmshQualityMeasure4Triangle &cr)
+double qmTriangle(const BDS_Point *p1, const BDS_Point *p2, const BDS_Point *p3,
+ const gmshQualityMeasure4Triangle &cr)
{
- return qmTriangle (p1->X,p1->Y,p1->Z,p2->X,p2->Y,p2->Z,p3->X,p3->Y,p3->Z,cr);
+ return qmTriangle(p1->X, p1->Y, p1->Z, p2->X, p2->Y, p2->Z, p3->X, p3->Y, p3->Z, cr);
}
-double qmTriangle ( BDS_Face *t, const gmshQualityMeasure4Triangle &cr)
+double qmTriangle(BDS_Face *t, const gmshQualityMeasure4Triangle &cr)
{
BDS_Point *n[4];
t->getNodes(n);
- return qmTriangle (n[0],n[1],n[2],cr);
+ return qmTriangle(n[0], n[1], n[2], cr);
}
-double qmTriangle ( MTriangle*t, const gmshQualityMeasure4Triangle &cr)
+double qmTriangle(MTriangle*t, const gmshQualityMeasure4Triangle &cr)
{
- return qmTriangle (t->getVertex(0),t->getVertex(1),t->getVertex(2),cr);
+ return qmTriangle(t->getVertex(0), t->getVertex(1), t->getVertex(2), cr);
}
-double qmTriangle ( const MVertex *v1, const MVertex *v2, const MVertex *v3, const gmshQualityMeasure4Triangle &cr)
+double qmTriangle(const MVertex *v1, const MVertex *v2, const MVertex *v3,
+ const gmshQualityMeasure4Triangle &cr)
{
- return qmTriangle (v1->x(),v1->y(),v1->z(),v2->x(),v2->y(),v2->z(),v3->x(),v3->y(),v3->z(),cr);
+ return qmTriangle(v1->x(), v1->y(), v1->z(), v2->x(), v2->y(), v2->z(),
+ v3->x(), v3->y(), v3->z(), cr);
}
// Triangle abc
// quality is between 0 and 1
-double qmTriangle ( const double &xa, const double &ya, const double &za,
- const double &xb, const double &yb, const double &zb,
- const double &xc, const double &yc, const double &zc,
- const gmshQualityMeasure4Triangle &cr)
+double qmTriangle(const double &xa, const double &ya, const double &za,
+ const double &xb, const double &yb, const double &zb,
+ const double &xc, const double &yc, const double &zc,
+ const gmshQualityMeasure4Triangle &cr)
{
double quality;
- switch(cr)
+ switch(cr){
+ case QMTRI_RHO:
{
- case QMTRI_RHO:
- {
- // quality = rho / R = 2 * inscribed radius / circumradius
-
- double a [3] = {xc-xb,yc-yb,zc-zb};
- double b [3] = {xa-xc,ya-yc,za-zc};
- double c [3] = {xb-xa,yb-ya,zb-za};
-
- const double la = norme (a);
- const double lb = norme (b);
- const double lc = norme (c);
-
- double pva [3]; prodve (b,c,pva); const double sina = norm3 (pva);
- double pvb [3]; prodve (c,a,pvb); const double sinb = norm3 (pvb);
- double pvc [3]; prodve (a,b,pvc); const double sinc = norm3 (pvc);
-
- if (sina == 0.0 && sinb == 0.0 && sinc == 0.0) quality = 0.0;
- else quality = 2 * (2*sina*sinb*sinc/(sina + sinb + sinc) );
- }
- break;
- default:
- throw;
- }
-
+ // quality = rho / R = 2 * inscribed radius / circumradius
+ double a [3] = {xc-xb,yc-yb,zc-zb};
+ double b [3] = {xa-xc,ya-yc,za-zc};
+ double c [3] = {xb-xa,yb-ya,zb-za};
+
+ const double la = norme (a);
+ const double lb = norme (b);
+ const double lc = norme (c);
+
+ double pva [3]; prodve (b,c,pva); const double sina = norm3 (pva);
+ double pvb [3]; prodve (c,a,pvb); const double sinb = norm3 (pvb);
+ double pvc [3]; prodve (a,b,pvc); const double sinc = norm3 (pvc);
+
+ if (sina == 0.0 && sinb == 0.0 && sinc == 0.0) quality = 0.0;
+ else quality = 2 * (2*sina*sinb*sinc/(sina + sinb + sinc) );
+ }
+ break;
+ default:
+ throw;
+ }
return quality;
}
-double qmTet ( MTetrahedron*t, const gmshQualityMeasure4Tet &cr, double *volume)
+double qmTet(MTetrahedron *t, const gmshQualityMeasure4Tet &cr, double *volume)
{
- return qmTet (t->getVertex(0),t->getVertex(1),t->getVertex(2),t->getVertex(3),cr,volume);
+ return qmTet(t->getVertex(0), t->getVertex(1), t->getVertex(2), t->getVertex(3),
+ cr, volume);
}
-double qmTet ( const MVertex *v1, const MVertex *v2, const MVertex *v3, const MVertex *v4, const gmshQualityMeasure4Tet &cr, double *volume)
+double qmTet(const MVertex *v1, const MVertex *v2, const MVertex *v3,
+ const MVertex *v4, const gmshQualityMeasure4Tet &cr, double *volume)
{
- return qmTet (v1->x(),v1->y(),v1->z(),v2->x(),v2->y(),v2->z(),v3->x(),v3->y(),v3->z(),v4->x(),v4->y(),v4->z(),cr,volume);
+ return qmTet(v1->x(), v1->y(), v1->z(), v2->x(), v2->y(), v2->z(),
+ v3->x(), v3->y(), v3->z(), v4->x(), v4->y(), v4->z(), cr, volume);
}
-double qmTet ( 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,
- const gmshQualityMeasure4Tet &cr, double *volume ){
+double qmTet(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,
+ const gmshQualityMeasure4Tet &cr, double *volume)
+{
double quality;
- switch(cr)
+ switch(cr){
+ case QMTET_ONE:
+ return 1.0;
+ case QMTET_3:
+ {
+ double mat[3][3];
+ mat[0][0] = x2 - x1;
+ mat[0][1] = x3 - x1;
+ mat[0][2] = x4 - x1;
+ mat[1][0] = y2 - y1;
+ mat[1][1] = y3 - y1;
+ mat[1][2] = y4 - y1;
+ mat[2][0] = z2 - z1;
+ mat[2][1] = z3 - z1;
+ mat[2][2] = z4 - z1;
+ *volume = fabs(det3x3(mat)) / 6.;
+ double l = ((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
+ l += ((x3-x1)*(x3-x1)+(y3-y1)*(y3-y1)+(z3-z1)*(z3-z1));
+ l += ((x4-x1)*(x4-x1)+(y4-y1)*(y4-y1)+(z4-z1)*(z4-z1));
+ l += ((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2)+(z3-z2)*(z3-z2));
+ l += ((x4-x2)*(x4-x2)+(y4-y2)*(y4-y2)+(z4-z2)*(z4-z2));
+ l += ((x3-x4)*(x3-x4)+(y3-y4)*(y3-y4)+(z3-z4)*(z3-z4));
+ return 12. * pow(3*fabs(*volume),2./3.)/l;
+ }
+ case QMTET_2:
{
- case QMTET_ONE:
- return 1.0;
- case QMTET_3:
- {
- double mat[3][3];
- mat[0][0] = x2 - x1;
- mat[0][1] = x3 - x1;
- mat[0][2] = x4 - x1;
- mat[1][0] = y2 - y1;
- mat[1][1] = y3 - y1;
- mat[1][2] = y4 - y1;
- mat[2][0] = z2 - z1;
- mat[2][1] = z3 - z1;
- mat[2][2] = z4 - z1;
- *volume = fabs(det3x3(mat)) / 6.;
- double l = ((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
- l += ((x3-x1)*(x3-x1)+(y3-y1)*(y3-y1)+(z3-z1)*(z3-z1));
- l +=((x4-x1)*(x4-x1)+(y4-y1)*(y4-y1)+(z4-z1)*(z4-z1));
- l +=((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2)+(z3-z2)*(z3-z2));
- l +=((x4-x2)*(x4-x2)+(y4-y2)*(y4-y2)+(z4-z2)*(z4-z2));
- l +=((x3-x4)*(x3-x4)+(y3-y4)*(y3-y4)+(z3-z4)*(z3-z4));
- return 12. * pow(3*fabs(*volume),2./3.)/l;
- }
- case QMTET_2:
- {
- double mat[3][3];
- mat[0][0] = x2 - x1;
- mat[0][1] = x3 - x1;
- mat[0][2] = x4 - x1;
- mat[1][0] = y2 - y1;
- mat[1][1] = y3 - y1;
- mat[1][2] = y4 - y1;
- mat[2][0] = z2 - z1;
- mat[2][1] = z3 - z1;
- mat[2][2] = z4 - z1;
- *volume = fabs(det3x3(mat)) / 6.;
- double p0[3] = { x1,y1,z1};
- double p1[3] = { x2,y2,z2};
- double p2[3] = { x3,y3,z3};
- double p3[3] = { x4,y4,z4};
- double s1 = fabs(triangle_area(p0, p1, p2));
- double s2 = fabs(triangle_area(p0, p2, p3));
- double s3 = fabs(triangle_area(p0, p1, p3));
- double s4 = fabs(triangle_area(p1, p2, p3));
- double rhoin = 3. * fabs(*volume) / (s1 + s2 + s3 + s4);
- double l = sqrt ((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
- l = std::max(l,sqrt ((x3-x1)*(x3-x1)+(y3-y1)*(y3-y1)+(z3-z1)*(z3-z1)));
- l = std::max(l,sqrt ((x4-x1)*(x4-x1)+(y4-y1)*(y4-y1)+(z4-z1)*(z4-z1)));
- l = std::max(l,sqrt ((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2)+(z3-z2)*(z3-z2)));
- l = std::max(l,sqrt ((x4-x2)*(x4-x2)+(y4-y2)*(y4-y2)+(z4-z2)*(z4-z2)));
- l = std::max(l,sqrt ((x3-x4)*(x3-x4)+(y3-y4)*(y3-y4)+(z3-z4)*(z3-z4)));
- return 2. * sqrt(6.) * rhoin / l;
- }
- break;
- default:
- throw;
+ double mat[3][3];
+ mat[0][0] = x2 - x1;
+ mat[0][1] = x3 - x1;
+ mat[0][2] = x4 - x1;
+ mat[1][0] = y2 - y1;
+ mat[1][1] = y3 - y1;
+ mat[1][2] = y4 - y1;
+ mat[2][0] = z2 - z1;
+ mat[2][1] = z3 - z1;
+ mat[2][2] = z4 - z1;
+ *volume = fabs(det3x3(mat)) / 6.;
+ double p0[3] = {x1,y1,z1};
+ double p1[3] = {x2,y2,z2};
+ double p2[3] = {x3,y3,z3};
+ double p3[3] = {x4,y4,z4};
+ double s1 = fabs(triangle_area(p0, p1, p2));
+ double s2 = fabs(triangle_area(p0, p2, p3));
+ double s3 = fabs(triangle_area(p0, p1, p3));
+ double s4 = fabs(triangle_area(p1, p2, p3));
+ double rhoin = 3. * fabs(*volume) / (s1 + s2 + s3 + s4);
+ double l = sqrt ((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
+ l = std::max(l,sqrt ((x3-x1)*(x3-x1)+(y3-y1)*(y3-y1)+(z3-z1)*(z3-z1)));
+ l = std::max(l,sqrt ((x4-x1)*(x4-x1)+(y4-y1)*(y4-y1)+(z4-z1)*(z4-z1)));
+ l = std::max(l,sqrt ((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2)+(z3-z2)*(z3-z2)));
+ l = std::max(l,sqrt ((x4-x2)*(x4-x2)+(y4-y2)*(y4-y2)+(z4-z2)*(z4-z2)));
+ l = std::max(l,sqrt ((x3-x4)*(x3-x4)+(y3-y4)*(y3-y4)+(z3-z4)*(z3-z4)));
+ return 2. * sqrt(6.) * rhoin / l;
}
+ break;
+ default:
+ throw;
+ }
}
diff --git a/Mesh/qualityMeasures.h b/Mesh/qualityMeasures.h
index 96d9078bd0..335a74437d 100644
--- a/Mesh/qualityMeasures.h
+++ b/Mesh/qualityMeasures.h
@@ -9,21 +9,25 @@ class MTetrahedron;
enum gmshQualityMeasure4Triangle {QMTRI_RHO};
enum gmshQualityMeasure4Tet {QMTET_1,QMTET_2,QMTET_3,QMTET_ONE};
-double qmTriangle ( MTriangle *f, const gmshQualityMeasure4Triangle &cr);
-double qmTriangle ( BDS_Face *f, const gmshQualityMeasure4Triangle &cr);
-double qmTriangle ( const BDS_Point *p1, const BDS_Point *p2, const BDS_Point *p3, const gmshQualityMeasure4Triangle &cr);
-double qmTriangle ( const MVertex *v1, const MVertex *v2, const MVertex *v3, const gmshQualityMeasure4Triangle &cr);
-double qmTriangle ( const double *d1, const double *d2, const double *d3, const gmshQualityMeasure4Triangle &cr);
-double qmTriangle ( 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 gmshQualityMeasure4Triangle &cr);
-double qmTet ( MTetrahedron *t, const gmshQualityMeasure4Tet &cr, double *volume = 0);
-double qmTet ( const MVertex *v1, const MVertex *v2, const MVertex *v3, const MVertex *v4, const gmshQualityMeasure4Tet &cr, double *volume = 0);
-double qmTet ( 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,
- const gmshQualityMeasure4Tet &cr, double *volume = 0);
+double qmTriangle(MTriangle *f, const gmshQualityMeasure4Triangle &cr);
+double qmTriangle(BDS_Face *f, const gmshQualityMeasure4Triangle &cr);
+double qmTriangle(const BDS_Point *p1, const BDS_Point *p2, const BDS_Point *p3,
+ const gmshQualityMeasure4Triangle &cr);
+double qmTriangle(const MVertex *v1, const MVertex *v2, const MVertex *v3,
+ const gmshQualityMeasure4Triangle &cr);
+double qmTriangle(const double *d1, const double *d2, const double *d3,
+ const gmshQualityMeasure4Triangle &cr);
+double qmTriangle(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 gmshQualityMeasure4Triangle &cr);
+double qmTet(MTetrahedron *t, const gmshQualityMeasure4Tet &cr, double *volume = 0);
+double qmTet(const MVertex *v1, const MVertex *v2, const MVertex *v3,
+ const MVertex *v4, const gmshQualityMeasure4Tet &cr, double *volume = 0);
+double qmTet(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,
+ const gmshQualityMeasure4Tet &cr, double *volume = 0);
#endif
--
GitLab