From 102fb5c72d77a6f01679b3fe6ad7d33c94de2a84 Mon Sep 17 00:00:00 2001
From: ageay <ageay>
Date: Tue, 26 Mar 2013 13:29:45 +0000
Subject: [PATCH] Barycenter on polyhedron with volume == 0...

---
 src/INTERP_KERNEL/VolSurfFormulae.hxx         | 81 ++++++++++++++-----
 src/MEDCoupling_Swig/MEDCouplingBasicsTest.py |  2 +
 2 files changed, 62 insertions(+), 21 deletions(-)

diff --git a/src/INTERP_KERNEL/VolSurfFormulae.hxx b/src/INTERP_KERNEL/VolSurfFormulae.hxx
index 1b2c56573..a96914a6e 100644
--- a/src/INTERP_KERNEL/VolSurfFormulae.hxx
+++ b/src/INTERP_KERNEL/VolSurfFormulae.hxx
@@ -501,31 +501,70 @@ namespace INTERP_KERNEL
     double area[3];
     areaVectorOfPolygon<ConnType,numPol>(connec,lgth,coords,area);
     double norm=sqrt(area[0]*area[0]+area[1]*area[1]+area[2]*area[2]);
-    area[0]/=norm; area[1]/=norm; area[2]/=norm;
-    res[0]=0.; res[1]=0.; res[2]=0.;
-    for(int i=1;i<lgth-1;i++)
+    if(norm>std::numeric_limits<double>::min())
+      {
+        area[0]/=norm; area[1]/=norm; area[2]/=norm;
+        res[0]=0.; res[1]=0.; res[2]=0.;
+        for(int i=1;i<lgth-1;i++)
+          {
+            double v[3];
+            double tmpArea[3];
+            v[0]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[0])]+
+                  coords[3*OTT<ConnType,numPol>::coo2C(connec[i])]+
+                  coords[3*OTT<ConnType,numPol>::coo2C(connec[i+1])])/3.;
+            v[1]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[0])+1]+
+                  coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+1]+
+                  coords[3*OTT<ConnType,numPol>::coo2C(connec[i+1])+1])/3.;
+            v[2]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[0])+2]+
+                  coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+2]+
+                  coords[3*OTT<ConnType,numPol>::coo2C(connec[i+1])+2])/3.;
+            ConnType tmpConn[3]={connec[0],connec[i],connec[i+1]};
+            areaVectorOfPolygon<ConnType,numPol>(tmpConn,3,coords,tmpArea);
+            double norm2=sqrt(tmpArea[0]*tmpArea[0]+tmpArea[1]*tmpArea[1]+tmpArea[2]*tmpArea[2]);
+            if(norm2>1e-12)
+              {
+                tmpArea[0]/=norm2; tmpArea[1]/=norm2; tmpArea[2]/=norm2;
+                double signOfArea=area[0]*tmpArea[0]+area[1]*tmpArea[1]+area[2]*tmpArea[2];
+                res[0]+=signOfArea*norm2*v[0]/norm; res[1]+=signOfArea*norm2*v[1]/norm; res[2]+=signOfArea*norm2*v[2]/norm;
+              }
+          }
+      }
+    else
       {
+        res[0]=0.; res[1]=0.; res[2]=0.;
+        if(lgth<1)
+          throw INTERP_KERNEL::Exception("computePolygonBarycenter3D : lgth of polygon is < 1 !");
+        norm=0.;
         double v[3];
-        double tmpArea[3];
-        v[0]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[0])]+
-              coords[3*OTT<ConnType,numPol>::coo2C(connec[i])]+
-              coords[3*OTT<ConnType,numPol>::coo2C(connec[i+1])])/3.;
-        v[1]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[0])+1]+
-              coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+1]+
-              coords[3*OTT<ConnType,numPol>::coo2C(connec[i+1])+1])/3.;
-        v[2]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[0])+2]+
-              coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+2]+
-              coords[3*OTT<ConnType,numPol>::coo2C(connec[i+1])+2])/3.;
-        ConnType tmpConn[3]={connec[0],connec[i],connec[i+1]};
-        areaVectorOfPolygon<ConnType,numPol>(tmpConn,3,coords,tmpArea);
-        double norm2=sqrt(tmpArea[0]*tmpArea[0]+tmpArea[1]*tmpArea[1]+tmpArea[2]*tmpArea[2]);
-        if(norm2>1e-12)
+        for(int i=0;i<lgth;i++)
+          {
+            v[0]=coords[3*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])]-coords[3*OTT<ConnType,numPol>::coo2C(connec[i])];
+            v[1]=coords[3*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])+1]-coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+1];
+            v[2]=coords[3*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])+2]-coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+2];
+            double norm2=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
+            res[0]+=(coords[3*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])]+coords[3*OTT<ConnType,numPol>::coo2C(connec[i])])/2.*norm2;
+            res[1]+=(coords[3*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])+1]+coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+1])/2.*norm2;
+            res[2]+=(coords[3*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])+2]+coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+2])/2.*norm2;
+            norm+=norm2;
+          }
+        if(norm>std::numeric_limits<double>::min())
           {
-            tmpArea[0]/=norm2; tmpArea[1]/=norm2; tmpArea[2]/=norm2;
-            double signOfArea=area[0]*tmpArea[0]+area[1]*tmpArea[1]+area[2]*tmpArea[2];
-            res[0]+=signOfArea*norm2*v[0]/norm; res[1]+=signOfArea*norm2*v[1]/norm; res[2]+=signOfArea*norm2*v[2]/norm;
+            res[0]/=norm; res[1]/=norm; res[2]/=norm;
+            return;
           }
-      }   
+        else
+          {
+            res[0]=0.; res[1]=0.; res[2]=0.;
+            for(int i=0;i<lgth;i++)
+              {
+                res[0]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[i])];
+                res[1]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+1];
+                res[2]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+2];
+              }
+            res[0]/=lgth; res[1]/=lgth; res[2]/=lgth;
+            return;
+          }
+      }
   }
 
   inline double integrationOverA3DLine(double u1, double v1, double u2, double v2, double A, double B, double C)
diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
index eec0bc36a..29a1816ae 100644
--- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
+++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
@@ -11778,6 +11778,8 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         m3=m2.buildExtrudedMesh(m1,0)
         m3.insertNextCell(NORM_POLYHED,[3,4,5,-1,8,7,6,-1,4,3,6,7,-1,5,4,7,8,-1,5,4,-1,3,5,8,6])# addition of face #4 with null surface
         self.assertTrue(m3.getBarycenterAndOwner().isEqual(DataArrayDouble([0.3333333333333333,0.3333333333333333,0.,0.3333333333333333,0.3333333333333333,1.,0.3333333333333333,0.3333333333333333,1.],3,3),1e-13))
+        m4,a,b,c,d=m3.buildDescendingConnectivity()
+        self.assertTrue(m4.getBarycenterAndOwner().isEqual(DataArrayDouble([0.3333333333333333,0.3333333333333333,0.,0.3333333333333333,0.3333333333333333,0.,0.5,0.,0.,0.5,0.5,0.,0.,0.5,0.,0.3333333333333333,0.3333333333333333,2.,0.5,0.,1.,0.5,0.5,1.,0.,0.5,1.,0.5,0.5,0.],10,3),1e-13))
         pass
 
     def setUp(self):
-- 
2.39.2