From 6e2dbc81b99759f68c16ac06fb3383f1f003a806 Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Fri, 15 Apr 2016 07:23:06 +0000
Subject: [PATCH] clean up

---
 Mesh/qualityMeasuresJacobian.cpp | 110 ++++---------
 Mesh/qualityMeasuresJacobian.h   |  21 ++-
 Numeric/MetricBasis.cpp          | 271 -------------------------------
 Numeric/MetricBasis.h            |  12 --
 4 files changed, 39 insertions(+), 375 deletions(-)

diff --git a/Mesh/qualityMeasuresJacobian.cpp b/Mesh/qualityMeasuresJacobian.cpp
index 92950cddae..bad7e09e03 100644
--- a/Mesh/qualityMeasuresJacobian.cpp
+++ b/Mesh/qualityMeasuresJacobian.cpp
@@ -118,11 +118,7 @@ double minAnisotropyMeasure(MElement *el, int n)
 {
   if (_CoeffDataAnisotropy::noMoreComputationPlease()) return -9;
 
-//  double minSampled = minSampledAnisotropyMeasure(el, n);
-
-//  double values[3];
-//  Msg::Info("R %g", el->getEigenvaluesMetric(-1,-1,-1,values));
-//  Msg::Info("R %g", el->getEigenvaluesMetric(-1,1,1,values));
+//  double minSampled = minSampledAnisotropyMeasure(el, n); //fordebug
 
   fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
   el->getNodesCoord(nodesXYZ);
@@ -183,22 +179,17 @@ double minAnisotropyMeasure(MElement *el, int n)
     }
   }
 
-  //
-//  TMBEG(0, _CoeffDataAnisotropy)
   std::vector<_CoeffData*> domains;
   domains.push_back(
       new _CoeffDataAnisotropy(coeffJacDetBez, coeffMetricBez,
                                jacBasis ? jacBasis->getBezier() : NULL,
                                gradBasis->getBezier(), 0, el)
   );
-//  TMEND(0)
 
-  static int i = 0;
-  ++i;
   _subdivideDomains(domains);
-  if (domains.size()/7+1 > 20) {
-    Msg::Info("%d: %d", el->getNum(), domains.size()/7+1);
-  }
+//  if (domains.size()/7+1 > 20) {//fordebug
+//    Msg::Info("%d: %d", el->getNum(), domains.size()/7+1);
+//  }
 
   double min = domains[0]->minB();
   delete domains[0];
@@ -207,30 +198,30 @@ double minAnisotropyMeasure(MElement *el, int n)
     delete domains[i];
   }
 
-//  if (min - minSampled > 1e-13) {
+//  if (min - minSampled > 1e-13) {//fordebug
 //    Msg::Error("no, works no el %d (%g vs %g, diff %g)", el->getNum(), min, minSampled, min-minSampled);
 //  }
 
   return min;
 }
 
-double minSampledAnisotropyMeasure(MElement *el, int deg, bool writeInFile)
-{
-  fullMatrix<double> samplingPoints;
-  bool serendip = false;
-  gmshGeneratePoints(FuncSpaceData(el, deg, &serendip), samplingPoints);
-
-  double values[3];
-  double min = std::numeric_limits<double>::infinity();
-  for (int i = 0; i < samplingPoints.size1(); ++i) {
-    min = std::min(min, el->getEigenvaluesMetric(samplingPoints(i, 0),
-                                                 samplingPoints(i, 1),
-                                                 samplingPoints(i, 2),
-                                                 values));
-  }
-
-  return min;
-}
+//double minSampledAnisotropyMeasure(MElement *el, int deg, bool writeInFile)//fordebug
+//{
+//  fullMatrix<double> samplingPoints;
+//  bool serendip = false;
+//  gmshGeneratePoints(FuncSpaceData(el, deg, &serendip), samplingPoints);
+//
+//  double values[3];
+//  double min = std::numeric_limits<double>::infinity();
+//  for (int i = 0; i < samplingPoints.size1(); ++i) {
+//    min = std::min(min, el->getEigenvaluesMetric(samplingPoints(i, 0),
+//                                                 samplingPoints(i, 1),
+//                                                 samplingPoints(i, 2),
+//                                                 values));
+//  }
+//
+//  return min;
+//}
 
 // Virtual class _CoeffData
 bool _lessMinB::operator()(_CoeffData *cd1, _CoeffData *cd2) const
@@ -406,7 +397,6 @@ _CoeffDataAnisotropy::_CoeffDataAnisotropy(fullVector<double> &det,
   _coeffsMetric(metric.getDataPtr(), metric.size1(), metric.size2()),
   _bfsDet(bfsDet), _bfsMet(bfsMet), _elForDebug(el), _numForDebug(num)
 {
-//  Msg::Info("new");
   if (!det.getOwnData() || !metric.getOwnData()) {
     Msg::Fatal("Cannot create an instance of _CoeffDataScaledJac from a "
                "fullVector or a fullMatrix that does not own its data.");
@@ -419,14 +409,12 @@ _CoeffDataAnisotropy::_CoeffDataAnisotropy(fullVector<double> &det,
   const_cast<fullVector<double>&>(_coeffsJacDet).setOwnData(true);
   const_cast<fullMatrix<double>&>(_coeffsMetric).setOwnData(true);
 
-  _computeAtCorner(_minL, _maxL);;
+  _computeAtCorner(_minL, _maxL);
 
-//  TMBEG(0, _computeLowerBound)
   if (_minL < 1e-3) _minB = 0;
   else _minB = _computeLowerBound();
-//  TMEND(0)
   // computation of _maxB not implemented for now
-  //statsForMatlab(_elForDebug, _numForDebug);
+  //statsForMatlab(_elForDebug, _numForDebug);//fordebug
 }
 
 void _CoeffDataAnisotropy::_computeAtCorner(double &min, double &max) const
@@ -442,14 +430,10 @@ void _CoeffDataAnisotropy::_computeAtCorner(double &min, double &max) const
     }
     p = std::sqrt(p);
     double qualSqrd;
-//      Msg::Info("qpJ: (%g, %g, %g) => aK: (%g, %g)", q, p, _coeffsJacDet(i), q/p, _coeffsJacDet(i)*_coeffsJacDet(i)/p/p/p);
     if (_coeffsMetric.size2() == 3)
       qualSqrd = _R2Dsafe(q, p);
-    else {
+    else
       qualSqrd = _R3Dsafe(q, p, _coeffsJacDet(i));
-//      Msg::Info("qpJ:  (%g, %g, %g) => aK: (%g, %g) => R = %g", q, p, _coeffsJacDet(i), q/p, _coeffsJacDet(i)*_coeffsJacDet(i)/p/p/p, std::sqrt(_R3Dsafe(q, p, _coeffsJacDet(i))));
-//      Msg::Info("(q^2, p^2) = (%g, %g)", q*q, p*p);
-    }
     min = std::min(min, qualSqrd);
     max = std::max(max, qualSqrd);
   }
@@ -466,21 +450,13 @@ double _CoeffDataAnisotropy::_computeLowerBound() const
   }
 
   // 3D element
-//  TMBEG(0, _computeLowerBoundA)
   double mina = _computeLowerBoundA();
-//  TMEND(0)
-//  TMBEG(1, _computeLowerBoundK)
   double minK = _computeLowerBoundK();
-//  TMEND(1)
 
-//  TMBEG(2, _moveInsideDomain)
   _moveInsideDomain(mina, minK, true);
-//  TMEND(2)
 
   double p_dRda, p_dRdK;
-//  TMBEG(3, _computePseudoDerivatives)
   _computePseudoDerivatives(mina, minK, p_dRda, p_dRdK);
-//  TMEND(3)
 
   if (p_dRda < 0) {
     // cases 1 & 2 => compute a lower bounding curve K = beta * a^3 + c * a
@@ -490,10 +466,8 @@ double _CoeffDataAnisotropy::_computeLowerBound() const
     double gamma = .5*(3*minK/mina - slope);
     double beta;
 
-//    TMBEG(4, 12_computeBoundingCurve)
     _computeBoundingCurve(beta, gamma, true);
-//    TMEND(4)
-    /*{
+    /*{//fordebug
       double beta = (minK/mina-c)/mina/mina;
       if (std::abs(p_dRda + p_dRdK * (3*beta*mina*mina+c)) > 1e-12) {
         Msg::Error("Derivative along curve not zero %g", p_dRda + p_dRdK * (3*beta*mina*mina+c));
@@ -515,7 +489,7 @@ double _CoeffDataAnisotropy::_computeLowerBound() const
     else {
       double p_dRda, p_dRdK;
       _computePseudoDerivatives(a_int, K_int, p_dRda, p_dRdK);
-      /*if (p_dRda + p_dRdK * (3*beta*a_int*a_int+c) < -1e-5) {
+      /*if (p_dRda + p_dRdK * (3*beta*a_int*a_int+c) < -1e-5) {//fordebug
         Msg::Error("Derivative along curve not positive or zero %g", p_dRda + p_dRdK * (3*beta*a_int*a_int+c));
         Msg::Info("(%g, %g) vs (%g, %g), diff (%g, %g)", mina, minK, a_int, K_int, a_int-mina, K_int-minK);
         double beta2 = (minK/mina-c)/mina/mina;
@@ -539,9 +513,7 @@ double _CoeffDataAnisotropy::_computeLowerBound() const
     double slope = -p_dRda/p_dRdK;
     double gamma = .5*(3*minK/mina - slope);
     double beta;
-//    TMBEG(5, 45_computeBoundingCurve)
     _computeBoundingCurve(beta, gamma, false);
-//    TMEND(5)
 
     double a_int = mina, K_int = minK;
     if (_intersectionCurveLeftCorner(beta, gamma, a_int, K_int)) {
@@ -621,55 +593,37 @@ double _CoeffDataAnisotropy::_computeLowerBoundK() const
 void _CoeffDataAnisotropy::_computeBoundingCurve(double &beta, double gamma,
                                                  bool compLowBound) const
 {
-  // compute a lower/upper bounding curve K = beta * a^3 + c * a
-  // with c fixed.
+  // compute a lower/upper bounding curve K = beta * a^3 + gamma * a
+  // with gamma fixed.
 
-//  TMBEG(1, num)
-//  TMBEG(11, J2)
   fullVector<double> J2;
   _bfsDet->getRaiser()->computeCoeff(_coeffsJacDet, _coeffsJacDet, J2);
-//  TMEND(11)
 
-//  TMBEG(12, q)
   fullVector<double> q(_coeffsMetric.size1());
   for (int i = 0; i < q.size(); ++i) {
     q(i) = _coeffsMetric(i, 0);
   }
-//  TMEND(12)
 
-//  TMBEG(13, qp2)
   fullVector<double> qp2;
   {
-//    TMBEG(131, qp2A)
     fullMatrix<double> terms_p, terms_qp2;
     terms_p.setAsProxy(_coeffsMetric, 1, _coeffsMetric.size2()-1);
     _bfsMet->getRaiser()->computeCoeff(q, terms_p, terms_p, terms_qp2);
-//    TMEND(131)
-//    TMBEG(132, qp2B)
     qp2.resize(terms_qp2.size1(), true);
     for (int i = 0; i < terms_qp2.size1(); ++i) {
       for (int j = 0; j < terms_qp2.size2(); ++j) {
         qp2(i) += terms_qp2(i, j);
       }
     }
-//    TMEND(132)
   }
-//  TMEND(13)
 
-//  TMBEG(14, J2-gqp2)
   fullVector<double> &coeffNumerator = J2;
   coeffNumerator.axpy(qp2, -gamma);
-//  TMEND(14)
-//  TMEND(1)
 
-//  TMBEG(2, den)
   fullVector<double> coeffDenominator;
   _bfsMet->getRaiser()->computeCoeff(q, q, q, coeffDenominator);
-//  TMEND(2)
 
-//  TMBEG(3, _computeBoundRationalBoundingCurve)
   beta = _computeBoundRational(coeffNumerator, coeffDenominator, compLowBound);
-//  TMEND(3)
 }
 
 void _CoeffDataAnisotropy::fillMetricCoeff(const GradientBasis *gradients,
@@ -901,7 +855,6 @@ bool _CoeffDataAnisotropy::_intersectionCurveLeftCorner(double beta, double gamm
 
 double _CoeffDataAnisotropy::_computeAbscissaMinR(double aStart, double K)
 {
-//  Msg::Info("In _computeAbscissaMinR... (%g, %g)", aStart, K);
   double a = aStart;
   double a3 = (3-a*a)*a;
   double da3 = std::numeric_limits<double>::infinity();
@@ -915,7 +868,6 @@ double _CoeffDataAnisotropy::_computeAbscissaMinR(double aStart, double K)
     da3 -= a3;
     // Make sure we are not going outside the domain:
     while (std::abs(_w(a+da, K)) > 1) da /= 2;
-//    Msg::Info("%g + %g = %g", a, da, a + da);
     a += da;
   }
   return a;
@@ -1013,10 +965,6 @@ double _CoeffDataAnisotropy::_R3Dsafe(double a, double K)
     }
 
     const double w = _wSafe(a, K);
-//    if (std::abs(w) > 1-1e-5) {
-//      if (w > 0) return (a - 1) / (a + 2);
-//      else       return (a - 2) / (a + 1);
-//    }
 
     const double phi = std::acos(w) / 3;
     const double R = (a + 2*std::cos(phi + 2*M_PI/3)) / (a + 2*std::cos(phi));
diff --git a/Mesh/qualityMeasuresJacobian.h b/Mesh/qualityMeasuresJacobian.h
index 29c753707a..2d4b71f3a7 100644
--- a/Mesh/qualityMeasuresJacobian.h
+++ b/Mesh/qualityMeasuresJacobian.h
@@ -18,8 +18,8 @@ namespace jacobianBasedQuality {
 void minMaxJacobianDeterminant(MElement *el, double &min, double &max);
 double minScaledJacobian(MElement *el);
 double minAnisotropyMeasure(MElement *el, int n = 5);
-double minSampledAnisotropyMeasure(MElement *el, int order,
-                                   bool writeInFile = false);
+//double minSampledAnisotropyMeasure(MElement *el, int order,//fordebug
+//                                   bool writeInFile = false);
 
 class _CoeffData
 {
@@ -33,7 +33,6 @@ public:
                          _depth(depth) {}
   virtual ~_CoeffData() {}
 
-  //inline int depth() const {return _depthSub;}
   inline double minL() const {return _minL;}
   inline double maxL() const {return _maxL;}
   inline double minB() const {return _minB;}
@@ -95,12 +94,12 @@ private:
   const bezierBasis *_bfsDet, *_bfsMet;
   MElement *_elForDebug;
   int _numForDebug;
-  static bool hasError;
-  static std::fstream file;
+  static bool hasError;//fordebug
+  static std::fstream file;//fordebug
 
 public:
-  static std::vector<double> mytimes;
-  static std::vector<int> mynumbers;
+  static std::vector<double> mytimes;//fordebug
+  static std::vector<int> mynumbers;//fordebug
   _CoeffDataAnisotropy(fullVector<double> &det,
                        fullMatrix<double> &metric,
                        const bezierBasis *bfsDet,
@@ -117,8 +116,8 @@ public:
                               fullMatrix<double> &coeff,
                               int dim, bool ideal);
 
-  static bool noMoreComputationPlease() {return hasError;}
-  static void youReceivedAnError() {hasError = true;}
+  static bool noMoreComputationPlease() {return hasError;}//fordebug
+  static void youReceivedAnError() {hasError = true;}//fordebug
 
 private:
   void _computeAtCorner(double &min, double &max) const;
@@ -148,8 +147,8 @@ private:
   }
 
 public:
-  void interpolateForMatlab(const fullMatrix<double> &nodes) const;
-  void statsForMatlab(MElement *el, int) const;
+  void interpolateForMatlab(const fullMatrix<double> &nodes) const;//fordebug
+  void statsForMatlab(MElement *el, int) const;//fordebug
 };
 
 //inline bool isValid(MElement *el) {
diff --git a/Numeric/MetricBasis.cpp b/Numeric/MetricBasis.cpp
index 3b15bf71b1..7294c03275 100644
--- a/Numeric/MetricBasis.cpp
+++ b/Numeric/MetricBasis.cpp
@@ -1676,277 +1676,6 @@ void MetricBasis::_minKsharp(const fullMatrix<double> &coeff,
     minK = std::max(.0, upperBound);
 }
 
-/*void MetricBasis::_maxAstKpos(const fullMatrix<double> &coeff,
-    const fullVector<double> &jac, double minK, double beta, double &maxa) const
-{
-  Msg::Fatal("Verifie si cette fonction est ok");
-  fullVector<double> P(coeff.size1());
-  for (int i = 0; i < coeff.size1(); ++i) {
-    P(i) = 0;
-    for (int l = 1; l < 7; ++l) {
-      P(i) += coeff(i, l) * coeff(i, l);
-    }
-    P(i) = std::sqrt(P(i));
-  }
-
-  double min = 1e10;
-
-  std::map<int, std::vector<IneqData> >::const_iterator itJ, itP;
-  itJ = _ineqJ2.begin();
-  itP = _ineqP3.begin();
-
-  while (itJ != _ineqJ2.end() && itP != _ineqP3.end()) {
-    double num = 0, den = 0;
-    for (unsigned int l = 0; l < itJ->second.size(); ++l) {
-      const int i = itJ->second[l].i;
-      const int j = itJ->second[l].j;
-      num += itJ->second[l].val * jac(i) * jac(j);
-    }
-    num *= beta;
-    for (unsigned int l = 0; l < itP->second.size(); ++l) {
-      const int i = itP->second[l].i;
-      const int j = itP->second[l].j;
-      const int k = itP->second[l].k;
-      num -= itP->second[l].val * coeff(i, 0) * coeff(j, 0) * coeff(k, 0);
-      den += itP->second[l].val * P(i) * P(j) * P(k);
-    }
-    min = std::min(min, num/den);
-    ++itJ;
-    ++itP;
-  }
-
-  maxa = std::pow(beta*minK-min, 1/3.);
-}
-
-void MetricBasis::_maxAstKneg(const fullMatrix<double> &coeff,
-    const fullVector<double> &jac, double minK, double beta, double &maxa) const
-{
-  Msg::Fatal("Verifie si cette fonction est ok");
-  fullVector<double> P(coeff.size1());
-  fullMatrix<double> Q(coeff.size1(), coeff.size1());
-  for (int i = 0; i < coeff.size1(); ++i) {
-    P(i) = 0;
-    for (int l = 1; l < 7; ++l) {
-      P(i) += coeff(i, l) * coeff(i, l);
-    }
-    P(i) = std::sqrt(P(i));
-    for (int j = 0; j < coeff.size1(); ++j) {
-      Q(i, j) = 0;
-      for (int l = 1; l < 7; ++l) {
-        Q(i, j) += coeff(i, l) * coeff(j, l);
-      }
-    }
-  }
-
-  double min = 1e10;
-
-  std::map<int, std::vector<IneqData> >::const_iterator itJ, itP;
-  itJ = _ineqJ2.begin();
-  itP = _ineqP3.begin();
-
-  while (itJ != _ineqJ2.end() && itP != _ineqP3.end()) {
-    double num = 0, den = 0;
-    for (unsigned int l = 0; l < itJ->second.size(); ++l) {
-      const int i = itJ->second[l].i;
-      const int j = itJ->second[l].j;
-      num += itJ->second[l].val * jac(i) * jac(j);
-    }
-    num *= beta;
-    for (unsigned int l = 0; l < itP->second.size(); ++l) {
-      const int i = itP->second[l].i;
-      const int j = itP->second[l].j;
-      const int k = itP->second[l].k;
-      num -= itP->second[l].val * coeff(i, 0) * coeff(j, 0) * coeff(k, 0);
-      double tmp = P(i) * Q(j, k);
-      tmp = std::min(tmp, P(j) * Q(i, k));
-      tmp = std::min(tmp, P(k) * Q(i, j));
-      den += itP->second[l].val * tmp;
-    }
-    min = std::min(min, num/den);
-    ++itJ;
-    ++itP;
-  }
-
-  maxa = std::pow(beta*minK-min, 1/3.);
-}
-
-void MetricBasis::_maxKstAfast(const fullMatrix<double> &coeff,
-    const fullVector<double> &jac, double mina, double beta, double &maxK) const
-{
-  Msg::Fatal("Verifie si cette fonction est ok");
-  fullVector<double> r(coeff.size1());
-  for (int i = 0; i < coeff.size1(); ++i) {
-    r(i) = 0;
-    for (int l = 1; l < 7; ++l) {
-      r(i) += coeff(i, l) * coeff(i, l);
-    }
-    r(i) = std::sqrt(r(i));
-  }
-
-  double min = 1e10;
-
-  std::map<int, std::vector<IneqData> >::const_iterator itJ, itP;
-  itJ = _ineqJ2.begin();
-  itP = _ineqP3.begin();
-
-  while (itJ != _ineqJ2.end() && itP != _ineqP3.end()) {
-    double num = 0, den = 0;
-    for (unsigned int l = 0; l < itJ->second.size(); ++l) {
-      const int i = itJ->second[l].i;
-      const int j = itJ->second[l].j;
-      num -= itJ->second[l].val * jac(i) * jac(j);
-    }
-    num *= beta;
-    for (unsigned int l = 0; l < itP->second.size(); ++l) {
-      const int i = itP->second[l].i;
-      const int j = itP->second[l].j;
-      const int k = itP->second[l].k;
-      num += itP->second[l].val * coeff(i, 0) * coeff(j, 0) * coeff(k, 0);
-      den += itP->second[l].val * r(i) * r(j) * r(k);
-    }
-    min = std::min(min, num/den);
-    ++itJ;
-    ++itP;
-  }
-
-  maxK = 1/beta*(mina*mina*mina-min);
-}
-
-void MetricBasis::_maxKstAsharp(const fullMatrix<double> &coeff,
-    const fullVector<double> &jac, double mina, double beta, double &maxK) const
-{
-  Msg::Fatal("Verifie si cette fonction est ok");
-  fullVector<double> P(coeff.size1());
-  fullMatrix<double> Q(coeff.size1(), coeff.size1());
-  for (int i = 0; i < coeff.size1(); ++i) {
-    P(i) = 0;
-    for (int l = 1; l < 7; ++l) {
-      P(i) += coeff(i, l) * coeff(i, l);
-    }
-    P(i) = std::sqrt(P(i));
-    // fill only the half
-    for (int j = i; j < coeff.size1(); ++j) {
-      Q(i, j) = 0;
-      for (int l = 1; l < 7; ++l) {
-        Q(i, j) += coeff(i, l) * coeff(j, l);
-      }
-    }
-  }
-
-  double min = 1e10;
-
-  std::map<int, std::vector<IneqData> >::const_iterator itJ, itP;
-  itJ = _ineqJ2.begin();
-  itP = _ineqP3.begin();
-
-  while (itJ != _ineqJ2.end() && itP != _ineqP3.end()) {
-    double num = 0, den = 0;
-    for (unsigned int l = 0; l < itJ->second.size(); ++l) {
-      const int i = itJ->second[l].i;
-      const int j = itJ->second[l].j;
-      num -= itJ->second[l].val * jac(i) * jac(j);
-    }
-    num *= beta;
-    for (unsigned int l = 0; l < itP->second.size(); ++l) {
-      const int i = itP->second[l].i;
-      const int j = itP->second[l].j;
-      const int k = itP->second[l].k;
-      num += itP->second[l].val * coeff(i, 0) * coeff(j, 0) * coeff(k, 0);
-      if (j == k) // TODO verif ok
-        den += itP->second[l].val * Q(i,i) * P(i);
-      else
-        den += itP->second[l].val * 1/3*(Q(i,j)*P(k)+Q(i,k)*P(j)+Q(j,k)*P(i));
-    }
-    min = std::min(min, num/den);
-    ++itJ;
-    ++itP;
-  }
-
-  maxK = 1/beta*(mina*mina*mina-min);
-}
-
-void MetricBasis::_computeBoundBeta(const fullMatrix<double> &coeff,
-                                    const fullVector<double> &jac,
-                                    double &beta, bool compLowBound) const
-{
-  // compute a lower/upper bound for function beta = K/a^3
-
-  double upperBound = std::numeric_limits<double>::infinity();
-  double lowerBound = -upperBound;
-
-  // value in case of a failure
-  if (compLowBound) beta = 0;
-  else              beta = 1;
-
-  std::map<int, std::vector<IneqData> >::const_iterator itJ, itP;
-  itJ = _ineqJ2.begin();
-  itP = _ineqP3.begin();
-
-  // _ineqJ2 and _ineqP3 have the same size
-  while (itJ != _ineqJ2.end()) {
-    double num = 0;
-    for (unsigned int k = 0; k < itJ->second.size(); ++k) {
-      const int i = itJ->second[k].i;
-      const int j = itJ->second[k].j;
-      num += itJ->second[k].val * jac(i) * jac(j);
-    }
-
-    double den = 0;
-    for (unsigned int l = 0; l < itP->second.size(); ++l) {
-      const int i = itP->second[l].i;
-      const int j = itP->second[l].j;
-      const int k = itP->second[l].k;
-      den += itP->second[l].val * coeff(i, 0) * coeff(j, 0) * coeff(k, 0);
-    }
-
-    // beta has to satisfy beta * den <= num if compLowBound = true
-    //                  or beta * den >= num if compLowBound = false
-    if (den == 0) {
-      if (compLowBound) {
-        if (num >= 0) {
-          ++itJ;
-          ++itP;
-          continue;
-        }
-        else return;
-        //TODO c'est mieux de retourner min(J^2)/max(q^3) ??
-      }
-      else {
-        if (num <= 0) {
-          ++itJ;
-          ++itP;
-          continue;
-        }
-        else return;
-        //TODO c'est mieux de retourner max(J^2)/min(q^3) ??
-      }
-    }
-    if (den > 0) {
-      if (compLowBound)
-        upperBound = std::min(upperBound, num/den);
-      else
-        lowerBound = std::max(lowerBound, num/den);
-    }
-    else {
-      if (compLowBound)
-        lowerBound = std::max(lowerBound, num/den);
-      else
-        upperBound = std::min(upperBound, num/den);
-    }
-
-    ++itJ;
-    ++itP;
-  }
-
-  if (lowerBound <= upperBound) {
-    if (compLowBound)
-      beta = std::max(.0, upperBound);
-    else
-      beta = std::min(1., lowerBound);
-  }
-   //TODO sinon, c'est mieux de retourner m..(J^2)/m..(q^3) ??
-}*/
-
 void MetricBasis::_computeBoundingCurve(const fullMatrix<double> &coeff,
                                         const fullVector<double> &jac,
                                         double &beta, double c,
diff --git a/Numeric/MetricBasis.h b/Numeric/MetricBasis.h
index 968e5d49fb..5d0a389f96 100644
--- a/Numeric/MetricBasis.h
+++ b/Numeric/MetricBasis.h
@@ -112,7 +112,6 @@ public:
                                                 double toleranceTensor,
                                                 double tolerance);
   void statsForMatlab(MElement *el, int deg, MetricData *md) const;
-  //void interpolate(const MElement*, const MetricData*, const double *uvw, double *minmaxQ) const;
   void interpolate(const MElement*, const MetricData*,
                    const fullMatrix<double> &nodes, fullMatrix<double> &R) const;
   void interpolateAfterNSubdivisions(const MElement *el,
@@ -143,17 +142,6 @@ private:
   void _minAfast(const fullMatrix<double>&, double &min) const;
   void _minKfast(const fullMatrix<double>&, const fullVector<double>&, double &min) const;
   void _minKsharp(const fullMatrix<double>&, const fullVector<double>&, double &min) const;
-  /*void _maxAstKpos(const fullMatrix<double>&, const fullVector<double>&,
-                 double minK, double beta, double &maxa) const;
-  void _maxAstKneg(const fullMatrix<double>&, const fullVector<double>&,
-                 double minK, double beta, double &maxa) const;
-  void _maxKstAfast(const fullMatrix<double>&, const fullVector<double>&,
-                 double mina, double beta, double &maxK) const;
-  void _maxKstAsharp(const fullMatrix<double>&, const fullVector<double>&,
-                 double mina, double beta, double &maxK) const;
-  void _computeBoundBeta(const fullMatrix<double>&,
-                         const fullVector<double>&,
-                         double &beta, bool lowerBound) const;*/
   void _computeBoundingCurve(const fullMatrix<double>&,
                              const fullVector<double>&,
                              double &beta, double c, bool lowerBound) const;
-- 
GitLab