[EDF10718] : Diameter computation of non linear cells.

[tools/medcoupling.git] / src / INTERP_KERNEL / VolSurfFormulae.hxx

diff --git a/src/INTERP_KERNEL/VolSurfFormulae.hxx b/src/INTERP_KERNEL/VolSurfFormulae.hxx
index 6e68ad4d1528265a1e25572c92de025204eddd80..5562de5a318aac55e40f67b9f85f5e24e47a9865 100644 (file)
--- a/src/INTERP_KERNEL/VolSurfFormulae.hxx
+++ b/src/INTERP_KERNEL/VolSurfFormulae.hxx
@@ -1,9 +1,9 @@
-// Copyright (C) 2007-2013  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2015  CEA/DEN, EDF R&D
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 // License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// version 2.1 of the License, or (at your option) any later version.
 //
 // This library is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
@@ -23,6 +23,8 @@
 
 #include "InterpolationUtils.hxx"
 #include "InterpKernelException.hxx"
+#include "InterpKernelGeo2DEdgeLin.hxx"
+#include "InterpKernelGeo2DEdgeArcCircle.hxx"
 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
 
 #include <sstream>
@@ -52,6 +54,18 @@ namespace INTERP_KERNEL
         return sqrt(ret);
       }
   }
+  
+  inline double calculateLgthForSeg3(const double *begin, const double *end, const double *middle, int spaceDim)
+  {
+    if(spaceDim==2)
+      {
+        Edge *ed=Edge::BuildEdgeFrom3Points(begin,middle,end);
+        double ret=ed->getCurveLength(); ed->decrRef();
+        return ret;
+        }
+    else
+      return calculateLgthForSeg2(begin,end,spaceDim);
+  }
 
   // ===========================
   // Calculate Area for triangle
@@ -707,37 +721,37 @@ namespace INTERP_KERNEL
   }
 
   template<>
-  inline void calculateBarycenter<2,0>(const double **/*pts*/, double */*bary*/)
+  inline void calculateBarycenter<2,0>(const double ** /*pts*/, double * /*bary*/)
   {
   }
   
   template<>
-  inline void calculateBarycenter<3,0>(const double **/*pts*/, double */*bary*/)
+  inline void calculateBarycenter<3,0>(const double ** /*pts*/, double * /*bary*/)
   {
   }
   
   template<>
-  inline void calculateBarycenter<4,0>(const double **/*pts*/, double */*bary*/)
+  inline void calculateBarycenter<4,0>(const double ** /*pts*/, double * /*bary*/)
   {
   }
   
   template<>
-  inline void calculateBarycenter<5,0>(const double **/*pts*/, double */*bary*/)
+  inline void calculateBarycenter<5,0>(const double ** /*pts*/, double * /*bary*/)
   {
   }
   
   template<>
-  inline void calculateBarycenter<6,0>(const double **/*pts*/, double */*bary*/)
+  inline void calculateBarycenter<6,0>(const double ** /*pts*/, double * /*bary*/)
   {
   }
   
   template<>
-  inline void calculateBarycenter<7,0>(const double **/*pts*/, double */*bary*/)
+  inline void calculateBarycenter<7,0>(const double ** /*pts*/, double * /*bary*/)
   {
   }
   
   template<>
-  inline void calculateBarycenter<8,0>(const double **/*pts*/, double */*bary*/)
+  inline void calculateBarycenter<8,0>(const double ** /*pts*/, double * /*bary*/)
   {
   }
 
@@ -768,23 +782,54 @@ namespace INTERP_KERNEL
         bary[i]=temp/nbPts;
       }
   }
-  
-  template<class ConnType, NumberingPolicy numPol>
-  inline void computePolygonBarycenter2D(const ConnType *connec, int lgth, const double *coords, double *res)
+
+  inline void computePolygonBarycenter2DEngine(double **coords, int lgth, double *res)
   {
     double area=0.;
     res[0]=0.; res[1]=0.;
     for(int i=0;i<lgth;i++)
       {
-        double cp=coords[2*OTT<ConnType,numPol>::coo2C(connec[i])]*coords[2*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])+1]-
-          coords[2*OTT<ConnType,numPol>::coo2C(connec[i])+1]*coords[2*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])];
+        double cp=coords[i][0]*coords[(i+1)%lgth][1]-coords[i][1]*coords[(i+1)%lgth][0];
         area+=cp;
-        res[0]+=cp*(coords[2*OTT<ConnType,numPol>::coo2C(connec[i])]+coords[2*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])]);
-        res[1]+=cp*(coords[2*OTT<ConnType,numPol>::coo2C(connec[i])+1]+coords[2*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])+1]);
+        res[0]+=cp*(coords[i][0]+coords[(i+1)%lgth][0]);
+        res[1]+=cp*(coords[i][1]+coords[(i+1)%lgth][1]);
       }
     res[0]/=3.*area;
     res[1]/=3.*area;
   }
+
+  template<class ConnType, NumberingPolicy numPol>
+  inline void computePolygonBarycenter2D(const ConnType *connec, int lgth, const double *coords, double *res)
+  {
+    double **coords2=new double *[lgth];
+    for(int i=0;i<lgth;i++)
+      coords2[i]=const_cast<double *>(coords+2*OTT<ConnType,numPol>::coo2C(connec[i]));
+    computePolygonBarycenter2DEngine(coords2,lgth,res);
+    delete [] coords2;
+  }
+  
+  inline void computeQPolygonBarycenter2D(double **coords, int nbOfPtsInPolygs, int spaceDim, double *res)
+  {
+    if(nbOfPtsInPolygs%2==0)
+      {
+        if(spaceDim==2)
+          {
+            std::vector<Node *> nodes(nbOfPtsInPolygs);
+            for(int i=0;i<nbOfPtsInPolygs;i++)
+              nodes[i]=new Node(coords[i][0],coords[i][1]);
+            QuadraticPolygon *pol=QuadraticPolygon::BuildArcCirclePolygon(nodes);
+            pol->getBarycenter(res);
+            delete pol;
+          }
+        else
+          return computePolygonBarycenter2DEngine(coords,nbOfPtsInPolygs/2,res);
+      }
+    else
+      {
+        std::ostringstream oss; oss << "INTERP_KERNEL::computeQPolygonBarycenter2D : nb of points in quadratic polygon is " << nbOfPtsInPolygs << " should be even !";
+        throw INTERP_KERNEL::Exception(oss.str().c_str());
+      }
+  }
 }
 
 #endif