From fcbf93d378a54cd6d65a495c1c958aa70bc8baac Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Thu, 3 Mar 2011 14:05:49 +0000
Subject: [PATCH] dg : precompute coordinates

---
 Geo/MElement.cpp     | 20 ++++++++++++++++++++
 Geo/MElement.h       |  2 ++
 Numeric/fullMatrix.h |  2 +-
 3 files changed, 23 insertions(+), 1 deletion(-)

diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index e4b8571cf4..40ca858f71 100644
--- a/Geo/MElement.cpp
+++ b/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);
diff --git a/Geo/MElement.h b/Geo/MElement.h
index cc10d2cc8c..22133a09f4 100644
--- a/Geo/MElement.h
+++ b/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
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 85a250c70e..2c9e7a6122 100644
--- a/Numeric/fullMatrix.h
+++ b/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++)
-- 
GitLab