From c1f05e3883b6617084b6e28bf35bd8aae889ca55 Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Tue, 18 Jun 2013 11:16:13 +0000
Subject: [PATCH] revert monomials to be fullMatrix<double> + Add
 pow_int(double, double)

---
 FunctionSpace/BasisLagrange.h |   4 +-
 Numeric/nodalBasis.h          |   2 -
 Numeric/polynomialBasis.cpp   | 181 ++++++++++++++++++----------------
 Numeric/polynomialBasis.h     |  11 ++-
 Numeric/pyramidalBasis.h      |   1 +
 Post/PViewData.cpp            |  64 ------------
 Post/PViewData.h              |  13 ---
 7 files changed, 107 insertions(+), 169 deletions(-)

diff --git a/FunctionSpace/BasisLagrange.h b/FunctionSpace/BasisLagrange.h
index 7e340a9e96..37d85d31e3 100644
--- a/FunctionSpace/BasisLagrange.h
+++ b/FunctionSpace/BasisLagrange.h
@@ -58,7 +58,7 @@ class BasisLagrange: public BasisLocal{
     getPreEvaluatedDerivatives(unsigned int orientation) const;
 
   const fullMatrix<double>& getCoefficient(void) const;
-  const fullMatrix<int>& getMonomial(void) const;
+  const fullMatrix<double>& getMonomial(void) const;
 
   std::vector<double>
     project(const MElement& element,
@@ -126,7 +126,7 @@ getCoefficient(void) const{
   return lBasis->coefficients;
 }
 
-inline const fullMatrix<int>& BasisLagrange::
+inline const fullMatrix<double>& BasisLagrange::
 getMonomial(void) const{
   return lBasis->monomials;
 }
diff --git a/Numeric/nodalBasis.h b/Numeric/nodalBasis.h
index 8e5cf93ac8..34f7216a74 100644
--- a/Numeric/nodalBasis.h
+++ b/Numeric/nodalBasis.h
@@ -14,8 +14,6 @@ class nodalBasis {
   int type, parentType, order, dimension, numFaces;
   bool serendip;
   fullMatrix<double> points;
-  fullMatrix<int> monomials;
-  fullMatrix<double> coefficients;
 
   nodalBasis(int tag);
   virtual ~nodalBasis() {}
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 9c1c98fb1b..dadc221bca 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -46,7 +46,7 @@ static void printClosure(polynomialBasis::clCont &fullClosure,
 
 
 static fullMatrix<double> generateLagrangeMonomialCoefficients
-  (const fullMatrix<int>& monomial, const fullMatrix<double>& point)
+  (const fullMatrix<double>& monomial, const fullMatrix<double>& point)
 {
   if(monomial.size1() != point.size1() || monomial.size2() != point.size2()){
     Msg::Fatal("Wrong sizes for Lagrange coefficients generation %d %d -- %d %d",
@@ -75,34 +75,41 @@ static fullMatrix<double> generateLagrangeMonomialCoefficients
 
 polynomialBasis::polynomialBasis(int tag) : nodalBasis(tag)
 {
-
+  fullMatrix<int> monom;
   switch (parentType) {
   case TYPE_PNT :
-    monomials = gmshGenerateMonomialsLine(0);
+    monom = gmshGenerateMonomialsLine(0);
     break;
   case TYPE_LIN :
-    monomials = gmshGenerateMonomialsLine(order);
+    monom = gmshGenerateMonomialsLine(order);
     break;
   case TYPE_TRI :
-    monomials = gmshGenerateMonomialsTriangle(order, serendip);
+    monom = gmshGenerateMonomialsTriangle(order, serendip);
     break;
   case TYPE_QUA :
-    monomials = serendip ? gmshGenerateMonomialsQuadSerendipity(order) :
+    monom = serendip ? gmshGenerateMonomialsQuadSerendipity(order) :
     gmshGenerateMonomialsQuadrangle(order);
     break;
   case TYPE_TET :
-    monomials = gmshGenerateMonomialsTetrahedron(order, serendip);
+    monom = gmshGenerateMonomialsTetrahedron(order, serendip);
     break;
   case TYPE_PRI :
-    monomials = serendip ? gmshGenerateMonomialsPrismSerendipity(order) :
+    monom = serendip ? gmshGenerateMonomialsPrismSerendipity(order) :
     gmshGenerateMonomialsPrism(order);
     break;
   case TYPE_HEX :
-    monomials = serendip ? gmshGenerateMonomialsHexaSerendipity(order) :
+    monom = serendip ? gmshGenerateMonomialsHexaSerendipity(order) :
     gmshGenerateMonomialsHexahedron(order);
     break;
   }
 
+  monomials.resize(monom.size1(), monom.size2());
+  for (int i = 0; i < monom.size1(); ++i) {
+    for (int j = 0; j < monom.size2(); ++j) {
+      monomials(i, j) = static_cast<double>(monom(i, j));
+    }
+  }
+
   coefficients = generateLagrangeMonomialCoefficients(monomials, points);
 }
 
@@ -187,7 +194,7 @@ void polynomialBasis::df(double u, double v, double w, double grads[][3]) const
       for(int j = 0; j < coefficients.size2(); j++){
         if (monomials(j, 0) > 0)
           grads[i][0] += coefficients(i, j) *
-            pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0);
+            pow_int(u, (monomials(j, 0) - 1)) * monomials(j, 0);
       }
     }
     break;
@@ -199,12 +206,12 @@ void polynomialBasis::df(double u, double v, double w, double grads[][3]) const
       for(int j = 0; j < coefficients.size2(); j++){
         if (monomials(j, 0) > 0)
           grads[i][0] += coefficients(i, j) *
-            pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0) *
-            pow_int(v, (int)monomials(j, 1));
+            pow_int(u, (monomials(j, 0) - 1)) * monomials(j, 0) *
+            pow_int(v, monomials(j, 1));
         if (monomials(j, 1) > 0)
           grads[i][1] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0)) *
-            pow_int(v, (int)(monomials(j, 1) - 1)) * monomials(j, 1);
+            pow_int(u, monomials(j, 0)) *
+            pow_int(v, (monomials(j, 1) - 1)) * monomials(j, 1);
       }
     }
     break;
@@ -216,19 +223,19 @@ void polynomialBasis::df(double u, double v, double w, double grads[][3]) const
       for(int j = 0; j < coefficients.size2(); j++){
         if (monomials(j, 0) > 0)
           grads[i][0] += coefficients(i, j) *
-            pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0) *
-            pow_int(v, (int)monomials(j, 1)) *
-            pow_int(w, (int)monomials(j, 2));
+            pow_int(u, (monomials(j, 0) - 1)) * monomials(j, 0) *
+            pow_int(v, monomials(j, 1)) *
+            pow_int(w, monomials(j, 2));
         if (monomials(j, 1) > 0)
           grads[i][1] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0)) *
-            pow_int(v, (int)(monomials(j, 1) - 1)) * monomials(j, 1) *
-            pow_int(w, (int)monomials(j, 2));
+            pow_int(u, monomials(j, 0)) *
+            pow_int(v, (monomials(j, 1) - 1)) * monomials(j, 1) *
+            pow_int(w, monomials(j, 2));
         if (monomials(j, 2) > 0)
           grads[i][2] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0)) *
-            pow_int(v, (int)monomials(j, 1)) *
-            pow_int(w, (int)(monomials(j, 2) - 1)) * monomials(j, 2);
+            pow_int(u, monomials(j, 0)) *
+            pow_int(v, monomials(j, 1)) *
+            pow_int(w, (monomials(j, 2) - 1)) * monomials(j, 2);
       }
     }
     break;
@@ -248,7 +255,7 @@ void polynomialBasis::ddf(double u, double v, double w, double hess[][3][3]) con
 
       for(int j = 0; j < coefficients.size2(); j++){
         if (monomials(j, 0) > 1) // second derivative !=0
-          hess[i][0][0] += coefficients(i, j) * pow_int(u, (int)(monomials(j, 0) - 2)) *
+          hess[i][0][0] += coefficients(i, j) * pow_int(u, (monomials(j, 0) - 2)) *
             monomials(j, 0) * (monomials(j, 0) - 1);
       }
     }
@@ -260,16 +267,16 @@ void polynomialBasis::ddf(double u, double v, double w, double hess[][3][3]) con
       hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0;
       for(int j = 0; j < coefficients.size2(); j++){
         if (monomials(j, 0) > 1) // second derivative !=0
-          hess[i][0][0] += coefficients(i, j) * pow_int(u, (int)(monomials(j, 0) - 2)) *
-            monomials(j, 0) * (monomials(j, 0) - 1) * pow_int(v, (int)monomials(j, 1));
+          hess[i][0][0] += coefficients(i, j) * pow_int(u, (monomials(j, 0) - 2)) *
+            monomials(j, 0) * (monomials(j, 0) - 1) * pow_int(v, monomials(j, 1));
         if ((monomials(j, 1) > 0) && (monomials(j, 0) > 0))
           hess[i][0][1] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) *
-            pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1);
+            pow_int(u, monomials(j, 0) - 1) * monomials(j, 0) *
+            pow_int(v, monomials(j, 1) - 1) * monomials(j, 1);
         if (monomials(j, 1) > 1)
           hess[i][1][1] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0)) *
-            pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1);
+            pow_int(u, monomials(j, 0)) *
+            pow_int(v, monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1);
       }
       hess[i][1][0] = hess[i][0][1];
     }
@@ -282,35 +289,35 @@ void polynomialBasis::ddf(double u, double v, double w, double hess[][3][3]) con
       for(int j = 0; j < coefficients.size2(); j++){
         if (monomials(j, 0) > 1)
           hess[i][0][0] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0) * (monomials(j, 0)-1) *
-            pow_int(v, (int)monomials(j, 1)) *
-            pow_int(w, (int)monomials(j, 2));
+            pow_int(u, monomials(j, 0) - 2) * monomials(j, 0) * (monomials(j, 0)-1) *
+            pow_int(v, monomials(j, 1)) *
+            pow_int(w, monomials(j, 2));
 
         if ((monomials(j, 0) > 0) && (monomials(j, 1) > 0))
           hess[i][0][1] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) *
-            pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1) *
-            pow_int(w, (int)monomials(j, 2));
+            pow_int(u, monomials(j, 0) - 1) * monomials(j, 0) *
+            pow_int(v, monomials(j, 1) - 1) * monomials(j, 1) *
+            pow_int(w, monomials(j, 2));
         if ((monomials(j, 0) > 0) && (monomials(j, 2) > 0))
           hess[i][0][2] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) *
-            pow_int(v, (int)monomials(j, 1)) *
-            pow_int(w, (int)monomials(j, 2) - 1) * monomials(j, 2);
+            pow_int(u, monomials(j, 0) - 1) * monomials(j, 0) *
+            pow_int(v, monomials(j, 1)) *
+            pow_int(w, monomials(j, 2) - 1) * monomials(j, 2);
         if (monomials(j, 1) > 1)
           hess[i][1][1] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0)) *
-            pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1)-1) *
-            pow_int(w, (int)monomials(j, 2));
+            pow_int(u, monomials(j, 0)) *
+            pow_int(v, monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1)-1) *
+            pow_int(w, monomials(j, 2));
         if ((monomials(j, 1) > 0) && (monomials(j, 2) > 0))
           hess[i][1][2] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0)) *
-            pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1) *
-            pow_int(w, (int)monomials(j, 2) - 1) * monomials(j, 2);
+            pow_int(u, monomials(j, 0)) *
+            pow_int(v, monomials(j, 1) - 1) * monomials(j, 1) *
+            pow_int(w, monomials(j, 2) - 1) * monomials(j, 2);
         if (monomials(j, 2) > 1)
           hess[i][2][2] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0)) *
-            pow_int(v, (int)monomials(j, 1)) *
-            pow_int(w, (int)monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1);
+            pow_int(u, monomials(j, 0)) *
+            pow_int(v, monomials(j, 1)) *
+            pow_int(w, monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1);
       }
       hess[i][1][0] = hess[i][0][1];
       hess[i][2][0] = hess[i][0][2];
@@ -334,7 +341,7 @@ void polynomialBasis::dddf(double u, double v, double w, double third[][3][3][3]
 
       for(int j = 0; j < coefficients.size2(); j++){
         if (monomials(j, 0) > 2) // third derivative !=0
-          third[i][0][0][0] += coefficients(i, j) * pow_int(u, (int)(monomials(j, 0) - 3)) *
+          third[i][0][0][0] += coefficients(i, j) * pow_int(u, (monomials(j, 0) - 3)) *
             monomials(j, 0) * (monomials(j, 0) - 1)*(monomials(j, 0) - 2);
       }
     }
@@ -348,20 +355,20 @@ void polynomialBasis::dddf(double u, double v, double w, double third[][3][3][3]
       for(int j = 0; j < coefficients.size2(); j++){
         if (monomials(j, 0) > 2) // second derivative !=0
           third[i][0][0][0] += coefficients(i, j) *
-            pow_int(u, (int)(monomials(j, 0) - 3))*monomials(j, 0) * (monomials(j, 0) - 1) *(monomials(j, 0) - 2) *
-            pow_int(v, (int)monomials(j, 1));
+            pow_int(u, (monomials(j, 0) - 3))*monomials(j, 0) * (monomials(j, 0) - 1) *(monomials(j, 0) - 2) *
+            pow_int(v, monomials(j, 1));
         if ((monomials(j, 0) > 1) && (monomials(j, 1) > 0))
           third[i][0][0][1] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0)*(monomials(j, 0) - 1) *
-            pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1);
+            pow_int(u, monomials(j, 0) - 2) * monomials(j, 0)*(monomials(j, 0) - 1) *
+            pow_int(v, monomials(j, 1) - 1) * monomials(j, 1);
         if ((monomials(j, 0) > 0) && (monomials(j, 1) > 1))
           third[i][0][1][1] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0) - 1) *monomials(j, 0)*
-            pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1);
+            pow_int(u, monomials(j, 0) - 1) *monomials(j, 0)*
+            pow_int(v, monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1);
         if (monomials(j, 1) > 2)
           third[i][1][1][1] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0)) *
-            pow_int(v, (int)monomials(j, 1) - 3) * monomials(j, 1) * (monomials(j, 1) - 1)*(monomials(j, 1) - 2);
+            pow_int(u, monomials(j, 0)) *
+            pow_int(v, monomials(j, 1) - 3) * monomials(j, 1) * (monomials(j, 1) - 1)*(monomials(j, 1) - 2);
       }
       third[i][0][1][0] = third[i][0][0][1];
       third[i][1][0][0] = third[i][0][0][1];
@@ -378,62 +385,62 @@ void polynomialBasis::dddf(double u, double v, double w, double third[][3][3][3]
       for(int j = 0; j < coefficients.size2(); j++){
         if (monomials(j, 0) > 2) // second derivative !=0
           third[i][0][0][0] += coefficients(i, j) *
-            pow_int(u, (int)(monomials(j, 0) - 3)) *monomials(j, 0) * (monomials(j, 0) - 1)*(monomials(j, 0) - 2) *
-            pow_int(v, (int)monomials(j, 1))*
-            pow_int(w, (int)monomials(j, 2));
+            pow_int(u, (monomials(j, 0) - 3)) *monomials(j, 0) * (monomials(j, 0) - 1)*(monomials(j, 0) - 2) *
+            pow_int(v, monomials(j, 1))*
+            pow_int(w, monomials(j, 2));
 
         if ((monomials(j, 0) > 1) && (monomials(j, 1) > 0))
           third[i][0][0][1] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0)*(monomials(j, 0) - 1) *
-            pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1)*
-            pow_int(w, (int)monomials(j, 2));
+            pow_int(u, monomials(j, 0) - 2) * monomials(j, 0)*(monomials(j, 0) - 1) *
+            pow_int(v, monomials(j, 1) - 1) * monomials(j, 1)*
+            pow_int(w, monomials(j, 2));
 
         if ((monomials(j, 0) > 1) && (monomials(j, 2) > 0))
           third[i][0][0][2] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0)*(monomials(j, 0) - 1) *
-            pow_int(v, (int)monomials(j, 1)) *
-            pow_int(w, (int)monomials(j, 2) - 1)* monomials(j, 2);
+            pow_int(u, monomials(j, 0) - 2) * monomials(j, 0)*(monomials(j, 0) - 1) *
+            pow_int(v, monomials(j, 1)) *
+            pow_int(w, monomials(j, 2) - 1)* monomials(j, 2);
 
         if ((monomials(j, 0) > 0) && (monomials(j, 1) > 1))
           third[i][0][1][1] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0) - 1) *monomials(j, 0)*
-            pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1)*
-            pow_int(w, (int)monomials(j, 2));
+            pow_int(u, monomials(j, 0) - 1) *monomials(j, 0)*
+            pow_int(v, monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1)*
+            pow_int(w, monomials(j, 2));
 
         if ((monomials(j, 0) > 0) && (monomials(j, 1) > 0)&& (monomials(j, 2) > 0))
           third[i][0][1][2] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0) - 1) *monomials(j, 0)*
-            pow_int(v, (int)monomials(j, 1) - 1) *monomials(j, 1) *
-            pow_int(w, (int)monomials(j, 2) - 1) *monomials(j, 2);
+            pow_int(u, monomials(j, 0) - 1) *monomials(j, 0)*
+            pow_int(v, monomials(j, 1) - 1) *monomials(j, 1) *
+            pow_int(w, monomials(j, 2) - 1) *monomials(j, 2);
 
         if ((monomials(j, 0) > 0) && (monomials(j, 2) > 1))
           third[i][0][2][2] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0) - 1) *monomials(j, 0)*
-            pow_int(v, (int)monomials(j, 1))*
-            pow_int(w, (int)monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1);
+            pow_int(u, monomials(j, 0) - 1) *monomials(j, 0)*
+            pow_int(v, monomials(j, 1))*
+            pow_int(w, monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1);
         if ((monomials(j, 1) > 2))
           third[i][1][1][1] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0)) *
-            pow_int(v, (int)monomials(j, 1)-3) * monomials(j, 1) * (monomials(j, 1) - 1)*(monomials(j, 1) - 2)*
-            pow_int(w, (int)monomials(j, 2));
+            pow_int(u, monomials(j, 0)) *
+            pow_int(v, monomials(j, 1)-3) * monomials(j, 1) * (monomials(j, 1) - 1)*(monomials(j, 1) - 2)*
+            pow_int(w, monomials(j, 2));
 
         if ((monomials(j, 1) > 1) && (monomials(j, 2) > 0))
           third[i][1][1][2] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0)) *
-            pow_int(v, (int)monomials(j, 1)-2) * monomials(j, 1) * (monomials(j, 1) - 1)*
-            pow_int(w, (int)monomials(j, 2) - 1) *monomials(j, 2) ;
+            pow_int(u, monomials(j, 0)) *
+            pow_int(v, monomials(j, 1)-2) * monomials(j, 1) * (monomials(j, 1) - 1)*
+            pow_int(w, monomials(j, 2) - 1) *monomials(j, 2) ;
 
         if ((monomials(j, 1) > 0) && (monomials(j, 2) > 1))
           third[i][1][2][2] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0)) *
-            pow_int(v, (int)monomials(j, 1)-1) *monomials(j, 1)*
-            pow_int(w, (int)monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1) ;
+            pow_int(u, monomials(j, 0)) *
+            pow_int(v, monomials(j, 1)-1) *monomials(j, 1)*
+            pow_int(w, monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1) ;
 
         if ((monomials(j, 2) > 2))
           third[i][2][2][2] += coefficients(i, j) *
-            pow_int(u, (int)monomials(j, 0)) *
-            pow_int(v, (int)monomials(j, 1))*
-            pow_int(w, (int)monomials(j, 2) - 3) * monomials(j, 2) * (monomials(j, 2) - 1)*(monomials(j, 2) - 2);
+            pow_int(u, monomials(j, 0)) *
+            pow_int(v, monomials(j, 1))*
+            pow_int(w, monomials(j, 2) - 3) * monomials(j, 2) * (monomials(j, 2) - 1)*(monomials(j, 2) - 2);
 
 
       }
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 9868a9e045..66c7c27a9f 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -60,15 +60,24 @@ inline double pow_int(const double &a, const int &n)
       return a4*a4*a2;
     }
   default :
-    return pow_int(a,n-1)*a;
+    return pow_int(a, n-9) * pow_int(a, 9);
   }
 }
 
+inline double pow_int(const double &a, const double &d)
+{
+  // Round double !
+  int n = static_cast<int>(d + .5);
+  return pow_int(a, n);
+}
+
 class polynomialBasis : public nodalBasis
 {
  public:
   // for now the only implemented polynomial basis are nodal poly
   // basis, we use the type of the corresponding gmsh element as type
+  fullMatrix<double> monomials;
+  fullMatrix<double> coefficients;
 
   polynomialBasis(int tag);
   ~polynomialBasis();
diff --git a/Numeric/pyramidalBasis.h b/Numeric/pyramidalBasis.h
index 8b0982b02f..db87118fc5 100644
--- a/Numeric/pyramidalBasis.h
+++ b/Numeric/pyramidalBasis.h
@@ -18,6 +18,7 @@ class pyramidalBasis: public nodalBasis
  private:
   // Orthogonal basis for the pyramid
   BergotBasis *bergot;
+  fullMatrix<double> coefficients;
 
  public:
   pyramidalBasis(int tag);
diff --git a/Post/PViewData.cpp b/Post/PViewData.cpp
index 1ce16372cd..6bcc7c11af 100644
--- a/Post/PViewData.cpp
+++ b/Post/PViewData.cpp
@@ -122,23 +122,6 @@ void PViewData::setInterpolationMatrices(int type,
   _interpolation[type].push_back(new fullMatrix<double>(expVal));
 }
 
-void PViewData::setInterpolationMatrices(int type,
-                                         const fullMatrix<double> &coefVal,
-                                         const fullMatrix<int> &expVal)
-{
-  if(!type || _interpolation[type].size()) return;
-
-  fullMatrix<double> *dExpVal;
-  dExpVal = new fullMatrix<double>(expVal.size1(), expVal.size2());
-  for (int i = 0; i < expVal.size1(); ++i) {
-    for (int j = 0; i < expVal.size2(); ++j) {
-      (*dExpVal)(i, j) = static_cast<double>(expVal(i, j));
-    }
-  }
-  _interpolation[type].push_back(new fullMatrix<double>(coefVal));
-  _interpolation[type].push_back(dExpVal);
-}
-
 void PViewData::setInterpolationMatrices(int type,
                                          const fullMatrix<double> &coefVal,
                                          const fullMatrix<double> &expVal,
@@ -151,53 +134,6 @@ void PViewData::setInterpolationMatrices(int type,
   _interpolation[type].push_back(new fullMatrix<double>(coefGeo));
   _interpolation[type].push_back(new fullMatrix<double>(expGeo));
 }
-void PViewData::setInterpolationMatrices(int type,
-                                         const fullMatrix<double> &coefVal,
-                                         const fullMatrix<int> &expVal,
-                                         const fullMatrix<double> &coefGeo,
-                                         const fullMatrix<int> &expGeo)
-{
-  if(!type || _interpolation[type].size()) return;
-
-  fullMatrix<double> *dExpVal, *dExpGeo;
-  dExpVal = new fullMatrix<double>(expVal.size1(), expVal.size2());
-  dExpGeo = new fullMatrix<double>(expGeo.size1(), expGeo.size2());
-
-  for (int i = 0; i < expVal.size1(); ++i) {
-    for (int j = 0; i < expVal.size2(); ++j) {
-      (*dExpVal)(i, j) = static_cast<double>(expVal(i, j));
-    }
-  }
-  for (int i = 0; i < expGeo.size1(); ++i) {
-    for (int j = 0; i < expGeo.size2(); ++j) {
-      (*dExpGeo)(i, j) = static_cast<double>(expGeo(i, j));
-    }
-  }
-  _interpolation[type].push_back(new fullMatrix<double>(coefVal));
-  _interpolation[type].push_back(dExpVal);
-  _interpolation[type].push_back(new fullMatrix<double>(coefGeo));
-  _interpolation[type].push_back(dExpGeo);
-}
-void PViewData::setInterpolationMatrices(int type,
-                                         const fullMatrix<double> &coefVal,
-                                         const fullMatrix<double> &expVal,
-                                         const fullMatrix<double> &coefGeo,
-                                         const fullMatrix<int> &expGeo)
-{
-  if(!type || _interpolation[type].size()) return;
-
-  fullMatrix<double> *dExpGeo;
-  dExpGeo = new fullMatrix<double>(expGeo.size1(), expGeo.size2());
-  for (int i = 0; i < expGeo.size1(); ++i) {
-    for (int j = 0; i < expGeo.size2(); ++j) {
-      (*dExpGeo)(i, j) = static_cast<double>(expGeo(i, j));
-    }
-  }
-  _interpolation[type].push_back(new fullMatrix<double>(coefVal));
-  _interpolation[type].push_back(new fullMatrix<double>(expVal));
-  _interpolation[type].push_back(new fullMatrix<double>(coefGeo));
-  _interpolation[type].push_back(dExpGeo);
-}
 
 int PViewData::getInterpolationMatrices(int type, std::vector<fullMatrix<double>*> &p)
 {
diff --git a/Post/PViewData.h b/Post/PViewData.h
index 73c3ae5a55..6d6bd13f52 100644
--- a/Post/PViewData.h
+++ b/Post/PViewData.h
@@ -205,24 +205,11 @@ class PViewData {
   void setInterpolationMatrices(int type,
                                 const fullMatrix<double> &coefVal,
                                 const fullMatrix<double> &expVal);
-  void setInterpolationMatrices(int type,
-                                const fullMatrix<double> &coefVal,
-                                const fullMatrix<int> &expVal);
   void setInterpolationMatrices(int type,
                                 const fullMatrix<double> &coefVal,
                                 const fullMatrix<double> &expVal,
                                 const fullMatrix<double> &coefGeo,
                                 const fullMatrix<double> &expGeo);
-  void setInterpolationMatrices(int type,
-                                const fullMatrix<double> &coefVal,
-                                const fullMatrix<int> &expVal,
-                                const fullMatrix<double> &coefGeo,
-                                const fullMatrix<int> &expGeo);
-  void setInterpolationMatrices(int type,
-                                const fullMatrix<double> &coefVal,
-                                const fullMatrix<double> &expVal,
-                                const fullMatrix<double> &coefGeo,
-                                const fullMatrix<int> &expGeo);
   int getInterpolationMatrices(int type, std::vector<fullMatrix<double>*> &p);
   bool haveInterpolationMatrices(int type=0);
 
-- 
GitLab