From d0ac51dd8ab278b4524f360a489a50c3430656fb Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Fri, 22 Aug 2014 13:20:28 +0000
Subject: [PATCH] for simplexe p1, metricCorner sample only at 1 point

---
 Geo/MElement.cpp        |  5 +++++
 Geo/MElement.h          |  1 +
 Numeric/MetricBasis.cpp | 32 +++++++++++++++++++++++++++++---
 Numeric/MetricBasis.h   |  3 +++
 4 files changed, 38 insertions(+), 3 deletions(-)

diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 40b6164f79..428a11c208 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 90772a261f..480b468221 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 f484e683ee..3d534828ee 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 2202b903fc..63cf8f11c8 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);
-- 
GitLab