From 97b954cf5a7e526b142444a5252b9da544b957ad Mon Sep 17 00:00:00 2001
From: Adrien Bruneton <adrien.bruneton@cea.fr>
Date: Tue, 18 Feb 2014 10:22:05 +0100
Subject: [PATCH] Bug fix: cell orientation for NORM_QPOLY with two SEG3 (a
 circle) was not working.

---
 src/MEDCoupling/MEDCouplingUMesh.cxx          | 30 +++++++++++++++++--
 src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 21 +++++++++++++
 2 files changed, 49 insertions(+), 2 deletions(-)

diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx
index 29d78628e..fe5cb46f8 100644
--- a/src/MEDCoupling/MEDCouplingUMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.cxx
@@ -7863,17 +7863,43 @@ void MEDCouplingUMesh::AppendExtrudedCell(const int *connBg, const int *connEnd,
  */
 bool MEDCouplingUMesh::IsPolygonWellOriented(bool isQuadratic, const double *vec, const int *begin, const int *end, const double *coords)
 {
+  std::size_t i, ip1;
   double v[3]={0.,0.,0.};
   std::size_t sz=std::distance(begin,end);
   if(isQuadratic)
     sz/=2;
-  for(std::size_t i=0;i<sz;i++)
+  for(i=0;i<sz;i++)
     {
       v[0]+=coords[3*begin[i]+1]*coords[3*begin[(i+1)%sz]+2]-coords[3*begin[i]+2]*coords[3*begin[(i+1)%sz]+1];
       v[1]+=coords[3*begin[i]+2]*coords[3*begin[(i+1)%sz]]-coords[3*begin[i]]*coords[3*begin[(i+1)%sz]+2];
       v[2]+=coords[3*begin[i]]*coords[3*begin[(i+1)%sz]+1]-coords[3*begin[i]+1]*coords[3*begin[(i+1)%sz]];
     }
-  return vec[0]*v[0]+vec[1]*v[1]+vec[2]*v[2]>0.;
+  double ret = vec[0]*v[0]+vec[1]*v[1]+vec[2]*v[2];
+
+  // Try using quadratic points if standard points are degenerated (for example a QPOLYG with two
+  // SEG3 forming a circle):
+  if (fabs(ret) < INTERP_KERNEL::DEFAULT_ABS_TOL && isQuadratic)
+    {
+      v[0] = 0.0; v[1] = 0.0; v[2] = 0.0;
+      for(std::size_t j=0;j<sz;j++)
+        {
+          if (j%2)  // current point i is quadratic, next point i+1 is standard
+            {
+              i = sz+j;
+              ip1 = (j+1)%sz; // ip1 = "i+1"
+            }
+          else      // current point i is standard, next point i+1 is quadratic
+            {
+              i = j;
+              ip1 = j+sz;
+            }
+          v[0]+=coords[3*begin[i]+1]*coords[3*begin[ip1]+2]-coords[3*begin[i]+2]*coords[3*begin[ip1]+1];
+          v[1]+=coords[3*begin[i]+2]*coords[3*begin[ip1]]-coords[3*begin[i]]*coords[3*begin[ip1]+2];
+          v[2]+=coords[3*begin[i]]*coords[3*begin[ip1]+1]-coords[3*begin[i]+1]*coords[3*begin[ip1]];
+        }
+      ret = vec[0]*v[0]+vec[1]*v[1]+vec[2]*v[2];
+    }
+  return (ret>0.);
 }
 
 /*!
diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
index 734391fa9..95f19142f 100644
--- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
+++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
@@ -2183,6 +2183,27 @@ class MEDCouplingBasicsTest(unittest.TestCase):
             pass
         pass
 
+    def testCellOrientation3(self):
+        from cmath import rect  
+
+        c = [rect(1.0, i*pi/4.0) for i in range(8)]
+        coords = [c[-1].real,c[-1].imag,  c[3].real,c[3].imag,
+                   c[5].real,c[5].imag,  c[1].real,c[1].imag]
+        connec = [0,1,2,3] 
+        baseMesh = MEDCouplingUMesh.New("circle", 2)  
+        baseMesh.allocateCells(1)
+        meshCoords = DataArrayDouble.New(coords, 4, 2)
+        baseMesh.setCoords(meshCoords)
+        baseMesh.insertNextCell(NORM_QPOLYG, connec)  # a circle
+        baseMesh.finishInsertingCells()  
+        baseMesh.changeSpaceDimension(3)
+        Oz = [0.0, 0.0, -1.0] 
+        cell_lst = baseMesh.are2DCellsNotCorrectlyOriented(Oz, False)
+        self.assertEqual(cell_lst.getNumberOfTuples(), 0)
+        Oz[2] = 1.0
+        cell_lst = baseMesh.are2DCellsNotCorrectlyOriented(Oz, False)
+        self.assertEqual(cell_lst.getNumberOfTuples(), 1)
+
     def testPolyhedronBarycenter(self):
         connN=[0,3,2,1, -1, 4,5,6,7, -1, 0,4,7,3, -1, 3,7,6,2, -1, 2,6,5,1, -1, 1,5,4,0];
         coords=[0.,0.,0., 1.,0.,0., 1.,1.,0., 0.,1.,0., 0.,0.,1., 1.,0.,1., 1.,1.,1., 0.,1.,1., 0.5, 0.5, 0.5];
-- 
2.39.2