Barycenter on polyhedron with volume == 0...

diff --git a/src/INTERP_KERNEL/VolSurfFormulae.hxx b/src/INTERP_KERNEL/VolSurfFormulae.hxx
index c5973706805361c56cf960d17a2f1c4e1c007b72..1b2c5657359b1be9db30b5369b0f35666e426e98 100644 (file)
--- a/src/INTERP_KERNEL/VolSurfFormulae.hxx
+++ b/src/INTERP_KERNEL/VolSurfFormulae.hxx
@@ -495,6 +495,39 @@ namespace INTERP_KERNEL
       }
   }
 
+  template<class ConnType, NumberingPolicy numPol>
+  inline void computePolygonBarycenter3D(const ConnType *connec, int lgth, const double *coords, double *res)
+  {
+    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++)
+      {
+        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;
+          }
+      }   
+  }
+
   inline double integrationOverA3DLine(double u1, double v1, double u2, double v2, double A, double B, double C)
   {
     return (u1-u2)*(6.*C*C*(v1+v2)+B*B*(v1*v1*v1+v1*v1*v2+v1*v2*v2+v2*v2*v2)+A*A*(2.*u1*u2*(v1+v2)+u1*u1*(3.*v1+v2)+u2*u2*(v1+3.*v2))+ 
@@ -515,6 +548,8 @@ namespace INTERP_KERNEL
         double normal[3];
         areaVectorOfPolygon<ConnType,numPol>(work,nbOfNodesOfCurFace,coords,normal);
         double normOfNormal=sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
+        if(normOfNormal<std::numeric_limits<double>::min())
+          continue;
         normal[0]/=normOfNormal; normal[1]/=normOfNormal; normal[2]/=normOfNormal;
         double u[2]={normal[1],-normal[0]};
         double s=sqrt(u[0]*u[0]+u[1]*u[1]);
@@ -561,7 +596,32 @@ namespace INTERP_KERNEL
         work=work2+1;
       }
     double vol=calculateVolumeForPolyh2<ConnType,numPol>(connec,lgth,coords);
-    res[0]/=vol; res[1]/=vol; res[2]/=vol;
+    if(fabs(vol)>std::numeric_limits<double>::min())
+      {
+        res[0]/=vol; res[1]/=vol; res[2]/=vol;
+      }
+    else
+      {
+        double sum=0.;
+        res[0]=0.; res[1]=0.; res[2]=0.;
+        work=connec;
+        for(std::size_t i=0;i<nbOfFaces;i++)
+          {
+            const int *work2=std::find(work+1,connec+lgth,-1);
+            int nbOfNodesOfCurFace=(int)std::distance(work,work2);
+            double normal[3];
+            areaVectorOfPolygon<ConnType,numPol>(work,nbOfNodesOfCurFace,coords,normal);
+            double normOfNormal=sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
+            if(normOfNormal<std::numeric_limits<double>::min())
+              continue;
+            sum+=normOfNormal;
+            double tmpBary[3];
+            computePolygonBarycenter3D<ConnType,numPol>(work,nbOfNodesOfCurFace,coords,tmpBary);
+            res[0]+=normOfNormal*tmpBary[0]; res[1]+=normOfNormal*tmpBary[1]; res[2]+=normOfNormal*tmpBary[2];
+            work=work2+1;
+          }
+        res[0]/=sum; res[1]/=sum; res[2]/=sum;
+      }
   }
 
   // ============================================================================================================================================
@@ -686,39 +746,6 @@ namespace INTERP_KERNEL
     res[0]/=3.*area;
     res[1]/=3.*area;
   }
-
-  template<class ConnType, NumberingPolicy numPol>
-  inline void computePolygonBarycenter3D(const ConnType *connec, int lgth, const double *coords, double *res)
-  {
-    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++)
-      {
-        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;
-          }
-      }   
-  }
 }
 
 #endif

diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
index e20306b2ad8a602047717d0ffae99c97b808c929..eec0bc36a8d387c4720a54642e681a6b1d21ddd5 100644 (file)
--- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
+++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
@@ -11714,7 +11714,7 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         self.assertTrue(d.isEqual(DataArrayInt([6,16,5,15,4,14,3,13,2,12,1,11,0,10],7,2)))
         pass
 
-    def testSwig2DAPow1(self):
+    def testSwigDAPow1(self):
         d=DataArrayInt(10)
         d.iota(0)
         d1=d.deepCpy()
@@ -11758,6 +11758,27 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         d4**=d3
         self.assertTrue(d4.isEqual(DataArrayDouble([1.,sqrt(2.),1.4422495703074083,sqrt(2.)]),1e-14))
         pass
+    
+    def testSwig2Baryenter3DForCellsWithVolumeZero1(self):
+        coo=DataArrayDouble([0.,0.,0.,1.,0.,0.,0.,1.,0.],3,3)
+        m2=MEDCouplingUMesh("mesh",2)
+        m2.allocateCells(0)
+        m2.insertNextCell(NORM_POLYGON,[0,1,2])
+        m2.setCoords(coo)
+        m2.checkCoherency1()
+        #
+        coo2=DataArrayDouble([0.,0.,0.,0.,0.,0.,0.,0.,2.],3,3)
+        m1=MEDCouplingUMesh("mesh",1)
+        m1.allocateCells(0)
+        m1.insertNextCell(NORM_SEG2,[0,1])
+        m1.insertNextCell(NORM_SEG2,[1,2])
+        m1.setCoords(coo2)
+        m1.checkCoherency1()
+        #
+        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))
+        pass
 
     def setUp(self):
         pass