src/MEDMEM/MEDMEM_Formulae.hxx

   1 #ifndef MEDMEM_FORMULAE
   2 #define MEDMEM_FORMULAE
   3
   4 #include <math.h>
   5
   6 inline void CalculateBarycenterDyn(const double **pts, int nbPts, int dim, double *bary);
   7
   8 inline double CalculateAreaForPolyg(const double **coords, int nbOfPtsInPolygs, int spaceDim);
   9
  10 inline double CalculateAreaForTria(const double *p1, const double *p2, const double *p3,int spaceDim)
  11 {
  12   if (spaceDim==2)
  13     return (-((p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1]))/2.0);
  14   else
  15     return (sqrt(((p2[1]-p1[1])*(p3[2]-p1[2]) - (p3[1]-p1[1])*(p2[2]-p1[2]))*
  16                  ((p2[1]-p1[1])*(p3[2]-p1[2]) - (p3[1]-p1[1])*(p2[2]-p1[2])) +
  17                  ((p3[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p3[2]-p1[2]))*
  18                  ((p3[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p3[2]-p1[2])) +
  19                  ((p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1]))*
  20                  ((p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1])))/2.0);
  21 }
  22
  23 inline double CalculateAreaForQuad(const double *p1, const double *p2, const double *p3, const double *p4, int spaceDim)
  24 {
  25   if (spaceDim==2)
  26     {
  27       double a1 = (p2[0]-p1[0])/4.0, a2 = (p2[1]-p1[1])/4.0;
  28       double b1 = (p3[0]-p4[0])/4.0, b2 = (p3[1]-p4[1])/4.0;
  29       double c1 = (p3[0]-p2[0])/4.0, c2 = (p3[1]-p2[1])/4.0;
  30       double d1 = (p4[0]-p1[0])/4.0, d2 = (p4[1]-p1[1])/4.0;
  31
  32       return  (- 4.0*(b1*c2 - c1*b2 + a1*c2 - c1*a2 + b1*d2 -
  33                       d1*b2 + a1*d2 - d1*a2));
  34     }
  35   else
  36     {
  37       return ((sqrt(((p2[1]-p1[1])*(p4[2]-p1[2]) - (p4[1]-p1[1])*(p2[2]-p1[2]))*
  38                    ((p2[1]-p1[1])*(p4[2]-p1[2]) - (p4[1]-p1[1])*(p2[2]-p1[2])) +
  39                    ((p4[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p4[2]-p1[2]))*
  40                    ((p4[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p4[2]-p1[2])) +
  41                    ((p2[0]-p1[0])*(p4[1]-p1[1]) - (p4[0]-p1[0])*(p2[1]-p1[1]))*
  42                    ((p2[0]-p1[0])*(p4[1]-p1[1]) - (p4[0]-p1[0])*(p2[1]-p1[1]))) +
  43               sqrt(((p4[1]-p3[1])*(p2[2]-p3[2]) - (p2[1]-p3[1])*(p4[2]-p3[2]))*
  44                    ((p4[1]-p3[1])*(p2[2]-p3[2]) - (p2[1]-p3[1])*(p4[2]-p3[2])) +
  45                    ((p2[0]-p3[0])*(p4[2]-p3[2]) - (p4[0]-p3[0])*(p2[2]-p3[2]))*
  46                    ((p2[0]-p3[0])*(p4[2]-p3[2]) - (p4[0]-p3[0])*(p2[2]-p3[2])) +
  47                    ((p4[0]-p3[0])*(p2[1]-p3[1]) - (p2[0]-p3[0])*(p4[1]-p3[1]))*
  48                    ((p4[0]-p3[0])*(p2[1]-p3[1]) - (p2[0]-p3[0])*(p4[1]-p3[1]))))/2.0);
  49     }
  50 }
  51
  52 inline void CalculateNormalForTria(const double *p1, const double *p2, const double *p3, double *normal)
  53 {
  54   normal[0] = ((p2[1]-p1[1])*(p3[2]-p1[2]) - (p3[1]-p1[1])*(p2[2]-p1[2]))/2.0;
  55   normal[1] = ((p3[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p3[2]-p1[2]))/2.0;
  56   normal[2] = ((p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1]))/2.0;
  57 }
  58
  59 inline void CalculateNormalForQuad(const double *p1, const double *p2, const double *p3, const double *p4, double *normal)
  60 {
  61   double xnormal1 = (p2[1]-p1[1])*(p4[2]-p1[2]) - (p4[1]-p1[1])*(p2[2]-p1[2]);
  62   double xnormal2 = (p4[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p4[2]-p1[2]);
  63   double xnormal3 = (p2[0]-p1[0])*(p4[1]-p1[1]) - (p4[0]-p1[0])*(p2[1]-p1[1]);
  64   double xarea = sqrt(xnormal1*xnormal1 + xnormal2*xnormal2 + xnormal3*xnormal3);
  65   xnormal1 = xnormal1/xarea;
  66   xnormal2 = xnormal2/xarea;
  67   xnormal3 = xnormal3/xarea;
  68   xarea = CalculateAreaForQuad(p1,p2,p3,p4,3);
  69   normal[0] = xnormal1*xarea ;
  70   normal[1] = xnormal2*xarea ;
  71   normal[2] = xnormal3*xarea ;
  72 }
  73
  74 inline void CalculateNormalForPolyg(const double **coords, int nbOfPtsInPolygs, double *normal)
  75 {
  76   double coordOfBary[3];
  77   CalculateBarycenterDyn(coords,nbOfPtsInPolygs,3,coordOfBary);
  78   double xnormal1 = (coords[0][1]-coords[1][1])*(coordOfBary[2]-coords[1][2])-(coords[0][2]-coords[1][2])*(coordOfBary[1]-coords[1][1]);
  79   double xnormal2 = (coords[0][2]-coords[1][2])*(coordOfBary[0]-coords[1][0])-(coords[0][0]-coords[1][0])*(coordOfBary[2]-coords[1][2]);
  80   double xnormal3 = (coords[0][0]-coords[1][0])*(coordOfBary[1]-coords[1][1])-(coords[0][1]-coords[1][1])*(coordOfBary[0]-coords[1][0]);
  81   double xarea = sqrt(xnormal1*xnormal1 + xnormal2*xnormal2 + xnormal3*xnormal3);
  82   xnormal1 = xnormal1/xarea;
  83   xnormal2 = xnormal2/xarea;
  84   xnormal3 = xnormal3/xarea;
  85   xarea = CalculateAreaForPolyg(coords,nbOfPtsInPolygs,3);
  86   normal[0] = xnormal1*xarea ;
  87   normal[1] = xnormal2*xarea ;
  88   normal[2] = xnormal3*xarea ;
  89 }
  90
  91 inline double CalculateAreaForPolyg(const double **coords, int nbOfPtsInPolygs, int spaceDim)
  92 {
  93   double coordOfBary[3];
  94   CalculateBarycenterDyn(coords,nbOfPtsInPolygs,spaceDim,coordOfBary);
  95   double ret=0.;
  96   for(int i=0;i<nbOfPtsInPolygs;i++)
  97     {
  98       double tmp=CalculateAreaForTria(coords[i],coords[(i+1)%nbOfPtsInPolygs],coordOfBary,spaceDim);
  99       ret+=tmp;
 100     }
 101   return ret;
 102 }
 103
 104 inline double CalculateVolumeForTetra(const double *p1, const double *p2, const double *p3, const double *p4)
 105 {
 106   return  ((p3[0]-p1[0])*((p2[1]-p1[1])*(p4[2]-p1[2]) - (p2[2]-p1[2])*(p4[1]-p1[1])) -
 107            (p2[0]-p1[0])*((p3[1]-p1[1])*(p4[2]-p1[2]) - (p3[2]-p1[2])*(p4[1]-p1[1])) +
 108            (p4[0]-p1[0])*((p3[1]-p1[1])*(p2[2]-p1[2]) - (p3[2]-p1[2])*(p2[1]-p1[1])))/6.0;
 109 }
 110
 111 inline double CalculateVolumeForPyra(const double *p1, const double *p2, const double *p3, const double *p4, const double *p5)
 112 {
 113   return (((p3[0]-p1[0])*((p2[1]-p1[1])*(p5[2]-p1[2]) - (p2[2]-p1[2])*(p5[1]-p1[1])) -
 114            (p2[0]-p1[0])*((p3[1]-p1[1])*(p5[2]-p1[2]) - (p3[2]-p1[2])*(p5[1]-p1[1])) +
 115            (p5[0]-p1[0])*((p3[1]-p1[1])*(p2[2]-p1[2]) - (p3[2]-p1[2])*(p2[1]-p1[1]))) +
 116           ((p4[0]-p1[0])*((p3[1]-p1[1])*(p5[2]-p1[2]) - (p3[2]-p1[2])*(p5[1]-p1[1])) -
 117            (p3[0]-p1[0])*((p4[1]-p1[1])*(p5[2]-p1[2]) - (p4[2]-p1[2])*(p5[1]-p1[1])) +
 118            (p5[0]-p1[0])*((p4[1]-p1[1])*(p3[2]-p1[2]) - (p4[2]-p1[2])*(p3[1]-p1[1])))
 119           )/6.0;
 120 }
 121
 122 inline double CalculateVolumeForPenta(const double *p1, const double *p2, const double *p3, const double *p4, const double *p5, const double *p6)
 123 {
 124   double a1 = (p2[0]-p3[0])/2.0, a2 = (p2[1]-p3[1])/2.0, a3 = (p2[2]-p3[2])/2.0;
 125   double b1 = (p5[0]-p6[0])/2.0, b2 = (p5[1]-p6[1])/2.0, b3 = (p5[2]-p6[2])/2.0;
 126   double c1 = (p4[0]-p1[0])/2.0, c2 = (p4[1]-p1[1])/2.0, c3 = (p4[2]-p1[2])/2.0;
 127   double d1 = (p5[0]-p2[0])/2.0, d2 = (p5[1]-p2[1])/2.0, d3 = (p5[2]-p2[2])/2.0;
 128   double e1 = (p6[0]-p3[0])/2.0, e2 = (p6[1]-p3[1])/2.0, e3 = (p6[2]-p3[2])/2.0;
 129   double f1 = (p1[0]-p3[0])/2.0, f2 = (p1[1]-p3[1])/2.0, f3 = (p1[2]-p3[2])/2.0;
 130   double h1 = (p4[0]-p6[0])/2.0, h2 = (p4[1]-p6[1])/2.0, h3 = (p4[2]-p6[2])/2.0;
 131
 132   double A = a1*c2*f3 - a1*c3*f2 - a2*c1*f3 + a2*c3*f1 +
 133     a3*c1*f2 - a3*c2*f1;
 134   double B = b1*c2*h3 - b1*c3*h2 - b2*c1*h3 + b2*c3*h1 +
 135     b3*c1*h2 - b3*c2*h1;
 136   double C = (a1*c2*h3 + b1*c2*f3) - (a1*c3*h2 + b1*c3*f2) -
 137     (a2*c1*h3 + b2*c1*f3) + (a2*c3*h1 + b2*c3*f1) +
 138     (a3*c1*h2 + b3*c1*f2) - (a3*c2*h1 + b3*c2*f1);
 139   double D = a1*d2*f3 - a1*d3*f2 - a2*d1*f3 + a2*d3*f1 +
 140     a3*d1*f2 - a3*d2*f1;
 141   double E = b1*d2*h3 - b1*d3*h2 - b2*d1*h3 + b2*d3*h1 +
 142     b3*d1*h2 - b3*d2*h1;
 143   double F = (a1*d2*h3 + b1*d2*f3) - (a1*d3*h2 + b1*d3*f2) -
 144     (a2*d1*h3 + b2*d1*f3) + (a2*d3*h1 + b2*d3*f1) +
 145     (a3*d1*h2 + b3*d1*f2) - (a3*d2*h1 + b3*d2*f1);
 146   double G = a1*e2*f3 - a1*e3*f2 - a2*e1*f3 + a2*e3*f1 +
 147     a3*e1*f2 - a3*e2*f1;
 148   double H = b1*e2*h3 - b1*e3*h2 - b2*e1*h3 + b2*e3*h1 +
 149     b3*e1*h2 - b3*e2*h1;
 150   double P = (a1*e2*h3 + b1*e2*f3) - (a1*e3*h2 + b1*e3*f2) -
 151     (a2*e1*h3 + b2*e1*f3) + (a2*e3*h1 + b2*e3*f1) +
 152     (a3*e1*h2 + b3*e1*f2) - (a3*e2*h1 + b3*e2*f1);
 153
 154   return (-2.0*(2.0*(A + B + D + E + G + H) + C + F + P)/9.0);
 155 }
 156
 157 inline double CalculateVolumeForHexa(const double *pt1, const double *pt2, const double *pt3, const double *pt4, const double *pt5, const double *pt6, const double *pt7, const double *pt8)
 158 {
 159   double a1 = (pt3[0]-pt4[0])/8.0, a2 = (pt3[1]-pt4[1])/8.0, a3 = (pt3[2]-pt4[2])/8.0;
 160   double b1 = (pt2[0]-pt1[0])/8.0, b2 = (pt2[1]-pt1[1])/8.0, b3 = (pt2[2]-pt1[2])/8.0;
 161   double c1 = (pt7[0]-pt8[0])/8.0, c2 = (pt7[1]-pt8[1])/8.0, c3 = (pt7[2]-pt8[2])/8.0;
 162   double d1 = (pt6[0]-pt5[0])/8.0, d2 = (pt6[1]-pt5[1])/8.0, d3 = (pt6[2]-pt5[2])/8.0;
 163   double e1 = (pt3[0]-pt2[0])/8.0, e2 = (pt3[1]-pt2[1])/8.0, e3 = (pt3[2]-pt2[2])/8.0;
 164   double f1 = (pt4[0]-pt1[0])/8.0, f2 = (pt4[1]-pt1[1])/8.0, f3 = (pt4[2]-pt1[2])/8.0;
 165   double h1 = (pt7[0]-pt6[0])/8.0, h2 = (pt7[1]-pt6[1])/8.0, h3 = (pt7[2]-pt6[2])/8.0;
 166   double p1 = (pt8[0]-pt5[0])/8.0, p2 = (pt8[1]-pt5[1])/8.0, p3 = (pt8[2]-pt5[2])/8.0;
 167   double q1 = (pt3[0]-pt7[0])/8.0, q2 = (pt3[1]-pt7[1])/8.0, q3 = (pt3[2]-pt7[2])/8.0;
 168   double r1 = (pt4[0]-pt8[0])/8.0, r2 = (pt4[1]-pt8[1])/8.0, r3 = (pt4[2]-pt8[2])/8.0;
 169   double s1 = (pt2[0]-pt6[0])/8.0, s2 = (pt2[1]-pt6[1])/8.0, s3 = (pt2[2]-pt6[2])/8.0;
 170   double t1 = (pt1[0]-pt5[0])/8.0, t2 = (pt1[1]-pt5[1])/8.0, t3 = (pt1[2]-pt5[2])/8.0;
 171
 172   double A = a1*e2*q3 - a1*e3*q2 - a2*e1*q3 + a2*e3*q1 +
 173     a3*e1*q2 - a3*e2*q1;
 174   double B = c1*h2*q3 - c1*h3*q2 - c2*h1*q3 + c2*h3*q1 +
 175     c3*h1*q2 - c3*h2*q1;
 176   double C = (a1*h2 + c1*e2)*q3 - (a1*h3 + c1*e3)*q2 -
 177     (a2*h1 + c2*e1)*q3 + (a2*h3 + c2*e3)*q1 +
 178     (a3*h1 + c3*e1)*q2 - (a3*h2 + c3*e2)*q1;
 179   double D = b1*e2*s3 - b1*e3*s2 - b2*e1*s3 + b2*e3*s1 +
 180     b3*e1*s2 - b3*e2*s1;
 181   double E = d1*h2*s3 - d1*h3*s2 - d2*h1*s3 + d2*h3*s1 +
 182     d3*h1*s2 - d3*h2*s1;
 183   double F = (b1*h2 + d1*e2)*s3 - (b1*h3 + d1*e3)*s2 -
 184     (b2*h1 + d2*e1)*s3 + (b2*h3 + d2*e3)*s1 +
 185     (b3*h1 + d3*e1)*s2 - (b3*h2 + d3*e2)*s1;
 186   double G = (a1*e2*s3 + b1*e2*q3) - (a1*e3*s2 + b1*e3*q2) -
 187     (a2*e1*s3 + b2*e1*q3) + (a2*e3*s1 + b2*e3*q1) +
 188     (a3*e1*s2 + b3*e1*q2) - (a3*e2*s1 + b3*e2*q1);
 189   double H = (c1*h2*s3 + d1*h2*q3) - (c1*h3*s2 + d1*h3*q2) -
 190     (c2*h1*s3 + d2*h1*q3) + (c2*h3*s1 + d2*h3*q1) +
 191     (c3*h1*s2 + d3*h1*q2) - (c3*h2*s1 + d3*h2*q1);
 192   double I = ((a1*h2 + c1*e2)*s3 + (b1*h2 + d1*e2)*q3) -
 193     ((a1*h3 + c1*e3)*s2 + (b1*h3 + d1*e3)*q2) -
 194     ((a2*h1 + c2*e1)*s3 + (b2*h1 + d2*e1)*q3) +
 195     ((a2*h3 + c2*e3)*s1 + (b2*h3 + d2*e3)*q1) +
 196     ((a3*h1 + c3*e1)*s2 + (b3*h1 + d3*e1)*q2) -
 197     ((a3*h2 + c3*e2)*s1 + (b3*h2 + d3*e2)*q1);
 198   double J = a1*f2*r3 - a1*f3*r2 - a2*f1*r3 + a2*f3*r1 +
 199     a3*f1*r2 - a3*f2*r1;
 200   double K = c1*p2*r3 - c1*p3*r2 - c2*p1*r3 + c2*p3*r1 +
 201     c3*p1*r2 - c3*p2*r1;
 202   double L = (a1*p2 + c1*f2)*r3 - (a1*p3 + c1*f3)*r2 -
 203     (a2*p1 + c2*f1)*r3 + (a2*p3 + c2*f3)*r1 +
 204     (a3*p1 + c3*f1)*r2 - (a3*p2 + c3*f2)*r1;
 205   double M = b1*f2*t3 - b1*f3*t2 - b2*f1*t3 + b2*f3*t1 +
 206     b3*f1*t2 - b3*f2*t1;
 207   double N = d1*p2*t3 - d1*p3*t2 - d2*p1*t3 + d2*p3*t1 +
 208     d3*p1*t2 - d3*p2*t1;
 209   double O = (b1*p2 + d1*f2)*t3 - (b1*p3 + d1*f3)*t2 -
 210     (b2*p1 + d2*f1)*t3 + (b2*p3 + d2*f3)*t1 +
 211     (b3*p1 + d3*f1)*t2 - (b3*p2 + d3*f2)*t1;
 212   double P = (a1*f2*t3 + b1*f2*r3) - (a1*f3*t2 + b1*f3*r2) -
 213     (a2*f1*t3 + b2*f1*r3) + (a2*f3*t1 + b2*f3*r1) +
 214     (a3*f1*t2 + b3*f1*r2) - (a3*f2*t1 + b3*f2*r1);
 215   double Q = (c1*p2*t3 + d1*p2*r3) - (c1*p3*t2 + d1*p3*r2) -
 216     (c2*p1*t3 + d2*p1*r3) + (c2*p3*t1 + d2*p3*r1) +
 217     (c3*p1*t2 + d3*p1*r2) - (c3*p2*t1 + d3*p2*r1);
 218   double R = ((a1*p2 + c1*f2)*t3 + (b1*p2 + d1*f2)*r3) -
 219     ((a1*p3 + c1*f3)*t2 + (b1*p3 + d1*f3)*r2) -
 220     ((a2*p1 + c2*f1)*t3 + (b2*p1 + d2*f1)*r3) +
 221     ((a2*p3 + c2*f3)*t1 + (b2*p3 + d2*f3)*r1) +
 222     ((a3*p1 + c3*f1)*t2 + (b3*p1 + d3*f1)*r2) -
 223     ((a3*p2 + c3*f2)*t1 + (b3*p2 + d3*f2)*r1);
 224   double S = (a1*e2*r3 + a1*f2*q3) - (a1*e3*r2 + a1*f3*q2) -
 225     (a2*e1*r3 + a2*f1*q3) + (a2*e3*r1 + a2*f3*q1) +
 226     (a3*e1*r2 + a3*f1*q2) - (a3*e2*r1 + a3*f2*q1);
 227   double T = (c1*h2*r3 + c1*p2*q3) - (c1*h3*r2 + c1*p3*q2) -
 228     (c2*h1*r3 + c2*p1*q3) + (c2*h3*r1 + c2*p3*q1) +
 229     (c3*h1*r2 + c3*p1*q2) - (c3*h2*r1 + c3*p2*q1);
 230   double U = ((a1*h2 + c1*e2)*r3 + (a1*p2 + c1*f2)*q3) -
 231     ((a1*h3 + c1*e3)*r2 + (a1*p3 + c1*f3)*q2) -
 232     ((a2*h1 + c2*e1)*r3 + (a2*p1 + c2*f1)*q3) +
 233     ((a2*h3 + c2*e3)*r1 + (a2*p3 + c2*f3)*q1) +
 234     ((a3*h1 + c3*e1)*r2 + (a3*p1 + c3*f1)*q2) -
 235     ((a3*h2 + c3*e2)*r1 + (a3*p2 + c3*f2)*q1);
 236   double V = (b1*e2*t3 + b1*f2*s3) - (b1*e3*t2 + b1*f3*s2) -
 237     (b2*e1*t3 + b2*f1*s3) + (b2*e3*t1 + b2*f3*s1) +
 238     (b3*e1*t2 + b3*f1*s2) - (b3*e2*t1 + b3*f2*s1);
 239   double W = (d1*h2*t3 + d1*p2*s3) - (d1*h3*t2 + d1*p3*s2) -
 240     (d2*h1*t3 + d2*p1*s3) + (d2*h3*t1 + d2*p3*s1) +
 241     (d3*h1*t2 + d3*p1*s2) - (d3*h2*t1 + d3*p2*s1);
 242   double X = ((b1*h2 + d1*e2)*t3 + (b1*p2 + d1*f2)*s3) -
 243     ((b1*h3 + d1*e3)*t2 + (b1*p3 + d1*f3)*s2) -
 244     ((b2*h1 + d2*e1)*t3 + (b2*p1 + d2*f1)*s3) +
 245     ((b2*h3 + d2*e3)*t1 + (b2*p3 + d2*f3)*s1) +
 246     ((b3*h1 + d3*e1)*t2 + (b3*p1 + d3*f1)*s2) -
 247     ((b3*h2 + d3*e2)*t1 + (b3*p2 + d3*f2)*s1);
 248   double Y = (a1*e2*t3 + a1*f2*s3 + b1*e2*r3 + b1*f2*q3) -
 249     (a1*e3*t2 + a1*f3*s2 + b1*e3*r2 + b1*f3*q2) -
 250     (a2*e1*t3 + a2*f1*s3 + b2*e1*r3 + b2*f1*q3) +
 251     (a2*e3*t1 + a2*f3*s1 + b2*e3*r1 + b2*f3*q1) +
 252     (a3*e1*t2 + a3*f1*s2 + b3*e1*r2 + b3*f1*q2) -
 253     (a3*e2*t1 + a3*f2*s1 + b3*e2*r1 + b3*f2*q1);
 254   double Z = (c1*h2*t3 + c1*p2*s3 + d1*h2*r3 + d1*p2*q3) -
 255     (c1*h3*t2 + c1*p3*s2 + d1*h3*r2 + d1*p3*q2) -
 256     (c2*h1*t3 + c2*p1*s3 + d2*h1*r3 + d2*p1*q3) +
 257     (c2*h3*t1 + c2*p3*s1 + d2*h3*r1 + d2*p3*q1) +
 258     (c3*h1*t2 + c3*p1*s2 + d3*h1*r2 + d3*p1*q2) -
 259     (c3*h2*t1 + c3*p2*s1 + d3*h2*r1 + d3*p2*q1);
 260   double AA = ((a1*h2 + c1*e2)*t3 + (a1*p2 + c1*f2)*s3 +
 261                (b1*h2 + d1*e2)*r3 + (b1*p2 + d1*f2)*q3) -
 262     ((a1*h3 + c1*e3)*t2 + (a1*p3 + c1*f3)*s2 +
 263      (b1*h3 + d1*e3)*r2 + (b1*p3 + d1*f3)*q2) -
 264     ((a2*h1 + c2*e1)*t3 + (a2*p1 + c2*f1)*s3 +
 265      (b2*h1 + d2*e1)*r3 + (b2*p1 + d2*f1)*q3) +
 266     ((a2*h3 + c2*e3)*t1 + (a2*p3 + c2*f3)*s1 +
 267      (b2*h3 + d2*e3)*r1 + (b2*p3 + d2*f3)*q1) +
 268     ((a3*h1 + c3*e1)*t2 + (a3*p1 + c3*f1)*s2 +
 269      (b3*h1 + d3*e1)*r2 + (b3*p1 + d3*f1)*q2) -
 270     ((a3*h2 + c3*e2)*t1 + (a3*p2 + c3*f2)*s1 +
 271      (b3*h2 + d3*e2)*r1 + (b3*p2 + d3*f2)*q1);
 272
 273   return  64.0*(8.0*(A + B + D + E + J + K + M + N) +
 274                 4.0*(C + F + G + H + L + O + P + Q + S + T +
 275                      V + W) + 2.0*(I + R + U + X + Y + Z) +
 276                 AA)/27.0;
 277 }
 278
 279 inline double CalculateVolumeForPolyh(const double ***pts, const int *nbOfNodesPerFaces, int nbOfFaces, const double *bary)
 280 {
 281   double volume=0.;
 282   for(int i=0;i<nbOfFaces;i++)
 283     {
 284       double normal[3];
 285       CalculateNormalForPolyg(pts[i],nbOfNodesPerFaces[i],normal);
 286       double vecForAlt[3];
 287       vecForAlt[0]=bary[0]-pts[i][0][0];
 288       vecForAlt[1]=bary[1]-pts[i][0][1];
 289       vecForAlt[2]=bary[2]-pts[i][0][2];
 290       volume+=fabs(vecForAlt[0]*normal[0]+vecForAlt[1]*normal[1]+vecForAlt[2]*normal[2]);
 291     }
 292   return volume/3.;
 293 }
 294
 295 template<int N>
 296 inline double addComponentsOfVec(const double **pts, int rk)
 297 {
 298   return pts[N-1][rk]+addComponentsOfVec<N-1>(pts,rk);
 299 }
 300
 301 template<>
 302 inline double addComponentsOfVec<1>(const double **pts, int rk)
 303 {
 304   return pts[0][rk];
 305 }
 306
 307 template<int N, int DIM>
 308 inline void CalculateBarycenter(const double **pts, double *bary)
 309 {
 310   bary[DIM-1]=addComponentsOfVec<N>(pts,DIM-1)/N;
 311   CalculateBarycenter<N,DIM-1>(pts,bary);
 312 }
 313
 314 template<>
 315 inline void CalculateBarycenter<2,0>(const double **pts, double *bary) { }
 316
 317 template<>
 318 inline void CalculateBarycenter<3,0>(const double **pts, double *bary) { }
 319
 320 template<>
 321 inline void CalculateBarycenter<4,0>(const double **pts, double *bary) { }
 322
 323 template<>
 324 inline void CalculateBarycenter<5,0>(const double **pts, double *bary) { }
 325
 326 template<>
 327 inline void CalculateBarycenter<6,0>(const double **pts, double *bary) { }
 328
 329 template<>
 330 inline void CalculateBarycenter<7,0>(const double **pts, double *bary) { }
 331
 332 template<>
 333 inline void CalculateBarycenter<8,0>(const double **pts, double *bary) { }
 334
 335 inline void CalculateBarycenterDyn(const double **pts, int nbPts, int dim, double *bary)
 336 {
 337   for(int i=0;i<dim;i++)
 338     {
 339       double temp=0.;
 340       for(int j=0;j<nbPts;j++)
 341         {
 342           temp+=pts[j][i];
 343         }
 344       bary[i]=temp/nbPts;
 345     }
 346 }
 347
 348 #endif