src/INTERP_KERNEL/VolSurfFormulae.hxx

   1 // Copyright (C) 2007-2024  CEA, EDF
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #ifndef __VOLSURFFORMULAE_HXX__
  22 #define __VOLSURFFORMULAE_HXX__
  23
  24 #include "InterpolationUtils.hxx"
  25 #include "InterpKernelException.hxx"
  26 #include "InterpKernelGeo2DEdgeLin.hxx"
  27 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
  28 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
  29 #include "MCIdType.hxx"
  30
  31 #include <sstream>
  32 #include <cmath>
  33
  34 namespace INTERP_KERNEL
  35 {
  36   inline void calculateBarycenterDyn(const double **pts, mcIdType nbPts,
  37                                      int dim, double *bary);
  38
  39   inline double calculateAreaForPolyg(const double **coords, mcIdType nbOfPtsInPolygs,
  40                                       int spaceDim);
  41
  42
  43   inline double calculateAreaForQPolyg(const double **coords, mcIdType nbOfPtsInPolygs,
  44                                        int spaceDim);
  45
  46   inline double calculateLgthForSeg2(const double *p1, const double *p2, int spaceDim)
  47   {
  48     if(spaceDim==1)
  49       return *p2-*p1;
  50     else
  51       {
  52         double ret=0;
  53         for(int i=0;i<spaceDim;i++)
  54           ret+=(p2[i]-p1[i])*(p2[i]-p1[i]);
  55         return sqrt(ret);
  56       }
  57   }
  58
  59   inline double calculateLgthForSeg3(const double *begin, const double *end, const double *middle, int spaceDim)
  60   {
  61     if(spaceDim==2)
  62       {
  63         Edge *ed=Edge::BuildEdgeFrom3Points(begin,middle,end);
  64         double ret=ed->getCurveLength(); ed->decrRef();
  65         return ret;
  66         }
  67     else
  68       return calculateLgthForSeg2(begin,end,spaceDim);
  69   }
  70
  71   // ===========================
  72   // Calculate Area for triangle
  73   // ===========================
  74   inline double calculateAreaForTria(const double *p1, const double *p2,
  75                                      const double *p3, int spaceDim)
  76   {
  77     double area ;
  78
  79     if ( spaceDim == 2 )
  80       {
  81         area = -((p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1]))/2.0;
  82       }
  83     else
  84       {
  85         area = sqrt(((p2[1]-p1[1])*(p3[2]-p1[2]) - (p3[1]-p1[1])*(p2[2]-p1[2]))*
  86                     ((p2[1]-p1[1])*(p3[2]-p1[2]) - (p3[1]-p1[1])*(p2[2]-p1[2]))
  87                     +
  88                     ((p3[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p3[2]-p1[2]))*
  89                     ((p3[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p3[2]-p1[2]))
  90                     +
  91                     ((p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1]))*
  92                     ((p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1])))/2.0;
  93       }
  94
  95     return area ;
  96   }
  97
  98   // =============================
  99   // Calculate Area for quadrangle
 100   // =============================
 101   inline double calculateAreaForQuad(const double *p1, const double *p2,
 102                                      const double *p3, const double *p4,
 103                                      int spaceDim)
 104   {
 105     double area ;
 106
 107     if (spaceDim==2)
 108       {
 109         double a1 = (p2[0]-p1[0])/4.0, a2 = (p2[1]-p1[1])/4.0;
 110         double b1 = (p3[0]-p4[0])/4.0, b2 = (p3[1]-p4[1])/4.0;
 111         double c1 = (p3[0]-p2[0])/4.0, c2 = (p3[1]-p2[1])/4.0;
 112         double d1 = (p4[0]-p1[0])/4.0, d2 = (p4[1]-p1[1])/4.0;
 113
 114         area = - 4.0*(  b1*c2 - c1*b2 + a1*c2 - c1*a2
 115                         + b1*d2 - d1*b2 + a1*d2 - d1*a2);
 116       }
 117     else
 118       {
 119         area = (sqrt(((p2[1]-p1[1])*(p4[2]-p1[2]) - (p4[1]-p1[1])*(p2[2]-p1[2]))*
 120                      ((p2[1]-p1[1])*(p4[2]-p1[2]) - (p4[1]-p1[1])*(p2[2]-p1[2]))
 121                      + ((p4[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p4[2]-p1[2]))*
 122                      ((p4[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p4[2]-p1[2]))
 123                      + ((p2[0]-p1[0])*(p4[1]-p1[1]) - (p4[0]-p1[0])*(p2[1]-p1[1]))*
 124                      ((p2[0]-p1[0])*(p4[1]-p1[1]) - (p4[0]-p1[0])*(p2[1]-p1[1])))
 125                 +
 126                 sqrt(((p4[1]-p3[1])*(p2[2]-p3[2]) - (p2[1]-p3[1])*(p4[2]-p3[2]))*
 127                      ((p4[1]-p3[1])*(p2[2]-p3[2]) - (p2[1]-p3[1])*(p4[2]-p3[2]))
 128                      + ((p2[0]-p3[0])*(p4[2]-p3[2]) - (p4[0]-p3[0])*(p2[2]-p3[2]))*
 129                      ((p2[0]-p3[0])*(p4[2]-p3[2]) - (p4[0]-p3[0])*(p2[2]-p3[2]))
 130                      + ((p4[0]-p3[0])*(p2[1]-p3[1]) - (p2[0]-p3[0])*(p4[1]-p3[1]))*
 131                      ((p4[0]-p3[0])*(p2[1]-p3[1]) - (p2[0]-p3[0])*(p4[1]-p3[1])))
 132                 )/2.0;
 133       }
 134
 135     return area ;
 136   }
 137
 138   // ====================================
 139   // Calculate Normal Vector for Triangle
 140   // ====================================
 141   inline void calculateNormalForTria(const double *p1, const double *p2,
 142                                      const double *p3, double *normal)
 143   {
 144     normal[0] = ((p2[1]-p1[1])*(p3[2]-p1[2]) - (p3[1]-p1[1])*(p2[2]-p1[2]))/2.0;
 145     normal[1] = ((p3[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p3[2]-p1[2]))/2.0;
 146     normal[2] = ((p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1]))/2.0;
 147   }
 148
 149   // ======================================
 150   // Calculate Normal Vector for Quadrangle
 151   // ======================================
 152   inline void calculateNormalForQuad(const double *p1, const double *p2,
 153                                      const double *p3, const double *p4,
 154                                      double *normal)
 155   {
 156     double xnormal1 = (p2[1]-p1[1])*(p4[2]-p1[2]) - (p4[1]-p1[1])*(p2[2]-p1[2]);
 157     double xnormal2 = (p4[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p4[2]-p1[2]);
 158     double xnormal3 = (p2[0]-p1[0])*(p4[1]-p1[1]) - (p4[0]-p1[0])*(p2[1]-p1[1]);
 159     double xarea = sqrt(xnormal1*xnormal1 + xnormal2*xnormal2 + xnormal3*xnormal3);
 160     xnormal1 = xnormal1/xarea;
 161     xnormal2 = xnormal2/xarea;
 162     xnormal3 = xnormal3/xarea;
 163     xarea = calculateAreaForQuad(p1,p2,p3,p4,3);
 164     normal[0] = xnormal1*xarea ;
 165     normal[1] = xnormal2*xarea ;
 166     normal[2] = xnormal3*xarea ;
 167   }
 168
 169   // ===================================
 170   // Calculate Normal Vector for Polygon
 171   // ===================================
 172   inline void calculateNormalForPolyg(const double **coords, mcIdType nbOfPtsInPolygs,
 173                                       double *normal)
 174   {
 175     double coordOfBary[3];
 176
 177     calculateBarycenterDyn(coords,nbOfPtsInPolygs,3,coordOfBary);
 178     double xnormal1 = (coords[0][1]-coords[1][1]) * (coordOfBary[2]-coords[1][2])
 179       - (coords[0][2]-coords[1][2]) * (coordOfBary[1]-coords[1][1]);
 180
 181     double xnormal2 = (coords[0][2]-coords[1][2]) * (coordOfBary[0]-coords[1][0])
 182       - (coords[0][0]-coords[1][0]) * (coordOfBary[2]-coords[1][2]);
 183
 184     double xnormal3 = (coords[0][0]-coords[1][0]) * (coordOfBary[1]-coords[1][1])
 185       - (coords[0][1]-coords[1][1]) * (coordOfBary[0]-coords[1][0]);
 186
 187     double xarea=sqrt(xnormal1*xnormal1 + xnormal2*xnormal2 + xnormal3*xnormal3);
 188
 189     if ( xarea < 1.e-6 )
 190       {
 191         //std::string diagnosis"Can not calculate normal - polygon is singular";
 192         throw std::exception();
 193       }
 194
 195     xnormal1 = -xnormal1/xarea;
 196     xnormal2 = -xnormal2/xarea;
 197     xnormal3 = -xnormal3/xarea;
 198     xarea = calculateAreaForPolyg(coords,nbOfPtsInPolygs,3);
 199     normal[0] = xnormal1*xarea ;
 200     normal[1] = xnormal2*xarea ;
 201     normal[2] = xnormal3*xarea ;
 202   }
 203
 204   // ==========================
 205   // Calculate Area for Polygon
 206   // ==========================
 207   inline double calculateAreaForPolyg(const double **coords, mcIdType nbOfPtsInPolygs,
 208                                       int spaceDim)
 209   {
 210     double ret=0.;
 211     double coordOfBary[3];
 212
 213     calculateBarycenterDyn(coords,nbOfPtsInPolygs,spaceDim,coordOfBary);
 214     for ( mcIdType i=0; i<nbOfPtsInPolygs; i++ )
 215       {
 216         double tmp = calculateAreaForTria(coords[i],coords[(i+1)%nbOfPtsInPolygs],
 217                                           coordOfBary,spaceDim);
 218         ret+=tmp;
 219       }
 220     return ret;
 221   }
 222
 223   double calculateAreaForQPolyg(const double **coords, mcIdType nbOfPtsInPolygs, int spaceDim)
 224   {
 225
 226     if(nbOfPtsInPolygs%2==0)
 227       {
 228         if(spaceDim==2)
 229           {
 230             std::vector<Node *> nodes(nbOfPtsInPolygs);
 231             for(mcIdType i=0;i<nbOfPtsInPolygs;i++)
 232               nodes[i]=new Node(coords[i][0],coords[i][1]);
 233             QuadraticPolygon *pol=QuadraticPolygon::BuildArcCirclePolygon(nodes);
 234             double ret=pol->getArea();
 235             delete pol;
 236             return -ret;
 237           }
 238         else
 239           return calculateAreaForPolyg(coords,nbOfPtsInPolygs/2,spaceDim);
 240       }
 241     else
 242       {
 243         std::ostringstream oss; oss << "INTERP_KERNEL::calculateAreaForQPolyg : nb of points in quadratic polygon is " << nbOfPtsInPolygs << " should be even !";
 244         throw INTERP_KERNEL::Exception(oss.str().c_str());
 245       }
 246   }
 247
 248   // ==========================
 249   // Calculate Volume for Tetra
 250   // ==========================
 251   inline double calculateVolumeForTetra(const double *p1, const double *p2,
 252                                         const double *p3, const double *p4)
 253   {
 254     return (  (p3[0]-p1[0])*(  (p2[1]-p1[1])*(p4[2]-p1[2])
 255                                - (p2[2]-p1[2])*(p4[1]-p1[1]) )
 256               - (p2[0]-p1[0])*(  (p3[1]-p1[1])*(p4[2]-p1[2])
 257                                  - (p3[2]-p1[2])*(p4[1]-p1[1]) )
 258               + (p4[0]-p1[0])*(  (p3[1]-p1[1])*(p2[2]-p1[2])
 259                                  - (p3[2]-p1[2])*(p2[1]-p1[1]) )
 260               ) / 6.0 ;
 261   }
 262
 263   // =========================
 264   // Calculate Volume for Pyra
 265   // =========================
 266   inline double calculateVolumeForPyra(const double *p1, const double *p2,
 267                                        const double *p3, const double *p4,
 268                                        const double *p5)
 269   {
 270     return ( ((p3[0]-p1[0])*(  (p2[1]-p1[1])*(p5[2]-p1[2])
 271                                - (p2[2]-p1[2])*(p5[1]-p1[1]) )
 272               -(p2[0]-p1[0])*(  (p3[1]-p1[1])*(p5[2]-p1[2])
 273                                 - (p3[2]-p1[2])*(p5[1]-p1[1]) )
 274               +(p5[0]-p1[0])*(  (p3[1]-p1[1])*(p2[2]-p1[2])
 275                                 - (p3[2]-p1[2])*(p2[1]-p1[1]) ))
 276              +
 277              ((p4[0]-p1[0])*(  (p3[1]-p1[1])*(p5[2]-p1[2])
 278                                - (p3[2]-p1[2])*(p5[1]-p1[1]) )
 279               -(p3[0]-p1[0])*(  (p4[1]-p1[1])*(p5[2]-p1[2])
 280                                 - (p4[2]-p1[2])*(p5[1]-p1[1]))
 281               +(p5[0]-p1[0])*(  (p4[1]-p1[1])*(p3[2]-p1[2])
 282                                 - (p4[2]-p1[2])*(p3[1]-p1[1]) ))
 283              ) / 6.0 ;
 284   }
 285
 286   // ==========================
 287   // Calculate Volume for Penta
 288   // ==========================
 289   inline double calculateVolumeForPenta(const double *p1, const double *p2,
 290                                         const double *p3, const double *p4,
 291                                         const double *p5, const double *p6)
 292   {
 293     double a1 = (p2[0]-p3[0])/2.0, a2 = (p2[1]-p3[1])/2.0, a3 = (p2[2]-p3[2])/2.0;
 294     double b1 = (p5[0]-p6[0])/2.0, b2 = (p5[1]-p6[1])/2.0, b3 = (p5[2]-p6[2])/2.0;
 295     double c1 = (p4[0]-p1[0])/2.0, c2 = (p4[1]-p1[1])/2.0, c3 = (p4[2]-p1[2])/2.0;
 296     double d1 = (p5[0]-p2[0])/2.0, d2 = (p5[1]-p2[1])/2.0, d3 = (p5[2]-p2[2])/2.0;
 297     double e1 = (p6[0]-p3[0])/2.0, e2 = (p6[1]-p3[1])/2.0, e3 = (p6[2]-p3[2])/2.0;
 298     double f1 = (p1[0]-p3[0])/2.0, f2 = (p1[1]-p3[1])/2.0, f3 = (p1[2]-p3[2])/2.0;
 299     double h1 = (p4[0]-p6[0])/2.0, h2 = (p4[1]-p6[1])/2.0, h3 = (p4[2]-p6[2])/2.0;
 300
 301     double A = a1*c2*f3 - a1*c3*f2 - a2*c1*f3 + a2*c3*f1 + a3*c1*f2 - a3*c2*f1;
 302     double B = b1*c2*h3 - b1*c3*h2 - b2*c1*h3 + b2*c3*h1 + b3*c1*h2 - b3*c2*h1;
 303     double C = (a1*c2*h3 + b1*c2*f3) - (a1*c3*h2 + b1*c3*f2)
 304       - (a2*c1*h3 + b2*c1*f3) + (a2*c3*h1 + b2*c3*f1)
 305       + (a3*c1*h2 + b3*c1*f2) - (a3*c2*h1 + b3*c2*f1);
 306     double D = a1*d2*f3 - a1*d3*f2 - a2*d1*f3 + a2*d3*f1 + a3*d1*f2 - a3*d2*f1;
 307     double E = b1*d2*h3 - b1*d3*h2 - b2*d1*h3 + b2*d3*h1 + b3*d1*h2 - b3*d2*h1;
 308     double F = (a1*d2*h3 + b1*d2*f3) - (a1*d3*h2 + b1*d3*f2)
 309       - (a2*d1*h3 + b2*d1*f3) + (a2*d3*h1 + b2*d3*f1)
 310       + (a3*d1*h2 + b3*d1*f2) - (a3*d2*h1 + b3*d2*f1);
 311     double G = a1*e2*f3 - a1*e3*f2 - a2*e1*f3 + a2*e3*f1 + a3*e1*f2 - a3*e2*f1;
 312     double H = b1*e2*h3 - b1*e3*h2 - b2*e1*h3 + b2*e3*h1 + b3*e1*h2 - b3*e2*h1;
 313     double P = (a1*e2*h3 + b1*e2*f3) - (a1*e3*h2 + b1*e3*f2)
 314       - (a2*e1*h3 + b2*e1*f3) + (a2*e3*h1 + b2*e3*f1)
 315       + (a3*e1*h2 + b3*e1*f2) - (a3*e2*h1 + b3*e2*f1);
 316
 317     return (-2.0*(2.0*(A + B + D + E + G + H) + C + F + P)/9.0);
 318   }
 319
 320   // =========================
 321   // Calculate Volume for Hexa
 322   // =========================
 323   inline double calculateVolumeForHexa(const double *pt1, const double *pt2,
 324                                        const double *pt3, const double *pt4,
 325                                        const double *pt5, const double *pt6,
 326                                        const double *pt7, const double *pt8)
 327   {
 328     double a1=(pt3[0]-pt4[0])/8.0, a2=(pt3[1]-pt4[1])/8.0, a3=(pt3[2]-pt4[2])/8.0;
 329     double b1=(pt2[0]-pt1[0])/8.0, b2=(pt2[1]-pt1[1])/8.0, b3=(pt2[2]-pt1[2])/8.0;
 330     double c1=(pt7[0]-pt8[0])/8.0, c2=(pt7[1]-pt8[1])/8.0, c3=(pt7[2]-pt8[2])/8.0;
 331     double d1=(pt6[0]-pt5[0])/8.0, d2=(pt6[1]-pt5[1])/8.0, d3=(pt6[2]-pt5[2])/8.0;
 332     double e1=(pt3[0]-pt2[0])/8.0, e2=(pt3[1]-pt2[1])/8.0, e3=(pt3[2]-pt2[2])/8.0;
 333     double f1=(pt4[0]-pt1[0])/8.0, f2=(pt4[1]-pt1[1])/8.0, f3=(pt4[2]-pt1[2])/8.0;
 334     double h1=(pt7[0]-pt6[0])/8.0, h2=(pt7[1]-pt6[1])/8.0, h3=(pt7[2]-pt6[2])/8.0;
 335     double p1=(pt8[0]-pt5[0])/8.0, p2=(pt8[1]-pt5[1])/8.0, p3=(pt8[2]-pt5[2])/8.0;
 336     double q1=(pt3[0]-pt7[0])/8.0, q2=(pt3[1]-pt7[1])/8.0, q3=(pt3[2]-pt7[2])/8.0;
 337     double r1=(pt4[0]-pt8[0])/8.0, r2=(pt4[1]-pt8[1])/8.0, r3=(pt4[2]-pt8[2])/8.0;
 338     double s1=(pt2[0]-pt6[0])/8.0, s2=(pt2[1]-pt6[1])/8.0, s3=(pt2[2]-pt6[2])/8.0;
 339     double t1=(pt1[0]-pt5[0])/8.0, t2=(pt1[1]-pt5[1])/8.0, t3=(pt1[2]-pt5[2])/8.0;
 340
 341     double A = a1*e2*q3 - a1*e3*q2 - a2*e1*q3 + a2*e3*q1 + a3*e1*q2 - a3*e2*q1;
 342     double B = c1*h2*q3 - c1*h3*q2 - c2*h1*q3 + c2*h3*q1 + c3*h1*q2 - c3*h2*q1;
 343     double C = (a1*h2 + c1*e2)*q3 - (a1*h3 + c1*e3)*q2
 344       - (a2*h1 + c2*e1)*q3 + (a2*h3 + c2*e3)*q1
 345       + (a3*h1 + c3*e1)*q2 - (a3*h2 + c3*e2)*q1;
 346     double D = b1*e2*s3 - b1*e3*s2 - b2*e1*s3 + b2*e3*s1 + b3*e1*s2 - b3*e2*s1;
 347     double E = d1*h2*s3 - d1*h3*s2 - d2*h1*s3 + d2*h3*s1 + d3*h1*s2 - d3*h2*s1;
 348     double F = (b1*h2 + d1*e2)*s3 - (b1*h3 + d1*e3)*s2
 349       - (b2*h1 + d2*e1)*s3 + (b2*h3 + d2*e3)*s1
 350       + (b3*h1 + d3*e1)*s2 - (b3*h2 + d3*e2)*s1;
 351     double G = (a1*e2*s3 + b1*e2*q3) - (a1*e3*s2 + b1*e3*q2)
 352       - (a2*e1*s3 + b2*e1*q3) + (a2*e3*s1 + b2*e3*q1)
 353       + (a3*e1*s2 + b3*e1*q2) - (a3*e2*s1 + b3*e2*q1);
 354     double H = (c1*h2*s3 + d1*h2*q3) - (c1*h3*s2 + d1*h3*q2)
 355       - (c2*h1*s3 + d2*h1*q3) + (c2*h3*s1 + d2*h3*q1)
 356       + (c3*h1*s2 + d3*h1*q2) - (c3*h2*s1 + d3*h2*q1);
 357     double I = ((a1*h2 + c1*e2)*s3 + (b1*h2 + d1*e2)*q3)
 358       - ((a1*h3 + c1*e3)*s2 + (b1*h3 + d1*e3)*q2)
 359       - ((a2*h1 + c2*e1)*s3 + (b2*h1 + d2*e1)*q3)
 360       + ((a2*h3 + c2*e3)*s1 + (b2*h3 + d2*e3)*q1)
 361       + ((a3*h1 + c3*e1)*s2 + (b3*h1 + d3*e1)*q2)
 362       - ((a3*h2 + c3*e2)*s1 + (b3*h2 + d3*e2)*q1);
 363     double J = a1*f2*r3 - a1*f3*r2 - a2*f1*r3 + a2*f3*r1 + a3*f1*r2 - a3*f2*r1;
 364     double K = c1*p2*r3 - c1*p3*r2 - c2*p1*r3 + c2*p3*r1 + c3*p1*r2 - c3*p2*r1;
 365     double L = (a1*p2 + c1*f2)*r3 - (a1*p3 + c1*f3)*r2
 366       - (a2*p1 + c2*f1)*r3 + (a2*p3 + c2*f3)*r1
 367       + (a3*p1 + c3*f1)*r2 - (a3*p2 + c3*f2)*r1;
 368     double M = b1*f2*t3 - b1*f3*t2 - b2*f1*t3 + b2*f3*t1 + b3*f1*t2 - b3*f2*t1;
 369     double N = d1*p2*t3 - d1*p3*t2 - d2*p1*t3 + d2*p3*t1 + d3*p1*t2 - d3*p2*t1;
 370     double O = (b1*p2 + d1*f2)*t3 - (b1*p3 + d1*f3)*t2
 371       - (b2*p1 + d2*f1)*t3 + (b2*p3 + d2*f3)*t1
 372       + (b3*p1 + d3*f1)*t2 - (b3*p2 + d3*f2)*t1;
 373     double P = (a1*f2*t3 + b1*f2*r3) - (a1*f3*t2 + b1*f3*r2)
 374       - (a2*f1*t3 + b2*f1*r3) + (a2*f3*t1 + b2*f3*r1)
 375       + (a3*f1*t2 + b3*f1*r2) - (a3*f2*t1 + b3*f2*r1);
 376     double Q = (c1*p2*t3 + d1*p2*r3) - (c1*p3*t2 + d1*p3*r2)
 377       - (c2*p1*t3 + d2*p1*r3) + (c2*p3*t1 + d2*p3*r1)
 378       + (c3*p1*t2 + d3*p1*r2) - (c3*p2*t1 + d3*p2*r1);
 379     double R = ((a1*p2 + c1*f2)*t3 + (b1*p2 + d1*f2)*r3)
 380       - ((a1*p3 + c1*f3)*t2 + (b1*p3 + d1*f3)*r2)
 381       - ((a2*p1 + c2*f1)*t3 + (b2*p1 + d2*f1)*r3)
 382       + ((a2*p3 + c2*f3)*t1 + (b2*p3 + d2*f3)*r1)
 383       + ((a3*p1 + c3*f1)*t2 + (b3*p1 + d3*f1)*r2)
 384       - ((a3*p2 + c3*f2)*t1 + (b3*p2 + d3*f2)*r1);
 385     double S = (a1*e2*r3 + a1*f2*q3) - (a1*e3*r2 + a1*f3*q2)
 386       - (a2*e1*r3 + a2*f1*q3) + (a2*e3*r1 + a2*f3*q1)
 387       + (a3*e1*r2 + a3*f1*q2) - (a3*e2*r1 + a3*f2*q1);
 388     double T = (c1*h2*r3 + c1*p2*q3) - (c1*h3*r2 + c1*p3*q2)
 389       - (c2*h1*r3 + c2*p1*q3) + (c2*h3*r1 + c2*p3*q1)
 390       + (c3*h1*r2 + c3*p1*q2) - (c3*h2*r1 + c3*p2*q1);
 391     double U = ((a1*h2 + c1*e2)*r3 + (a1*p2 + c1*f2)*q3)
 392       - ((a1*h3 + c1*e3)*r2 + (a1*p3 + c1*f3)*q2)
 393       - ((a2*h1 + c2*e1)*r3 + (a2*p1 + c2*f1)*q3)
 394       + ((a2*h3 + c2*e3)*r1 + (a2*p3 + c2*f3)*q1)
 395       + ((a3*h1 + c3*e1)*r2 + (a3*p1 + c3*f1)*q2)
 396       - ((a3*h2 + c3*e2)*r1 + (a3*p2 + c3*f2)*q1);
 397     double V = (b1*e2*t3 + b1*f2*s3) - (b1*e3*t2 + b1*f3*s2)
 398       - (b2*e1*t3 + b2*f1*s3) + (b2*e3*t1 + b2*f3*s1)
 399       + (b3*e1*t2 + b3*f1*s2) - (b3*e2*t1 + b3*f2*s1);
 400     double W = (d1*h2*t3 + d1*p2*s3) - (d1*h3*t2 + d1*p3*s2)
 401       - (d2*h1*t3 + d2*p1*s3) + (d2*h3*t1 + d2*p3*s1)
 402       + (d3*h1*t2 + d3*p1*s2) - (d3*h2*t1 + d3*p2*s1);
 403     double X = ((b1*h2 + d1*e2)*t3 + (b1*p2 + d1*f2)*s3)
 404       - ((b1*h3 + d1*e3)*t2 + (b1*p3 + d1*f3)*s2)
 405       - ((b2*h1 + d2*e1)*t3 + (b2*p1 + d2*f1)*s3)
 406       + ((b2*h3 + d2*e3)*t1 + (b2*p3 + d2*f3)*s1)
 407       + ((b3*h1 + d3*e1)*t2 + (b3*p1 + d3*f1)*s2)
 408       - ((b3*h2 + d3*e2)*t1 + (b3*p2 + d3*f2)*s1);
 409     double Y = (a1*e2*t3 + a1*f2*s3 + b1*e2*r3 + b1*f2*q3)
 410       - (a1*e3*t2 + a1*f3*s2 + b1*e3*r2 + b1*f3*q2)
 411       - (a2*e1*t3 + a2*f1*s3 + b2*e1*r3 + b2*f1*q3)
 412       + (a2*e3*t1 + a2*f3*s1 + b2*e3*r1 + b2*f3*q1)
 413       + (a3*e1*t2 + a3*f1*s2 + b3*e1*r2 + b3*f1*q2)
 414       - (a3*e2*t1 + a3*f2*s1 + b3*e2*r1 + b3*f2*q1);
 415     double Z = (c1*h2*t3 + c1*p2*s3 + d1*h2*r3 + d1*p2*q3)
 416       - (c1*h3*t2 + c1*p3*s2 + d1*h3*r2 + d1*p3*q2)
 417       - (c2*h1*t3 + c2*p1*s3 + d2*h1*r3 + d2*p1*q3)
 418       + (c2*h3*t1 + c2*p3*s1 + d2*h3*r1 + d2*p3*q1)
 419       + (c3*h1*t2 + c3*p1*s2 + d3*h1*r2 + d3*p1*q2)
 420       - (c3*h2*t1 + c3*p2*s1 + d3*h2*r1 + d3*p2*q1);
 421     double AA = ((a1*h2 + c1*e2)*t3 + (a1*p2 + c1*f2)*s3
 422                  +(b1*h2 + d1*e2)*r3 + (b1*p2 + d1*f2)*q3)
 423       - ((a1*h3 + c1*e3)*t2 + (a1*p3 + c1*f3)*s2
 424          +(b1*h3 + d1*e3)*r2 + (b1*p3 + d1*f3)*q2)
 425       - ((a2*h1 + c2*e1)*t3 + (a2*p1 + c2*f1)*s3
 426          +(b2*h1 + d2*e1)*r3 + (b2*p1 + d2*f1)*q3)
 427       + ((a2*h3 + c2*e3)*t1 + (a2*p3 + c2*f3)*s1
 428          +(b2*h3 + d2*e3)*r1 + (b2*p3 + d2*f3)*q1)
 429       + ((a3*h1 + c3*e1)*t2 + (a3*p1 + c3*f1)*s2
 430          +(b3*h1 + d3*e1)*r2 + (b3*p1 + d3*f1)*q2)
 431       - ((a3*h2 + c3*e2)*t1 + (a3*p2 + c3*f2)*s1
 432          + (b3*h2 + d3*e2)*r1 + (b3*p2 + d3*f2)*q1);
 433
 434     return  64.0*(  8.0*(A + B + D + E + J + K + M + N)
 435                     + 4.0*(C + F + G + H + L + O + P + Q + S + T + V + W)
 436                     + 2.0*(I + R + U + X + Y + Z) + AA ) / 27.0 ;
 437   }
 438
 439   // =========================================================================================================================
 440   // Calculate Volume for Generic Polyedron, even not convex one, WARNING !!! The polyedron's faces must be correctly ordered
 441   // =========================================================================================================================
 442   inline double calculateVolumeForPolyh(const double ***pts,
 443                                         const int *nbOfNodesPerFaces,
 444                                         int nbOfFaces,
 445                                         const double *bary)
 446   {
 447     double volume=0.;
 448
 449     for ( int i=0; i<nbOfFaces; i++ )
 450       {
 451         double normal[3];
 452         double vecForAlt[3];
 453
 454         calculateNormalForPolyg(pts[i],nbOfNodesPerFaces[i],normal);
 455         vecForAlt[0]=bary[0]-pts[i][0][0];
 456         vecForAlt[1]=bary[1]-pts[i][0][1];
 457         vecForAlt[2]=bary[2]-pts[i][0][2];
 458         volume+=vecForAlt[0]*normal[0]+vecForAlt[1]*normal[1]+vecForAlt[2]*normal[2];
 459       }
 460     return -volume/3.;
 461   }
 462
 463   /*!
 464    * Calculate Volume for Generic Polyedron, even not convex one, WARNING !!! The polyedron's faces must be correctly ordered.
 465    * 2nd API avoiding to create double** arrays. The returned value could be negative if polyhedrons faces are not oriented with normal outside of the
 466    * polyhedron
 467    */
 468   template<class ConnType, NumberingPolicy numPol>
 469   inline double calculateVolumeForPolyh2(const ConnType *connec, mcIdType lgth, const double *coords)
 470   {
 471     std::size_t nbOfFaces=std::count(connec,connec+lgth,-1)+1;
 472     double volume=0.;
 473     const ConnType *work=connec;
 474     for(std::size_t iFace=0;iFace<nbOfFaces;iFace++)
 475       {
 476         const ConnType *work2=std::find(work+1,connec+lgth,-1);
 477         std::size_t nbOfNodesOfCurFace=std::distance(work,work2);
 478         double areaVector[3]={0.,0.,0.};
 479         for(std::size_t ptId=0;ptId<nbOfNodesOfCurFace;ptId++)
 480           {
 481             const double *pti=coords+3*OTT<ConnType,numPol>::coo2C(work[ptId]);
 482             const double *pti1=coords+3*OTT<ConnType,numPol>::coo2C(work[(ptId+1)%nbOfNodesOfCurFace]);
 483             areaVector[0]+=pti[1]*pti1[2]-pti[2]*pti1[1];
 484             areaVector[1]+=pti[2]*pti1[0]-pti[0]*pti1[2];
 485             areaVector[2]+=pti[0]*pti1[1]-pti[1]*pti1[0];
 486           }
 487         const double *pt=coords+3*work[0];
 488         volume+=pt[0]*areaVector[0]+pt[1]*areaVector[1]+pt[2]*areaVector[2];
 489         work=work2+1;
 490       }
 491     return volume/6.;
 492   }
 493
 494   /*!
 495    * This method returns the area oriented vector of a polygon. This method is useful for normal computation without
 496    * any troubles if several edges are colinears.
 497    * @param res must be of size at least 3 to store the result.
 498    */
 499   template<class ConnType, NumberingPolicy numPol>
 500   inline void areaVectorOfPolygon(const ConnType *connec, mcIdType lgth, const double *coords, double *res)
 501   {
 502     res[0]=0.; res[1]=0.; res[2]=0.;
 503     for(mcIdType ptId=0;ptId<lgth;ptId++)
 504       {
 505         const double *pti=coords+3*OTT<ConnType,numPol>::coo2C(connec[ptId]);
 506         const double *pti1=coords+3*OTT<ConnType,numPol>::coo2C(connec[(ptId+1)%lgth]);
 507         res[0]+=pti[1]*pti1[2]-pti[2]*pti1[1];
 508         res[1]+=pti[2]*pti1[0]-pti[0]*pti1[2];
 509         res[2]+=pti[0]*pti1[1]-pti[1]*pti1[0];
 510       }
 511   }
 512
 513   template<class ConnType, NumberingPolicy numPol>
 514   inline void computePolygonBarycenter3D(const ConnType *connec, mcIdType lgth, const double *coords, double *res)
 515   {
 516     double area[3];
 517     areaVectorOfPolygon<ConnType,numPol>(connec,lgth,coords,area);
 518     double norm=sqrt(area[0]*area[0]+area[1]*area[1]+area[2]*area[2]);
 519     if(norm>std::numeric_limits<double>::min())
 520       {
 521         area[0]/=norm; area[1]/=norm; area[2]/=norm;
 522         res[0]=0.; res[1]=0.; res[2]=0.;
 523         for(mcIdType i=1;i<lgth-1;i++)
 524           {
 525             double v[3];
 526             double tmpArea[3];
 527             v[0]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[0])]+
 528                   coords[3*OTT<ConnType,numPol>::coo2C(connec[i])]+
 529                   coords[3*OTT<ConnType,numPol>::coo2C(connec[i+1])])/3.;
 530             v[1]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[0])+1]+
 531                   coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+1]+
 532                   coords[3*OTT<ConnType,numPol>::coo2C(connec[i+1])+1])/3.;
 533             v[2]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[0])+2]+
 534                   coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+2]+
 535                   coords[3*OTT<ConnType,numPol>::coo2C(connec[i+1])+2])/3.;
 536             ConnType tmpConn[3]={connec[0],connec[i],connec[i+1]};
 537             areaVectorOfPolygon<ConnType,numPol>(tmpConn,3,coords,tmpArea);
 538             double norm2=sqrt(tmpArea[0]*tmpArea[0]+tmpArea[1]*tmpArea[1]+tmpArea[2]*tmpArea[2]);
 539             if(norm2>1e-12)
 540               {
 541                 tmpArea[0]/=norm2; tmpArea[1]/=norm2; tmpArea[2]/=norm2;
 542                 double signOfArea=area[0]*tmpArea[0]+area[1]*tmpArea[1]+area[2]*tmpArea[2];
 543                 res[0]+=signOfArea*norm2*v[0]/norm; res[1]+=signOfArea*norm2*v[1]/norm; res[2]+=signOfArea*norm2*v[2]/norm;
 544               }
 545           }
 546       }
 547     else
 548       {
 549         res[0]=0.; res[1]=0.; res[2]=0.;
 550         if(lgth<1)
 551           throw INTERP_KERNEL::Exception("computePolygonBarycenter3D : lgth of polygon is < 1 !");
 552         norm=0.;
 553         double v[3];
 554         for(mcIdType i=0;i<lgth;i++)
 555           {
 556             v[0]=coords[3*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])]-coords[3*OTT<ConnType,numPol>::coo2C(connec[i])];
 557             v[1]=coords[3*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])+1]-coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+1];
 558             v[2]=coords[3*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])+2]-coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+2];
 559             double norm2=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
 560             res[0]+=(coords[3*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])]+coords[3*OTT<ConnType,numPol>::coo2C(connec[i])])/2.*norm2;
 561             res[1]+=(coords[3*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])+1]+coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+1])/2.*norm2;
 562             res[2]+=(coords[3*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])+2]+coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+2])/2.*norm2;
 563             norm+=norm2;
 564           }
 565         if(norm>std::numeric_limits<double>::min())
 566           {
 567             res[0]/=norm; res[1]/=norm; res[2]/=norm;
 568             return;
 569           }
 570         else
 571           {
 572             res[0]=0.; res[1]=0.; res[2]=0.;
 573             for(mcIdType i=0;i<lgth;i++)
 574               {
 575                 res[0]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[i])];
 576                 res[1]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+1];
 577                 res[2]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+2];
 578               }
 579             res[0]/=FromIdType<double>(lgth); res[1]/=FromIdType<double>(lgth); res[2]/=FromIdType<double>(lgth);
 580             return;
 581           }
 582       }
 583   }
 584
 585   inline double integrationOverA3DLine(double u1, double v1, double u2, double v2, double A, double B, double C)
 586   {
 587     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))+
 588                     4.*C*(A*(2*u1*v1+u2*v1+u1*v2+2.*u2*v2)+B*(v1*v1+v1*v2+v2*v2))+A*B*(u1*(3.*v1*v1+2.*v1*v2+v2*v2)+u2*(v1*v1+2.*v1*v2+3.*v2*v2)))/24.;
 589   }
 590
 591   template<class ConnType, NumberingPolicy numPol>
 592   inline void barycenterOfPolyhedron(const ConnType *connec, mcIdType lgth, const double *coords, double *res)
 593   {
 594     std::size_t nbOfFaces=std::count(connec,connec+lgth,-1)+1;
 595     res[0]=0.; res[1]=0.; res[2]=0.;
 596     const ConnType *work=connec;
 597     for(std::size_t i=0;i<nbOfFaces;i++)
 598       {
 599         const ConnType *work2=std::find(work+1,connec+lgth,-1);
 600         int nbOfNodesOfCurFace=(int)std::distance(work,work2);
 601         // projection to (u,v) of each faces of polyh to compute integral(x^2/2) on each faces.
 602         double normal[3];
 603         areaVectorOfPolygon<ConnType,numPol>(work,nbOfNodesOfCurFace,coords,normal);
 604         double normOfNormal=sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
 605         if(normOfNormal<std::numeric_limits<double>::min())
 606           continue;
 607         normal[0]/=normOfNormal; normal[1]/=normOfNormal; normal[2]/=normOfNormal;
 608         double u[2]={normal[1],-normal[0]};
 609         double s=sqrt(u[0]*u[0]+u[1]*u[1]);
 610         double c=normal[2];
 611         if(fabs(s)>1e-12)
 612           {
 613             u[0]/=std::abs(s); u[1]/=std::abs(s);
 614           }
 615         else
 616           { u[0]=1.; u[1]=0.; }
 617         //C : high in plane of polyhedron face : always constant
 618         double w=normal[0]*coords[3*OTT<ConnType,numPol>::coo2C(work[0])]+
 619           normal[1]*coords[3*OTT<ConnType,numPol>::coo2C(work[0])+1]+
 620           normal[2]*coords[3*OTT<ConnType,numPol>::coo2C(work[0])+2];
 621         // A,B,D,F,G,H,L,M,N coeffs of rotation matrix defined by (u,c,s)
 622         double A=u[0]*u[0]*(1-c)+c;
 623         double B=u[0]*u[1]*(1-c);
 624         double D=u[1]*s;
 625         double F=B;
 626         double G=u[1]*u[1]*(1-c)+c;
 627         double H=-u[0]*s;
 628         double L=-u[1]*s;
 629         double M=u[0]*s;
 630         double N=c;
 631         double CX=-w*D;
 632         double CY=-w*H;
 633         double CZ=-w*N;
 634         for(int j=0;j<nbOfNodesOfCurFace;j++)
 635           {
 636             const double *p1=coords+3*OTT<ConnType,numPol>::coo2C(work[j]);
 637             const double *p2=coords+3*OTT<ConnType,numPol>::coo2C(work[(j+1)%nbOfNodesOfCurFace]);
 638             double u1=A*p1[0]+B*p1[1]+D*p1[2];
 639             double u2=A*p2[0]+B*p2[1]+D*p2[2];
 640             double v1=F*p1[0]+G*p1[1]+H*p1[2];
 641             double v2=F*p2[0]+G*p2[1]+H*p2[2];
 642             //
 643             double gx=integrationOverA3DLine(u1,v1,u2,v2,A,B,CX);
 644             double gy=integrationOverA3DLine(u1,v1,u2,v2,F,G,CY);
 645             double gz=integrationOverA3DLine(u1,v1,u2,v2,L,M,CZ);
 646             res[0]+=gx*normal[0];
 647             res[1]+=gy*normal[1];
 648             res[2]+=gz*normal[2];
 649           }
 650         work=work2+1;
 651       }
 652     double vol=calculateVolumeForPolyh2<ConnType,numPol>(connec,lgth,coords);
 653     if(fabs(vol)>std::numeric_limits<double>::min())
 654       {
 655         res[0]/=vol; res[1]/=vol; res[2]/=vol;
 656       }
 657     else
 658       {
 659         double sum=0.;
 660         res[0]=0.; res[1]=0.; res[2]=0.;
 661         work=connec;
 662         for(std::size_t i=0;i<nbOfFaces;i++)
 663           {
 664             const ConnType *work2=std::find(work+1,connec+lgth,-1);
 665             int nbOfNodesOfCurFace=(int)std::distance(work,work2);
 666             double normal[3];
 667             areaVectorOfPolygon<ConnType,numPol>(work,nbOfNodesOfCurFace,coords,normal);
 668             double normOfNormal=sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
 669             if(normOfNormal<std::numeric_limits<double>::min())
 670               continue;
 671             sum+=normOfNormal;
 672             double tmpBary[3];
 673             computePolygonBarycenter3D<ConnType,numPol>(work,nbOfNodesOfCurFace,coords,tmpBary);
 674             res[0]+=normOfNormal*tmpBary[0]; res[1]+=normOfNormal*tmpBary[1]; res[2]+=normOfNormal*tmpBary[2];
 675             work=work2+1;
 676           }
 677         res[0]/=sum; res[1]/=sum; res[2]/=sum;
 678       }
 679   }
 680
 681   // ============================================================================================================================================
 682   // Calculate Volume for NON Generic Polyedron. Only polydrons with bary included in pts is supported by this method. Result is always positive.
 683   // ============================================================================================================================================
 684   inline double calculateVolumeForPolyhAbs(const double ***pts,
 685                                            const int *nbOfNodesPerFaces,
 686                                            int nbOfFaces,
 687                                            const double *bary)
 688   {
 689     double volume=0.;
 690
 691     for ( int i=0; i<nbOfFaces; i++ )
 692       {
 693         double normal[3];
 694         double vecForAlt[3];
 695
 696         calculateNormalForPolyg(pts[i],nbOfNodesPerFaces[i],normal);
 697         vecForAlt[0]=bary[0]-pts[i][0][0];
 698         vecForAlt[1]=bary[1]-pts[i][0][1];
 699         vecForAlt[2]=bary[2]-pts[i][0][2];
 700         volume+=fabs(vecForAlt[0]*normal[0]+vecForAlt[1]*normal[1]+vecForAlt[2]*normal[2]);
 701       }
 702     return volume/3.;
 703   }
 704
 705   template<int N>
 706   inline double addComponentsOfVec(const double **pts, int rk)
 707   {
 708     return pts[N-1][rk]+addComponentsOfVec<N-1>(pts,rk);
 709   }
 710
 711   template<>
 712   inline double addComponentsOfVec<1>(const double **pts, int rk)
 713   {
 714     return pts[0][rk];
 715   }
 716
 717   template<int N, int DIM>
 718   inline void calculateBarycenter(const double **pts, double *bary)
 719   {
 720     bary[DIM-1]=addComponentsOfVec<N>(pts,DIM-1)/N;
 721     calculateBarycenter<N,DIM-1>(pts,bary);
 722   }
 723
 724   template<>
 725   inline void calculateBarycenter<2,0>(const double ** /*pts*/, double * /*bary*/)
 726   {
 727   }
 728
 729   template<>
 730   inline void calculateBarycenter<3,0>(const double ** /*pts*/, double * /*bary*/)
 731   {
 732   }
 733
 734   template<>
 735   inline void calculateBarycenter<4,0>(const double ** /*pts*/, double * /*bary*/)
 736   {
 737   }
 738
 739   template<>
 740   inline void calculateBarycenter<5,0>(const double ** /*pts*/, double * /*bary*/)
 741   {
 742   }
 743
 744   template<>
 745   inline void calculateBarycenter<6,0>(const double ** /*pts*/, double * /*bary*/)
 746   {
 747   }
 748
 749   template<>
 750   inline void calculateBarycenter<7,0>(const double ** /*pts*/, double * /*bary*/)
 751   {
 752   }
 753
 754   template<>
 755   inline void calculateBarycenter<8,0>(const double ** /*pts*/, double * /*bary*/)
 756   {
 757   }
 758
 759   inline void calculateBarycenterDyn(const double **pts, mcIdType nbPts,
 760                                      int dim, double *bary)
 761   {
 762     for(int i=0;i<dim;i++)
 763       {
 764         double temp=0.;
 765         for(mcIdType j=0;j<nbPts;j++)
 766           {
 767             temp+=pts[j][i];
 768           }
 769         bary[i]=temp/FromIdType<double>(nbPts);
 770       }
 771   }
 772
 773   template<int SPACEDIM>
 774   inline void calculateBarycenterDyn2(const double *pts, mcIdType nbPts, double *bary)
 775   {
 776     for(int i=0;i<SPACEDIM;i++)
 777       {
 778         double temp=0.;
 779         for(mcIdType j=0;j<nbPts;j++)
 780           {
 781             temp+=pts[j*SPACEDIM+i];
 782           }
 783         bary[i]=temp/FromIdType<double>(nbPts);
 784       }
 785   }
 786
 787   inline void computePolygonBarycenter2DEngine(double **coords, mcIdType lgth, double *res)
 788   {
 789     double area=0.;
 790     res[0]=0.; res[1]=0.;
 791     for(mcIdType i=0;i<lgth;i++)
 792       {
 793         double cp=coords[i][0]*coords[(i+1)%lgth][1]-coords[i][1]*coords[(i+1)%lgth][0];
 794         area+=cp;
 795         res[0]+=cp*(coords[i][0]+coords[(i+1)%lgth][0]);
 796         res[1]+=cp*(coords[i][1]+coords[(i+1)%lgth][1]);
 797       }
 798     res[0]/=3.*area;
 799     res[1]/=3.*area;
 800   }
 801
 802   template<class ConnType, NumberingPolicy numPol>
 803   inline void computePolygonBarycenter2D(const ConnType *connec, mcIdType lgth, const double *coords, double *res)
 804   {
 805     double **coords2=new double *[lgth];
 806     for(mcIdType i=0;i<lgth;i++)
 807       coords2[i]=const_cast<double *>(coords+2*OTT<ConnType,numPol>::coo2C(connec[i]));
 808     computePolygonBarycenter2DEngine(coords2,lgth,res);
 809     delete [] coords2;
 810   }
 811
 812   inline void computeQPolygonBarycenter2D(double **coords, mcIdType nbOfPtsInPolygs, int spaceDim, double *res)
 813   {
 814     if(nbOfPtsInPolygs%2==0)
 815       {
 816         if(spaceDim==2)
 817           {
 818             std::vector<Node *> nodes(nbOfPtsInPolygs);
 819             for(mcIdType i=0;i<nbOfPtsInPolygs;i++)
 820               nodes[i]=new Node(coords[i][0],coords[i][1]);
 821             QuadraticPolygon *pol=QuadraticPolygon::BuildArcCirclePolygon(nodes);
 822             pol->getBarycenter(res);
 823             delete pol;
 824           }
 825         else
 826           return computePolygonBarycenter2DEngine(coords,nbOfPtsInPolygs/2,res);
 827       }
 828     else
 829       {
 830         std::ostringstream oss; oss << "INTERP_KERNEL::computeQPolygonBarycenter2D : nb of points in quadratic polygon is " << nbOfPtsInPolygs << " should be even !";
 831         throw INTERP_KERNEL::Exception(oss.str().c_str());
 832       }
 833   }
 834 }
 835
 836 #endif