MetricBasis: better bounds

Numeric/JacobianBasis.cpp

@@ -88,11 +88,15 @@ void GradientBasis::getGradientsFromNodes(const fullMatrix<double> &nodes,
   if (dxyzdZ) gradShapeMatZ.mult(nodes, *dxyzdZ);
 }
 
-JacobianBasis::JacobianBasis(int tag)
+JacobianBasis::JacobianBasis(int tag, int jacOrder)
 {
   const int parentType = ElementType::ParentTypeFromTag(tag);
   const int order = ElementType::OrderFromTag(tag);
-  const int jacobianOrder = JacobianBasis::jacobianOrder(parentType, order);
+  int jacobianOrder;
+  if (jacOrder < 0)
+    jacobianOrder = JacobianBasis::jacobianOrder(parentType, order);
+  else
+    jacobianOrder = jacOrder;
   const int primJacobianOrder = JacobianBasis::jacobianOrder(parentType, 1);
 
   fullMatrix<double> lagPoints;                                  // Sampling points

Numeric/JacobianBasis.h

@@ -44,7 +44,7 @@ class JacobianBasis {
  public :
   const bezierBasis *bezier;
 
-  JacobianBasis(int tag);
+  JacobianBasis(int tag, int jacOrder = -1);
 
   // Get methods
   inline int getNumJacNodes() const { return numJacNodes; }

Numeric/MetricBasis.h

@@ -11,25 +11,82 @@
 #include "fullMatrix.h"
 
 class MetricCoefficient {
+private:
+  class MetricData {
+   public:
+    fullMatrix<double> *_metcoeffs;
+    fullVector<double> *_jaccoeffs;
+    double _RminBez;
+    int _depth;
+
+   public:
+    MetricData(fullMatrix<double> *m, fullVector<double> *j, double r, int d) :
+      _metcoeffs(m), _jaccoeffs(j), _RminBez(r), _depth(d) {}
+    ~MetricData() {
+      delete _metcoeffs;
+      delete _jaccoeffs;
+    }
+  };
+  struct lessMinB {
+    bool operator()(const MetricData *md1, const MetricData *md2) const {
+      return md1->_RminBez > md2->_RminBez;
+    }
+  };
+
+  class IneqData {
+   public:
+    int i, j, k;
+    double val;
+
+   public:
+    IneqData(double val, int i, int j, int k = -1) : i(i), j(j), k(k), val(val) {}
+  };
+  class gterIneq {
+   public:
+    bool operator()(const IneqData &id1, const IneqData &id2) const {
+      return id1.val > id2.val;
+    }
+  };
+
  private:
   MElement *_element;
   const JacobianBasis *_jacobian;
   const GradientBasis *_gradients;
   const bezierBasis *_bezier;
   fullMatrix<double> _coefficientsLag, _coefficientsBez;
+  fullMatrix<double> _inequations;
+  std::vector<std::map<std::pair<int, int>, double> > _ineq2; //TODO : pas top
+  std::map<int, std::vector<IneqData> > _ineqJ2, _ineqP3;
+  int __maxdepth, __numSubdivision;
+  std::vector<int> __numSub;
 
  public:
   MetricCoefficient(MElement*);
 
   void getCoefficients(fullMatrix<double>&, bool bezier = true);
   void interpolate(const double *uvw, double *minmaxQ);
-  double getBoundRmin(double tol);
+  double getBoundRmin(double tol, int which);
 
   static int metricOrder(int tag);
 
  private:
   void _computeBezCoeff();
+  void _fillInequalities(int order);
+  int _lightenInequalities(double, double, int&, int&); //TODO change
   void _interpolateBezierPyramid(const double *uvw, double *minmaxQ);
+  void _computeRmin1(fullMatrix<double>&, fullVector<double>&,
+                    double &RminLag, double &RminBez, int) const;
+  void _computeRmin2(fullMatrix<double>&, fullVector<double>&,
+                    double &RminLag, double &RminBez, int depth, int which) const;
+  double _subdivideForRmin(MetricData*, double RminLag, double tol, int which) const;
+  double _minp(const fullMatrix<double>&) const;
+  double _minq(const fullMatrix<double>&) const;
+  double _maxp(const fullMatrix<double>&) const;
+  double _maxq(const fullMatrix<double>&) const;
+  void _minMaxA(const fullMatrix<double>&, double &min, double &max) const;
+  void _minMaxA2(const fullMatrix<double>&, double &min, double &max) const;
+  void _minMaxJacobianSqr(const fullVector<double>&, double &min, double &max) const;
+  void _minJ2P3(const fullMatrix<double>&, const fullVector<double>&, double &min) const;
 };
 
 #endif

Numeric/bezierBasis.cpp

@@ -551,59 +551,58 @@ void bezierBasis::_construct(int parentType, int p)
     return;
   }
 
-  int dimSimplex;
   switch (parentType) {
     case TYPE_PNT :
       dim = 0;
       numLagCoeff = 1;
-      dimSimplex = 0;
+      _dimSimplex = 0;
       _exponents = gmshGenerateMonomialsLine(0);
       subPoints.push_back(gmshGeneratePointsLine(0));
       break;
     case TYPE_LIN : {
       dim = 1;
-      numLagCoeff = 2;
-      dimSimplex = 0;
+      numLagCoeff = order == 0 ? 1 : 2;
+      _dimSimplex = 0;
       _exponents = gmshGenerateMonomialsLine(order);
       subPoints = generateSubPointsLine(order);
       break;
     }
     case TYPE_TRI : {
       dim = 2;
-      numLagCoeff = 3;
-      dimSimplex = 2;
+      numLagCoeff = order == 0 ? 1 : 3;
+      _dimSimplex = 2;
       _exponents = gmshGenerateMonomialsTriangle(order);
       subPoints = generateSubPointsTriangle(order);
       break;
     }
     case TYPE_QUA : {
       dim = 2;
-      numLagCoeff = 4;
-      dimSimplex = 0;
+      numLagCoeff = order == 0 ? 1 : 4;
+      _dimSimplex = 0;
       _exponents = gmshGenerateMonomialsQuadrangle(order);
       subPoints = generateSubPointsQuad(order);
       break;
     }
     case TYPE_TET : {
       dim = 3;
-      numLagCoeff = 4;
-      dimSimplex = 3;
+      numLagCoeff = order == 0 ? 1 : 4;
+      _dimSimplex = 3;
       _exponents = gmshGenerateMonomialsTetrahedron(order);
       subPoints = generateSubPointsTetrahedron(order);
       break;
     }
     case TYPE_PRI : {
       dim = 3;
-      numLagCoeff = 6;
-      dimSimplex = 2;
+      numLagCoeff = order == 0 ? 1 : 6;
+      _dimSimplex = 2;
       _exponents = gmshGenerateMonomialsPrism(order);
       subPoints = generateSubPointsPrism(order);
       break;
     }
     case TYPE_HEX : {
       dim = 3;
-      numLagCoeff = 8;
-      dimSimplex = 0;
+      numLagCoeff = order == 0 ? 1 : 8;
+      _dimSimplex = 0;
       _exponents = gmshGenerateMonomialsHexahedron(order);
       subPoints = generateSubPointsHex(order);
       break;
@@ -614,7 +613,7 @@ void bezierBasis::_construct(int parentType, int p)
       dim = 3;
       order = 0;
       numLagCoeff = 4;
-      dimSimplex = 3;
+      _dimSimplex = 3;
       _exponents = gmshGenerateMonomialsTetrahedron(order);
       subPoints = generateSubPointsTetrahedron(order);
       break;
@@ -625,7 +624,7 @@ void bezierBasis::_construct(int parentType, int p)
   fullMatrix<double> bezierPoints = _exponents;
   bezierPoints.scale(1./order);
 
-  matrixBez2Lag = generateBez2LagMatrix(_exponents, bezierPoints, order, dimSimplex);
+  matrixBez2Lag = generateBez2LagMatrix(_exponents, bezierPoints, order, _dimSimplex);
   matrixBez2Lag.invert(matrixLag2Bez);
-  subDivisor = generateSubDivisor(_exponents, subPoints, matrixLag2Bez, order, dimSimplex);
+  subDivisor = generateSubDivisor(_exponents, subPoints, matrixLag2Bez, order, _dimSimplex);
 }

Numeric/bezierBasis.h

@@ -19,6 +19,8 @@ class bezierBasis {
   int numLagCoeff;
   int dim, type, order;
   int numDivisions;
+ public:
+  int _dimSimplex;
   fullMatrix<double> _exponents;
 
  public :