diff --git a/Geo/MHexahedron.cpp b/Geo/MHexahedron.cpp
index 31f5e74d7184e8778d4aa5280d7c06a7a0f97662..39aa47f33d3c4fc36d9187169516db21ee041b1c 100644
--- a/Geo/MHexahedron.cpp
+++ b/Geo/MHexahedron.cpp
@@ -10,7 +10,10 @@
#include "nodalBasis.h"
#include "polynomialBasis.h"
#include "MQuadrangle.h"
+
+#if defined(HAVE_MESH)
#include "qualityMeasures.h"
+#endif
int MHexahedron::getVolumeSign()
{
@@ -38,32 +41,8 @@ void MHexahedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
double MHexahedron::angleShapeMeasure()
{
-
#if defined(HAVE_MESH)
- double angleMax = 0.0;
- double angleMin = M_PI;
- double zeta = 0.0;
- for (int i=0; i<getNumFaces(); i++){
- std::vector<MVertex*> vv;
- vv.push_back(getFace(i).getVertex(0));
- vv.push_back(getFace(i).getVertex(1));
- vv.push_back(getFace(i).getVertex(2));
- vv.push_back(getFace(i).getVertex(3));
- // MVertex *v0 = new MVertex(0, 0, 0); vv.push_back(v0);
- // MVertex *v1 = new MVertex(1., 0, 0);vv.push_back(v1);
- // MVertex *v2 = new MVertex(2., 1., 0);vv.push_back(v2);
- // MVertex *v3 = new MVertex(1, 1., 0);vv.push_back(v3);
- for (int j=0; j<4; j++){
- SVector3 a(vv[(j+2)%4]->x()-vv[(j+1)%4]->x(),vv[(j+2)%4]->y()-vv[(j+1)%4]->y(),vv[(j+2)%4]->z()-vv[(j+1)%4]->z() );
- SVector3 b(vv[(j+1)%4]->x()-vv[(j)%4]->x(), vv[(j+1)%4]->y()-vv[(j)%4]->y(), vv[(j+1)%4]->z()-vv[(j)%4]->z() );
- double angle = acos( dot(a,b)/(norm(a)*norm(b))); //*180/M_PI;
- angleMax = std::max(angleMax, angle);
- angleMin = std::min(angleMin, angle);
- }
- //printf("angle max =%g min =%g \n", angleMax*180/M_PI, angleMin*180/M_PI);
- }
- zeta = 1.-std::max((angleMax-0.5*M_PI)/(0.5*M_PI),(0.5*M_PI-angleMin)/(0.5*M_PI));
- return zeta;
+ return qmHexahedron::angles(this);
#else
return 1.;
#endif
diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index 41c2b69ef4bb9de236a45a7c94b8c3198a45d6e1..e1dcf182b6291fc57f74ba2af2614329b448e001 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -8,6 +8,10 @@
#include "BasisFactory.h"
#include "Context.h"
+#if defined(HAVE_MESH)
+#include "qualityMeasures.h"
+#endif
+
int MPrism::getVolumeSign()
{
double mat[3][3];
@@ -517,35 +521,8 @@ void MPrismN::getFaceVertices(const int num, std::vector<MVertex*> &v) const
}
-static double scaled_jacobian(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
- double val;
- double l1,l2,l3;
- SVector3 vec1,vec2,vec3;
-
- vec1 = SVector3(b->x()-a->x(),b->y()-a->y(),b->z()-a->z());
- vec2 = SVector3(c->x()-a->x(),c->y()-a->y(),c->z()-a->z());
- vec3 = SVector3(d->x()-a->x(),d->y()-a->y(),d->z()-a->z());
-
- l1 = vec1.norm();
- l2 = vec2.norm();
- l3 = vec3.norm();
-
- val = dot(vec1,crossprod(vec2,vec3));
- return fabs(val)/(l1*l2*l3);
-}
-
-double MPrism::gammaShapeMeasure() {
- MVertex *a = getVertex(0),*b= getVertex(1),*c= getVertex(2);
- MVertex *d = getVertex(3),*e= getVertex(4),*f= getVertex(5);
- const double j [6] = {
- scaled_jacobian(a,b,c,d),
- scaled_jacobian(b,a,c,e),
- scaled_jacobian(c,a,b,f),
- scaled_jacobian(d,a,e,f),
- scaled_jacobian(e,b,d,f),
- scaled_jacobian(f,c,d,e)};
- const double result = *std::min_element(j,j+6);
- printf("%12.5E\n",result);
- return result;
+double MPrism::gammaShapeMeasure()
+{
+ return qmPrism::minNCJ(this);
}
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 56aa05f5dc09f9fdbebd5ae6078bac686caa3611..6d5e98e26c212d194bc5cb23c89943f05ff5afe4 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -186,56 +186,26 @@ void MQuadrangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
double MQuadrangle::etaShapeMeasure()
{
- double AR = 1;//(minEdge()/maxEdge());
-
- SVector3 v01 (_v[1]->x()-_v[0]->x(),_v[1]->y()-_v[0]->y(),_v[1]->z()-_v[0]->z());
- SVector3 v12 (_v[2]->x()-_v[1]->x(),_v[2]->y()-_v[1]->y(),_v[2]->z()-_v[1]->z());
- SVector3 v23 (_v[3]->x()-_v[2]->x(),_v[3]->y()-_v[2]->y(),_v[3]->z()-_v[2]->z());
- SVector3 v30 (_v[0]->x()-_v[3]->x(),_v[0]->y()-_v[3]->y(),_v[0]->z()-_v[3]->z());
-
- SVector3 a = crossprod(v01,v12);
- SVector3 b = crossprod(v12,v23);
- SVector3 c = crossprod(v23,v30);
- SVector3 d = crossprod(v30,v01);
-
- double sign = 1.0;
- if (dot(a,b) < 0 || dot(a,c) < 0 || dot(a,d) < 0 )sign = -1;
- // FIXME ...
- // if (a.z() > 0 || b.z() > 0 || c.z() > 0 || d.z() > 0) sign = -1;
-
- double a1 = 180 * angle3Vertices(_v[0], _v[1], _v[2]) / M_PI;
- double a2 = 180 * angle3Vertices(_v[1], _v[2], _v[3]) / M_PI;
- double a3 = 180 * angle3Vertices(_v[2], _v[3], _v[0]) / M_PI;
- double a4 = 180 * angle3Vertices(_v[3], _v[0], _v[1]) / M_PI;
-
- a1 = std::min(180.,a1);
- a2 = std::min(180.,a2);
- a3 = std::min(180.,a3);
- a4 = std::min(180.,a4);
- double angle = fabs(90. - a1);
- angle = std::max(fabs(90. - a2),angle);
- angle = std::max(fabs(90. - a3),angle);
- angle = std::max(fabs(90. - a4),angle);
-
- return sign*(1.-angle/90) * AR;
+#if defined(HAVE_MESH)
+ return qmQuadrangle::eta(this);
+#else
+ return 0.;
+#endif
}
-/// a shape measure for quadrangles
-/// assume (for now) 2D elements --
-/// sf = (1 \pm xi) (1 \pm eta) / 4
-/// dsf_xi = \pm (1 \pm eta) / 4
-/// 1 + eta , -(1+eta) , -(1-eta), 1-eta
-/// dsf_eta = \pm (1 \pm xi) / 4
-/// 1 + xi , 1 - xi , -(1-xi), -(1+xi)
double MQuadrangle::gammaShapeMeasure(){
- return etaShapeMeasure();
+#if defined(HAVE_MESH)
+ return qmQuadrangle::gamma(this);
+#else
+ return 0.;
+#endif
}
double MQuadrangle::angleShapeMeasure()
{
#if defined(HAVE_MESH)
- return qmQuadrangleAngles(this);
+ return qmQuadrangle::angles(this);
#else
return 1.;
#endif
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index 43e5587f546a76db4baf1b673e7f30ce6c19b9ae..ea2c49892c59a2aea82b971e65dc1c77f9f2afe3 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -67,7 +67,7 @@ double MTetrahedron::gammaShapeMeasure()
{
#if defined(HAVE_MESH)
double vol;
- return qmTet(this, QMTET_2, &vol);
+ return qmTetrahedron::qm(this, qmTetrahedron::QMTET_GAMMA, &vol);
#else
return 0.;
#endif
@@ -77,7 +77,7 @@ double MTetrahedron::etaShapeMeasure()
{
#if defined(HAVE_MESH)
double vol;
- return qmTet(this, QMTET_3, &vol);
+ return qmTetrahedron::qm(this, qmTetrahedron::QMTET_ETA, &vol);
#else
return 0.;
#endif
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index 587d7a1a520771ad273ac22aecfd4795d94e1169..5c21e2e54045f4bc683162440c1f6ec2d1517806 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -65,7 +65,7 @@ double MTriangle::getOuterRadius()
double MTriangle::angleShapeMeasure()
{
#if defined(HAVE_MESH)
- return qmTriangleAngles(this);
+ return qmTriangle::angles(this);
#else
return 0.;
#endif
@@ -73,20 +73,18 @@ double MTriangle::angleShapeMeasure()
double MTriangle::etaShapeMeasure()
{
- double a1 = 180 * angle3Vertices(_v[0], _v[1], _v[2]) / M_PI;
- double a2 = 180 * angle3Vertices(_v[1], _v[2], _v[0]) / M_PI;
- double a3 = 180 * angle3Vertices(_v[2], _v[0], _v[1]) / M_PI;
-
- double amin = std::min(std::min(a1,a2),a3);
- double angle = fabs(60. - amin);
- return 1.-angle/60;
+#if defined(HAVE_MESH)
+ return qmTriangle::eta(this);
+#else
+ return 0.;
+#endif
}
double MTriangle::gammaShapeMeasure()
{
#if defined(HAVE_MESH)
- return qmTriangle(this, QMTRI_RHO);
+ return qmTriangle::gamma(this);
#else
return 0.;
#endif
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 357759af15e5d97c7b864e2eab9787b7eb44b1e9..af1c2137eb51ff25b623eb286a140e0eb623dda4 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -779,10 +779,10 @@ bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1, BDS_Point *_p2, BDS_P
double cosnq; prosca(n, q, &cosnq);
double cosonq; prosca(on, oq, &cosonq);
- double qa1 = qmTriangle(_p1, _p2, _p3, QMTRI_RHO);
- double qa2 = qmTriangle(_q1, _q2, _q3, QMTRI_RHO);
- double qb1 = qmTriangle(_op1, _op2, _op3, QMTRI_RHO);
- double qb2 = qmTriangle(_oq1, _oq2, _oq3, QMTRI_RHO);
+ double qa1 = qmTriangle::gamma(_p1, _p2, _p3);
+ double qa2 = qmTriangle::gamma(_q1, _q2, _q3);
+ double qb1 = qmTriangle::gamma(_op1, _op2, _op3);
+ double qb2 = qmTriangle::gamma(_oq1, _oq2, _oq3);
// we swap for a better configuration
double mina = std::min(qa1,qa2);
@@ -1113,8 +1113,8 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge *e, BDS_Point *p)
pt[1][nt] = (pts[1] == p) ? o : pts[1];
pt[2][nt] = (pts[2] == p) ? o : pts[2];
-// double qnew = qmTriangle(pt[0][nt], pt[1][nt], pt[2][nt], QMTRI_RHO);
-// double qold = qmTriangle(pts[0], pts[1], pts[2], QMTRI_RHO);
+// double qnew = qmTriangle::gamma(pt[0][nt], pt[1][nt], pt[2][nt]);
+// double qold = qmTriangle::gamma(pts[0], pts[1], pts[2]);
// if(qold > 1.e-4 && qnew < 1.e-4) return false;
nt++;
// pt[0][nt] = (pts[0] == p) ? o->iD : pts[0]->iD;
@@ -1263,14 +1263,14 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
p->X = gp.x();
p->Y = gp.y();
p->Z = gp.z();
- newWorst = std::min(newWorst, qmTriangle(*it, QMTRI_RHO));
+ newWorst = std::min(newWorst, qmTriangle::gamma(*it));
double norm1[3],norm2[3];
normal_triangle(n[0], n[1], n[2], norm1);
p->X = oldX;
p->Y = oldY;
p->Z = oldZ;
normal_triangle(n[0], n[1], n[2], norm2);
- oldWorst = std::min(oldWorst, qmTriangle(*it, QMTRI_RHO));
+ oldWorst = std::min(oldWorst, qmTriangle::gamma(*it));
double ps;
prosca(norm1, norm2, &ps);
double threshold = (isSphere ? 0.95 : 0.5);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 5dd68634fa5ac85e30ed3697da7c7de45c1f1323..cac622ed26719fa285bb638f6d3e50e560a7e5ff 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -58,7 +58,7 @@ static void computeElementShapes(GFace *gf, double &worst, double &avg,
nT = 0;
greaterThan = 0;
for(unsigned int i = 0; i < gf->triangles.size(); i++){
- double q = qmTriangle(gf->triangles[i], QMTRI_RHO);
+ double q = qmTriangle::gamma(gf->triangles[i]);
if(q > .9) greaterThan++;
avg += q;
worst = std::min(worst, q);
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index c96846ea13412dea12be2460bf360be6ca5be156..0103985d009b6c32a3f6d6eb42be4d2e7895e1a0 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -234,10 +234,10 @@ bool evalSwapForOptimize(BDS_Edge *e, GFace *gf, BDS_Mesh &m)
// First, evaluate what we gain in element quality if the
// swap is performed
- double qa1 = qmTriangle(p11, p12, p13, QMTRI_RHO);
- double qa2 = qmTriangle(p21, p22, p23, QMTRI_RHO);
- double qb1 = qmTriangle(p31, p32, p33, QMTRI_RHO);
- double qb2 = qmTriangle(p41, p42, p43, QMTRI_RHO);
+ double qa1 = qmTriangle::gamma(p11, p12, p13);
+ double qa2 = qmTriangle::gamma(p21, p22, p23);
+ double qb1 = qmTriangle::gamma(p31, p32, p33);
+ double qb2 = qmTriangle::gamma(p41, p42, p43);
double qa = std::min(qa1, qa2);
double qb = std::min(qb1, qb2);
double qualIndicator = qb - qa;
@@ -357,10 +357,10 @@ bool evalSwap(BDS_Edge *e, double &qa, double &qb)
if(e->numfaces() != 2) return false;
e->oppositeof (op);
- double qa1 = qmTriangle(e->p1, e->p2, op[0], QMTRI_RHO);
- double qa2 = qmTriangle(e->p1, e->p2, op[1], QMTRI_RHO);
- double qb1 = qmTriangle(e->p1, op[0], op[1], QMTRI_RHO);
- double qb2 = qmTriangle(e->p2, op[0], op[1], QMTRI_RHO);
+ double qa1 = qmTriangle::gamma(e->p1, e->p2, op[0]);
+ double qa2 = qmTriangle::gamma(e->p1, e->p2, op[1]);
+ double qb1 = qmTriangle::gamma(e->p1, op[0], op[1]);
+ double qb2 = qmTriangle::gamma(e->p2, op[0], op[1]);
qa = std::min(qa1, qa2);
qb = std::min(qb1, qb2);
return true;
@@ -381,10 +381,10 @@ int edgeSwapTestQuality(BDS_Edge *e, double fact=1.1, bool force=false)
if (!edgeSwapTestAngle(e, cos(CTX::instance()->mesh.allowSwapEdgeAngle * M_PI / 180.)))
return -1;
- double qa1 = qmTriangle(e->p1, e->p2, op[0], QMTRI_RHO);
- double qa2 = qmTriangle(e->p1, e->p2, op[1], QMTRI_RHO);
- double qb1 = qmTriangle(e->p1, op[0], op[1], QMTRI_RHO);
- double qb2 = qmTriangle(e->p2, op[0], op[1], QMTRI_RHO);
+ double qa1 = qmTriangle::gamma(e->p1, e->p2, op[0]);
+ double qa2 = qmTriangle::gamma(e->p1, e->p2, op[1]);
+ double qb1 = qmTriangle::gamma(e->p1, op[0], op[1]);
+ double qb2 = qmTriangle::gamma(e->p2, op[0], op[1]);
double qa = std::min(qa1, qa2);
double qb = std::min(qb1, qb2);
if(qb > fact * qa) return 1;
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index d778a94a2a84233f179a0c20ad3bd7b1a0669523..b5ec7972b131000561831dcd05f9ee53e2f0f627 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -537,7 +537,7 @@ int _removeFourTrianglesNodes(GFace *gf,bool replace_by_quads)
for(int i=0;i<4;i++){
newt[i] = new MTriangle(edges[i][0],edges[(i+1)%4][0],edges[(i+2)%4][0]);
surf[i] = surfaceFaceUV(newt[i],gf);
- qual[i] = qmTriangle(newt[i],QMTRI_RHO);
+ qual[i] = qmTriangle::gamma(newt[i]);
}
double q02=(fabs((surf[0]+surf[2]-surfaceRef)/surfaceRef)<1e-8) ?
std::min(qual[0],qual[2]) : -1;
@@ -2851,10 +2851,10 @@ bool edgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalEdge,
switch(cr){
case SWCR_QUAL:
{
- const double triQualityRef = std::min(qmTriangle(t1->tri(), QMTRI_RHO),
- qmTriangle(t2->tri(), QMTRI_RHO));
- const double triQuality = std::min(qmTriangle(t1b, QMTRI_RHO),
- qmTriangle(t2b, QMTRI_RHO));
+ const double triQualityRef = std::min(qmTriangle::gamma(t1->tri()),
+ qmTriangle::gamma(t2->tri()));
+ const double triQuality = std::min(qmTriangle::gamma(t1b),
+ qmTriangle::gamma(t2b));
if (!edgeSwapDelProj(v1,v2,v3,v4)){
if(triQuality < triQualityRef){
delete t1b;
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 7bbe110baf3505918daa44a2e726411ce3906f0e..9113a1a535fb1746bb986e430159a461da5b9568 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -1778,7 +1778,7 @@ void optimizeMeshGRegionGmsh::operator() (GRegion *gr)
if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == EXTRUDED_ENTITY) return;
Msg::Info("Optimizing volume %d", gr->tag());
- optimizeMesh(gr, QMTET_2);
+ optimizeMesh(gr, qmTetrahedron::QMTET_GAMMA);
}
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 29d6c584dcde3598f33cc7dfe191c7af9baea1cf..fa04be33182491dad64feade6ebf94547dfdd6aa 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -734,7 +734,7 @@ void non_recursive_classify(MTet4 *t, std::list<MTet4*> &theRegion,
void adaptMeshGRegion::operator () (GRegion *gr)
{
- const qualityMeasure4Tet qm = QMTET_2;
+ const qmTetrahedron::Measures qm = qmTetrahedron::QMTET_GAMMA;
typedef std::list<MTet4 *> CONTAINER ;
CONTAINER allTets;
@@ -793,7 +793,7 @@ void adaptMeshGRegion::operator () (GRegion *gr)
if (!(*it)->isDeleted()){
for (int i = 0; i < 4; i++){
for (int j = 0; j < 4; j++){
- if (collapseVertex(newTets, *it, i, j, QMTET_2)){
+ if (collapseVertex(newTets, *it, i, j, qmTetrahedron::QMTET_GAMMA)){
nbCollapse++; i = j = 10;
}
}
@@ -917,7 +917,7 @@ void adaptMeshGRegion::operator () (GRegion *gr)
}
//template <class CONTAINER, class DATA>
-void optimizeMesh(GRegion *gr, const qualityMeasure4Tet &qm)
+void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
{
typedef std::list<MTet4 *> CONTAINER ;
@@ -1443,7 +1443,7 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
double qq = (*it)->getQuality();
if (qq < .4)
for (int i = 0; i < 4; i++){
- if (smoothVertex(*it, i, QMTET_2)) nbReloc++;
+ if (smoothVertex(*it, i, qmTetrahedron::QMTET_GAMMA)) nbReloc++;
}
}
}
@@ -1656,7 +1656,7 @@ void bowyerWatsonFrontalLayers(GRegion *gr, bool hex)
double qq = (*it)->getQuality();
if (qq < .4)
for (int i = 0; i < 4; i++){
- if (smoothVertex(*it, i, QMTET_2)) nbReloc++;
+ if (smoothVertex(*it, i, qmTetrahedron::QMTET_GAMMA)) nbReloc++;
}
}
}
diff --git a/Mesh/meshGRegionDelaunayInsertion.h b/Mesh/meshGRegionDelaunayInsertion.h
index 996c8c75e28b693589234388ae61d8b41140293c..888bd9336fea6cd3a5b81b923af11295207d1a2f 100644
--- a/Mesh/meshGRegionDelaunayInsertion.h
+++ b/Mesh/meshGRegionDelaunayInsertion.h
@@ -72,12 +72,12 @@ class MTet4
{
neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
}
- MTet4(MTetrahedron *t, const qualityMeasure4Tet &qm)
+ MTet4(MTetrahedron *t, const qmTetrahedron::Measures &qm)
: deleted(false), base(t), gr(0)
{
neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
double vol;
- circum_radius = qmTet(t, qm, &vol);
+ circum_radius = qmTetrahedron::qm(t, qm, &vol);
}
void circumcenter(double *res)
{
@@ -266,6 +266,6 @@ class MTet4Factory
container &getAllTets(){ return allTets; }
};
-void optimizeMesh(GRegion *gr, const qualityMeasure4Tet &qm);
+void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm);
#endif
diff --git a/Mesh/meshGRegionLocalMeshMod.cpp b/Mesh/meshGRegionLocalMeshMod.cpp
index a89b892eba464bd300f741b5691b1d60e79eabdc..131c9ccf70f5f94e1c84c5129f0cf05e5b93bdf2 100644
--- a/Mesh/meshGRegionLocalMeshMod.cpp
+++ b/Mesh/meshGRegionLocalMeshMod.cpp
@@ -210,7 +210,7 @@ void BuildSwapPattern7(SwapPattern *sc)
bool edgeSwap(std::vector<MTet4 *> &newTets,
MTet4 *tet,
int iLocalEdge,
- const qualityMeasure4Tet &cr)
+ const qmTetrahedron::Measures &cr)
{
std::vector<MTet4*> cavity;
std::vector<MTet4*> outside;
@@ -250,8 +250,8 @@ bool edgeSwap(std::vector<MTet4 *> &newTets,
int p1 = sp.triangles[i][0];
int p2 = sp.triangles[i][1];
int p3 = sp.triangles[i][2];
- tetQuality1[i] = qmTet(ring[p1], ring[p2], ring[p3], v1, cr, &(volume1[i]));
- tetQuality2[i] = qmTet(ring[p1], ring[p2], ring[p3], v2, cr, &(volume2[i]));
+ tetQuality1[i] = qmTetrahedron::qm(ring[p1], ring[p2], ring[p3], v1, cr, &(volume1[i]));
+ tetQuality2[i] = qmTetrahedron::qm(ring[p1], ring[p2], ring[p3], v2, cr, &(volume2[i]));
}
// look for the best triangulation, i.e. the one that maximize the
@@ -318,7 +318,7 @@ bool edgeSwap(std::vector<MTet4 *> &newTets,
}
bool edgeSplit(std::vector<MTet4 *> &newTets, MTet4 *tet, MVertex *newVertex,
- int iLocalEdge, const qualityMeasure4Tet &cr)
+ int iLocalEdge, const qmTetrahedron::Measures &cr)
{
std::vector<MTet4*> cavity;
std::vector<MTet4*> outside;
@@ -352,7 +352,7 @@ bool edgeSplit(std::vector<MTet4 *> &newTets, MTet4 *tet, MVertex *newVertex,
// swap a face i.e. remove a face shared by 2 tets
bool faceSwap(std::vector<MTet4 *> &newTets, MTet4 *t1, int iLocalFace,
- const qualityMeasure4Tet &cr)
+ const qmTetrahedron::Measures &cr)
{
MTet4 *t2 = t1->getNeigh(iLocalFace);
if (!t2) return false;
@@ -383,11 +383,11 @@ bool faceSwap(std::vector<MTet4 *> &newTets, MTet4 *t1, int iLocalFace,
double q2 = t2->getQuality();
double vol3;
- double q3 = qmTet(f1, f2, v1, v2, cr, &vol3);
+ double q3 = qmTetrahedron::qm(f1, f2, v1, v2, cr, &vol3);
double vol4;
- double q4 = qmTet(f2, f3, v1, v2, cr, &vol4);
+ double q4 = qmTetrahedron::qm(f2, f3, v1, v2, cr, &vol4);
double vol5;
- double q5 = qmTet(f3, f1, v1, v2, cr, &vol5);
+ double q5 = qmTetrahedron::qm(f3, f1, v1, v2, cr, &vol5);
if (fabs(vol1 + vol2 - vol3 - vol4 - vol5) > 1.e-10 * (vol1 + vol2)) return false;
@@ -481,7 +481,7 @@ void buildVertexCavity_recur(MTet4 *t, MVertex *v, std::vector<MTet4*> &cavity)
// after that crap, the sliver is trashed
bool sliverRemoval(std::vector<MTet4*> &newTets, MTet4 *t,
- const qualityMeasure4Tet &cr)
+ const qmTetrahedron::Measures &cr)
{
std::vector<MTet4*> cavity;
std::vector<MTet4*> outside;
@@ -507,7 +507,7 @@ bool sliverRemoval(std::vector<MTet4*> &newTets, MTet4 *t,
else if (nbSwappable == 1){
// classical case, the sliver has 5 edges on the boundary
// try to swap first
- if (edgeSwap(newTets,t,iSwappable,QMTET_3))return true;
+ if (edgeSwap(newTets,t,iSwappable,qmTetrahedron::QMTET_ETA))return true;
// if it does not work, split, smooth and collapse
MVertex *v1 = t->tet()->getVertex(edges[iSwappable][0]);
MVertex *v2 = t->tet()->getVertex(edges[iSwappable][1]);
@@ -516,7 +516,7 @@ bool sliverRemoval(std::vector<MTet4*> &newTets, MTet4 *t,
0.5 * (v1->z() + v2->z()), t->onWhat());
t->onWhat()->mesh_vertices.push_back(newv);
- if (!edgeSplit(newTets, t, newv, iSwappable, QMTET_ONE)) return false;
+ if (!edgeSplit(newTets, t, newv, iSwappable, qmTetrahedron::QMTET_ONE)) return false;
for (int i = 0; i < 4; i++){
if (newTets[newTets.size() - 1]->tet()->getVertex(i) == newv){
smoothVertex(newTets[newTets.size() - 1], i, cr);
@@ -556,7 +556,7 @@ bool collapseVertex(std::vector<MTet4 *> &newTets,
MTet4 *t,
int iVertex,
int iTarget,
- const qualityMeasure4Tet &cr,
+ const qmTetrahedron::Measures &cr,
const localMeshModAction action,
double *minQual)
{
@@ -609,7 +609,7 @@ bool collapseVertex(std::vector<MTet4 *> &newTets,
}
for (unsigned int i = 0; i < toUpdate.size(); i++){
double vv;
- newQuals[i] = qmTet(toUpdate[i]->tet(),cr,&vv);
+ newQuals[i] = qmTetrahedron::qm(toUpdate[i]->tet(),cr,&vv);
worstAfter = std::min(worstAfter,newQuals[i]);
volume_update += vv;
}
@@ -648,7 +648,7 @@ bool collapseVertex(std::vector<MTet4 *> &newTets,
return true;
}
-bool smoothVertex(MTet4 *t, int iVertex, const qualityMeasure4Tet &cr)
+bool smoothVertex(MTet4 *t, int iVertex, const qmTetrahedron::Measures &cr)
{
if(t->isDeleted()){
Msg::Error("Impossible to collapse vertex");
@@ -702,7 +702,7 @@ bool smoothVertex(MTet4 *t, int iVertex, const qualityMeasure4Tet &cr)
}
for (unsigned int i = 0; i < cavity.size(); i++){
double volume;
- newQuals[i] = qmTet(cavity[i]->tet(),cr,&volume);
+ newQuals[i] = qmTetrahedron::qm(cavity[i]->tet(),cr,&volume);
volumeAfter += volume;
worstAfter = std::min(worstAfter,newQuals[i]);
}
@@ -742,7 +742,7 @@ double smoothing_objective_function_3D(double X, double Y, double Z,
std::vector < MTet4 * >::iterator ite = ts.end();
double qMin = 1, vol;
while(it != ite){
- qMin = std::min(qmTet((*it)->tet(), QMTET_2, &vol), qMin);
+ qMin = std::min(qmTetrahedron::qm((*it)->tet(), qmTetrahedron::QMTET_GAMMA, &vol), qMin);
++it;
}
v->x() = oldX;
@@ -776,7 +776,7 @@ double smooth_obj_3D(double *XYZ, void *data)
return smoothing_objective_function_3D(XYZ[0], XYZ[1], XYZ[2], svd->v, svd->ts);
}
-bool smoothVertexOptimize(MTet4 *t, int iVertex, const qualityMeasure4Tet &cr)
+bool smoothVertexOptimize(MTet4 *t, int iVertex, const qmTetrahedron::Measures &cr)
{
if(t->tet()->getVertex(iVertex)->onWhat()->dim() < 3) return false;
@@ -817,7 +817,7 @@ bool smoothVertexOptimize(MTet4 *t, int iVertex, const qualityMeasure4Tet &cr)
}
for(unsigned int i = 0; i < vd.ts.size(); i++){
double volume;
- newQuals[i] = qmTet(vd.ts[i]->tet(), cr, &volume);
+ newQuals[i] = qmTetrahedron::qm(vd.ts[i]->tet(), cr, &volume);
volumeAfter += volume;
}
diff --git a/Mesh/meshGRegionLocalMeshMod.h b/Mesh/meshGRegionLocalMeshMod.h
index a838150ceda5144ca6c753b030bd921252404892..8cbde9c9c8a6465d104a411ff61ee486be17ff0c 100644
--- a/Mesh/meshGRegionLocalMeshMod.h
+++ b/Mesh/meshGRegionLocalMeshMod.h
@@ -16,28 +16,28 @@
enum localMeshModAction {GMSH_DOIT, GMSH_EVALONLY};
bool edgeSwap(std::vector<MTet4*> &newTets, MTet4 *tet,
- int iLocalEdge, const qualityMeasure4Tet &cr);
+ int iLocalEdge, const qmTetrahedron::Measures &cr);
bool faceSwap(std::vector<MTet4*> &newTets, MTet4 *tet,
- int iLocalFace, const qualityMeasure4Tet &cr);
+ int iLocalFace, const qmTetrahedron::Measures &cr);
bool smoothVertex(MTet4 *t, int iLocalVertex,
- const qualityMeasure4Tet &cr);
+ const qmTetrahedron::Measures &cr);
bool smoothVertexOptimize(MTet4 *t, int iVertex,
- const qualityMeasure4Tet &cr);
+ const qmTetrahedron::Measures &cr);
bool collapseVertex(std::vector<MTet4*> &newTets, MTet4 *t,
int iVertex, int iTarget,
- const qualityMeasure4Tet &cr,
+ const qmTetrahedron::Measures &cr,
const localMeshModAction = GMSH_DOIT,
double *result = 0);
bool egeSplit(std::vector<MTet4*> &newTets, MTet4 *tet,
MVertex *newVertex, int iLocalEdge,
- const qualityMeasure4Tet &cr);
+ const qmTetrahedron::Measures &cr);
bool sliverRemoval(std::vector<MTet4*> &newTets, MTet4 *t,
- const qualityMeasure4Tet &cr);
+ const qmTetrahedron::Measures &cr);
#endif
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 84af3697cd345e433f1950d7ab3c25a74a49a42c..ea9615bf1255aeeefb0af65b8eb4ff62c92f2d73 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -9,330 +9,14 @@
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MTetrahedron.h"
+#include "MPrism.h"
+#include "MHexahedron.h"
#include "Numeric.h"
#include "polynomialBasis.h"
#include "GmshMessage.h"
#include <limits>
#include <string.h>
-double qmTriangle(const BDS_Point *p1, const BDS_Point *p2, const BDS_Point *p3,
- const qualityMeasure4Triangle &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 qualityMeasure4Triangle &cr)
-{
- BDS_Point *n[4];
- t->getNodes(n);
- return qmTriangle(n[0], n[1], n[2], cr);
-}
-
-double qmTriangle(MTriangle*t, const qualityMeasure4Triangle &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 qualityMeasure4Triangle &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 qualityMeasure4Triangle &cr)
-{
- double quality;
- switch(cr){
- 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};
- norme(a);
- norme(b);
- 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;
- // condition number
- case QMTRI_COND:
- {
- /*
- double a [3] = {xc - xa, yc - ya, zc - za};
- double b [3] = {xb - xa, yb - ya, zb - za};
- double c [3] ; prodve(a, b, c); norme(c);
- double A[3][3] = {{a[0] , b[0] , c[0]} ,
- {a[1] , b[1] , c[1]} ,
- {a[2] , b[2] , c[2]}};
- */
- quality = -1;
- }
- break;
- default:
- Msg::Error("Unknown quality measure");
- return 0.;
- }
-
- return quality;
-}
-
-double qmTet(MTetrahedron *t, const qualityMeasure4Tet &cr, double *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 qualityMeasure4Tet &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);
-}
-
-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 qualityMeasure4Tet &cr, double *volume)
-{
- 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:
- {
- 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;
- case QMTET_COND:
- {
- /// condition number is defined as (see Knupp & Freitag in IJNME)
- double INVW[3][3] = {{1,-1./sqrt(3.),-1./sqrt(6.)},{0,2/sqrt(3.),-1./sqrt(6.)},{0,0,sqrt(1.5)}};
- double A[3][3] = {{x2-x1,y2-y1,z2-z1},{x3-x1,y3-y1,z3-z1},{x4-x1,y4-y1,z4-z1}};
- double S[3][3],INVS[3][3];
- matmat(A,INVW,S);
- *volume = inv3x3(S,INVS) * 0.70710678118654762;//2/sqrt(2);
- double normS = norm2 (S);
- double normINVS = norm2 (INVS);
- return normS * normINVS;
- }
- default:
- Msg::Error("Unknown quality measure");
- return 0.;
- }
-}
-
-/*
-double conditionNumberAndDerivativeOfTet(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 INVW[3][3] = {{1,-1./sqrt(3.),-1./sqrt(6.)},{0,2/sqrt(3.),-1./sqrt(6.)},{0,0,sqrt(1.5)}};
- double A[3][3] = {{x2-x1,y2-y1,z2-z1},{x3-x1,y3-y1,z3-z1},{x4-x1,y4-y1,z4-z1}};
- double S[3][3],INVS[3][3];
- matmat(A,INVW,S);
- double sigma = inv3x3(S,INVS);
- double normS = norm2 (S);
- double normINVS = norm2 (INVS);
- conditionNumber = normS * normINVS;
-
-}
-*/
-
-double qmTriangleAngles (MTriangle *e) {
- double a = 500;
- double worst_quality = std::numeric_limits<double>::max();
- double mat[3][3];
- double mat2[3][3];
- double den = atan(a*(M_PI/9)) + atan(a*(M_PI/9));
-
- // This matrix is used to "rotate" the triangle to get each vertex
- // as the "origin" of the mapping in turn
- double rot[3][3];
- rot[0][0]=-1; rot[0][1]=1; rot[0][2]=0;
- rot[1][0]=-1; rot[1][1]=0; rot[1][2]=0;
- rot[2][0]= 0; rot[2][1]=0; rot[2][2]=1;
- double tmp[3][3];
-
- //double minAngle = 120.0;
- for (int i = 0; i < e->getNumPrimaryVertices(); i++) {
- const double u = i == 1 ? 1 : 0;
- const double v = i == 2 ? 1 : 0;
- const double w = 0;
- e->getJacobian(u, v, w, mat);
- e->getPrimaryJacobian(u,v,w,mat2);
- for (int j = 0; j < i; j++) {
- matmat(rot,mat,tmp);
- memcpy(mat, tmp, sizeof(mat));
- }
- //get angle
- double v1[3] = {mat[0][0], mat[0][1], mat[0][2] };
- double v2[3] = {mat[1][0], mat[1][1], mat[1][2] };
- double v3[3] = {mat2[0][0], mat2[0][1], mat2[0][2] };
- double v4[3] = {mat2[1][0], mat2[1][1], mat2[1][2] };
- norme(v1);
- norme(v2);
- norme(v3);
- norme(v4);
- double v12[3], v34[3];
- prodve(v1,v2,v12);
- prodve(v3,v4,v34);
- norme(v12);
- norme(v34);
- double orientation;
- prosca(v12,v34,&orientation);
-
- // If the triangle is "flipped" it's no good
- if (orientation < 0)
- return -std::numeric_limits<double>::max();
-
- double c;
- prosca(v1,v2,&c);
- double x = acos(c)-M_PI/3;
- //double angle = (x+M_PI/3)/M_PI*180;
- double quality = (atan(a*(x+M_PI/9)) + atan(a*(M_PI/9-x)))/den;
- worst_quality = std::min(worst_quality, quality);
-
- //minAngle = std::min(angle, minAngle);
- //printf("Angle %g ", angle);
- // printf("Quality %g\n",quality);
- }
- //printf("MinAngle %g \n", minAngle);
- //return minAngle;
-
- return worst_quality;
-}
-
-
-double qmQuadrangleAngles (MQuadrangle *e) {
- double a = 100;
- double worst_quality = std::numeric_limits<double>::max();
- double mat[3][3];
- double mat2[3][3];
- double den = atan(a*(M_PI/4)) + atan(a*(2*M_PI/4 - (M_PI/4)));
-
- // This matrix is used to "rotate" the triangle to get each vertex
- // as the "origin" of the mapping in turn
- //double rot[3][3];
- //rot[0][0]=-1; rot[0][1]=1; rot[0][2]=0;
- //rot[1][0]=-1; rot[1][1]=0; rot[1][2]=0;
- //rot[2][0]= 0; rot[2][1]=0; rot[2][2]=1;
- //double tmp[3][3];
-
- const double u[9] = {-1,-1, 1, 1, 0,0,1,-1,0};
- const double v[9] = {-1, 1, 1,-1, -1,1,0,0,0};
-
- for (int i = 0; i < 9; i++) {
-
- e->getJacobian(u[i], v[i], 0, mat);
- e->getPrimaryJacobian(u[i],v[i],0,mat2);
- //for (int j = 0; j < i; j++) {
- // matmat(rot,mat,tmp);
- // memcpy(mat, tmp, sizeof(mat));
- //}
-
- //get angle
- double v1[3] = {mat[0][0], mat[0][1], mat[0][2] };
- double v2[3] = {mat[1][0], mat[1][1], mat[1][2] };
- double v3[3] = {mat2[0][0], mat2[0][1], mat2[0][2] };
- double v4[3] = {mat2[1][0], mat2[1][1], mat2[1][2] };
- norme(v1);
- norme(v2);
- norme(v3);
- norme(v4);
- double v12[3], v34[3];
- prodve(v1,v2,v12);
- prodve(v3,v4,v34);
- norme(v12);
- norme(v34);
- double orientation;
- prosca(v12,v34,&orientation);
-
- // If the if the triangle is "flipped" it's no good
- // if (orientation < 0)
- // return -std::numeric_limits<double>::max();
-
- double c;
- prosca(v1,v2,&c);
- double x = fabs(acos(c))-M_PI/2;
- //double angle = fabs(acos(c))*180/M_PI;
- double quality = (atan(a*(x+M_PI/4)) + atan(a*(2*M_PI/4 - (x+M_PI/4))))/den;
- worst_quality = std::min(worst_quality, quality);
- }
-
- return worst_quality;
-
-}
-
namespace {
@@ -471,6 +155,7 @@ inline void triAreaAndGrad(const double fact, const fullMatrix<double> &vga,
}
+// Revert vector and gradients
inline void revertVG(const fullMatrix<double> &vg, fullMatrix<double> &res)
{
res(0, 0) = -vg(0, 3); res(0, 1) = -vg(0, 4); res(0, 2) = -vg(0, 5); res(0, 6) = -vg(0, 6);
@@ -482,10 +167,138 @@ inline void revertVG(const fullMatrix<double> &vg, fullMatrix<double> &res)
}
-void measuresTriangle::getNCJ(const double &x0, const double &y0, const double &z0,
- const double &x1, const double &y1, const double &z1,
- const double &x2, const double &y2, const double &z2,
- fullVector<double> &ncj)
+double qmTriangle::gamma(const BDS_Point *p1, const BDS_Point *p2, const BDS_Point *p3)
+{
+ return gamma(p1->X, p1->Y, p1->Z, p2->X, p2->Y, p2->Z, p3->X, p3->Y, p3->Z);
+}
+
+
+double qmTriangle::gamma(BDS_Face *t)
+{
+ BDS_Point *n[4];
+ t->getNodes(n);
+ return gamma(n[0], n[1], n[2]);
+}
+
+
+double qmTriangle::gamma(MTriangle*t)
+{
+ return gamma(t->getVertex(0), t->getVertex(1), t->getVertex(2));
+}
+
+
+double qmTriangle::gamma(const MVertex *v1, const MVertex *v2, const MVertex *v3)
+{
+ return gamma(v1->x(), v1->y(), v1->z(), v2->x(), v2->y(), v2->z(), v3->x(), v3->y(), v3->z());
+}
+
+
+// Triangle abc
+// quality is between 0 and 1
+double qmTriangle::gamma(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)
+{
+ // 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};
+ norme(a);
+ norme(b);
+ 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) return 0.0;
+ else return 2 * (2 * sina * sinb * sinc / (sina + sinb + sinc));
+}
+
+
+double qmTriangle::eta(MTriangle *el)
+{
+ MVertex *_v[3] = {el->getVertex(0), el->getVertex(1), el->getVertex(2)};
+
+ double a1 = 180 * angle3Vertices(_v[0], _v[1], _v[2]) / M_PI;
+ double a2 = 180 * angle3Vertices(_v[1], _v[2], _v[0]) / M_PI;
+ double a3 = 180 * angle3Vertices(_v[2], _v[0], _v[1]) / M_PI;
+
+ double amin = std::min(std::min(a1,a2),a3);
+ double angle = fabs(60. - amin);
+ return 1.-angle/60;
+}
+
+
+double qmTriangle::angles(MTriangle *e)
+{
+ double a = 500;
+ double worst_quality = std::numeric_limits<double>::max();
+ double mat[3][3];
+ double mat2[3][3];
+ double den = atan(a*(M_PI/9)) + atan(a*(M_PI/9));
+
+ // This matrix is used to "rotate" the triangle to get each vertex
+ // as the "origin" of the mapping in turn
+ double rot[3][3];
+ rot[0][0]=-1; rot[0][1]=1; rot[0][2]=0;
+ rot[1][0]=-1; rot[1][1]=0; rot[1][2]=0;
+ rot[2][0]= 0; rot[2][1]=0; rot[2][2]=1;
+ double tmp[3][3];
+
+ //double minAngle = 120.0;
+ for (int i = 0; i < e->getNumPrimaryVertices(); i++) {
+ const double u = i == 1 ? 1 : 0;
+ const double v = i == 2 ? 1 : 0;
+ const double w = 0;
+ e->getJacobian(u, v, w, mat);
+ e->getPrimaryJacobian(u,v,w,mat2);
+ for (int j = 0; j < i; j++) {
+ matmat(rot,mat,tmp);
+ memcpy(mat, tmp, sizeof(mat));
+ }
+ //get angle
+ double v1[3] = {mat[0][0], mat[0][1], mat[0][2] };
+ double v2[3] = {mat[1][0], mat[1][1], mat[1][2] };
+ double v3[3] = {mat2[0][0], mat2[0][1], mat2[0][2] };
+ double v4[3] = {mat2[1][0], mat2[1][1], mat2[1][2] };
+ norme(v1);
+ norme(v2);
+ norme(v3);
+ norme(v4);
+ double v12[3], v34[3];
+ prodve(v1,v2,v12);
+ prodve(v3,v4,v34);
+ norme(v12);
+ norme(v34);
+ double orientation;
+ prosca(v12,v34,&orientation);
+
+ // If the triangle is "flipped" it's no good
+ if (orientation < 0)
+ return -std::numeric_limits<double>::max();
+
+ double c;
+ prosca(v1,v2,&c);
+ double x = acos(c)-M_PI/3;
+ //double angle = (x+M_PI/3)/M_PI*180;
+ double quality = (atan(a*(x+M_PI/9)) + atan(a*(M_PI/9-x)))/den;
+ worst_quality = std::min(worst_quality, quality);
+
+ //minAngle = std::min(angle, minAngle);
+ //printf("Angle %g ", angle);
+ // printf("Quality %g\n",quality);
+ }
+ //printf("MinAngle %g \n", minAngle);
+ //return minAngle;
+
+ return worst_quality;
+}
+
+
+void qmTriangle::NCJ(const double &x0, const double &y0, const double &z0,
+ const double &x1, const double &y1, const double &z1,
+ const double &x2, const double &y2, const double &z2,
+ fullVector<double> &ncj)
{
// Compute unit vectors for each edge
double x01n, y01n, z01n, x12n, y12n, z12n, x20n, y20n, z20n;
@@ -495,20 +308,20 @@ void measuresTriangle::getNCJ(const double &x0, const double &y0, const double &
// Compute NCJ at vertex from unit vectors a and b as 0.5*||a^b||/A_equilateral
// Factor = 2./sqrt(3.) = 0.5/A_equilateral
- static const double fact = 1.1547005383792517;
+ static const double fact = 2./sqrt(3.);
ncj(0) = triArea(fact, x01n, y01n, z01n, -x20n, -y20n, -z20n);
ncj(1) = triArea(fact, x12n, y12n, z12n, -x01n, -y01n, -z01n);
ncj(2) = triArea(fact, x20n, y20n, z20n, -x12n, -y12n, -z12n);
}
-void measuresTriangle::getNCJAndGradients(const double &x0, const double &y0, const double &z0,
- const double &x1, const double &y1, const double &z1,
- const double &x2, const double &y2, const double &z2,
- fullMatrix<double> &res)
+void qmTriangle::NCJAndGradients(const double &x0, const double &y0, const double &z0,
+ const double &x1, const double &y1, const double &z1,
+ const double &x2, const double &y2, const double &z2,
+ fullMatrix<double> &res)
{
// Factor = 2./sqrt(3.) = 0.5/A_equilateral
- static const double fact = 1.1547005383792517;
+ static const double fact = 2./sqrt(3.);
// Compute unit vector and its gradients for each edge
fullMatrix<double> vg01(3,7);
@@ -553,11 +366,111 @@ void measuresTriangle::getNCJAndGradients(const double &x0, const double &y0, co
}
-void measuresQuadrangle::getNCJ(const double &x0, const double &y0, const double &z0,
- 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,
- fullVector<double> &ncj)
+double qmQuadrangle::eta(MQuadrangle *el) {
+ double AR = 1;//(minEdge()/maxEdge());
+
+ MVertex *_v[4] = {el->getVertex(0), el->getVertex(1), el->getVertex(2), el->getVertex(3)};
+
+ SVector3 v01 (_v[1]->x()-_v[0]->x(),_v[1]->y()-_v[0]->y(),_v[1]->z()-_v[0]->z());
+ SVector3 v12 (_v[2]->x()-_v[1]->x(),_v[2]->y()-_v[1]->y(),_v[2]->z()-_v[1]->z());
+ SVector3 v23 (_v[3]->x()-_v[2]->x(),_v[3]->y()-_v[2]->y(),_v[3]->z()-_v[2]->z());
+ SVector3 v30 (_v[0]->x()-_v[3]->x(),_v[0]->y()-_v[3]->y(),_v[0]->z()-_v[3]->z());
+
+ SVector3 a = crossprod(v01,v12);
+ SVector3 b = crossprod(v12,v23);
+ SVector3 c = crossprod(v23,v30);
+ SVector3 d = crossprod(v30,v01);
+
+ double sign = 1.0;
+ if (dot(a,b) < 0 || dot(a,c) < 0 || dot(a,d) < 0 )sign = -1;
+ // FIXME ...
+ // if (a.z() > 0 || b.z() > 0 || c.z() > 0 || d.z() > 0) sign = -1;
+
+ double a1 = 180 * angle3Vertices(_v[0], _v[1], _v[2]) / M_PI;
+ double a2 = 180 * angle3Vertices(_v[1], _v[2], _v[3]) / M_PI;
+ double a3 = 180 * angle3Vertices(_v[2], _v[3], _v[0]) / M_PI;
+ double a4 = 180 * angle3Vertices(_v[3], _v[0], _v[1]) / M_PI;
+
+ a1 = std::min(180.,a1);
+ a2 = std::min(180.,a2);
+ a3 = std::min(180.,a3);
+ a4 = std::min(180.,a4);
+ double angle = fabs(90. - a1);
+ angle = std::max(fabs(90. - a2),angle);
+ angle = std::max(fabs(90. - a3),angle);
+ angle = std::max(fabs(90. - a4),angle);
+
+ return sign*(1.-angle/90) * AR;
+}
+
+
+double qmQuadrangle::angles(MQuadrangle *e)
+{
+ double a = 100;
+ double worst_quality = std::numeric_limits<double>::max();
+ double mat[3][3];
+ double mat2[3][3];
+ double den = atan(a*(M_PI/4)) + atan(a*(2*M_PI/4 - (M_PI/4)));
+
+ // This matrix is used to "rotate" the triangle to get each vertex
+ // as the "origin" of the mapping in turn
+ //double rot[3][3];
+ //rot[0][0]=-1; rot[0][1]=1; rot[0][2]=0;
+ //rot[1][0]=-1; rot[1][1]=0; rot[1][2]=0;
+ //rot[2][0]= 0; rot[2][1]=0; rot[2][2]=1;
+ //double tmp[3][3];
+
+ const double u[9] = {-1,-1, 1, 1, 0,0,1,-1,0};
+ const double v[9] = {-1, 1, 1,-1, -1,1,0,0,0};
+
+ for (int i = 0; i < 9; i++) {
+
+ e->getJacobian(u[i], v[i], 0, mat);
+ e->getPrimaryJacobian(u[i],v[i],0,mat2);
+ //for (int j = 0; j < i; j++) {
+ // matmat(rot,mat,tmp);
+ // memcpy(mat, tmp, sizeof(mat));
+ //}
+
+ //get angle
+ double v1[3] = {mat[0][0], mat[0][1], mat[0][2] };
+ double v2[3] = {mat[1][0], mat[1][1], mat[1][2] };
+ double v3[3] = {mat2[0][0], mat2[0][1], mat2[0][2] };
+ double v4[3] = {mat2[1][0], mat2[1][1], mat2[1][2] };
+ norme(v1);
+ norme(v2);
+ norme(v3);
+ norme(v4);
+ double v12[3], v34[3];
+ prodve(v1,v2,v12);
+ prodve(v3,v4,v34);
+ norme(v12);
+ norme(v34);
+ double orientation;
+ prosca(v12,v34,&orientation);
+
+ // If the if the triangle is "flipped" it's no good
+ // if (orientation < 0)
+ // return -std::numeric_limits<double>::max();
+
+ double c;
+ prosca(v1,v2,&c);
+ double x = fabs(acos(c))-M_PI/2;
+ //double angle = fabs(acos(c))*180/M_PI;
+ double quality = (atan(a*(x+M_PI/4)) + atan(a*(2*M_PI/4 - (x+M_PI/4))))/den;
+ worst_quality = std::min(worst_quality, quality);
+ }
+
+ return worst_quality;
+
+}
+
+
+void qmQuadrangle::NCJ(const double &x0, const double &y0, const double &z0,
+ 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,
+ fullVector<double> &ncj)
{
// Compute unit vectors for each edge
double x01n, y01n, z01n, x12n, y12n, z12n, x23n, y23n, z23n, x30n, y30n, z30n;
@@ -576,11 +489,11 @@ void measuresQuadrangle::getNCJ(const double &x0, const double &y0, const double
}
-void measuresQuadrangle::getNCJAndGradients(const double &x0, const double &y0, const double &z0,
- 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,
- fullMatrix<double> &res)
+void qmQuadrangle::NCJAndGradients(const double &x0, const double &y0, const double &z0,
+ 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,
+ fullMatrix<double> &res)
{
// Factor = 1. = 0.5/A_triRect
static const double fact = 1.;
@@ -643,3 +556,215 @@ void measuresQuadrangle::getNCJAndGradients(const double &x0, const double &y0,
res(3, 9) = res3(0); res(3, 10) = res3(1); res(3, 11) = res3(2);
res(3, 12) = res3(9);
}
+
+
+double qmTetrahedron::qm(MTetrahedron *t, const Measures &cr, double *volume)
+{
+ return qm(t->getVertex(0), t->getVertex(1), t->getVertex(2), t->getVertex(3), cr, volume);
+}
+
+
+double qmTetrahedron::qm(const MVertex *v1, const MVertex *v2, const MVertex *v3,
+ const MVertex *v4, const Measures &cr, double *volume)
+{
+ return qm(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 qmTetrahedron::qm(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 Measures &cr, double *volume)
+{
+ switch(cr){
+ case QMTET_ONE:
+ return 1.0;
+ case QMTET_ETA:
+ return eta(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, volume);
+ case QMTET_GAMMA:
+ return gamma(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, volume);
+ case QMTET_COND:
+ return cond(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, volume);
+ default:
+ Msg::Error("Unknown quality measure");
+ return 0.;
+ }
+}
+
+
+double qmTetrahedron::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)
+{
+ 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;
+}
+
+
+double qmTetrahedron::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)
+{
+ 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;
+}
+
+
+double qmTetrahedron::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)
+{
+ /// condition number is defined as (see Knupp & Freitag in IJNME)
+ double INVW[3][3] = {{1,-1./sqrt(3.),-1./sqrt(6.)},{0,2/sqrt(3.),-1./sqrt(6.)},{0,0,sqrt(1.5)}};
+ double A[3][3] = {{x2-x1,y2-y1,z2-z1},{x3-x1,y3-y1,z3-z1},{x4-x1,y4-y1,z4-z1}};
+ double S[3][3],INVS[3][3];
+ matmat(A,INVW,S);
+ *volume = inv3x3(S,INVS) * 0.70710678118654762;//2/sqrt(2);
+ double normS = norm2 (S);
+ double normINVS = norm2 (INVS);
+ return normS * normINVS;
+}
+
+
+// TODO: Replace this
+static double prismNCJ(const MVertex* a, const MVertex* b, const MVertex* c, const MVertex* d)
+{
+ static const double fact = 2./sqrt(3);
+
+ const SVector3 vec1 = SVector3(b->x()-a->x(),b->y()-a->y(),b->z()-a->z());
+ const SVector3 vec2 = SVector3(c->x()-a->x(),c->y()-a->y(),c->z()-a->z());
+ const SVector3 vec3 = SVector3(d->x()-a->x(),d->y()-a->y(),d->z()-a->z());
+
+ const double l1 = vec1.norm();
+ const double l2 = vec2.norm();
+ const double l3 = vec3.norm();
+
+ const double val = dot(vec1,crossprod(vec2,vec3));
+ return fact*fabs(val)/(l1*l2*l3);
+}
+
+
+double qmPrism::minNCJ(const MPrism *el) {
+ const MVertex *a = el->getVertex(0),*b= el->getVertex(1), *c= el->getVertex(2);
+ const MVertex *d = el->getVertex(3),*e= el->getVertex(4), *f= el->getVertex(5);
+ const double j[6] = {
+ prismNCJ(a,b,c,d),
+ prismNCJ(b,a,c,e),
+ prismNCJ(c,a,b,f),
+ prismNCJ(d,a,e,f),
+ prismNCJ(e,b,d,f),
+ prismNCJ(f,c,d,e)};
+ const double result = *std::min_element(j,j+6);
+ return result;
+}
+
+
+//void qmPrism::NCJ(const double &x0, const double &y0, const double &z0,
+// 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,
+// fullVector<double> &ncj)
+//{
+// // Compute unit vectors for each edge
+// double x01n, y01n, z01n, x12n, y12n, z12n, x23n, y23n, z23n, x30n, y30n, z30n;
+// unitVec(x0, y0, z0, x1, y1, z1, x01n, y01n, z01n);
+// unitVec(x1, y1, z1, x2, y2, z2, x12n, y12n, z12n);
+// unitVec(x2, y2, z2, x3, y3, z3, x23n, y23n, z23n);
+// unitVec(x3, y3, z3, x0, y0, z0, x30n, y30n, z30n);
+//
+// // Compute NCJ at vertex from unit vectors a and b as 0.5*||a^b||/A_equilateral
+// // Factor = 2./sqrt(3.) = 0.5/A_equilateral
+// static const double fact = 1.1547005383792517;
+// ncj(0) = triArea(fact, x01n, y01n, z01n, -x20n, -y20n, -z20n);
+// ncj(1) = triArea(fact, x12n, y12n, z12n, -x01n, -y01n, -z01n);
+// ncj(2) = triArea(fact, x20n, y20n, z20n, -x12n, -y12n, -z12n);
+//}
+
+
+// TODO: Remove this (useless as quality measure)
+double qmHexahedron::angles(MHexahedron *el)
+{
+ double angleMax = 0.0;
+ double angleMin = M_PI;
+ double zeta = 0.0;
+ for (int i=0; i<el->getNumFaces(); i++){
+ std::vector<MVertex*> vv;
+ vv.push_back(el->getFace(i).getVertex(0));
+ vv.push_back(el->getFace(i).getVertex(1));
+ vv.push_back(el->getFace(i).getVertex(2));
+ vv.push_back(el->getFace(i).getVertex(3));
+ // MVertex *v0 = new MVertex(0, 0, 0); vv.push_back(v0);
+ // MVertex *v1 = new MVertex(1., 0, 0);vv.push_back(v1);
+ // MVertex *v2 = new MVertex(2., 1., 0);vv.push_back(v2);
+ // MVertex *v3 = new MVertex(1, 1., 0);vv.push_back(v3);
+ for (int j=0; j<4; j++){
+ SVector3 a(vv[(j+2)%4]->x()-vv[(j+1)%4]->x(),vv[(j+2)%4]->y()-vv[(j+1)%4]->y(),vv[(j+2)%4]->z()-vv[(j+1)%4]->z() );
+ SVector3 b(vv[(j+1)%4]->x()-vv[(j)%4]->x(), vv[(j+1)%4]->y()-vv[(j)%4]->y(), vv[(j+1)%4]->z()-vv[(j)%4]->z() );
+ double angle = acos( dot(a,b)/(norm(a)*norm(b))); //*180/M_PI;
+ angleMax = std::max(angleMax, angle);
+ angleMin = std::min(angleMin, angle);
+ }
+ //printf("angle max =%g min =%g \n", angleMax*180/M_PI, angleMin*180/M_PI);
+ }
+ zeta = 1.-std::max((angleMax-0.5*M_PI)/(0.5*M_PI),(0.5*M_PI-angleMin)/(0.5*M_PI));
+ return zeta;
+}
diff --git a/Mesh/qualityMeasures.h b/Mesh/qualityMeasures.h
index 9637e25057be7394216ffcd0831dbad0a4c78689..8086d20856051492719b78ec61de90f5402f64dc 100644
--- a/Mesh/qualityMeasures.h
+++ b/Mesh/qualityMeasures.h
@@ -14,61 +14,121 @@ class MVertex;
class MTriangle;
class MQuadrangle;
class MTetrahedron;
+class MPrism;
+class MHexahedron;
class MElement;
-enum qualityMeasure4Triangle {QMTRI_RHO, QMTRI_COND};
-enum qualityMeasure4Tet{QMTET_1, QMTET_2, QMTET_3, QMTET_ONE, QMTET_COND};
-double qmTriangleAngles (MTriangle *e);
-double qmQuadrangleAngles (MQuadrangle *e);
+class qmTriangle
+{
+public:
+ static double gamma(MTriangle *f);
+ static double gamma(BDS_Face *f);
+ static double gamma(const BDS_Point *p1, const BDS_Point *p2, const BDS_Point *p3);
+ static double gamma(const MVertex *v1, const MVertex *v2, const MVertex *v3);
+ static double gamma(const double *d1, const double *d2, const double *d3);
+ 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);
+ static double eta(MTriangle *el);
+ static double angles(MTriangle *e);
+ static double minNCJ(const MTriangle *e);
+ static void NCJ(const double &x0, const double &y0, const double &z0,
+ const double &x1, const double &y1, const double &z1,
+ const double &x2, const double &y2, const double &z2,
+ fullVector<double> &ncj);
+ static void NCJAndGradients(const double &x0, const double &y0, const double &z0,
+ const double &x1, const double &y1, const double &z1,
+ const double &x2, const double &y2, const double &z2,
+ fullMatrix<double> &ncj);
+};
-double qmTriangle(MTriangle *f, const qualityMeasure4Triangle &cr);
-double qmTriangle(BDS_Face *f, const qualityMeasure4Triangle &cr);
-double qmTriangle(const BDS_Point *p1, const BDS_Point *p2, const BDS_Point *p3,
- const qualityMeasure4Triangle &cr);
-double qmTriangle(const MVertex *v1, const MVertex *v2, const MVertex *v3,
- const qualityMeasure4Triangle &cr);
-double qmTriangle(const double *d1, const double *d2, const double *d3,
- const qualityMeasure4Triangle &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 qualityMeasure4Triangle &cr);
-double qmTet(MTetrahedron *t, const qualityMeasure4Tet &cr, double *volume = 0);
-double qmTet(const MVertex *v1, const MVertex *v2, const MVertex *v3,
- const MVertex *v4, const qualityMeasure4Tet &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 qualityMeasure4Tet &cr, double *volume = 0);
-class measuresTriangle
+class qmQuadrangle
{
public:
- static void getNCJ(const double &x0, const double &y0, const double &z0,
- const double &x1, const double &y1, const double &z1,
- const double &x2, const double &y2, const double &z2,
- fullVector<double> &ncj);
- static void getNCJAndGradients(const double &x0, const double &y0, const double &z0,
- const double &x1, const double &y1, const double &z1,
- const double &x2, const double &y2, const double &z2,
- fullMatrix<double> &ncj);
+ static double gamma(MQuadrangle *el) { return eta(el); }
+ static double eta(MQuadrangle *el);
+ static double angles(MQuadrangle *e);
+ static double minNCJ(const MQuadrangle *e);
+ static void NCJ(const double &x0, const double &y0, const double &z0,
+ 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,
+ fullVector<double> &ncj);
+ static void NCJAndGradients(const double &x0, const double &y0, const double &z0,
+ 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,
+ fullMatrix<double> &ncj);
};
-class measuresQuadrangle
+
+class qmTetrahedron
{
public:
- static void getNCJ(const double &x0, const double &y0, const double &z0,
- const double &x1, const double &y1, const double &z1,
+ enum Measures {QMTET_GAMMA, QMTET_ETA, QMTET_ONE, QMTET_COND};
+ static double qm(MTetrahedron *t, const Measures &cr, double *volume = 0);
+ static double qm(const BDS_Point *p1, const BDS_Point *p2, const BDS_Point *p3);
+ static double qm(const MVertex *v1, const MVertex *v2, const MVertex *v3,
+ const MVertex *v4, const Measures &cr, double *volume = 0);
+ static double qm(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 Measures &cr, double *volume = 0);
+ 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 = 0);
+ 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 = 0);
+ 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,
- fullVector<double> &ncj);
- static void getNCJAndGradients(const double &x0, const double &y0, const double &z0,
- 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,
- fullMatrix<double> &ncj);
+ const double &x4, const double &y4, const double &z4,
+ double *volume = 0);
+ static double minNCJ(const MTetrahedron *e);
+ static void NCJ(const double &x0, const double &y0, const double &z0,
+ 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,
+ fullVector<double> &ncj);
+ static void NCJAndGradients(const double &x0, const double &y0, const double &z0,
+ 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,
+ fullMatrix<double> &ncj);
+};
+
+
+class qmPrism
+{
+public:
+ static double minNCJ(const MPrism *el);
+// static void NCJ(const double &x0, const double &y0, const double &z0,
+// 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,
+// fullVector<double> &ncj);
+};
+
+
+class qmHexahedron
+{
+public:
+ static double angles(MHexahedron *el);
+// static void NCJ(const double &x0, const double &y0, const double &z0,
+// 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,
+// fullVector<double> &ncj);
};
#endif