No commit message

No commit message

Numeric/BasisFactory.cpp

@@ -14,6 +14,7 @@
 
 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;
 
@@ -77,6 +78,17 @@ const JacobianBasis* BasisFactory::getJacobianBasis(int tag)
   return J;
 }
 
+const MetricBasis* BasisFactory::getMetricBasis(int tag)
+{
+  std::map<int, MetricBasis*>::const_iterator it = ms.find(tag);
+  if (it != ms.end())
+    return it->second;
+
+  MetricBasis* M = new MetricBasis(tag);
+  ms.insert(std::make_pair(tag, M));
+  return M;
+}
+
 const GradientBasis* BasisFactory::getGradientBasis(int tag, int order)
 {
   std::pair<int, int> key(tag, order);

Numeric/BasisFactory.h

@@ -10,6 +10,7 @@
 #include "MPyramid.h"
 #include "nodalBasis.h"
 #include "JacobianBasis.h"
+#include "MetricBasis.h"
 
 class BasisFactory
 {
@@ -19,6 +20,7 @@ class BasisFactory
  private:
   static std::map<int, nodalBasis*> fs;
   static std::map<int, JacobianBasis*> js;
+  static std::map<int, MetricBasis*> ms;
   static Cont_gradBasis gs;
   static Cont_bezierBasis bs;
   // store bezier bases by parentType and order (no serendipity..)
@@ -27,6 +29,8 @@ class BasisFactory
   // Caution: the returned pointer can be NULL
   static const nodalBasis* getNodalBasis(int tag);
   static const JacobianBasis* getJacobianBasis(int tag);
+  static const MetricBasis* getMetricBasis(int tag);
+
   static const GradientBasis* getGradientBasis(int tag, int order);
   static const bezierBasis* getBezierBasis(int parentTag, int order);
   static inline const bezierBasis* getBezierBasis(int tag) {

Numeric/JacobianBasis.cpp

@@ -659,12 +659,18 @@ void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
 
 void JacobianBasis::interpolate(const fullVector<double> &jacobian,
                                 const fullMatrix<double> &uvw,
-                                fullMatrix<double> &result) const
+                                fullMatrix<double> &result,
+                                bool areBezier) const
 {
   fullMatrix<double> bezM(jacobian.size(), 1);
   fullVector<double> bez;
   bez.setAsProxy(bezM, 0);
+
+  if (areBezier)
+    bez.setAll(jacobian);
+  else
     lag2Bez(jacobian, bez);
+
   bezier->interpolate(bezM, uvw, result);
 }
 

Numeric/JacobianBasis.h

@@ -113,7 +113,7 @@ class JacobianBasis {
   //
   void interpolate(const fullVector<double> &jacobian,
                    const fullMatrix<double> &uvw,
-                   fullMatrix<double> &result) const;
+                   fullMatrix<double> &result, bool areBezier = false) const;
 
   //
   static int jacobianOrder(int parentType, int order);

Numeric/MetricBasis.h

@@ -9,42 +9,138 @@
 #include "MElement.h"
 #include "JacobianBasis.h"
 #include "fullMatrix.h"
+#include <fstream>
 
-class MetricCoefficient {
+class MetricBasis {
+  friend class MetricCoefficient;
+  friend class GMSH_AnalyseCurvedMeshPlugin;
 private:
+  const JacobianBasis *_jacobian;
+  const GradientBasis *_gradients;
+  const bezierBasis *_bezier;
+  static double _tol;
+  static int _which;
+
+  int __maxdepth, __numSubdivision;
+  std::vector<int> __numSub;
+  MElement *__curElem;
+
+  std::fstream file;
+
+  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 MetricData {
    public:
     fullMatrix<double> *_metcoeffs;
     fullVector<double> *_jaccoeffs;
     double _RminBez;
-    int _depth;
+    int _depth, _num;
 
    public:
-    MetricData(fullMatrix<double> *m, fullVector<double> *j, double r, int d) :
-      _metcoeffs(m), _jaccoeffs(j), _RminBez(r), _depth(d) {}
+    MetricData(fullMatrix<double> *m, fullVector<double> *j, double r, int d, int num) :
+      _metcoeffs(m), _jaccoeffs(j), _RminBez(r), _depth(d), _num(num) {}
     ~MetricData() {
       delete _metcoeffs;
       delete _jaccoeffs;
     }
   };
+
+  std::map<int, std::vector<IneqData> > _ineqJ2, _ineqP3, _ineqA;
+
+public:
+  MetricBasis(int elementTag);
+
+  static void setTol(double tol) {_tol = tol;}
+  static double getTol() {return _tol;}
+  static void setWhich(int which) {_which = which;}
+
+  double getBoundRmin(const MElement*, MetricData*&, fullMatrix<double>&) const;
+  double getMinR(const MElement*, MetricData*&, int) const;
+  static double boundRmin(const MElement *el);
+  //double getBoundRmin(int, MElement**, double*);
+  //static double boundRmin(int, MElement**, double*, bool sameType = false);
+
+  void interpolate(const MElement*, const MetricData*, const double *uvw, double *minmaxQ, bool write = false) const;
+
+  static int metricOrder(int tag);
+
+private:
+  void _fillInequalities(int order);
+  void _lightenInequalities(int&, int&, int&); //TODO change
+
+  void _computeRmin1(const fullMatrix<double>&, const fullVector<double>&,
+                    double &RminLag, double &RminBez, int) const;
+  void _computeRmin2(const fullMatrix<double>&, const fullVector<double>&,
+                    double &RminLag, double &RminBez, int depth, int which) const;
+  void _computeRmin3(const fullMatrix<double>&, const fullVector<double>&,
+                    double &RminLag, double &RminBez, int depth, bool debug = false) const;
+
+  double _subdivideForRmin(MetricData*, double RminLag, double tol, int which) const;
+
+  double _minp(const fullMatrix<double>&) const;
+  double _minp2(const fullMatrix<double>&) const;
+  double _minq(const fullMatrix<double>&) const;
+  double _maxp(const fullMatrix<double>&) const;
+  double _maxq(const fullMatrix<double>&) const;
+  void _minMaxA2(const fullMatrix<double>&, double &min, double &max) const;
+  void _minJ2P3(const fullMatrix<double>&, const fullVector<double>&, double &min) const;
+  void _maxAstK(const fullMatrix<double>&, const fullVector<double>&,
+                double minK, double a1, double &maxa) const;                    //wrong
+  void _maxAstK2(const fullMatrix<double>&, const fullVector<double>&,
+                 double minK, double beta, double &maxa) const;                 //wrong
+  void _maxAstK3(const fullMatrix<double>&, const fullVector<double>&,
+                 double minK, double beta, double &maxa) const;                 //poor bound
+  void _maxAstK4(const fullMatrix<double>&, const fullVector<double>&,
+                 double minK, double beta, double &maxa) const;                 //better bound, WI positive ?
+  void _maxKstA(const fullMatrix<double>&, const fullVector<double>&,
+                 double mina, double &maxK) const;                              //wrong
+  void _maxKstA2(const fullMatrix<double>&, const fullVector<double>&,
+                 double mina, double beta, double &maxK) const;                 //faster
+  void _maxKstA3(const fullMatrix<double>&, const fullVector<double>&,
+                 double mina, double beta, double &maxK) const;                 //better bound
+  void _minMaxJacobianSqr(const fullVector<double>&, double &min, double &max) const;
+
+private:
+  class gterIneq {
+   public:
+    bool operator()(const IneqData &id1, const IneqData &id2) const {
+      return id1.val > id2.val;
+    }
+  };
   struct lessMinB {
     bool operator()(const MetricData *md1, const MetricData *md2) const {
       return md1->_RminBez > md2->_RminBez;
     }
   };
+};
 
-  class IneqData {
+class MetricCoefficient {
+private:
+  class MetricData {
    public:
-    int i, j, k;
-    double val;
+    fullMatrix<double> *_metcoeffs;
+    fullVector<double> *_jaccoeffs;
+    double _RminBez;
+    int _depth;
 
    public:
-    IneqData(double val, int i, int j, int k = -1) : i(i), j(j), k(k), val(val) {}
+    MetricData(fullMatrix<double> *m, fullVector<double> *j, double r, int d) :
+      _metcoeffs(m), _jaccoeffs(j), _RminBez(r), _depth(d) {}
+    ~MetricData() {
+      delete _metcoeffs;
+      delete _jaccoeffs;
+    }
   };
-  class gterIneq {
-   public:
-    bool operator()(const IneqData &id1, const IneqData &id2) const {
-      return id1.val > id2.val;
+  struct lessMinB {
+    bool operator()(const MetricData *md1, const MetricData *md2) const {
+      return md1->_RminBez > md2->_RminBez;
     }
   };
 
@@ -54,9 +150,9 @@ private:
   const GradientBasis *_gradients;
   const bezierBasis *_bezier;
   fullMatrix<double> _coefficientsLag, _coefficientsBez;
-  std::map<int, std::vector<IneqData> > _ineqJ2, _ineqP3, _ineqA;
   int __maxdepth, __numSubdivision;
   std::vector<int> __numSub;
+  MetricBasis *_basis;
 
  public:
   MetricCoefficient(MElement*);
@@ -65,25 +161,9 @@ private:
   void interpolate(const double *uvw, double *minmaxQ);
   double getBoundRmin(double tol, int which);
 
-  static int metricOrder(int tag);
-
  private:
   void _computeBezCoeff();
-  void _fillInequalities(int order);
-  void _lightenInequalities(double, int&, 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 _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