From e68b5cfe68310edf93ae670bdffaf061f73cd6b8 Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Wed, 1 Oct 2014 22:02:17 +0000
Subject: [PATCH] Fixed error in computation of inverse conditioning

---
 Numeric/JacobianBasis.cpp | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index b26a4b4f5a..828b8973c9 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -156,7 +156,7 @@ inline void calcIDI1D(double dxdX, double dxdY, double dxdZ,
                     + dydX*dydX + dydY*dydY + dydZ*dydZ
                     + dzdX*dzdX + dzdY*dzdY + dzdZ*dzdZ;
   const double invProd = 1./(nJSq*nInvJSq), sqrtInvProd = sqrt(invProd);
-  IDI(i, 3*numMapNodes) = sqrtInvProd;
+  IDI(i, 3*numMapNodes) = 3.*sqrtInvProd;
   for (int j = 0; j < numMapNodes; j++) {
     const double &dPhidX = gSMatX(i, j);
     const double &dJdxj = JDJ(i, j), &dJdyj = JDJ(i, j+numMapNodes),
@@ -181,6 +181,7 @@ inline void calcIDI1D(double dxdX, double dxdY, double dxdZ,
 
 // Compute inverse condition number and its gradients
 // w.r.t. node positions, at one location in a 2D element
+// FIXME: Precision issues?
 inline void calcIDI2D(double dxdX, double dxdY, double dxdZ,
                       double dydX, double dydY, double dydZ,
                       double dzdX, double dzdY, double dzdZ,
@@ -202,7 +203,7 @@ inline void calcIDI2D(double dxdX, double dxdY, double dxdZ,
                     + dydX*dydX + dydY*dydY + dydZ*dydZ
                     + dzdX*dzdX + dzdY*dzdY + dzdZ*dzdZ;
   const double invProd = 1./(nJSq*nInvJSq), sqrtInvProd = sqrt(invProd);
-  IDI(i, 3*numMapNodes) = 2.*sqrtInvProd;
+  IDI(i, 3*numMapNodes) = 3.*sqrtInvProd;
   for (int j = 0; j < numMapNodes; j++) {
     const double &dPhidX = gSMatX(i, j);
     const double &dPhidY = gSMatY(i, j);
@@ -963,8 +964,6 @@ inline void JacobianBasis::getInvCondGeneral(int nJacNodes, const fullMatrix<dou
                                                    dydX, dydY, dydZ,
                                                    dzdX, dzdY, dzdZ);
         invCond(i) = 3./sqrt(nJSq*nInvJSq);
-        std::cout << "DBGTT: normal = (" << dxdZ << ", " << dydZ<< ", " << dzdZ
-                  << "), nJSq = " << nJSq << ", nInvJSq = " << nInvJSq << "\n";
       }
       break;
     }
-- 
GitLab