From 4cda89344e96aa596d9704ac7a34da89dbc5e8b7 Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Wed, 15 Oct 2014 07:14:20 +0000
Subject: [PATCH] fix bug 2D element with zero Jacobian determinant

---
 Numeric/JacobianBasis.cpp | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index c443f09fea..2bf7bc4fe0 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -421,6 +421,10 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes,
     case 1 : {
       fullMatrix<double> normals(2,3);
       const double invScale = getPrimNormals1D(nodesXYZ,normals);
+      if (invScale == 0) {
+        for (int i = 0; i < nJacNodes; i++) jacobian(i) = 0;
+        return;
+      }
       if (scaling) {
         const double scale = 1./invScale;
         normals(0,0) *= scale; normals(0,1) *= scale; normals(0,2) *= scale;            // Faster to scale 1 normal than afterwards
@@ -441,6 +445,10 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes,
     case 2 : {
       fullMatrix<double> normal(1,3);
       const double invScale = getPrimNormal2D(nodesXYZ, normal, idealNorm);
+      if (invScale == 0) {
+        for (int i = 0; i < nJacNodes; i++) jacobian(i) = 0;
+        return;
+      }
       if (scaling) {
         const double scale = 1./invScale;
         normal(0,0) *= scale; normal(0,1) *= scale; normal(0,2) *= scale;               // Faster to scale normal than afterwards
@@ -547,6 +555,10 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes,
         }
         fullMatrix<double> normals(2,3);
         const double invScale = getPrimNormals1D(nodesXYZ,normals);
+        if (invScale == 0) {
+          for (int i = 0; i < nJacNodes; i++) jacobian(i,iEl) = 0;
+          continue;
+        }
         if (scaling) {
           const double scale = 1./invScale;
           normals(0,0) *= scale; normals(0,1) *= scale; normals(0,2) *= scale;                // Faster to scale 1 normal than afterwards
@@ -576,6 +588,10 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes,
         }
         fullMatrix<double> normal(1,3);
         const double invScale = getPrimNormal2D(nodesXYZ, normal, idealNorm);
+        if (invScale == 0) {
+          for (int i = 0; i < nJacNodes; i++) jacobian(i,iEl) = 0;
+          continue;
+        }
         if (scaling) {
           const double scale = 1./invScale;
           normal(0,0) *= scale; normal(0,1) *= scale; normal(0,2) *= scale;                   // Faster to scale normal than afterwards
-- 
GitLab