Bezier: Implementation of subdivision for pyramids

Numeric/JacobianBasis.cpp

@@ -676,7 +676,7 @@ void JacobianBasis::interpolate(const fullVector<double> &jacobian,
 
 fullMatrix<double> JacobianBasis::generateJacMonomialsPyramid(int order)
 {
-  int nbMonomials = (order+3)*((order+3)+1)*(2*(order+3)+1)/6 - 5;
+  const int nbMonomials = (order+3)*(order+3)*(order+1);
   fullMatrix<double> monomials(nbMonomials, 3);
 
   if (order == 0) {
@@ -706,28 +706,28 @@ fullMatrix<double> JacobianBasis::generateJacMonomialsPyramid(int order)
   monomials(4, 1) = 0;
   monomials(4, 2) = order;
 
-  monomials(5, 0) = 2;
+  monomials(5, 0) = order+2;
   monomials(5, 1) = 0;
   monomials(5, 2) = order;
 
-  monomials(6, 0) = 2;
-  monomials(6, 1) = 2;
+  monomials(6, 0) = order+2;
+  monomials(6, 1) = order+2;
   monomials(6, 2) = order;
 
   monomials(7, 0) = 0;
-  monomials(7, 1) = 2;
+  monomials(7, 1) = order+2;
   monomials(7, 2) = order;
 
   int index = 8;
 
-  static const int bottom_edges[8][2] = {
+  static const int bottom_edges[4][2] = {
     {0, 1},
     {1, 2},
     {2, 3},
     {3, 0}
   };
 
-  // bottom "edges"
+  // bottom & top "edges"
   for (int iedge = 0; iedge < 4; ++iedge) {
     int i0 = bottom_edges[iedge][0];
     int i1 = bottom_edges[iedge][1];
@@ -740,22 +740,14 @@ fullMatrix<double> JacobianBasis::generateJacMonomialsPyramid(int order)
       monomials(index, 1) = monomials(i0, 1) + i * u_2;
       monomials(index, 2) = 0;
     }
+    for (int i = 1; i < order + 2; ++i, ++index) {
+      monomials(index, 0) = monomials(i0, 0) + i * u_1;
+      monomials(index, 1) = monomials(i0, 1) + i * u_2;
+      monomials(index, 2) = order;
     }
-
-  // top "edges"
-  for (int iedge = 0; iedge < 4; ++iedge, ++index) {
-    int i0 = bottom_edges[iedge][0] + 4;
-    int i1 = bottom_edges[iedge][1] + 4;
-
-    int u_1 = (monomials(i1,0)-monomials(i0,0)) / 2;
-    int u_2 = (monomials(i1,1)-monomials(i0,1)) / 2;
-
-    monomials(index, 0) = monomials(i0, 0) + u_1;
-    monomials(index, 1) = monomials(i0, 1) + u_2;
-    monomials(index, 2) = monomials(i0, 2);
   }
 
-  // bottom "face"
+  // bottom & top "face"
   fullMatrix<double> uv = gmshGenerateMonomialsQuadrangle(order);
   uv.add(1);
   for (int i = 0; i < uv.size1(); ++i, ++index) {
@@ -763,16 +755,15 @@ fullMatrix<double> JacobianBasis::generateJacMonomialsPyramid(int order)
     monomials(index, 1) = uv(i, 1);
     monomials(index, 2) = 0;
   }
-
-  // top "face"
-  monomials(index, 0) = 1;
-  monomials(index, 1) = 1;
+  for (int i = 0; i < uv.size1(); ++i, ++index) {
+    monomials(index, 0) = uv(i, 0);
+    monomials(index, 1) = uv(i, 1);
     monomials(index, 2) = order;
-  ++index;
+  }
 
   // other monomials
+  uv = gmshGenerateMonomialsQuadrangle(order + 2);
   for (int k = 1; k < order; ++k) {
-    fullMatrix<double> uv = gmshGenerateMonomialsQuadrangle(order + 2 - k);
     for (int i = 0; i < uv.size1(); ++i, ++index) {
       monomials(index, 0) = uv(i, 0);
       monomials(index, 1) = uv(i, 1);
@@ -788,10 +779,10 @@ fullMatrix<double> JacobianBasis::generateJacPointsPyramid(int order)
 
   const double p = order + 2;
   for (int i = 0; i < points.size1(); ++i) {
-    points(i, 2) = points(i, 2) * 1. / p;
-    const double duv = -1. + points(i, 2);
-    points(i, 0) = duv + points(i, 0) * 2. / p;
-    points(i, 1) = duv + points(i, 1) * 2. / p;
+    points(i, 2) = points(i, 2) / p;
+    const double scale = 1. - points(i, 2);
+    points(i, 0) = (-1. + 2. * points(i, 0) / p) * scale;
+    points(i, 1) = (-1. + 2. * points(i, 1) / p) * scale;
   }
 
   return points;

Numeric/MetricBasis.cpp

@@ -169,7 +169,7 @@ double MetricBasis::getMinR(MElement *el, MetricData *&md, int deg) const
       else*/
       _computeRmin(*coeff, *jac, RminLag, RminBez, 0, false);
 
-      double betaOpt = beta, minaOpt, maxaOpt = 0., RminBezOpt;
+      double betaOpt = beta, minaOpt = mina, maxaOpt = maxa3, RminBezOpt;
       {
         /*const */double phi = std::acos(.5*(minK-maxa3*maxa3*maxa3+3*maxa3))/3;
         RminBezOpt = (maxa3+2*std::cos(phi+2*M_PI/3))/(maxa3+2*std::cos(phi));

Numeric/bezierBasis.cpp

@@ -229,7 +229,7 @@ namespace {
     std::vector< fullMatrix<double> > subPoints(4);
     fullMatrix<double> prox;
 
-      subPoints[0] = JacobianBasis::generateJacMonomialsPyramid(order);
+    subPoints[0] = JacobianBasis::generateJacMonomialsPyramid(0);
     subPoints[0].scale(.5/(order+2));
 
     subPoints[1].copy(subPoints[0]);
@@ -246,22 +246,12 @@ namespace {
 
     return subPoints;
   }
-
-    Msg::Error("Subpoints for pyramid p>1 not implemented !");
-    return std::vector< fullMatrix<double> >();
-
-    /*
+  else {
     std::vector< fullMatrix<double> > subPoints(8);
     fullMatrix<double> ref, prox;
 
     subPoints[0] = JacobianBasis::generateJacMonomialsPyramid(order);
-    int nPts = subPoints[0].size1();
-    for (int i = 0; i < nPts; ++i) {
-      const double factor = .5 / (order + 2 - subPoints[0](i, 2));
-      subPoints[0](i, 0) = subPoints[0](i, 0) * factor;
-      subPoints[0](i, 1) = subPoints[0](i, 1) * factor;
-      subPoints[0](i, 2) = subPoints[0](i, 2) * .5 / (order + 2);
-    }
+    subPoints[0].scale(.5/(order+2));
 
     subPoints[1].copy(subPoints[0]);
     prox.setAsProxy(subPoints[1], 0, 1);
@@ -291,6 +281,7 @@ namespace {
     prox.setAsProxy(subPoints[7], 2, 1);
     prox.add(.5);
 
+    const int nPts = subPoints[0].size1();
     for (int i = 0; i < 8; ++i) {
       for (int j = 0; j < nPts; ++j) {
         const double factor = (1. - subPoints[i](j, 2));
@@ -300,14 +291,14 @@ namespace {
     }
 
     return subPoints;
-    */
+  }
 }
 
 // Matrices generation
 int nChoosek(int n, int k)
 {
   if (n < k || k < 0) {
-      Msg::Error("Wrong argument for combination.");
+    Msg::Error("Wrong argument for combination. (%d, %d)", n, k);
     return 1;
   }
 
@@ -376,25 +367,22 @@ namespace {
     return fullMatrix<double>(1, 1);
   }
 
-    int ndofs = exponent.size1();
+  const int ndofs = exponent.size1();
 
   fullMatrix<double> bez2Lag(ndofs, ndofs);
   for (int i = 0; i < ndofs; i++) {
     for (int j = 0; j < ndofs; j++) {
-        bez2Lag(i, j) =
-            nChoosek(order, exponent(j, 2))
-            * pow_int(point(i, 2), exponent(j, 2))
-            * pow_int(1. - point(i, 2), order - exponent(j, 2));
-  
-        double p = order + 2 - exponent(j, 2);
       double denom = 1. - point(i, 2);
-        bez2Lag(i, j) *=
-              nChoosek(p, exponent(j, 0))
-            * nChoosek(p, exponent(j, 1))
+      bez2Lag(i, j) =
+            nChoosek(order + 2, exponent(j, 0))
+          * nChoosek(order + 2, exponent(j, 1))
+          * nChoosek(order, exponent(j, 2))
           * pow_int(point(i, 0) / denom, exponent(j, 0))
           * pow_int(point(i, 1) / denom, exponent(j, 1))
-            * pow_int(1. - point(i, 0) / denom, p - exponent(j, 0))
-            * pow_int(1. - point(i, 1) / denom, p - exponent(j, 1));
+          * pow_int(point(i, 2)        , exponent(j, 2))
+          * pow_int(1. - point(i, 0) / denom, order + 2 - exponent(j, 0))
+          * pow_int(1. - point(i, 1) / denom, order + 2 - exponent(j, 1))
+          * pow_int(1. - point(i, 2)        , order     - exponent(j, 2));
     }
   }
   return bez2Lag;
@@ -531,22 +519,22 @@ void bezierBasis::_construct(int parentType, int p)
     dim = 3;
     numLagCoeff = 8;
     _exponents = JacobianBasis::generateJacMonomialsPyramid(order);
-    if (order < 2) {
+
     subPoints = generateSubPointsPyr(order);
     numDivisions = static_cast<int>(subPoints.size());
-    }
-    else {
-      numDivisions = 8;
-      static int a = 0;
-      if (++a == 1) Msg::Warning("Subdivision not available for pyramid p>1");
-    }
 
-    fullMatrix<double> bezierPoints = _exponents;
-    bezierPoints.scale(1. / (order + 2));
+    fullMatrix<double> bezierPoints;
+    bezierPoints.resize(_exponents.size1(), _exponents.size2());
+    const double p = order + 2;
+    for (int i = 0; i < bezierPoints.size1(); ++i) {
+      bezierPoints(i, 2) = _exponents(i, 2) / p;
+      const double scale = 1. - bezierPoints(i, 2);
+      bezierPoints(i, 0) = _exponents(i, 0) / p * scale;
+      bezierPoints(i, 1) = _exponents(i, 1) / p * scale;
+    }
 
     matrixBez2Lag = generateBez2LagMatrixPyramid(_exponents, bezierPoints, order);
     matrixBez2Lag.invert(matrixLag2Bez);
-    if (order < 2)
     subDivisor = generateSubDivisorPyramid(_exponents, subPoints, matrixLag2Bez, order);
     return;
   }

Plugin/AnalyseCurvedMesh.cpp

@@ -110,7 +110,9 @@ std::string GMSH_AnalyseCurvedMeshPlugin::getHelp() const
     "\n"
     "- Recompute = {0,1}: If the mesh has changed, set to 1 to recompute the bounds.\n"
     "\n"
-    "- Tolerance = ]0, 1[: Tolerance on the computation of min(R) or min(J)";
+    "- Tolerance = ]0, 1[: Tolerance on the computation of min(R) or min(J). "
+    "It should be at most 0.01 but it can be set to 1 to just check the validity of "
+    "the mesh.";
 }
 
 // Execution