From f2dedb9eb2b5b98c7da0bd79d84cfd137fac37f7 Mon Sep 17 00:00:00 2001
From: Gauthier Becker <gauthierbecker@gmail.com>
Date: Thu, 2 Aug 2012 14:38:57 +0000
Subject: [PATCH] Cache f,gradf, hessf and thirdDevf in the function space for
 3D elements to be faster (up to 50%) the values are stored for each Gauss
 points

---
 Geo/MElement.cpp | 30 ++++++++++++++++++++++++++++++
 Geo/MElement.h   |  3 +++
 2 files changed, 33 insertions(+)

diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 2d54c47d4f..ca89017d99 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -369,6 +369,24 @@ double MElement::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
   return _computeDeterminantAndRegularize(this, jac);
 }
 
+double MElement::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])
+{
+  jac[0][0] = jac[0][1] = jac[0][2] = 0.;
+  jac[1][0] = jac[1][1] = jac[1][2] = 0.;
+  jac[2][0] = jac[2][1] = jac[2][2] = 0.;
+
+  for (int i = 0; i < getNumShapeFunctions(); i++) {
+    const MVertex *v = getShapeFunctionNode(i);
+    for (int j = 0; j < 3; j++) {
+      double mult = gsf[i][j];
+      jac[j][0] += v->x() * mult;
+      jac[j][1] += v->y() * mult;
+      jac[j][2] += v->z() * mult;
+    }
+  }
+  return _computeDeterminantAndRegularize(this, jac);
+}
+
 double MElement::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
 {
   jac[0][0] = jac[0][1] = jac[0][2] = 0.;
@@ -404,6 +422,18 @@ void MElement::pnt(double u, double v, double w, SPoint3 &p)
   p = SPoint3(x, y, z);
 }
 
+void MElement::pnt(const std::vector<double> &sf, SPoint3 &p)
+{
+  double x = 0., y = 0., z = 0.;
+  for (int j = 0; j < getNumShapeFunctions(); j++) {
+    const MVertex *v = getShapeFunctionNode(j);
+    x += sf[j] * v->x();
+    y += sf[j] * v->y();
+    z += sf[j] * v->z();
+  }
+  p = SPoint3(x, y, z);
+}
+
 void MElement::primaryPnt(double u, double v, double w, SPoint3 &p)
 {
   double x = 0., y = 0., z = 0.;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 5d04899653..f7b824e3b9 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -243,6 +243,8 @@ class MElement
   // return the Jacobian of the element evaluated at point (u,v,w) in
   // parametric coordinates
   double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
+  // To be compatible with _vgrads of functionSpace without having to put under fullMatrix form
+  double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]);
   double getJacobian(double u, double v, double w, double jac[3][3]);
   inline double getJacobian(double u, double v, double w, fullMatrix<double> &j){
     double JAC[3][3];
@@ -266,6 +268,7 @@ class MElement
   // get the point in cartesian coordinates corresponding to the point
   // (u,v,w) in parametric coordinates
   virtual void pnt(double u, double v, double w, SPoint3 &p);
+  virtual void pnt(const std::vector<double> &sf,SPoint3 &p); // To be compatible with functionSpace without changing form
   virtual void primaryPnt(double u, double v, double w, SPoint3 &p);
 
   // invert the parametrisation
-- 
GitLab