No commit message

No commit message

Geo/MPyramid.h

@@ -165,7 +165,6 @@ class MPyramid : public MElement {
   }
   static int faces_pyramid(const int face, const int vert)
   {
-    // only triangular faces
     static const int f[5][4] = {
       {0, 1, 4, -1},
       {3, 0, 4, -1},

Numeric/JacobianBasis.cpp

@@ -27,8 +27,8 @@ 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);
+  const int jacobianOrder = JacobianBasis::jacobianOrder(parentType, order);
+  const int primJacobianOrder = JacobianBasis::jacobianOrder(parentType, 1);
 
   fullMatrix<double> lagPoints;                                  // Sampling points
 
@@ -728,7 +728,8 @@ fullMatrix<double> JacobianBasis::generateJacPointsPyramid(int order)
   return points;
 }
 
-int JacobianBasis::jacobianOrder(int parentType, int order) {
+int JacobianBasis::jacobianOrder(int parentType, int order)
+{
   switch (parentType) {
     case TYPE_PNT : return 0;
     case TYPE_LIN : return order - 1;
@@ -745,3 +746,14 @@ int JacobianBasis::jacobianOrder(int parentType, int order) {
       return 0;
   }
 }
+
+
+void JacobianBasis::getGradientsFromNodes(const fullMatrix<double> &nodes,
+                                          fullMatrix<double> *dxyzdX,
+                                          fullMatrix<double> *dxyzdY,
+                                          fullMatrix<double> *dxyzdZ) const
+{
+  if (dxyzdX) gradShapeMatX.mult(nodes, *dxyzdX);
+  if (dxyzdY) gradShapeMatY.mult(nodes, *dxyzdY);
+  if (dxyzdZ) gradShapeMatZ.mult(nodes, *dxyzdZ);
+}

Numeric/JacobianBasis.h

@@ -14,8 +14,6 @@
 
 class JacobianBasis {
  private:
-  const bezierBasis *bezier;
-
   fullMatrix<double> gradShapeMatX, gradShapeMatY, gradShapeMatZ;
   fullMatrix<double> gradShapeMatXFast, gradShapeMatYFast, gradShapeMatZFast;
   fullVector<double> primGradShapeBarycenterX, primGradShapeBarycenterY, primGradShapeBarycenterZ;
@@ -41,6 +39,8 @@ class JacobianBasis {
                                        fullMatrix<double> &JDJ) const;
 
  public :
+  const bezierBasis *bezier;
+
   JacobianBasis(int tag);
 
   // Get methods
@@ -106,6 +106,12 @@ class JacobianBasis {
   }
 
   // Jacobian basis order and pyramidal basis
+  void getGradientsFromNodes(const fullMatrix<double> &nodes,
+                             fullMatrix<double> *dxyzdX,
+                             fullMatrix<double> *dxyzdY,
+                             fullMatrix<double> *dxyzdZ) const;
+
+  //
   static int jacobianOrder(int parentType, int order);
   static fullMatrix<double> generateJacMonomialsPyramid(int order);
   static fullMatrix<double> generateJacPointsPyramid(int order);

Numeric/bezierBasis.cpp

@@ -247,6 +247,10 @@ namespace {
       return subPoints;
     }
 
+    Msg::Error("Subpoints for pyramid p>1 not implemented !");
+    return std::vector< fullMatrix<double> >();
+
+    /*
     std::vector< fullMatrix<double> > subPoints(8);
     fullMatrix<double> ref, prox;
 
@@ -296,6 +300,7 @@ namespace {
     }
 
     return subPoints;
+    */
   }
 
   // Matrices generation
@@ -378,18 +383,18 @@ namespace {
       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));
-
+            * 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))
-            * 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));
+            * 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));
       }
     }
     return bez2Lag;
@@ -458,16 +463,23 @@ void bezierBasis::_construct(int parentType, int p)
     dim = 3;
     numLagCoeff = 8;
     exponents = JacobianBasis::generateJacMonomialsPyramid(order);
-    subPoints = generateSubPointsPyr(order);
-
-    numDivisions = static_cast<int>(subPoints.size());
+    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));
 
     matrixBez2Lag = generateBez2LagMatrixPyramid(exponents, bezierPoints, order);
     matrixBez2Lag.invert(matrixLag2Bez);
-    subDivisor = generateSubDivisorPyramid(exponents, subPoints, matrixLag2Bez, order);
+    if (order < 2)
+      subDivisor = generateSubDivisorPyramid(exponents, subPoints, matrixLag2Bez, order);
     return;
   }