add Pyramids to bezierBasis & corr JacobianBasis::generateJacPointsPyramid(..)

+ add jacobianOrder(..)

Numeric/JacobianBasis.cpp

@@ -18,53 +18,37 @@ JacobianBasis::JacobianBasis(int tag)
 {
   const int parentType = ElementType::ParentTypeFromTag(tag);
   const int order = ElementType::OrderFromTag(tag);
+  const int jacobianOrder =  JacobianBasis::jacobianOrder(parentType, order);
+  const int primJacobianOrder =  JacobianBasis::jacobianOrder(parentType, 1);
 
-  int jacobianOrder, primJacobianOrder;
   switch (parentType) {
     case TYPE_PNT :
-      primJacobianOrder = 0;
-      jacobianOrder = 0;
       lagPoints = gmshGeneratePointsLine(0);
       break;
     case TYPE_LIN :
-      primJacobianOrder = 0;
-      jacobianOrder = order - 1;
       lagPoints = gmshGeneratePointsLine(jacobianOrder);
       break;
     case TYPE_TRI :
-      primJacobianOrder = 0;
-      jacobianOrder = 2*order - 2;
       lagPoints = gmshGeneratePointsTriangle(jacobianOrder);
       break;
     case TYPE_QUA :
-      primJacobianOrder = 1;
-      jacobianOrder = 2*order - 1;
       lagPoints = gmshGeneratePointsQuadrangle(jacobianOrder);
       break;
     case TYPE_TET :
-      primJacobianOrder = 0;
-      jacobianOrder = 3*order - 3;
       lagPoints = gmshGeneratePointsTetrahedron(jacobianOrder);
       break;
     case TYPE_PRI :
-      primJacobianOrder = 2;
-      jacobianOrder = 3*order - 1;
       lagPoints = gmshGeneratePointsPrism(jacobianOrder);
       break;
     case TYPE_HEX :
-      primJacobianOrder = 2;
-      jacobianOrder = 3*order - 1;
       lagPoints = gmshGeneratePointsHexahedron(jacobianOrder);
       break;
     case TYPE_PYR :
-      primJacobianOrder = 0;
-      jacobianOrder = 3*order - 3;
       lagPoints = generateJacPointsPyramid(jacobianOrder);
       break;
     default :
       Msg::Error("Unknown Jacobian function space for element type %d", tag);
-      jacobianOrder = 0;
-      break;
+      return;
   }
   numJacNodes = lagPoints.size1();
 
@@ -680,3 +664,36 @@ fullMatrix<double> JacobianBasis::generateJacMonomialsPyramid(int order)
   }
   return monomials;
 }
+
+fullMatrix<double> JacobianBasis::generateJacPointsPyramid(int order)
+{
+  fullMatrix<double> points = generateJacMonomialsPyramid(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;
+  }
+
+  return points;
+}
+
+int JacobianBasis::jacobianOrder(int parentType, int order) {
+  switch (parentType) {
+    case TYPE_PNT : return 0;
+    case TYPE_LIN : return order - 1;
+    case TYPE_TRI : return 2*order - 2;
+    case TYPE_QUA : return 2*order - 1;
+    case TYPE_TET : return 3*order - 3;
+    case TYPE_PRI : return 3*order - 1;
+    case TYPE_HEX : return 3*order - 1;
+    case TYPE_PYR : return 3*order - 3;
+    // note : for the pyramid, the jacobian space is
+    // different from the space of the mapping
+    default :
+      Msg::Error("Unknown element type %d, return order 0", parentType);
+      return 0;
+  }
+}

Numeric/JacobianBasis.h

@@ -67,15 +67,9 @@ class JacobianBasis {
   }
 
   //
+  static int jacobianOrder(int parentType, int order);
   static fullMatrix<double> generateJacMonomialsPyramid(int order);
-  static inline fullMatrix<double> generateJacPointsPyramid(int order) {
-    fullMatrix<double> points = generateJacMonomialsPyramid(order);
-
-    if (order == 0) return points;
-
-    points.scale(1. / (order+2.));
-    return points;
-  }
+  static fullMatrix<double> generateJacPointsPyramid(int order);
 };
 
 #endif

Numeric/bezierBasis.cpp

@@ -13,7 +13,6 @@
 #include "BasisFactory.h"
 #include "Numeric.h"
 
-
 // Sub Control Points
 static std::vector< fullMatrix<double> > generateSubPointsLine(int order)
 {
@@ -223,6 +222,81 @@ static std::vector< fullMatrix<double> > generateSubPointsHex(int order)
   return subPoints;
 }
 
+static std::vector< fullMatrix<double> > generateSubPointsPyr(int order)
+{
+  if (order == 0) {
+    std::vector< fullMatrix<double> > subPoints(4);
+    fullMatrix<double> prox;
+
+    subPoints[0] = JacobianBasis::generateJacMonomialsPyramid(order);
+    subPoints[0].scale(.5/(order+2));
+
+    subPoints[1].copy(subPoints[0]);
+    prox.setAsProxy(subPoints[1], 0, 1);
+    prox.add(.5);
+
+    subPoints[2].copy(subPoints[0]);
+    prox.setAsProxy(subPoints[2], 1, 1);
+    prox.add(.5);
+
+    subPoints[3].copy(subPoints[1]);
+    prox.setAsProxy(subPoints[3], 1, 1);
+    prox.add(.5);
+
+    return subPoints;
+  }
+
+  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[1].copy(subPoints[0]);
+  prox.setAsProxy(subPoints[1], 0, 1);
+  prox.add(.5);
+
+  subPoints[2].copy(subPoints[0]);
+  prox.setAsProxy(subPoints[2], 1, 1);
+  prox.add(.5);
+
+  subPoints[3].copy(subPoints[1]);
+  prox.setAsProxy(subPoints[3], 1, 1);
+  prox.add(.5);
+
+  subPoints[4].copy(subPoints[0]);
+  prox.setAsProxy(subPoints[4], 2, 1);
+  prox.add(.5);
+
+  subPoints[5].copy(subPoints[1]);
+  prox.setAsProxy(subPoints[5], 2, 1);
+  prox.add(.5);
+
+  subPoints[6].copy(subPoints[2]);
+  prox.setAsProxy(subPoints[6], 2, 1);
+  prox.add(.5);
+
+  subPoints[7].copy(subPoints[3]);
+  prox.setAsProxy(subPoints[7], 2, 1);
+  prox.add(.5);
+
+  for (int i = 0; i < 8; ++i) {
+    for (int j = 0; j < nPts; ++j) {
+      const double factor = (1. - subPoints[i](j, 2));
+      subPoints[i](j, 0) = subPoints[i](j, 0) * factor;
+      subPoints[i](j, 1) = subPoints[i](j, 1) * factor;
+    }
+  }
+
+  return subPoints;
+}
+
 // Matrices generation
 static int nChoosek(int n, int k)
 {
@@ -246,9 +320,8 @@ static fullMatrix<double> generateBez2LagMatrix
   (const fullMatrix<double> &exponent, const fullMatrix<double> &point,
    int order, int dimSimplex)
 {
-
   if(exponent.size1() != point.size1() || exponent.size2() != point.size2()){
-    Msg::Fatal("Wrong sizes for Bezier coefficients generation %d %d -- %d %d",
+    Msg::Fatal("Wrong sizes for bez2lag matrix generation %d %d -- %d %d",
       exponent.size1(),point.size1(),
       exponent.size2(),point.size2());
     return fullMatrix<double>(1, 1);
@@ -285,6 +358,42 @@ static fullMatrix<double> generateBez2LagMatrix
   return bez2Lag;
 }
 
+static fullMatrix<double> generateBez2LagMatrixPyramid
+  (const fullMatrix<double> &exponent, const fullMatrix<double> &point,
+   int order)
+{
+  if(exponent.size1() != point.size1() || exponent.size2() != point.size2() ||
+      exponent.size2() != 3){
+    Msg::Fatal("Wrong sizes for bez2lag matrix generation %d %d -- %d %d",
+      exponent.size1(),point.size1(),
+      exponent.size2(),point.size2());
+    return fullMatrix<double>(1, 1);
+  }
+
+  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(point(i, 2), exponent(j, 2)) //FIXME use pow_int
+          * pow(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))
+          * pow(point(i, 0) / denom, exponent(j, 0))
+          * pow(point(i, 1) / denom, exponent(j, 1))
+          * pow(1. - point(i, 0) / denom, p - exponent(j, 0))
+          * pow(1. - point(i, 1) / denom, p - exponent(j, 1));
+    }
+  }
+  return bez2Lag;
+}
+
 static fullMatrix<double> generateSubDivisor
   (const fullMatrix<double> &exponents, const std::vector< fullMatrix<double> > &subPoints,
    const fullMatrix<double> &lag2Bez, int order, int dimSimplex)
@@ -311,102 +420,130 @@ static fullMatrix<double> generateSubDivisor
   return subDivisor;
 }
 
+static fullMatrix<double> generateSubDivisorPyramid
+  (const fullMatrix<double> &exponents, const std::vector< fullMatrix<double> > &subPoints,
+   const fullMatrix<double> &lag2Bez, int order)
+{
+  if (exponents.size1() != lag2Bez.size1() || exponents.size1() != lag2Bez.size2()){
+    Msg::Fatal("Wrong sizes for Bezier Divisor %d %d -- %d %d",
+      exponents.size1(), lag2Bez.size1(),
+      exponents.size1(), lag2Bez.size2());
+    return fullMatrix<double>(1, 1);
+  }
+
+  int nbPts = lag2Bez.size1();
+  int nbSubPts = nbPts * subPoints.size();
+
+  fullMatrix<double> intermediate2(nbPts, nbPts);
+  fullMatrix<double> subDivisor(nbSubPts, nbPts);
+
+  for (unsigned int i = 0; i < subPoints.size(); i++) {
+    fullMatrix<double> intermediate1 =
+      generateBez2LagMatrixPyramid(exponents, subPoints[i], order);
+    lag2Bez.mult(intermediate1, intermediate2);
+    subDivisor.copy(intermediate2, 0, nbPts, 0, nbPts, i*nbPts, 0);
+  }
+  return subDivisor;
+}
+
 void bezierBasis::_construct(int parentType, int p)
 {
   order = p;
+  fullMatrix<double> exponents;
+  std::vector< fullMatrix<double> > subPoints;
 
-  /*if (parentType == TYPE_PYR) {
+  if (parentType == TYPE_PYR) {
     dim = 3;
-    numLagPts = 5;
-    bezierPoints = gmshGeneratePointsPyramid(order,false);
-    lagPoints = bezierPoints;
-    matrixLag2Bez.resize(bezierPoints.size1(), bezierPoints.size1(), false);
-    matrixBez2Lag.resize(bezierPoints.size1(), bezierPoints.size1(), false);
-    for (int i=0; i<matrixBez2Lag.size1(); ++i) {
-      matrixLag2Bez(i,i) = 1.;
-      matrixBez2Lag(i,i) = 1.;
-    }
-    // TODO: Implement subdidivsor
-  }
-  else {*/
-    int dimSimplex;
-    fullMatrix<double> exponents;
-    std::vector< fullMatrix<double> > subPoints;
-
-    switch (parentType) {
-      case TYPE_PNT :
-        dim = 0;
-        numLagCoeff = 1;
-        dimSimplex = 0;
-        exponents = gmshGenerateMonomialsLine(0);
-        break;
-      case TYPE_LIN : {
-        dim = 1;
-        numLagCoeff = 2;
-        dimSimplex = 0;
-        exponents = gmshGenerateMonomialsLine(order);
-        subPoints = generateSubPointsLine(order);
-        break;
-      }
-      case TYPE_TRI : {
-        dim = 2;
-        numLagCoeff = 3;
-        dimSimplex = 2;
-        exponents = gmshGenerateMonomialsTriangle(order);
-        subPoints = generateSubPointsTriangle(order);
-        break;
-      }
-      case TYPE_QUA : {
-        dim = 2;
-        numLagCoeff = 4;
-        dimSimplex = 0;
-        exponents = gmshGenerateMonomialsQuadrangle(order);
-        subPoints = generateSubPointsQuad(order);
-        break;
-      }
-      case TYPE_TET : {
-        dim = 3;
-        numLagCoeff = 4;
-        dimSimplex = 3;
-        exponents = gmshGenerateMonomialsTetrahedron(order);
-        subPoints = generateSubPointsTetrahedron(order);
-        break;
-      }
-      case TYPE_PRI : {
-        dim = 3;
-        numLagCoeff = 6;
-        dimSimplex = 2;
-        exponents = gmshGenerateMonomialsPrism(order);
-        subPoints = generateSubPointsPrism(order);
-        break;
-      }
-      case TYPE_HEX : {
-        dim = 3;
-        numLagCoeff = 8;
-        dimSimplex = 0;
-        exponents = gmshGenerateMonomialsHexahedron(order);
-        subPoints = generateSubPointsHex(order);
-        break;
-      }
-      default : {
-        Msg::Error("Unknown function space of parentType %d : "
-            "reverting to TET_1", parentType);
-        dim = 3;
-        order = 0;
-        numLagCoeff = 4;
-        dimSimplex = 3;
-        exponents = gmshGenerateMonomialsTetrahedron(order);
-        subPoints = generateSubPointsTetrahedron(order);
-        break;
-      }
-    }
+    numLagCoeff = 8;
+    exponents = JacobianBasis::generateJacMonomialsPyramid(order);
+    subPoints = generateSubPointsPyr(order);
+
+    numDivisions = static_cast<int>(subPoints.size());
 
     fullMatrix<double> bezierPoints = exponents;
-    bezierPoints.scale(1./order);
+    bezierPoints.scale(1. / (order + 2));
 
-    numDivisions = static_cast<int>(subPoints.size());
-    matrixBez2Lag = generateBez2LagMatrix(exponents, bezierPoints, order, dimSimplex);
+    matrixBez2Lag = generateBez2LagMatrixPyramid(exponents, bezierPoints, order);
     matrixBez2Lag.invert(matrixLag2Bez);
-    subDivisor = generateSubDivisor(exponents, subPoints, matrixLag2Bez, order, dimSimplex);
-  //}
+    subDivisor = generateSubDivisorPyramid(exponents, subPoints, matrixLag2Bez, order);
+    return;
+  }
+
+  int dimSimplex;
+  switch (parentType) {
+    case TYPE_PNT :
+      dim = 0;
+      numLagCoeff = 1;
+      dimSimplex = 0;
+      exponents = gmshGenerateMonomialsLine(0);
+      subPoints.push_back(gmshGeneratePointsLine(0));
+      break;
+    case TYPE_LIN : {
+      dim = 1;
+      numLagCoeff = 2;
+      dimSimplex = 0;
+      exponents = gmshGenerateMonomialsLine(order);
+      subPoints = generateSubPointsLine(order);
+      break;
+    }
+    case TYPE_TRI : {
+      dim = 2;
+      numLagCoeff = 3;
+      dimSimplex = 2;
+      exponents = gmshGenerateMonomialsTriangle(order);
+      subPoints = generateSubPointsTriangle(order);
+      break;
+    }
+    case TYPE_QUA : {
+      dim = 2;
+      numLagCoeff = 4;
+      dimSimplex = 0;
+      exponents = gmshGenerateMonomialsQuadrangle(order);
+      subPoints = generateSubPointsQuad(order);
+      break;
+    }
+    case TYPE_TET : {
+      dim = 3;
+      numLagCoeff = 4;
+      dimSimplex = 3;
+      exponents = gmshGenerateMonomialsTetrahedron(order);
+      subPoints = generateSubPointsTetrahedron(order);
+      break;
+    }
+    case TYPE_PRI : {
+      dim = 3;
+      numLagCoeff = 6;
+      dimSimplex = 2;
+      exponents = gmshGenerateMonomialsPrism(order);
+      subPoints = generateSubPointsPrism(order);
+      break;
+    }
+    case TYPE_HEX : {
+      dim = 3;
+      numLagCoeff = 8;
+      dimSimplex = 0;
+      exponents = gmshGenerateMonomialsHexahedron(order);
+      subPoints = generateSubPointsHex(order);
+      break;
+    }
+    default : {
+      Msg::Error("Unknown function space of parentType %d : "
+          "reverting to TET_1", parentType);
+      dim = 3;
+      order = 0;
+      numLagCoeff = 4;
+      dimSimplex = 3;
+      exponents = gmshGenerateMonomialsTetrahedron(order);
+      subPoints = generateSubPointsTetrahedron(order);
+      break;
+    }
+  }
+  numDivisions = static_cast<int>(subPoints.size());
+
+  fullMatrix<double> bezierPoints = exponents;
+  bezierPoints.scale(1./order);
+
+  matrixBez2Lag = generateBez2LagMatrix(exponents, bezierPoints, order, dimSimplex);
+  matrixBez2Lag.invert(matrixLag2Bez);
+  subDivisor = generateSubDivisor(exponents, subPoints, matrixLag2Bez, order, dimSimplex);
 }