remove useless things in GradientBasisx

Numeric/JacobianBasis.cpp

@@ -25,7 +25,6 @@ namespace {
 
 GradientBasis::GradientBasis(int tag, int order)
 {
-  _bezier = BasisFactory::getBezierBasis(ElementType::ParentTypeFromTag(tag), order);
   const int type = ElementType::ParentTypeFromTag(tag);
 
   fullMatrix<double> samplingPoints;
@@ -57,7 +56,7 @@ GradientBasis::GradientBasis(int tag, int order)
       //samplingPoints = JacobianBasis::generateJacPointsPyramid(order);
       break;
     default :
-      Msg::Error("Unknown Jacobian function space for element type %d", type);
+      Msg::Error("Unknown Jacobian function space for element tag %d", tag);
       return;
   }
   const int numSampPnts = samplingPoints.size1();
@@ -125,7 +124,7 @@ JacobianBasis::JacobianBasis(int tag)
       lagPoints = generateJacPointsPyramid(jacobianOrder);
       break;
     default :
-      Msg::Error("Unknown Jacobian function space for element type %d", parentType);
+      Msg::Error("Unknown Jacobian function space for element tag %d", tag);
       return;
   }
   numJacNodes = lagPoints.size1();
@@ -134,7 +133,7 @@ JacobianBasis::JacobianBasis(int tag)
   bezier = BasisFactory::getBezierBasis(parentType, jacobianOrder);
 
   // Store shape function gradients of mapping at Jacobian nodes
-  _gradBasis = new GradientBasis(tag, jacobianOrder);
+  _gradBasis = BasisFactory::getGradientBasis(tag, jacobianOrder);
 
   // Compute matrix for lifting from primary Jacobian basis to Jacobian basis
   int primJacType = ElementType::getTag(parentType, primJacobianOrder, false);
@@ -626,8 +625,8 @@ void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
   for (int l = 0; l < numJacNodes; l++) {
     double jac[2][2] = {{0., 0.}, {0., 0.}};
     for (int i = 0; i < numMapNodes; i++) {
-      const double &dPhidX = _gradBasis->getgradShapeMatX(l,i);
-      const double &dPhidY = _gradBasis->getgradShapeMatY(l,i);
+      const double &dPhidX = _gradBasis->gradShapeMatX(l,i);
+      const double &dPhidY = _gradBasis->gradShapeMatY(l,i);
       const double dpsidx = dPhidX * invJaci[0][0] + dPhidY * invJaci[1][0];
       const double dpsidy = dPhidX * invJaci[0][1] + dPhidY * invJaci[1][1];
       jac[0][0] += nodesXYZ(i,0) * dpsidx;

Numeric/JacobianBasis.h

@@ -16,30 +16,17 @@ class GradientBasis {
  public:
   fullMatrix<double> gradShapeMatX, gradShapeMatY, gradShapeMatZ;
 
- private:
-  const bezierBasis *_bezier;
-
  public :
 
   GradientBasis(int tag, int order);
 
   int getNumSamplingPoints() const {return gradShapeMatX.size1();}
   int getNumMapNodes() const {return gradShapeMatX.size2();}
-  const bezierBasis* getBezierBasis() const {return _bezier;}
 
   void getGradientsFromNodes(const fullMatrix<double> &nodes,
                              fullMatrix<double> *dxyzdX,
                              fullMatrix<double> *dxyzdY,
                              fullMatrix<double> *dxyzdZ) const;
-  inline double getgradShapeMatX(int i, int j) const {return gradShapeMatX(i, j);}
-  inline double getgradShapeMatY(int i, int j) const {return gradShapeMatY(i, j);}
-  inline double getgradShapeMatZ(int i, int j) const {return gradShapeMatZ(i, j);}
-  inline void lag2Bez(const fullVector<double> &jac, fullVector<double> &bez) const {
-    _bezier->matrixLag2Bez.mult(jac,bez);
-  }
-  inline void lag2Bez(const fullMatrix<double> &jac, fullMatrix<double> &bez) const {
-    _bezier->matrixLag2Bez.mult(jac,bez);
-  }
 };
 
 

Numeric/MetricBasis.cpp

@@ -12,9 +12,11 @@
 MetricCoefficient::MetricCoefficient(MElement *el) : _element(el)
 {
   const int tag = el->getTypeForMSH();
+  const int type = ElementType::ParentTypeFromTag(tag);
   const int metricOrder = MetricCoefficient::metricOrder(tag);
   _jacobian = BasisFactory::getJacobianBasis(tag);
   _gradients = BasisFactory::getGradientBasis(tag, metricOrder);
+  _bezier = BasisFactory::getBezierBasis(type, metricOrder);
 
   int nSampPnts = _gradients->getNumSamplingPoints();
   int nMapping = _gradients->getNumMapNodes();
@@ -143,7 +145,7 @@ void MetricCoefficient::interpolate(const double *uvw, double *minmaxQ)
     return;
   }
 
-  int order = _gradients->getBezierBasis()->getOrder();
+  int order = _bezier->getOrder();
 
   int dimSimplex;
   fullMatrix<double> exponents;
@@ -277,6 +279,69 @@ void MetricCoefficient::interpolate(const double *uvw, double *minmaxQ)
   }
 }
 
+double MetricCoefficient::getBoundRmin(double tol)
+{
+  fullMatrix<double> nodes(_element->getNumShapeFunctions(),3);
+  fullVector<double> jac(_jacobian->getNumJacNodes());
+  _element->getNodesCoord(nodes);
+  _jacobian->getSignedJacobian(nodes, jac);
+
+  if(jac.size() != _coefficientsBez.size1())
+    Msg::Fatal("not same size, jac %d, coeff %d", jac.size(), _coefficientsBez.size1());
+
+  // compute RminLag
+  double RminLag = 1.;
+  const double factor1 = std::sqrt(6) / 3;
+  const double factor2 = std::sqrt(6) * 3;
+
+  for (int i = 0; i < _bezier->getNumLagCoeff(); ++i) {
+    double q = _coefficientsBez(i, 0);
+    double tmp = pow_int(_coefficientsBez(i, 1), 2);
+    tmp += pow_int(_coefficientsBez(i, 2), 2);
+    tmp += pow_int(_coefficientsBez(i, 3), 2);
+    tmp += 2 * pow_int(_coefficientsBez(i, 4), 2);
+    tmp += 2 * pow_int(_coefficientsBez(i, 5), 2);
+    tmp += 2 * pow_int(_coefficientsBez(i, 6), 2);
+    tmp = tmp*factor1;
+    if (tmp < 1e-3 * q) {
+      tmp = (q - tmp) / (q + tmp);
+      RminLag = std::min(RminLag, tmp);
+    }
+    else {
+      double phi = jac(i)*jac(i);
+      phi -= q * q * q;
+      phi += .5 * q * tmp * tmp;
+      phi /= tmp * tmp * tmp;
+      phi *= factor2;
+      if (phi >  1.1 || phi < -1.1) Msg::Warning("phi %g", phi);
+      if (phi >  1) phi =  1;
+      if (phi < -1) phi = -1;
+      phi = std::acos(phi)/3;
+      tmp = (q + tmp * std::cos(phi + 2*M_PI/3)) / (q + tmp * std::cos(phi));
+      RminLag = std::min(RminLag, tmp);
+    }
+  }
+
+  // compute RminBez
+  /*double RminBez = 1.;
+  double minq = minq();
+  double minp = minp();
+  double maxq = maxq();
+  double maxp = maxp();
+  if (minq < 1e-3 * maxp) {
+    double tmp = (minq - maxp) / (minq + maxp);
+    RminLag = std::min(RminBez, tmp);
+  }
+  else {
+
+  }*/
+
+
+  /*compute RminBez
+  if RminLag-RminBez < tol : return RminBez
+  return subdivide*/
+
+}
 
 int MetricCoefficient::metricOrder(int tag)
 {
@@ -355,5 +420,5 @@ void MetricCoefficient::_computeBezCoeff()
 {
   _coefficientsBez.resize(_coefficientsLag.size1(),
                           _coefficientsLag.size2());
-  _gradients->lag2Bez(_coefficientsLag, _coefficientsBez);
+  _bezier->matrixLag2Bez.mult(_coefficientsLag, _coefficientsBez);
 }

Numeric/MetricBasis.h

@@ -15,6 +15,7 @@ class MetricCoefficient {
   MElement *_element;
   const JacobianBasis *_jacobian;
   const GradientBasis *_gradients;
+  const bezierBasis *_bezier;
   fullMatrix<double> _coefficientsLag, _coefficientsBez;
 
  public:
@@ -22,7 +23,7 @@ class MetricCoefficient {
 
   void getCoefficients(fullMatrix<double>&, bool bezier = true);
   void interpolate(const double *uvw, double *minmaxQ);
-  void getBoundRmin(double tol);
+  double getBoundRmin(double tol);
 
   static int metricOrder(int tag);