dg : precompute coordinates

Geo/MElement.cpp

@@ -1207,6 +1207,26 @@ const fullMatrix<double> &MElement::getGradShapeFunctionsAtIntegrationPoints (in
   }
   return mat;
 }
+
+const fullMatrix<double> &MElement::getShapeFunctionsAtIntegrationPoints (int integrationOrder, int functionSpaceOrder) {
+  static std::map <std::pair<int,int>, fullMatrix<double> > F;
+  const polynomialBasis *fs = getFunctionSpace (functionSpaceOrder);
+  fullMatrix<double> &mat = F [std::make_pair(fs->type, integrationOrder)];
+  if (mat.size1()!=0) return mat;
+  int npts;
+  IntPt *pts;
+  getIntegrationPoints (integrationOrder, &npts, &pts);
+  mat.resize (fs->points.size1(), npts*3);
+  double f[512];
+  for (int i = 0; i < npts; i++) {
+    fs->f (pts[i].pt[0], pts[i].pt[1], pts[i].pt[2], f);
+    for (int j = 0; j < fs->points.size1(); j++) {
+      mat (j, i) = f[j];
+    }
+  }
+  return mat;
+}
+
 const fullMatrix<double> &MElement::getGradShapeFunctionsAtNodes (int functionSpaceOrder) {
   static std::map <int, fullMatrix<double> > DF;
   const polynomialBasis *fs = getFunctionSpace (functionSpaceOrder);

Geo/MElement.h

@@ -239,6 +239,8 @@ class MElement
                                      int order=-1);
   const fullMatrix<double> &getGradShapeFunctionsAtIntegrationPoints
     (int integrationOrder, int functionSpaceOrder=-1);
+  const fullMatrix<double> &getShapeFunctionsAtIntegrationPoints
+    (int integrationOrder, int functionSpaceOrder=-1);
   const fullMatrix<double> &getGradShapeFunctionsAtNodes (int functionSpaceOrder=-1);
 
   // return the Jacobian of the element evaluated at point (u,v,w) in

Numeric/fullMatrix.h

@@ -420,7 +420,7 @@ class fullMatrix
   }
 #endif
   ;
-  inline fullMatrix<scalar> transpose()
+  inline fullMatrix<scalar> transpose() const
   {
     fullMatrix<scalar> T(size2(), size1());
     for(int i = 0; i < size1(); i++)