1 // Copyright (C) 2007-2008 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 #ifndef VOLSURFFORMULAE
20 #define VOLSURFFORMULAE
24 namespace INTERP_KERNEL
26 inline void calculateBarycenterDyn(const double **pts, int nbPts,
27 int dim, double *bary);
29 inline double calculateAreaForPolyg(const double **coords, int nbOfPtsInPolygs,
33 // ===========================
34 // Calculate Area for triangle
35 // ===========================
36 inline double calculateAreaForTria(const double *p1, const double *p2,
37 const double *p3, int spaceDim)
43 area = -((p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1]))/2.0;
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]))
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]))
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;
60 // =============================
61 // Calculate Area for quadrangle
62 // =============================
63 inline double calculateAreaForQuad(const double *p1, const double *p2,
64 const double *p3, const double *p4,
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;
76 area = - 4.0*( b1*c2 - c1*b2 + a1*c2 - c1*a2
77 + b1*d2 - d1*b2 + a1*d2 - d1*a2);
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])))
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])))
100 // ====================================
101 // Calculate Normal Vector for Triangle
102 // ====================================
103 inline void calculateNormalForTria(const double *p1, const double *p2,
104 const double *p3, double *normal)
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;
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,
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 ;
131 // ===================================
132 // Calculate Normal Vector for Polygon
133 // ===================================
134 inline void calculateNormalForPolyg(const double **coords, int nbOfPtsInPolygs,
137 double coordOfBary[3];
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]);
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]);
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]);
149 double xarea=sqrt(xnormal1*xnormal1 + xnormal2*xnormal2 + xnormal3*xnormal3);
153 //std::string diagnosis"Can not calculate normal - polygon is singular";
154 throw std::exception();
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 ;
166 // ==========================
167 // Calculate Area for Polygon
168 // ==========================
169 inline double calculateAreaForPolyg(const double **coords, int nbOfPtsInPolygs,
173 double coordOfBary[3];
175 calculateBarycenterDyn(coords,nbOfPtsInPolygs,spaceDim,coordOfBary);
176 for ( int i=0; i<nbOfPtsInPolygs; i++ )
178 double tmp = calculateAreaForTria(coords[i],coords[(i+1)%nbOfPtsInPolygs],
179 coordOfBary,spaceDim);
185 // ==========================
186 // Calculate Volume for Tetra
187 // ==========================
188 inline double calculateVolumeForTetra(const double *p1, const double *p2,
189 const double *p3, const double *p4)
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]) )
200 // =========================
201 // Calculate Volume for Pyra
202 // =========================
203 inline double calculateVolumeForPyra(const double *p1, const double *p2,
204 const double *p3, const double *p4,
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]) ))
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]) ))
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)
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;
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);
254 return (-2.0*(2.0*(A + B + D + E + G + H) + C + F + P)/9.0);
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)
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;
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);
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 ;
376 // =========================
377 // Calculate Volume for Poly
378 // =========================
379 inline double calculateVolumeForPolyh(const double ***pts,
380 const int *nbOfNodesPerFaces,
386 for ( int i=0; i<nbOfFaces; i++ )
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];
401 inline double addComponentsOfVec(const double **pts, int rk)
403 return pts[N-1][rk]+addComponentsOfVec<N-1>(pts,rk);
407 inline double addComponentsOfVec<1>(const double **pts, int rk)
412 template<int N, int DIM>
413 inline void calculateBarycenter(const double **pts, double *bary)
415 bary[DIM-1]=addComponentsOfVec<N>(pts,DIM-1)/N;
416 calculateBarycenter<N,DIM-1>(pts,bary);
420 inline void calculateBarycenter<2,0>(const double **pts, double *bary)
425 inline void calculateBarycenter<3,0>(const double **pts, double *bary)
430 inline void calculateBarycenter<4,0>(const double **pts, double *bary)
435 inline void calculateBarycenter<5,0>(const double **pts, double *bary)
440 inline void calculateBarycenter<6,0>(const double **pts, double *bary)
445 inline void calculateBarycenter<7,0>(const double **pts, double *bary)
450 inline void calculateBarycenter<8,0>(const double **pts, double *bary)
454 inline void calculateBarycenterDyn(const double **pts, int nbPts,
455 int dim, double *bary)
457 for(int i=0;i<dim;i++)
460 for(int j=0;j<nbPts;j++)