From bec8701af421652115cb59d5bb71c7eb5b159a77 Mon Sep 17 00:00:00 2001
From: Amaury Johnen <amaury.johnen@uclouvain.be>
Date: Fri, 16 Jun 2017 11:23:53 +0200
Subject: [PATCH] fix bugs
---
Mesh/Generator.cpp | 2 +-
Mesh/qualityMeasuresJacobian.cpp | 71 ++++++++++++++++++++++++++++----
Mesh/qualityMeasuresJacobian.h | 4 +-
3 files changed, 65 insertions(+), 12 deletions(-)
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index b34fba5a0e..a956ebf86d 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -144,7 +144,7 @@ static void GetQualityMeasure(std::vector<T*> &ele,
if(s > (2*j-100) / 100. && s <= (2*j-98) / 100.) quality[0][j]++;
if(g > j / 100. && g <= (j + 1) / 100.) quality[1][j]++;
if(r > j / 100. && r <= (j + 1) / 100.) quality[2][j]++;
- if(e > (2*j-100) / 100. && e <= (2*j-98) / 100.) quality[4][j]++;
+ if(e > (2*j-100) / 100. && e <= (2*j-98) / 100.) quality[3][j]++;
}
}
}
diff --git a/Mesh/qualityMeasuresJacobian.cpp b/Mesh/qualityMeasuresJacobian.cpp
index fa960a6ed2..6d82b43ba1 100644
--- a/Mesh/qualityMeasuresJacobian.cpp
+++ b/Mesh/qualityMeasuresJacobian.cpp
@@ -20,9 +20,9 @@ static const double cTri = 2/std::sqrt(3);
static const double cTet = std::sqrt(2);
static const double cPyr = std::sqrt(2);
-static void computeCoeffLengthVectors_(const fullMatrix<double> &mat,
- fullMatrix<double> &coeff,
- int type, int numCoeff = -1)
+static inline void computeCoeffLengthVectors_(const fullMatrix<double> &mat,
+ fullMatrix<double> &coeff,
+ int type, int numCoeff = -1)
{
int sz1 = numCoeff > -1 ? numCoeff : mat.size1();
@@ -106,10 +106,10 @@ static void computeCoeffLengthVectors_(const fullMatrix<double> &mat,
}
}
-static void computeIGE_(const fullVector<double> &det,
- const fullMatrix<double> &v,
- fullVector<double> &ige,
- int type)
+static inline void computeIGE_(const fullVector<double> &det,
+ const fullMatrix<double> &v,
+ fullVector<double> &ige,
+ int type)
{
int sz = std::min(det.size(), v.size1());
ige.resize(sz);
@@ -124,18 +124,21 @@ static void computeIGE_(const fullVector<double> &det,
for (int i = 0; i < sz; i++) {
ige(i) = det(i)/v(i, 0)/v(i, 1)/v(i, 2);
}
+ break;
case TYPE_TRI:
for (int i = 0; i < sz; i++) {
ige(i) = cTri * det(i) * (1/v(i,0)/v(i,1) +
1/v(i,0)/v(i,2) +
1/v(i,1)/v(i,2)) / 3;
}
+ break;
case TYPE_PRI:
for (int i = 0; i < sz; i++) {
ige(i) = cTri * det(i) * (1/v(i,0)/v(i,1)/v(i,2) +
1/v(i,0)/v(i,3)/v(i,2) +
1/v(i,1)/v(i,3)/v(i,2)) / 3;
}
+ break;
case TYPE_TET:
for (int i = 0; i < sz; i++) {
ige(i) = cTet * det(i) * (1/v(i,0)/v(i,5)/v(i,1) +
@@ -149,8 +152,9 @@ static void computeIGE_(const fullVector<double> &det,
1/v(i,2)/v(i,3)/v(i,0) +
1/v(i,2)/v(i,3)/v(i,1) +
1/v(i,2)/v(i,3)/v(i,4) +
- 1/v(i,2)/v(i,3)/v(i,5)) / 12;
+ 1/v(i,2)/v(i,3)/v(i,5))/ 12;
}
+ break;
case TYPE_PYR:
for (int i = 0; i < sz; i++) {
ige(i) = cPyr * det(i) * (1/v(i,0)/v(i,1)/v(i,2) +
@@ -162,7 +166,56 @@ static void computeIGE_(const fullVector<double> &det,
1/v(i,4)/v(i,5)/v(i,2) +
1/v(i,4)/v(i,5)/v(i,3) ) / 8;
}
+ break;
}
+// for (int i = 0; i < sz; i++) {
+// double sJ;
+// switch (type) {
+// case TYPE_QUA:
+// ige(i) = det(i) / v(i, 0) / v(i, 1);
+// break;
+// case TYPE_HEX:
+// ige(i) = det(i) / v(i, 0) / v(i, 1) / v(i, 2);
+// break;
+// case TYPE_TRI:
+// ige(i) = cTri * det(i) * (1 / v(i, 0) / v(i, 1) +
+// 1 / v(i, 0) / v(i, 2) +
+// 1 / v(i, 1) / v(i, 2)) / 3;
+// break;
+// case TYPE_TET:
+// ige(i) = cTet * det(i) * (1 / v(i, 0) / v(i, 5) / v(i, 1) +
+// 1 / v(i, 0) / v(i, 5) / v(i, 2) +
+// 1 / v(i, 0) / v(i, 5) / v(i, 3) +
+// 1 / v(i, 0) / v(i, 5) / v(i, 4) +
+// 1 / v(i, 1) / v(i, 4) / v(i, 0) +
+// 1 / v(i, 1) / v(i, 4) / v(i, 2) +
+// 1 / v(i, 1) / v(i, 4) / v(i, 3) +
+// 1 / v(i, 1) / v(i, 4) / v(i, 5) +
+// 1 / v(i, 2) / v(i, 3) / v(i, 0) +
+// 1 / v(i, 2) / v(i, 3) / v(i, 1) +
+// 1 / v(i, 2) / v(i, 3) / v(i, 4) +
+// 1 / v(i, 2) / v(i, 3) / v(i, 5)) / 12;
+// break;
+// case TYPE_PRI:
+// ige(i) = cTri * det(i) * (1 / v(i, 0) / v(i, 1) / v(i, 2) +
+// 1 / v(i, 0) / v(i, 3) / v(i, 2) +
+// 1 / v(i, 1) / v(i, 3) / v(i, 2)) / 3;
+// break;
+// case TYPE_PYR:
+// ige(i) = cPyr * det(i) * (1 / v(i, 0) / v(i, 1) / v(i, 2) +
+// 1 / v(i, 0) / v(i, 1) / v(i, 3) +
+// 1 / v(i, 0) / v(i, 1) / v(i, 4) +
+// 1 / v(i, 0) / v(i, 1) / v(i, 5) +
+// 1 / v(i, 2) / v(i, 3) / v(i, 4) +
+// 1 / v(i, 2) / v(i, 3) / v(i, 5) +
+// 1 / v(i, 4) / v(i, 5) / v(i, 2) +
+// 1 / v(i, 4) / v(i, 5) / v(i, 3)) / 8;
+// break;
+// default:
+// Msg::Error("Unkown type for IGE computation");
+// return;
+// }
+// }
}
namespace jacobianBasedQuality {
@@ -727,7 +780,7 @@ double _CoeffDataIGE::_computeLowerBound() const
}
fullMatrix<double> v;
- _getCoeffLengthVectors(v, false);
+ computeCoeffLengthVectors_(_coeffsJacMat, v, _type);
fullVector<double> prox[6];
for (int i = 0; i < v.size2(); ++i) {
diff --git a/Mesh/qualityMeasuresJacobian.h b/Mesh/qualityMeasuresJacobian.h
index 5466765b0b..6c3b62992a 100644
--- a/Mesh/qualityMeasuresJacobian.h
+++ b/Mesh/qualityMeasuresJacobian.h
@@ -32,7 +32,7 @@ class _CoeffData
protected:
double _minL, _maxL; //Extremum of Jac at corners
double _minB, _maxB; //Extremum of bezier coefficients
- int _depth;
+ const int _depth;
public:
_CoeffData(int depth) : _minL(0), _maxL(0), _minB(0), _maxB(0),
@@ -78,7 +78,7 @@ private:
const fullVector<double> _coeffsJacDet;
const fullMatrix<double> _coeffsJacMat;
const bezierBasis *_bfsDet, *_bfsMat;
- int _type;
+ const int _type;
public:
_CoeffDataIGE(fullVector<double> &det,
--
GitLab