From: Anthony Geay <anthony.geay@edf.fr>
Date: Mon, 19 Nov 2018 12:29:11 +0000 (+0100)
Subject: BUG correction : During P0P0 remapping on tetra. If tetra is very small the TetraAffi... 
X-Git-Tag: V9_2_0rc1^0
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=bf04b84b8d956d03fd96e2c88bcde46dd8bee769;p=tools%2Fmedcoupling.git

BUG correction : During P0P0 remapping on tetra. If tetra is very small the TetraAffine matrix could be falsely detected as not inversible.
---

diff --git a/src/INTERP_KERNEL/TetraAffineTransform.cxx b/src/INTERP_KERNEL/TetraAffineTransform.cxx
index be3b8090e..a06e5f16f 100644
--- a/src/INTERP_KERNEL/TetraAffineTransform.cxx
+++ b/src/INTERP_KERNEL/TetraAffineTransform.cxx
@@ -62,10 +62,14 @@ namespace INTERP_KERNEL
     calculateDeterminant();
     
     LOG(3, "determinant before inverse = " << _determinant);
+
+    double ni(1./INTERP_KERNEL::normInf(_linear_transform));
+    ni = ni*ni*ni;
     
     // check that tetra is non-planar -> determinant is not zero
+    // AGY : the check to 0. must integrate the infinite norm of _linear_transform matrix.
     // otherwise set _determinant to zero to signal caller that transformation did not work
-    if(epsilonEqual(_determinant, 0.0))
+    if(epsilonEqual(ni*_determinant, 0.0))
       {
         _determinant = 0.0;
         return;
diff --git a/src/INTERP_KERNEL/VectorUtils.hxx b/src/INTERP_KERNEL/VectorUtils.hxx
index 00b2bf042..21aaa723f 100644
--- a/src/INTERP_KERNEL/VectorUtils.hxx
+++ b/src/INTERP_KERNEL/VectorUtils.hxx
@@ -185,6 +185,23 @@ namespace INTERP_KERNEL
     return relError < relTol;
   }
 
+  inline double sumOfAbsoluteValues(const double row[3])
+  {
+    double ret(0.);
+    std::for_each(row,row+3,[&ret](double v) { ret += std::abs(v); });
+    return ret;
+  }
+
+  /*!
+   * Returns the infinite norm of a 3x3 input matrix \a mat.
+   * The max of absolute value of row sum.
+   */
+  inline double normInf(const double mat[9])
+  {
+    double ret(std::max(sumOfAbsoluteValues(mat),sumOfAbsoluteValues(mat+3)));
+    return std::max(ret,sumOfAbsoluteValues(mat+6));
+  }
+
 }
 
 #endif
diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py
index 679d8fe82..d6481d805 100644
--- a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py
+++ b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py
@@ -1278,6 +1278,27 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         rem2.setCrudeMatrix(ms,mt,"P0P0",mat)
         self.assertTrue(rem2.transferField(fs,1e300).getArray().isEqual(DataArrayDouble([54.25,11.75,79.25,16.75]),1e-12))
         pass
+
+    def testSmallTetraCell(self):
+        """This test is a non regression test. When using tetra/tetra P0P0 interpolation on very small cells the
+        3x3 matrix in the TetraAffine contains very small values and so the determinant is small (cubic).
+        So the tetra was detected as flat. Now the infinite norm of matrix is considered to establish if matrix is inversible or not."""
+        coords = [(-0.019866666666666668, 0.02, 0.002), (-0.020000073463967143, 0.019999926535763005, 0.0018666666666666673), (-0.020000073463967143, 0.019999926535763005, 0.002), (-0.020000072974206463, 0.019866593202430387, 0.002)]
+        m=MEDCouplingUMesh("mesh",3)
+        m.allocateCells()
+        m.insertNextCell(NORM_TETRA4,[0,1,2,3])
+        m.setCoords(DataArrayDouble(coords))
+        rem=MEDCouplingRemapper()
+        rem.setPrecision(1e-12)
+        rem.prepare(m,m,"P0P0")
+        mat=rem.getCrudeMatrix()
+        self.assertTrue(len(mat)==1)
+        self.assertTrue(len(mat[0])==1)
+        self.assertTrue(list(mat[0].keys())==[0])
+        res=list(mat[0].values())[0]
+        ref=float(m.getMeasureField(True).getArray())
+        self.assertTrue(abs(res-ref)/ref<1e-12)
+        pass
         
     def checkMatrix(self,mat1,mat2,nbCols,eps):
         self.assertEqual(len(mat1),len(mat2))