src/INTERP_KERNEL/VolSurfFormulae.hxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D
   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.
   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 #ifndef VOLSURFFORMULAE
  20 #define VOLSURFFORMULAE
  21
  22 #include <math.h>
  23
  24 namespace INTERP_KERNEL
  25 {
  26   inline void calculateBarycenterDyn(const double **pts, int nbPts,
  27                                      int dim, double *bary);
  28
  29   inline double calculateAreaForPolyg(const double **coords, int nbOfPtsInPolygs,
  30                                       int spaceDim);
  31
  32
  33   // ===========================
  34   // Calculate Area for triangle
  35   // ===========================
  36   inline double calculateAreaForTria(const double *p1, const double *p2,
  37                                      const double *p3, int spaceDim)
  38   {
  39     double area ;
  40
  41     if ( spaceDim == 2 )
  42       {
  43         area = -((p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1]))/2.0;
  44       }
  45     else
  46       {
  47         area = sqrt(((p2[1]-p1[1])*(p3[2]-p1[2]) - (p3[1]-p1[1])*(p2[2]-p1[2]))*
  48                     ((p2[1]-p1[1])*(p3[2]-p1[2]) - (p3[1]-p1[1])*(p2[2]-p1[2]))
  49                     +
  50                     ((p3[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p3[2]-p1[2]))*
  51                     ((p3[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p3[2]-p1[2]))
  52                     +
  53                     ((p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1]))*
  54                     ((p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1])))/2.0;
  55       }
  56
  57     return area ;
  58   }
  59
  60   // =============================
  61   // Calculate Area for quadrangle
  62   // =============================
  63   inline double calculateAreaForQuad(const double *p1, const double *p2,
  64                                      const double *p3, const double *p4,
  65                                      int spaceDim)
  66   {
  67     double area ;
  68
  69     if (spaceDim==2)
  70       {
  71         double a1 = (p2[0]-p1[0])/4.0, a2 = (p2[1]-p1[1])/4.0;
  72         double b1 = (p3[0]-p4[0])/4.0, b2 = (p3[1]-p4[1])/4.0;
  73         double c1 = (p3[0]-p2[0])/4.0, c2 = (p3[1]-p2[1])/4.0;
  74         double d1 = (p4[0]-p1[0])/4.0, d2 = (p4[1]-p1[1])/4.0;
  75
  76         area = - 4.0*(  b1*c2 - c1*b2 + a1*c2 - c1*a2
  77                         + b1*d2 - d1*b2 + a1*d2 - d1*a2);
  78       }
  79     else
  80       {
  81         area = (sqrt(((p2[1]-p1[1])*(p4[2]-p1[2]) - (p4[1]-p1[1])*(p2[2]-p1[2]))*
  82                      ((p2[1]-p1[1])*(p4[2]-p1[2]) - (p4[1]-p1[1])*(p2[2]-p1[2]))
  83                      + ((p4[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p4[2]-p1[2]))*
  84                      ((p4[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p4[2]-p1[2]))
  85                      + ((p2[0]-p1[0])*(p4[1]-p1[1]) - (p4[0]-p1[0])*(p2[1]-p1[1]))*
  86                      ((p2[0]-p1[0])*(p4[1]-p1[1]) - (p4[0]-p1[0])*(p2[1]-p1[1])))
  87                 +
  88                 sqrt(((p4[1]-p3[1])*(p2[2]-p3[2]) - (p2[1]-p3[1])*(p4[2]-p3[2]))*
  89                      ((p4[1]-p3[1])*(p2[2]-p3[2]) - (p2[1]-p3[1])*(p4[2]-p3[2]))
  90                      + ((p2[0]-p3[0])*(p4[2]-p3[2]) - (p4[0]-p3[0])*(p2[2]-p3[2]))*
  91                      ((p2[0]-p3[0])*(p4[2]-p3[2]) - (p4[0]-p3[0])*(p2[2]-p3[2]))
  92                      + ((p4[0]-p3[0])*(p2[1]-p3[1]) - (p2[0]-p3[0])*(p4[1]-p3[1]))*
  93                      ((p4[0]-p3[0])*(p2[1]-p3[1]) - (p2[0]-p3[0])*(p4[1]-p3[1])))
  94                 )/2.0;
  95       }
  96
  97     return area ;
  98   }
  99
 100   // ====================================
 101   // Calculate Normal Vector for Triangle
 102   // ====================================
 103   inline void calculateNormalForTria(const double *p1, const double *p2,
 104                                      const double *p3, double *normal)
 105   {
 106     normal[0] = ((p2[1]-p1[1])*(p3[2]-p1[2]) - (p3[1]-p1[1])*(p2[2]-p1[2]))/2.0;
 107     normal[1] = ((p3[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p3[2]-p1[2]))/2.0;
 108     normal[2] = ((p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1]))/2.0;
 109   }
 110
 111   // ======================================
 112   // Calculate Normal Vector for Quadrangle
 113   // ======================================
 114   inline void calculateNormalForQuad(const double *p1, const double *p2,
 115                                      const double *p3, const double *p4,
 116                                      double *normal)
 117   {
 118     double xnormal1 = (p2[1]-p1[1])*(p4[2]-p1[2]) - (p4[1]-p1[1])*(p2[2]-p1[2]);
 119     double xnormal2 = (p4[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p4[2]-p1[2]);
 120     double xnormal3 = (p2[0]-p1[0])*(p4[1]-p1[1]) - (p4[0]-p1[0])*(p2[1]-p1[1]);
 121     double xarea = sqrt(xnormal1*xnormal1 + xnormal2*xnormal2 + xnormal3*xnormal3);
 122     xnormal1 = xnormal1/xarea;
 123     xnormal2 = xnormal2/xarea;
 124     xnormal3 = xnormal3/xarea;
 125     xarea = calculateAreaForQuad(p1,p2,p3,p4,3);
 126     normal[0] = xnormal1*xarea ;
 127     normal[1] = xnormal2*xarea ;
 128     normal[2] = xnormal3*xarea ;
 129   }
 130
 131   // ===================================
 132   // Calculate Normal Vector for Polygon
 133   // ===================================
 134   inline void calculateNormalForPolyg(const double **coords, int nbOfPtsInPolygs,
 135                                       double *normal)
 136   {
 137     double coordOfBary[3];
 138
 139     calculateBarycenterDyn(coords,nbOfPtsInPolygs,3,coordOfBary);
 140     double xnormal1 = (coords[0][1]-coords[1][1]) * (coordOfBary[2]-coords[1][2])
 141       - (coords[0][2]-coords[1][2]) * (coordOfBary[1]-coords[1][1]);
 142
 143     double xnormal2 = (coords[0][2]-coords[1][2]) * (coordOfBary[0]-coords[1][0])
 144       - (coords[0][0]-coords[1][0]) * (coordOfBary[2]-coords[1][2]);
 145
 146     double xnormal3 = (coords[0][0]-coords[1][0]) * (coordOfBary[1]-coords[1][1])
 147       - (coords[0][1]-coords[1][1]) * (coordOfBary[0]-coords[1][0]);
 148
 149     double xarea=sqrt(xnormal1*xnormal1 + xnormal2*xnormal2 + xnormal3*xnormal3);
 150
 151     if ( xarea < 1.e-6 )
 152       {
 153         //std::string diagnosis"Can not calculate normal - polygon is singular";
 154         throw std::exception();
 155       }
 156
 157     xnormal1 = -xnormal1/xarea;
 158     xnormal2 = -xnormal2/xarea;
 159     xnormal3 = -xnormal3/xarea;
 160     xarea = calculateAreaForPolyg(coords,nbOfPtsInPolygs,3);
 161     normal[0] = xnormal1*xarea ;
 162     normal[1] = xnormal2*xarea ;
 163     normal[2] = xnormal3*xarea ;
 164   }
 165
 166   // ==========================
 167   // Calculate Area for Polygon
 168   // ==========================
 169   inline double calculateAreaForPolyg(const double **coords, int nbOfPtsInPolygs,
 170                                       int spaceDim)
 171   {
 172     double ret=0.;
 173     double coordOfBary[3];
 174
 175     calculateBarycenterDyn(coords,nbOfPtsInPolygs,spaceDim,coordOfBary);
 176     for ( int i=0; i<nbOfPtsInPolygs; i++ )
 177       {
 178         double tmp = calculateAreaForTria(coords[i],coords[(i+1)%nbOfPtsInPolygs],
 179                                           coordOfBary,spaceDim);
 180         ret+=tmp;
 181       }
 182     return ret;
 183   }
 184
 185   // ==========================
 186   // Calculate Volume for Tetra
 187   // ==========================
 188   inline double calculateVolumeForTetra(const double *p1, const double *p2,
 189                                         const double *p3, const double *p4)
 190   {
 191     return (  (p3[0]-p1[0])*(  (p2[1]-p1[1])*(p4[2]-p1[2])
 192                                - (p2[2]-p1[2])*(p4[1]-p1[1]) )
 193               - (p2[0]-p1[0])*(  (p3[1]-p1[1])*(p4[2]-p1[2])
 194                                  - (p3[2]-p1[2])*(p4[1]-p1[1]) )
 195               + (p4[0]-p1[0])*(  (p3[1]-p1[1])*(p2[2]-p1[2])
 196                                  - (p3[2]-p1[2])*(p2[1]-p1[1]) )
 197               ) / 6.0 ;
 198   }
 199
 200   // =========================
 201   // Calculate Volume for Pyra
 202   // =========================
 203   inline double calculateVolumeForPyra(const double *p1, const double *p2,
 204                                        const double *p3, const double *p4,
 205                                        const double *p5)
 206   {
 207     return ( ((p3[0]-p1[0])*(  (p2[1]-p1[1])*(p5[2]-p1[2])
 208                                - (p2[2]-p1[2])*(p5[1]-p1[1]) )
 209               -(p2[0]-p1[0])*(  (p3[1]-p1[1])*(p5[2]-p1[2])
 210                                 - (p3[2]-p1[2])*(p5[1]-p1[1]) )
 211               +(p5[0]-p1[0])*(  (p3[1]-p1[1])*(p2[2]-p1[2])
 212                                 - (p3[2]-p1[2])*(p2[1]-p1[1]) ))
 213              +
 214              ((p4[0]-p1[0])*(  (p3[1]-p1[1])*(p5[2]-p1[2])
 215                                - (p3[2]-p1[2])*(p5[1]-p1[1]) )
 216               -(p3[0]-p1[0])*(  (p4[1]-p1[1])*(p5[2]-p1[2])
 217                                 - (p4[2]-p1[2])*(p5[1]-p1[1]))
 218               +(p5[0]-p1[0])*(  (p4[1]-p1[1])*(p3[2]-p1[2])
 219                                 - (p4[2]-p1[2])*(p3[1]-p1[1]) ))
 220              ) / 6.0 ;
 221   }
 222
 223   // ==========================
 224   // Calculate Volume for Penta
 225   // ==========================
 226   inline double calculateVolumeForPenta(const double *p1, const double *p2,
 227                                         const double *p3, const double *p4,
 228                                         const double *p5, const double *p6)
 229   {
 230     double a1 = (p2[0]-p3[0])/2.0, a2 = (p2[1]-p3[1])/2.0, a3 = (p2[2]-p3[2])/2.0;
 231     double b1 = (p5[0]-p6[0])/2.0, b2 = (p5[1]-p6[1])/2.0, b3 = (p5[2]-p6[2])/2.0;
 232     double c1 = (p4[0]-p1[0])/2.0, c2 = (p4[1]-p1[1])/2.0, c3 = (p4[2]-p1[2])/2.0;
 233     double d1 = (p5[0]-p2[0])/2.0, d2 = (p5[1]-p2[1])/2.0, d3 = (p5[2]-p2[2])/2.0;
 234     double e1 = (p6[0]-p3[0])/2.0, e2 = (p6[1]-p3[1])/2.0, e3 = (p6[2]-p3[2])/2.0;
 235     double f1 = (p1[0]-p3[0])/2.0, f2 = (p1[1]-p3[1])/2.0, f3 = (p1[2]-p3[2])/2.0;
 236     double h1 = (p4[0]-p6[0])/2.0, h2 = (p4[1]-p6[1])/2.0, h3 = (p4[2]-p6[2])/2.0;
 237
 238     double A = a1*c2*f3 - a1*c3*f2 - a2*c1*f3 + a2*c3*f1 + a3*c1*f2 - a3*c2*f1;
 239     double B = b1*c2*h3 - b1*c3*h2 - b2*c1*h3 + b2*c3*h1 + b3*c1*h2 - b3*c2*h1;
 240     double C = (a1*c2*h3 + b1*c2*f3) - (a1*c3*h2 + b1*c3*f2)
 241       - (a2*c1*h3 + b2*c1*f3) + (a2*c3*h1 + b2*c3*f1)
 242       + (a3*c1*h2 + b3*c1*f2) - (a3*c2*h1 + b3*c2*f1);
 243     double D = a1*d2*f3 - a1*d3*f2 - a2*d1*f3 + a2*d3*f1 + a3*d1*f2 - a3*d2*f1;
 244     double E = b1*d2*h3 - b1*d3*h2 - b2*d1*h3 + b2*d3*h1 + b3*d1*h2 - b3*d2*h1;
 245     double F = (a1*d2*h3 + b1*d2*f3) - (a1*d3*h2 + b1*d3*f2)
 246       - (a2*d1*h3 + b2*d1*f3) + (a2*d3*h1 + b2*d3*f1)
 247       + (a3*d1*h2 + b3*d1*f2) - (a3*d2*h1 + b3*d2*f1);
 248     double G = a1*e2*f3 - a1*e3*f2 - a2*e1*f3 + a2*e3*f1 + a3*e1*f2 - a3*e2*f1;
 249     double H = b1*e2*h3 - b1*e3*h2 - b2*e1*h3 + b2*e3*h1 + b3*e1*h2 - b3*e2*h1;
 250     double P = (a1*e2*h3 + b1*e2*f3) - (a1*e3*h2 + b1*e3*f2)
 251       - (a2*e1*h3 + b2*e1*f3) + (a2*e3*h1 + b2*e3*f1)
 252       + (a3*e1*h2 + b3*e1*f2) - (a3*e2*h1 + b3*e2*f1);
 253
 254     return (-2.0*(2.0*(A + B + D + E + G + H) + C + F + P)/9.0);
 255   }
 256
 257   // =========================
 258   // Calculate Volume for Hexa
 259   // =========================
 260   inline double calculateVolumeForHexa(const double *pt1, const double *pt2,
 261                                        const double *pt3, const double *pt4,
 262                                        const double *pt5, const double *pt6,
 263                                        const double *pt7, const double *pt8)
 264   {
 265     double a1=(pt3[0]-pt4[0])/8.0, a2=(pt3[1]-pt4[1])/8.0, a3=(pt3[2]-pt4[2])/8.0;
 266     double b1=(pt2[0]-pt1[0])/8.0, b2=(pt2[1]-pt1[1])/8.0, b3=(pt2[2]-pt1[2])/8.0;
 267     double c1=(pt7[0]-pt8[0])/8.0, c2=(pt7[1]-pt8[1])/8.0, c3=(pt7[2]-pt8[2])/8.0;
 268     double d1=(pt6[0]-pt5[0])/8.0, d2=(pt6[1]-pt5[1])/8.0, d3=(pt6[2]-pt5[2])/8.0;
 269     double e1=(pt3[0]-pt2[0])/8.0, e2=(pt3[1]-pt2[1])/8.0, e3=(pt3[2]-pt2[2])/8.0;
 270     double f1=(pt4[0]-pt1[0])/8.0, f2=(pt4[1]-pt1[1])/8.0, f3=(pt4[2]-pt1[2])/8.0;
 271     double h1=(pt7[0]-pt6[0])/8.0, h2=(pt7[1]-pt6[1])/8.0, h3=(pt7[2]-pt6[2])/8.0;
 272     double p1=(pt8[0]-pt5[0])/8.0, p2=(pt8[1]-pt5[1])/8.0, p3=(pt8[2]-pt5[2])/8.0;
 273     double q1=(pt3[0]-pt7[0])/8.0, q2=(pt3[1]-pt7[1])/8.0, q3=(pt3[2]-pt7[2])/8.0;
 274     double r1=(pt4[0]-pt8[0])/8.0, r2=(pt4[1]-pt8[1])/8.0, r3=(pt4[2]-pt8[2])/8.0;
 275     double s1=(pt2[0]-pt6[0])/8.0, s2=(pt2[1]-pt6[1])/8.0, s3=(pt2[2]-pt6[2])/8.0;
 276     double t1=(pt1[0]-pt5[0])/8.0, t2=(pt1[1]-pt5[1])/8.0, t3=(pt1[2]-pt5[2])/8.0;
 277
 278     double A = a1*e2*q3 - a1*e3*q2 - a2*e1*q3 + a2*e3*q1 + a3*e1*q2 - a3*e2*q1;
 279     double B = c1*h2*q3 - c1*h3*q2 - c2*h1*q3 + c2*h3*q1 + c3*h1*q2 - c3*h2*q1;
 280     double C = (a1*h2 + c1*e2)*q3 - (a1*h3 + c1*e3)*q2
 281       - (a2*h1 + c2*e1)*q3 + (a2*h3 + c2*e3)*q1
 282       + (a3*h1 + c3*e1)*q2 - (a3*h2 + c3*e2)*q1;
 283     double D = b1*e2*s3 - b1*e3*s2 - b2*e1*s3 + b2*e3*s1 + b3*e1*s2 - b3*e2*s1;
 284     double E = d1*h2*s3 - d1*h3*s2 - d2*h1*s3 + d2*h3*s1 + d3*h1*s2 - d3*h2*s1;
 285     double F = (b1*h2 + d1*e2)*s3 - (b1*h3 + d1*e3)*s2
 286       - (b2*h1 + d2*e1)*s3 + (b2*h3 + d2*e3)*s1
 287       + (b3*h1 + d3*e1)*s2 - (b3*h2 + d3*e2)*s1;
 288     double G = (a1*e2*s3 + b1*e2*q3) - (a1*e3*s2 + b1*e3*q2)
 289       - (a2*e1*s3 + b2*e1*q3) + (a2*e3*s1 + b2*e3*q1)
 290       + (a3*e1*s2 + b3*e1*q2) - (a3*e2*s1 + b3*e2*q1);
 291     double H = (c1*h2*s3 + d1*h2*q3) - (c1*h3*s2 + d1*h3*q2)
 292       - (c2*h1*s3 + d2*h1*q3) + (c2*h3*s1 + d2*h3*q1)
 293       + (c3*h1*s2 + d3*h1*q2) - (c3*h2*s1 + d3*h2*q1);
 294     double I = ((a1*h2 + c1*e2)*s3 + (b1*h2 + d1*e2)*q3)
 295       - ((a1*h3 + c1*e3)*s2 + (b1*h3 + d1*e3)*q2)
 296       - ((a2*h1 + c2*e1)*s3 + (b2*h1 + d2*e1)*q3)
 297       + ((a2*h3 + c2*e3)*s1 + (b2*h3 + d2*e3)*q1)
 298       + ((a3*h1 + c3*e1)*s2 + (b3*h1 + d3*e1)*q2)
 299       - ((a3*h2 + c3*e2)*s1 + (b3*h2 + d3*e2)*q1);
 300     double J = a1*f2*r3 - a1*f3*r2 - a2*f1*r3 + a2*f3*r1 + a3*f1*r2 - a3*f2*r1;
 301     double K = c1*p2*r3 - c1*p3*r2 - c2*p1*r3 + c2*p3*r1 + c3*p1*r2 - c3*p2*r1;
 302     double L = (a1*p2 + c1*f2)*r3 - (a1*p3 + c1*f3)*r2
 303       - (a2*p1 + c2*f1)*r3 + (a2*p3 + c2*f3)*r1
 304       + (a3*p1 + c3*f1)*r2 - (a3*p2 + c3*f2)*r1;
 305     double M = b1*f2*t3 - b1*f3*t2 - b2*f1*t3 + b2*f3*t1 + b3*f1*t2 - b3*f2*t1;
 306     double N = d1*p2*t3 - d1*p3*t2 - d2*p1*t3 + d2*p3*t1 + d3*p1*t2 - d3*p2*t1;
 307     double O = (b1*p2 + d1*f2)*t3 - (b1*p3 + d1*f3)*t2
 308       - (b2*p1 + d2*f1)*t3 + (b2*p3 + d2*f3)*t1
 309       + (b3*p1 + d3*f1)*t2 - (b3*p2 + d3*f2)*t1;
 310     double P = (a1*f2*t3 + b1*f2*r3) - (a1*f3*t2 + b1*f3*r2)
 311       - (a2*f1*t3 + b2*f1*r3) + (a2*f3*t1 + b2*f3*r1)
 312       + (a3*f1*t2 + b3*f1*r2) - (a3*f2*t1 + b3*f2*r1);
 313     double Q = (c1*p2*t3 + d1*p2*r3) - (c1*p3*t2 + d1*p3*r2)
 314       - (c2*p1*t3 + d2*p1*r3) + (c2*p3*t1 + d2*p3*r1)
 315       + (c3*p1*t2 + d3*p1*r2) - (c3*p2*t1 + d3*p2*r1);
 316     double R = ((a1*p2 + c1*f2)*t3 + (b1*p2 + d1*f2)*r3)
 317       - ((a1*p3 + c1*f3)*t2 + (b1*p3 + d1*f3)*r2)
 318       - ((a2*p1 + c2*f1)*t3 + (b2*p1 + d2*f1)*r3)
 319       + ((a2*p3 + c2*f3)*t1 + (b2*p3 + d2*f3)*r1)
 320       + ((a3*p1 + c3*f1)*t2 + (b3*p1 + d3*f1)*r2)
 321       - ((a3*p2 + c3*f2)*t1 + (b3*p2 + d3*f2)*r1);
 322     double S = (a1*e2*r3 + a1*f2*q3) - (a1*e3*r2 + a1*f3*q2)
 323       - (a2*e1*r3 + a2*f1*q3) + (a2*e3*r1 + a2*f3*q1)
 324       + (a3*e1*r2 + a3*f1*q2) - (a3*e2*r1 + a3*f2*q1);
 325     double T = (c1*h2*r3 + c1*p2*q3) - (c1*h3*r2 + c1*p3*q2)
 326       - (c2*h1*r3 + c2*p1*q3) + (c2*h3*r1 + c2*p3*q1)
 327       + (c3*h1*r2 + c3*p1*q2) - (c3*h2*r1 + c3*p2*q1);
 328     double U = ((a1*h2 + c1*e2)*r3 + (a1*p2 + c1*f2)*q3)
 329       - ((a1*h3 + c1*e3)*r2 + (a1*p3 + c1*f3)*q2)
 330       - ((a2*h1 + c2*e1)*r3 + (a2*p1 + c2*f1)*q3)
 331       + ((a2*h3 + c2*e3)*r1 + (a2*p3 + c2*f3)*q1)
 332       + ((a3*h1 + c3*e1)*r2 + (a3*p1 + c3*f1)*q2)
 333       - ((a3*h2 + c3*e2)*r1 + (a3*p2 + c3*f2)*q1);
 334     double V = (b1*e2*t3 + b1*f2*s3) - (b1*e3*t2 + b1*f3*s2)
 335       - (b2*e1*t3 + b2*f1*s3) + (b2*e3*t1 + b2*f3*s1)
 336       + (b3*e1*t2 + b3*f1*s2) - (b3*e2*t1 + b3*f2*s1);
 337     double W = (d1*h2*t3 + d1*p2*s3) - (d1*h3*t2 + d1*p3*s2)
 338       - (d2*h1*t3 + d2*p1*s3) + (d2*h3*t1 + d2*p3*s1)
 339       + (d3*h1*t2 + d3*p1*s2) - (d3*h2*t1 + d3*p2*s1);
 340     double X = ((b1*h2 + d1*e2)*t3 + (b1*p2 + d1*f2)*s3)
 341       - ((b1*h3 + d1*e3)*t2 + (b1*p3 + d1*f3)*s2)
 342       - ((b2*h1 + d2*e1)*t3 + (b2*p1 + d2*f1)*s3)
 343       + ((b2*h3 + d2*e3)*t1 + (b2*p3 + d2*f3)*s1)
 344       + ((b3*h1 + d3*e1)*t2 + (b3*p1 + d3*f1)*s2)
 345       - ((b3*h2 + d3*e2)*t1 + (b3*p2 + d3*f2)*s1);
 346     double Y = (a1*e2*t3 + a1*f2*s3 + b1*e2*r3 + b1*f2*q3)
 347       - (a1*e3*t2 + a1*f3*s2 + b1*e3*r2 + b1*f3*q2)
 348       - (a2*e1*t3 + a2*f1*s3 + b2*e1*r3 + b2*f1*q3)
 349       + (a2*e3*t1 + a2*f3*s1 + b2*e3*r1 + b2*f3*q1)
 350       + (a3*e1*t2 + a3*f1*s2 + b3*e1*r2 + b3*f1*q2)
 351       - (a3*e2*t1 + a3*f2*s1 + b3*e2*r1 + b3*f2*q1);
 352     double Z = (c1*h2*t3 + c1*p2*s3 + d1*h2*r3 + d1*p2*q3)
 353       - (c1*h3*t2 + c1*p3*s2 + d1*h3*r2 + d1*p3*q2)
 354       - (c2*h1*t3 + c2*p1*s3 + d2*h1*r3 + d2*p1*q3)
 355       + (c2*h3*t1 + c2*p3*s1 + d2*h3*r1 + d2*p3*q1)
 356       + (c3*h1*t2 + c3*p1*s2 + d3*h1*r2 + d3*p1*q2)
 357       - (c3*h2*t1 + c3*p2*s1 + d3*h2*r1 + d3*p2*q1);
 358     double AA = ((a1*h2 + c1*e2)*t3 + (a1*p2 + c1*f2)*s3
 359                  +(b1*h2 + d1*e2)*r3 + (b1*p2 + d1*f2)*q3)
 360       - ((a1*h3 + c1*e3)*t2 + (a1*p3 + c1*f3)*s2
 361          +(b1*h3 + d1*e3)*r2 + (b1*p3 + d1*f3)*q2)
 362       - ((a2*h1 + c2*e1)*t3 + (a2*p1 + c2*f1)*s3
 363          +(b2*h1 + d2*e1)*r3 + (b2*p1 + d2*f1)*q3)
 364       + ((a2*h3 + c2*e3)*t1 + (a2*p3 + c2*f3)*s1
 365          +(b2*h3 + d2*e3)*r1 + (b2*p3 + d2*f3)*q1)
 366       + ((a3*h1 + c3*e1)*t2 + (a3*p1 + c3*f1)*s2
 367          +(b3*h1 + d3*e1)*r2 + (b3*p1 + d3*f1)*q2)
 368       - ((a3*h2 + c3*e2)*t1 + (a3*p2 + c3*f2)*s1
 369          + (b3*h2 + d3*e2)*r1 + (b3*p2 + d3*f2)*q1);
 370
 371     return  64.0*(  8.0*(A + B + D + E + J + K + M + N)
 372                     + 4.0*(C + F + G + H + L + O + P + Q + S + T + V + W)
 373                     + 2.0*(I + R + U + X + Y + Z) + AA ) / 27.0 ;
 374   }
 375
 376   // =========================
 377   // Calculate Volume for Poly
 378   // =========================
 379   inline double calculateVolumeForPolyh(const double ***pts,
 380                                         const int *nbOfNodesPerFaces,
 381                                         int nbOfFaces,
 382                                         const double *bary)
 383   {
 384     double volume=0.;
 385
 386     for ( int i=0; i<nbOfFaces; i++ )
 387       {
 388         double normal[3];
 389         double vecForAlt[3];
 390
 391         calculateNormalForPolyg(pts[i],nbOfNodesPerFaces[i],normal);
 392         vecForAlt[0]=bary[0]-pts[i][0][0];
 393         vecForAlt[1]=bary[1]-pts[i][0][1];
 394         vecForAlt[2]=bary[2]-pts[i][0][2];
 395         volume+=vecForAlt[0]*normal[0]+vecForAlt[1]*normal[1]+vecForAlt[2]*normal[2];
 396       }
 397     return -volume/3.;
 398   }
 399
 400   template<int N>
 401   inline double addComponentsOfVec(const double **pts, int rk)
 402   {
 403     return pts[N-1][rk]+addComponentsOfVec<N-1>(pts,rk);
 404   }
 405
 406   template<>
 407   inline double addComponentsOfVec<1>(const double **pts, int rk)
 408   {
 409     return pts[0][rk];
 410   }
 411
 412   template<int N, int DIM>
 413   inline void calculateBarycenter(const double **pts, double *bary)
 414   {
 415     bary[DIM-1]=addComponentsOfVec<N>(pts,DIM-1)/N;
 416     calculateBarycenter<N,DIM-1>(pts,bary);
 417   }
 418
 419   template<>
 420   inline void calculateBarycenter<2,0>(const double **pts, double *bary)
 421   {
 422   }
 423
 424   template<>
 425   inline void calculateBarycenter<3,0>(const double **pts, double *bary)
 426   {
 427   }
 428
 429   template<>
 430   inline void calculateBarycenter<4,0>(const double **pts, double *bary)
 431   {
 432   }
 433
 434   template<>
 435   inline void calculateBarycenter<5,0>(const double **pts, double *bary)
 436   {
 437   }
 438
 439   template<>
 440   inline void calculateBarycenter<6,0>(const double **pts, double *bary)
 441   {
 442   }
 443
 444   template<>
 445   inline void calculateBarycenter<7,0>(const double **pts, double *bary)
 446   {
 447   }
 448
 449   template<>
 450   inline void calculateBarycenter<8,0>(const double **pts, double *bary)
 451   {
 452   }
 453
 454   inline void calculateBarycenterDyn(const double **pts, int nbPts,
 455                                      int dim, double *bary)
 456   {
 457     for(int i=0;i<dim;i++)
 458       {
 459         double temp=0.;
 460         for(int j=0;j<nbPts;j++)
 461           {
 462             temp+=pts[j][i];
 463           }
 464         bary[i]=temp/nbPts;
 465       }
 466   }
 467 }
 468
 469 #endif