diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index a4e360daaa272678f1c3363ce209a9e85db74e7c..c85fe0d7908dee82fd193808c17c51733e33ae86 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -6,6 +6,8 @@
#ifndef _GMSH_DEFINES_H_
#define _GMSH_DEFINES_H_
+//#define METRICSHAPEMEASURE
+
// IO file formats (numbers should not be changed)
#define FORMAT_MSH 1
#define FORMAT_UNV 2
diff --git a/Fltk/graphicWindow.cpp b/Fltk/graphicWindow.cpp
index 9ec0a454626f20e7f043129bdf827e97255284c4..3909284927e03dca7900734da84211f3be175764 100644
--- a/Fltk/graphicWindow.cpp
+++ b/Fltk/graphicWindow.cpp
@@ -1624,7 +1624,11 @@ static std::vector<std::string> getInfoStrings(MElement *ele)
{
std::ostringstream sstream;
sstream << " Quality: "
+#ifdef METRICSHAPEMEASURE
+ << "metric = " << ele->metricShapeMeasure() << " "
+#else
<< "rho = " << ele->rhoShapeMeasure() << " "
+#endif
<< "gamma = " << ele->gammaShapeMeasure() << " "
<< "eta = " << ele->etaShapeMeasure();
info.push_back(sstream.str());
diff --git a/Fltk/statisticsWindow.cpp b/Fltk/statisticsWindow.cpp
index 455227cd197ae6c2ba3879bc120e278094844750..0feddfb744c59b17dea9ac307d526297369281c6 100644
--- a/Fltk/statisticsWindow.cpp
+++ b/Fltk/statisticsWindow.cpp
@@ -75,7 +75,14 @@ static void statistics_histogram_cb(Fl_Widget *w, void *data)
else if(name == "Eta3D")
d[e->getNum()].push_back(e->etaShapeMeasure());
else if(name == "Rho3D")
- d[e->getNum()].push_back(e->rhoShapeMeasure());
+#ifdef METRICSHAPEMEASURE
+ if (e->getDim() == 3)
+ d[e->getNum()].push_back(e->metricShapeMeasure());
+ else
+ d[e->getNum()].push_back(1);
+#else
+ d[e->getNum()].push_back(e->rhoShapeMeasure());
+#endif
else
d[e->getNum()].push_back(e->distoShapeMeasure());
}
diff --git a/Geo/GModelVertexArrays.cpp b/Geo/GModelVertexArrays.cpp
index 52f7600f90e25ae0710046877bbc3b5ee27dffff..bae77544602b4d41370fe0d9db7929ca109ecaf4 100644
--- a/Geo/GModelVertexArrays.cpp
+++ b/Geo/GModelVertexArrays.cpp
@@ -104,7 +104,11 @@ bool isElementVisible(MElement *ele)
if(CTX::instance()->mesh.qualityType == 3)
q = ele->distoShapeMeasure();
else if(CTX::instance()->mesh.qualityType == 2)
+#ifdef METRICSHAPEMEASURE
+ q = ele->metricShapeMeasure();
+#else
q = ele->rhoShapeMeasure();
+#endif
else if(CTX::instance()->mesh.qualityType == 1)
q = ele->etaShapeMeasure();
else
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 108be656e24437a6fade04c53840970a65c5cb49..40b6164f7987041bb510727cf865a07b4c614d46 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -22,6 +22,7 @@
#include "GEntity.h"
#include "StringUtils.h"
#include "Numeric.h"
+#include "MetricBasis.h"
#include "Context.h"
#define SQU(a) ((a)*(a))
@@ -112,6 +113,11 @@ double MElement::rhoShapeMeasure()
return 0.;
}
+double MElement::metricShapeMeasure()
+{
+ return MetricBasis::minRCorner(this);
+}
+
double MElement::maxDistToStraight()
{
const nodalBasis *lagBasis = getFunctionSpace();
@@ -924,7 +930,11 @@ void MElement::writePOS(FILE *fp, bool printElementary, bool printElementNumber,
}
}
if(printRho){
+#ifdef METRICSHAPEMEASURE
+ double rho = metricShapeMeasure();
+#else
double rho = rhoShapeMeasure();
+#endif
for(int i = 0; i < n; i++){
if(first) first = false; else fprintf(fp, ",");
fprintf(fp, "%g", rho);
diff --git a/Geo/MElement.h b/Geo/MElement.h
index db684dbd13b84267e0771d43e201ff9183e12303..90772a261fa47851a4fc5c25b7d53dc97aaf5ca6 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -191,6 +191,7 @@ class MElement
}
virtual double angleShapeMeasure() { return 1.0; }
virtual void scaledJacRange(double &jmin, double &jmax, GEntity *ge = 0);
+ virtual double metricShapeMeasure();
// get the radius of the inscribed circle/sphere if it exists,
// otherwise get the minimum radius of all the circles/spheres
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index c1fc939494ed6d2386b90273709ed74a070dfa55..6c3999989ccd48cf172497b41999c3f41e19d6a4 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -64,10 +64,19 @@ static void GetQualityMeasure(std::vector<T*> &ele,
eta += e;
etaMin = std::min(etaMin, e);
etaMax = std::max(etaMax, e);
+#ifdef METRICSHAPEMEASURE
+ static int aa = 0;
+ if (++aa == 1) Msg::Warning("Computing metric shape measure instead of rho shape measure !");
+ double r = ele[i]->metricShapeMeasure();
+ rho += r;
+ rhoMin = std::min(rhoMin, r);
+ rhoMax = std::max(rhoMax, r);
+#else
double r = ele[i]->rhoShapeMeasure();
rho += r;
rhoMin = std::min(rhoMin, r);
rhoMax = std::max(rhoMax, r);
+#endif
double jmin,jmax; ele[i]->scaledJacRange(jmin,jmax);
double d = jmin;
disto += d;
diff --git a/Numeric/BasisFactory.cpp b/Numeric/BasisFactory.cpp
index d39b63f57cc0050abd138748d70da6c3da6ded64..b2906da4c8888cafd4297ab46bc9f61f8762c3a6 100644
--- a/Numeric/BasisFactory.cpp
+++ b/Numeric/BasisFactory.cpp
@@ -13,10 +13,10 @@
#include "MElement.h"
std::map<int, nodalBasis*> BasisFactory::fs;
-std::map<int, JacobianBasis*> BasisFactory::js;
std::map<int, MetricBasis*> BasisFactory::ms;
BasisFactory::Cont_bezierBasis BasisFactory::bs;
BasisFactory::Cont_gradBasis BasisFactory::gs;
+BasisFactory::Cont_jacBasis BasisFactory::js;
const nodalBasis* BasisFactory::getNodalBasis(int tag)
{
@@ -67,14 +67,15 @@ const nodalBasis* BasisFactory::getNodalBasis(int tag)
return inserted.first->second;
}
-const JacobianBasis* BasisFactory::getJacobianBasis(int tag)
+const JacobianBasis* BasisFactory::getJacobianBasis(int tag, int order)
{
- std::map<int, JacobianBasis*>::const_iterator it = js.find(tag);
+ std::pair<int, int> key(tag, order);
+ Cont_jacBasis::iterator it = js.find(key);
if (it != js.end())
return it->second;
- JacobianBasis* J = new JacobianBasis(tag);
- js.insert(std::make_pair(tag, J));
+ JacobianBasis* J = new JacobianBasis(tag, order);
+ js.insert(std::make_pair(key, J));
return J;
}
@@ -122,13 +123,27 @@ void BasisFactory::clearAll()
}
fs.clear();
- std::map<int, JacobianBasis*>::iterator itJ = js.begin();
+ std::map<int, MetricBasis*>::iterator itM = ms.begin();
+ while (itM != ms.end()) {
+ delete itM->second;
+ itM++;
+ }
+ ms.clear();
+
+ Cont_jacBasis::iterator itJ = js.begin();
while (itJ != js.end()) {
delete itJ->second;
itJ++;
}
js.clear();
+ Cont_gradBasis::iterator itG = gs.begin();
+ while (itG != gs.end()) {
+ delete itG->second;
+ itG++;
+ }
+ gs.clear();
+
Cont_bezierBasis::iterator itB = bs.begin();
while (itB != bs.end()) {
delete itB->second;
diff --git a/Numeric/BasisFactory.h b/Numeric/BasisFactory.h
index 82b5c3731aecc9ea79abafe0e01ef9a5c6c1b975..f2f4e41ee97f647c315a347a2e6d5655af6482dc 100644
--- a/Numeric/BasisFactory.h
+++ b/Numeric/BasisFactory.h
@@ -15,13 +15,14 @@
class BasisFactory
{
+ typedef std::map<std::pair<int, int>, JacobianBasis*> Cont_jacBasis;
typedef std::map<std::pair<int, int>, bezierBasis*> Cont_bezierBasis;
typedef std::map<std::pair<int, int>, GradientBasis*> Cont_gradBasis;
private:
static std::map<int, nodalBasis*> fs;
- static std::map<int, JacobianBasis*> js;
static std::map<int, MetricBasis*> ms;
+ static Cont_jacBasis js;
static Cont_gradBasis gs;
static Cont_bezierBasis bs;
// store bezier bases by parentType and order (no serendipity..)
@@ -29,7 +30,11 @@ class BasisFactory
public:
// Caution: the returned pointer can be NULL
static const nodalBasis* getNodalBasis(int tag);
- static const JacobianBasis* getJacobianBasis(int tag);
+ static const JacobianBasis* getJacobianBasis(int tag, int order);
+ static const JacobianBasis* getJacobianBasis(int tag) {
+ return getJacobianBasis(ElementType::ParentTypeFromTag(tag),
+ ElementType::OrderFromTag(tag) );
+ }
static const MetricBasis* getMetricBasis(int tag);
static const GradientBasis* getGradientBasis(int tag, int order);
diff --git a/Numeric/MetricBasis.cpp b/Numeric/MetricBasis.cpp
index 093a1d3d86f8270795734a65aa1481c92872b6f9..f484e683ee4e080901981474129b35de07a6161d 100644
--- a/Numeric/MetricBasis.cpp
+++ b/Numeric/MetricBasis.cpp
@@ -75,6 +75,107 @@ double MetricBasis::boundMinR(MElement *el)
return metric->getBoundRmin(el, md, dummy);
}
+double MetricBasis::minRCorner(MElement *el)
+{
+ int tag = el->getTypeForMSH();
+ int order = 1;
+ if (el->getType() == TYPE_TRI || el->getType() == TYPE_TET) order = 0;
+
+ const GradientBasis *gradients = BasisFactory::getGradientBasis(tag, 1);
+ const JacobianBasis *jacobian = BasisFactory::getJacobianBasis(tag, 1);
+
+ int nSampPnts = jacobian->getNumJacNodes();
+ if (el->getType() == TYPE_PYR) nSampPnts = 4;
+
+ int nMapping = gradients->getNumMapNodes();
+ fullMatrix<double> nodes(nMapping, 3);
+ el->getNodesCoord(nodes);
+
+ // Metric coefficients
+ fullMatrix<double> metCoeffLag;
+
+ switch (el->getDim()) {
+ case 0 :
+ return -1.;
+ case 1 :
+ case 2 :
+ Msg::Fatal("not implemented");
+ break;
+
+ case 3 :
+ {
+ fullMatrix<double> dxyzdX(nSampPnts,3), dxyzdY(nSampPnts,3), dxyzdZ(nSampPnts,3);
+ gradients->getGradientsFromNodes(nodes, &dxyzdX, &dxyzdY, &dxyzdZ);
+
+ metCoeffLag.resize(nSampPnts, 7);
+ for (int i = 0; i < nSampPnts; i++) {
+ const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
+ const double &dxdY = dxyzdY(i,0), &dydY = dxyzdY(i,1), &dzdY = dxyzdY(i,2);
+ const double &dxdZ = dxyzdZ(i,0), &dydZ = dxyzdZ(i,1), &dzdZ = dxyzdZ(i,2);
+ const double dvxdX = dxdX*dxdX + dydX*dydX + dzdX*dzdX;
+ const double dvxdY = dxdY*dxdY + dydY*dydY + dzdY*dzdY;
+ const double dvxdZ = dxdZ*dxdZ + dydZ*dydZ + dzdZ*dzdZ;
+ metCoeffLag(i, 0) = (dvxdX + dvxdY + dvxdZ) / 3;
+ metCoeffLag(i, 1) = dvxdX - metCoeffLag(i, 0);
+ metCoeffLag(i, 2) = dvxdY - metCoeffLag(i, 0);
+ metCoeffLag(i, 3) = dvxdZ - metCoeffLag(i, 0);
+ const double fact = std::sqrt(2);
+ metCoeffLag(i, 4) = fact * (dxdX*dxdY + dydX*dydY + dzdX*dzdY);
+ metCoeffLag(i, 5) = fact * (dxdZ*dxdY + dydZ*dydY + dzdZ*dzdY);
+ metCoeffLag(i, 6) = fact * (dxdX*dxdZ + dydX*dydZ + dzdX*dzdZ);
+ }
+ }
+ break;
+ }
+
+ // Jacobian coefficients
+ fullVector<double> jacLag(jacobian->getNumJacNodes());
+ jacobian->getSignedJacobian(nodes, jacLag);
+
+ // Compute min_corner(R)
+ double Rmin = 1.;
+
+ for (int i = 0; i < nSampPnts; ++i) {
+ if (jacLag(i) <= 0.) {
+ Rmin = std::min(Rmin, 0.);
+ }
+ else {
+ double q = metCoeffLag(i, 0);
+ double p = 0;
+ for (int k = 1; k < 7; ++k) {
+ p += pow_int(metCoeffLag(i, k), 2);
+ }
+ p = std::sqrt(p/6);
+ const double a = q/p;
+ if (a > 1e4) {
+ Rmin = std::min(Rmin, std::sqrt((a - std::sqrt(3)) / (a + std::sqrt(3))));
+ }
+ else {
+ const double x = .5 * (jacLag(i)/p/p*jacLag(i)/p - a*a*a + 3*a);
+ if (x > 1.1 || x < -1.1) {
+ Msg::Error("x much smaller than -1 or much greater than 1");
+ }
+
+ double tmpR;
+ if (x >= 1)
+ tmpR = (a - 1) / (a + 2);
+ else if (x <= -1)
+ tmpR = (a - 2) / (a + 1);
+ else {
+ const double phi = std::acos(x)/3;
+ tmpR = (a + 2*std::cos(phi + 2*M_PI/3)) / (a + 2*std::cos(phi));
+ }
+ if (tmpR < 0) {
+ if (tmpR < -1e-7) Msg::Error("negative R !!!");
+ else tmpR = 0;
+ }
+ Rmin = std::min(Rmin, std::sqrt(tmpR));
+ }
+ }
+ }
+ return Rmin;
+}
+
double MetricBasis::sampleR(MElement *el, int order)
{
MetricBasis *metric = (MetricBasis*)BasisFactory::getMetricBasis(el->getTypeForMSH());
diff --git a/Numeric/MetricBasis.h b/Numeric/MetricBasis.h
index f1cf8e061caa22b7e18c4c471520954738d2bcd9..2202b903fc36af447734a3c2d684001ea808c65c 100644
--- a/Numeric/MetricBasis.h
+++ b/Numeric/MetricBasis.h
@@ -66,6 +66,7 @@ public:
double getMinR(MElement*, MetricData*&, int) const;
bool notStraight(MElement*, double &metric, int order) const;
static double boundMinR(MElement *el);
+ static double minRCorner(MElement *el);
static double sampleR(MElement *el, int order);
//double getBoundRmin(int, MElement**, double*);
//static double boundRmin(int, MElement**, double*, bool sameType = false);