diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp index 40b6164f7987041bb510727cf865a07b4c614d46..428a11c2082209b1d0e44922b5787f54bcdc6f8e 100644 --- a/Geo/MElement.cpp +++ b/Geo/MElement.cpp @@ -118,6 +118,11 @@ double MElement::metricShapeMeasure() return MetricBasis::minRCorner(this); } +double MElement::metricShapeMeasure2() +{ + return MetricBasis::boundMinR(this); +} + double MElement::maxDistToStraight() { const nodalBasis *lagBasis = getFunctionSpace(); diff --git a/Geo/MElement.h b/Geo/MElement.h index 90772a261fa47851a4fc5c25b7d53dc97aaf5ca6..480b468221d6b8a7a663f7bf98976f26ffb95c66 100644 --- a/Geo/MElement.h +++ b/Geo/MElement.h @@ -192,6 +192,7 @@ class MElement virtual double angleShapeMeasure() { return 1.0; } virtual void scaledJacRange(double &jmin, double &jmax, GEntity *ge = 0); virtual double metricShapeMeasure(); + virtual double metricShapeMeasure2(); // get the radius of the inscribed circle/sphere if it exists, // otherwise get the minimum radius of all the circles/spheres diff --git a/Numeric/MetricBasis.cpp b/Numeric/MetricBasis.cpp index f484e683ee4e080901981474129b35de07a6161d..3d534828ee5b6f2c98c3b5f45e73360bf37b3d86 100644 --- a/Numeric/MetricBasis.cpp +++ b/Numeric/MetricBasis.cpp @@ -11,7 +11,12 @@ #include "OS.h" #include <sstream> -double MetricBasis::_tol = 1e-2; +double MetricBasis::_tol = 1e-3; +double MetricBasis::tm0 = 0; +double MetricBasis::tm1 = 0; +double MetricBasis::tm2 = 0; +double MetricBasis::tm3 = 0; +double MetricBasis::tm4 = 0; int MetricBasis::_which = 0; namespace { @@ -77,20 +82,26 @@ double MetricBasis::boundMinR(MElement *el) double MetricBasis::minRCorner(MElement *el) { + //double time = Cpu(); 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); + const GradientBasis *gradients = BasisFactory::getGradientBasis(tag, order); + const JacobianBasis *jacobian = BasisFactory::getJacobianBasis(tag, order); int nSampPnts = jacobian->getNumJacNodes(); if (el->getType() == TYPE_PYR) nSampPnts = 4; int nMapping = gradients->getNumMapNodes(); fullMatrix<double> nodes(nMapping, 3); + //tm0 += Cpu() - time; + //time = Cpu(); el->getNodesCoord(nodes); + //tm1 += Cpu() - time; + //time = Cpu(); + // Metric coefficients fullMatrix<double> metCoeffLag; @@ -128,10 +139,16 @@ double MetricBasis::minRCorner(MElement *el) break; } + //tm2 += Cpu() - time; + //time = Cpu(); + // Jacobian coefficients fullVector<double> jacLag(jacobian->getNumJacNodes()); jacobian->getSignedJacobian(nodes, jacLag); + //tm3 += Cpu() - time; + //time = Cpu(); + // Compute min_corner(R) double Rmin = 1.; @@ -173,6 +190,9 @@ double MetricBasis::minRCorner(MElement *el) } } } + + //tm4 += Cpu() - time; + //time = Cpu(); return Rmin; } @@ -459,7 +479,9 @@ double MetricBasis::getBoundRmin(MElement *el, MetricData *&md, fullMatrix<doubl break; } +#ifndef METRICSHAPEMEASURE lagCoeff = metCoeffLag; +#endif fullMatrix<double> *metCoeff; metCoeff = new fullMatrix<double>(nSampPnts, metCoeffLag.size2()); _bezier->matrixLag2Bez.mult(metCoeffLag, *metCoeff); @@ -489,14 +511,17 @@ double MetricBasis::getBoundRmin(MElement *el, MetricData *&md, fullMatrix<doubl //Msg::Info("el %d", el->getNum()); double mina, maxa; _minMaxA(*metCoeff, mina, maxa); +#ifndef METRICSHAPEMEASURE static int cntRight = 0, cntTOT = 0; ++cntTOT; if (maxa-mina < 1e-10) { ++cntRight; } +#endif //Msg::Info("right %d/%d", cntRight, cntTOT); +#ifndef METRICSHAPEMEASURE fullVector<double> *jjac = new fullVector<double>(*jac); fullMatrix<double> *mmet = new fullMatrix<double>(*metCoeff); /*for (int i = 0; i < jjac->size(); ++i) { @@ -507,6 +532,7 @@ double MetricBasis::getBoundRmin(MElement *el, MetricData *&md, fullMatrix<doubl }*/ md = new MetricData(mmet, jjac, RminBez, 0, 0); //Msg::Info("+1 %d", md); +#endif if (RminLag-RminBez < MetricBasis::_tol) { //Msg::Info("RETURNING %g", RminBez); diff --git a/Numeric/MetricBasis.h b/Numeric/MetricBasis.h index 2202b903fc36af447734a3c2d684001ea808c65c..63cf8f11c8677129a5d180eca6169af3b5d4fe1c 100644 --- a/Numeric/MetricBasis.h +++ b/Numeric/MetricBasis.h @@ -26,6 +26,8 @@ private: std::vector<int> __numSub; MElement *__curElem; + static double tm0, tm1, tm2, tm3, tm4; + std::fstream file; class IneqData { @@ -67,6 +69,7 @@ public: bool notStraight(MElement*, double &metric, int order) const; static double boundMinR(MElement *el); static double minRCorner(MElement *el); + static void printTm() {Msg::Info("%g %g %g %g %g", tm0, tm1, tm2, tm3, tm4);} static double sampleR(MElement *el, int order); //double getBoundRmin(int, MElement**, double*); //static double boundRmin(int, MElement**, double*, bool sameType = false);