add faster and more accurate function MetricBasis::minSampledRnew(..)

Numeric/BasisFactory.h

@@ -61,6 +61,9 @@ class BasisFactory
   static const GradientBasis* getGradientBasis(int tag, int order) {
     return getGradientBasis(FuncSpaceData(true, tag, order));
   }
+  static const GradientBasis* getGradientBasis(int tag) {
+    return getGradientBasis(FuncSpaceData(tag));
+  }
 
   // Bezier
   static const bezierBasis* getBezierBasis(FuncSpaceData);

Numeric/FuncSpaceData.cpp

@@ -6,19 +6,20 @@
 #include "FuncSpaceData.h"
 #include "MElement.h"
 
-FuncSpaceData::FuncSpaceData(MElement *el, int *serendip) :
+FuncSpaceData::FuncSpaceData(const MElement *el, const bool *serendip) :
   _tag(el->getTypeForMSH()), _spaceOrder(el->getPolynomialOrder()),
   _nij(0), _nk(_spaceOrder), _pyramidalSpace(el->getType() == TYPE_PYR),
   _serendipity(serendip ? *serendip : el->getIsOnlySerendipity())
 {}
 
-FuncSpaceData::FuncSpaceData(MElement *el, int order, int *serendip) :
+FuncSpaceData::FuncSpaceData(const MElement *el, int order, const bool *serendip) :
   _tag(el->getTypeForMSH()), _spaceOrder(order),
   _nij(0), _nk(_spaceOrder), _pyramidalSpace(el->getType() == TYPE_PYR),
   _serendipity(serendip ? *serendip : el->getIsOnlySerendipity())
 {}
 
-FuncSpaceData::FuncSpaceData(MElement *el, bool pyr, int nij, int nk, int *serendip) :
+FuncSpaceData::FuncSpaceData(const MElement *el, bool pyr, int nij, int nk,
+                             const bool *serendip) :
   _tag(el->getTypeForMSH()), _spaceOrder(pyr ? nij+nk : std::max(nij, nk)),
   _nij(nij), _nk(nk), _pyramidalSpace(pyr),
   _serendipity(serendip ? *serendip : el->getIsOnlySerendipity())
@@ -27,7 +28,7 @@ FuncSpaceData::FuncSpaceData(MElement *el, bool pyr, int nij, int nk, int *seren
     Msg::Error("Creation of pyramidal space data for a non-pyramid element !");
 }
 
-FuncSpaceData::FuncSpaceData(int tag, int *serendip) :
+FuncSpaceData::FuncSpaceData(int tag, const bool *serendip) :
   _tag(tag), _spaceOrder(ElementType::OrderFromTag(tag)),
   _nij(0), _nk(_spaceOrder),
   _pyramidalSpace(ElementType::ParentTypeFromTag(tag) == TYPE_PYR),

Numeric/FuncSpaceData.h

@@ -31,12 +31,14 @@ public:
     _pyramidalSpace(false), _serendipity(false) {}
 
   // Constructors using MElement*
-  FuncSpaceData(MElement *el, int *serendip = NULL);
-  FuncSpaceData(MElement *el, int order, int *serendip = NULL);
-  FuncSpaceData(MElement *el, bool pyr, int nij, int nk, int *serendip = NULL);
+  FuncSpaceData(const MElement *el, const bool *serendip = NULL);
+  FuncSpaceData(const MElement *el, int order, const bool *serendip = NULL);
+  FuncSpaceData(const MElement *el,
+                bool pyr, int nij, int nk,
+                const bool *serendip = NULL);
 
   // Constructors using element tag
-  FuncSpaceData(int tag, int *serendip = NULL);
+  FuncSpaceData(int tag, const bool *serendip = NULL);
 
   // constructors using element tag or element type
   FuncSpaceData(bool isTag, int tagOrType, int order,

Numeric/JacobianBasis.cpp

@@ -854,3 +854,35 @@ int JacobianBasis::jacobianOrder(int parentType, int order)
       return 0;
   }
 }
+
+FuncSpaceData JacobianBasis::jacobianMatrixSpace(int type, int order)
+{
+  if (type == TYPE_PYR) {
+    Msg::Fatal("jacobianMatrixSpace not yet implemented for pyramids");
+    return FuncSpaceData(false, type, false, 1, 0);
+  }
+  int jacOrder = -1;
+  switch (type) {
+    case TYPE_PNT :
+      jacOrder = 0;
+      break;
+
+    case TYPE_LIN :
+    case TYPE_TRI :
+    case TYPE_TET :
+      jacOrder = order - 1;
+      break;
+
+    case TYPE_QUA :
+    case TYPE_PRI :
+    case TYPE_HEX :
+      jacOrder = order;
+      break;
+
+    default :
+      Msg::Error("Unknown element type %d, return order 0", type);
+      return 0;
+  }
+
+  return FuncSpaceData(true, ElementType::getTag(type, order), jacOrder);
+}

Numeric/JacobianBasis.h

@@ -202,6 +202,7 @@ public:
   //
   static int jacobianOrder(int tag);
   static int jacobianOrder(int parentType, int order);
+  static FuncSpaceData jacobianMatrixSpace(int type, int order);
 
 
  private :

Numeric/MetricBasis.cpp

@@ -153,36 +153,8 @@ double MetricBasis::minSampledR(MElement *el, int order)
 double MetricBasis::getMinSampledR(MElement *el, int deg) const
 {
   fullMatrix<double> samplingPoints;
-
-  switch (el->getType()) {
-    case TYPE_PNT :
-      samplingPoints = gmshGeneratePointsLine(0);
-      break;
-    case TYPE_LIN :
-      samplingPoints = gmshGeneratePointsLine(deg);
-      break;
-    case TYPE_TRI :
-      samplingPoints = gmshGeneratePointsTriangle(deg,false);
-      break;
-    case TYPE_QUA :
-      samplingPoints = gmshGeneratePointsQuadrangle(deg,false);
-      break;
-    case TYPE_TET :
-      samplingPoints = gmshGeneratePointsTetrahedron(deg,false);
-      break;
-    case TYPE_PRI :
-      samplingPoints = gmshGeneratePointsPrism(deg,false);
-      break;
-    case TYPE_HEX :
-      samplingPoints = gmshGeneratePointsHexahedron(deg,false);
-      break;
-    case TYPE_PYR :
-      samplingPoints = gmshGeneratePointsPyramidGeneral(true, deg+1, deg);
-      break;
-    default :
-      Msg::Error("Unknown Jacobian function space for element type %d", el->getType());
-      return -1;
-  }
+  bool serendip = false;
+  gmshGeneratePoints(FuncSpaceData(el, deg, &serendip), samplingPoints);
 
   MetricData *md;
   _getMetricData(el, md);
@@ -199,6 +171,73 @@ double MetricBasis::getMinSampledR(MElement *el, int deg) const
   return min;
 }
 
+double MetricBasis::minSampledRnew(MElement *el, int deg)
+{
+  int dim = el->getDim();
+  if (dim < 2) return 1;
+
+  fullMatrix<double> nodes(el->getNumVertices(), 3);
+  el->getNodesCoord(nodes);
+
+  // if deg is small: sample at every point
+  // otherwise: compute Bezier coefficients and interpolate
+  //const bool BezierInterp = deg > 100;
+
+  FuncSpaceData jacMatSpace;
+  /*if (BezierInterp) {
+    jacMatSpace =
+        JacobianBasis::jacobianMatrixSpace(el->getType(), el->getPolynomialOrder());
+  }
+  else {*/
+    const bool serendip = false;
+    jacMatSpace = FuncSpaceData(el, deg, &serendip);
+  //}
+
+  const GradientBasis *gb = BasisFactory::getGradientBasis(jacMatSpace);
+  fullMatrix<double> coeffJacMat(gb->getNumSamplingPoints(), 3*dim);
+  fullMatrix<double> dxyzdX(coeffJacMat, 0, 3);
+  fullMatrix<double> dxyzdY(coeffJacMat, 3, 3);
+  if (dim == 2)
+    gb->getIdealGradientsFromNodes(nodes, &dxyzdX, &dxyzdY, NULL);
+  else {
+    fullMatrix<double> dxyzdZ(coeffJacMat, 6, 3);
+    gb->getIdealGradientsFromNodes(nodes, &dxyzdX, &dxyzdY, &dxyzdZ);
+  }
+
+  /*if (BezierInterp) {
+    const bezierBasis *bb = BasisFactory::getBezierBasis(jacMatSpace);
+    fullMatrix<double> points;
+    const bool serendip = false;
+    gmshGeneratePoints(FuncSpaceData(el, deg, &serendip), points);
+
+    fullMatrix<double> bezCoeffJacMat;
+    bb->lag2Bez(coeffJacMat, bezCoeffJacMat);
+    bb->interpolate(bezCoeffJacMat, points, coeffJacMat);
+  }*/
+
+  double minR = 1;
+  for (int k = 0; k < coeffJacMat.size1(); ++k) {
+    fullMatrix<double> metric(dim, dim);
+    for (int i = 0; i < dim; ++i) {
+      for (int j = i; j < dim; ++j) {
+        for (int n = 0; n < 3; ++n)
+          metric(i, j) += coeffJacMat(k, 3*i+n) * coeffJacMat(k, 3*j+n);
+      }
+    }
+    metric(1, 0) = metric(0, 1);
+    if (dim == 3) {
+      metric(2, 0) = metric(0, 2);
+      metric(2, 1) = metric(1, 2);
+    }
+    fullVector<double> valReal(dim), valImag(dim);
+    fullMatrix<double> vecLeft(dim, dim), vecRight(dim, dim);
+    metric.eig(valReal, valImag, vecLeft, vecRight, true);
+    minR = std::min(minR, std::sqrt(valReal(0) / valReal(dim-1)));
+  }
+
+  return minR;
+}
+
 double MetricBasis::getBoundMinR(MElement *el) const
 {
   int nSampPnts = _gradients->getNumSamplingPoints();
@@ -607,6 +646,30 @@ bool MetricBasis::validateBezierForMetricAndJacobian()
       metricBasis->interpolateAfterNSubdivisions(el, numSubdiv, numSampPnt,
                                                  isub, uvw, metric_Bez);
 
+      /*{
+        static int deg = 2;
+        MetricBasis::minSampledR(el, deg);
+        MetricBasis::minSampledRnew(el, deg);
+
+        double time = Cpu();
+        double minR1;
+        for (int cn = 0; cn < 100; ++cn)
+          minR1 = MetricBasis::minSampledR(el, deg);
+        double tm1 = Cpu()-time;
+
+        time = Cpu();
+        double minR2;
+        for (int cn = 0; cn < 100; ++cn)
+          minR2 = MetricBasis::minSampledRnew(el, deg);
+        double tm2 = Cpu()-time;
+
+        if (std::abs(minR1-minR2) < 1e-14)
+          Msg::Info("ok deg %d, times %g %g speedup %g", deg, tm1, tm2, tm1/tm2);
+        else
+          Msg::Error("not ok deg %d: %g vs %g, times %g %g speedup %g", deg, minR1, minR2, tm1, tm2, tm1/tm2);
+        //++deg;
+      }*/
+
       int numBadMatch = 0;
       int numBadMatchTensor = 0;
       double maxBadMatch = 0;

Numeric/MetricBasis.h

@@ -64,6 +64,7 @@ public:
 
   static double boundMinR(MElement *el);
   static double minSampledR(MElement *el, int order);
+  static double minSampledRnew(MElement *el, int order);
   double getBoundMinR(MElement*) const;
   double getMinSampledR(MElement*, int order) const;
 

Numeric/bezierBasis.cpp

@@ -499,7 +499,7 @@ void bezierBasis::interpolate(const fullMatrix<double> &coeffs,
                               fullMatrix<double> &result,
                               bool bezCoord) const
 {
-  if (result.size1() < uvw.size1() || result.size2() < coeffs.size2())
+  if (result.size1() != uvw.size1() || result.size2() != coeffs.size2())
     result.resize(uvw.size1(), coeffs.size2());
 
   fullMatrix<double> bezuvw = uvw;
@@ -535,11 +535,24 @@ void bezierBasis::interpolate(const fullMatrix<double> &coeffs,
   }
 }
 
+void bezierBasis::lag2Bez(const fullMatrix<double> &lag,
+                          fullMatrix<double> &bez) const
+{
+  if (lag.size1() != matrixLag2Bez.size1()) {
+    Msg::Error("matrix not the right size in lag2Bez function %d vs %d",
+        lag.size1(), matrixLag2Bez.size1());
+  }
+  if (bez.size1() != lag.size1() || bez.size2() != lag.size2()) {
+    bez.resize(lag.size1(), lag.size2());
+  }
+  matrixLag2Bez.mult(lag, bez);
+}
+
 void bezierBasis::subdivideBezCoeff(const fullMatrix<double> &coeff,
                                     fullMatrix<double> &subCoeff) const
 {
-  if (subCoeff.size1() < subDivisor.size1()
-      || subCoeff.size2() < coeff.size2()  ) {
+  if (subCoeff.size1() != subDivisor.size1()
+      || subCoeff.size2() != coeff.size2()  ) {
     subCoeff.resize(subDivisor.size1(), coeff.size2());
   }
   subDivisor.mult(coeff, subCoeff);
@@ -548,7 +561,7 @@ void bezierBasis::subdivideBezCoeff(const fullMatrix<double> &coeff,
 void bezierBasis::subdivideBezCoeff(const fullVector<double> &coeff,
                                     fullVector<double> &subCoeff) const
 {
-  if (subCoeff.size() < subDivisor.size1()) {
+  if (subCoeff.size() != subDivisor.size1()) {
     subCoeff.resize(subDivisor.size1());
   }
   subDivisor.mult(coeff, subCoeff);

Numeric/bezierBasis.h

@@ -48,6 +48,9 @@ class bezierBasis {
   // generate Bezier points
   void generateBezierPoints(fullMatrix<double>&) const;
 
+  // transform coeff Lagrange into Bezier coeff
+  void lag2Bez(const fullMatrix<double> &lag, fullMatrix<double> &bez) const;
+
   // Subdivide Bezier coefficients
   void subdivideBezCoeff(const fullMatrix<double> &coeff,
                          fullMatrix<double> &subCoeff) const;