From 3d767bbe20c2cc633425c864a713fb1d6bed5c29 Mon Sep 17 00:00:00 2001
From: abn <adrien.bruneton@cea.fr>
Date: Wed, 16 Jan 2019 16:12:00 +0100
Subject: [PATCH] orientCorreclty2DCells() -> bug fix: quadratic cell
 orientation not well handled.

---
 src/MEDCoupling/MEDCouplingUMesh.cxx          | 15 +++++++-----
 .../MEDCouplingBasicsTest1.py                 | 24 +++++++++++++++++++
 2 files changed, 33 insertions(+), 6 deletions(-)

diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx
index cce5fe0bf..d447ef6da 100644
--- a/src/MEDCoupling/MEDCouplingUMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.cxx
@@ -7010,26 +7010,29 @@ bool MEDCouplingUMesh::IsPolygonWellOriented(bool isQuadratic, const double *vec
   if(!isQuadratic)
     for(i=0;i<sz;i++)
       {
+        // Algorithm: sum in v the cross products of (e1, e2) where e_i it the vector between (0,0,0) and point i
+        // and e2 is linear point directly following e1 in the connectivity. All points are used.
         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]];
       }
   else
     {
-      // Use all points if quadratic (taking only linear points might lead to issues if the linearized version of the
+      // Same algorithm as above but also using intermediate quadratic points.
+      // (taking only linear points might lead to issues if the linearized version of the
       // polygon is not convex or self-intersecting ... see testCellOrientation4)
-      sz /= 2;
+      int hsz = sz/2;
       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"
+              i = hsz+(j-1)/2;
+              ip1 = ((j-1)/2 + 1)%hsz; // ip1 means "i+1", i.e. next point
             }
           else      // current point i is standard, next point i+1 is quadratic
             {
-              i = j;
-              ip1 = j+sz;
+              i = j/2;
+              ip1 = j/2+hsz;
             }
           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];
diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest1.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest1.py
index 92298fe33..0459b0376 100644
--- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest1.py
+++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest1.py
@@ -2287,6 +2287,30 @@ class MEDCouplingBasicsTest1(unittest.TestCase):
         self.assertEqual(mesh.getNodalConnectivityIndex().getValues(), cIRef)
         pass
 
+    def testCellOrientation5(self):
+        """ Non regression for NORM_QPOLYG  """
+        mesh = MEDCouplingUMesh('Mesh_3', 2)
+        coo = DataArrayDouble([(-34.3035,5.1),(-35.2018,4.59163),(-34.9509,6.21985),(-35.0858,5.4072),(-34.7527,4.84582),(-34.6641,5.63857)])
+        mesh.setCoords(coo)
+        c = DataArrayInt([6, 2, 1, 0, 3, 4, 5])
+        cI = DataArrayInt([0, 7])
+        mesh.setConnectivity(c, cI)
+        vec = [0., 0., -1.]
+        mesh.changeSpaceDimension(3)
+        mesh.orientCorrectly2DCells(vec, False)
+        mesh.changeSpaceDimension(2)
+        cRef = [6, 2, 0, 1, 5, 4, 3]
+        cIRef = [0, 7]
+        self.assertEqual(mesh.getNodalConnectivity().getValues(), cRef)
+        self.assertEqual(mesh.getNodalConnectivityIndex().getValues(), cIRef)
+        # Second call doest not change anything:
+        mesh.changeSpaceDimension(3)
+        mesh.orientCorrectly2DCells(vec, False)
+        mesh.changeSpaceDimension(2)
+        self.assertEqual(mesh.getNodalConnectivity().getValues(), cRef)
+        self.assertEqual(mesh.getNodalConnectivityIndex().getValues(), cIRef)
+        pass
+
     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