1 // Copyright (C) 2007-2016 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, or (at your option) any later version.
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 // Author : Anthony Geay (CEA/DEN)
21 #ifndef __VOLSURFFORMULAE_HXX__
22 #define __VOLSURFFORMULAE_HXX__
24 #include "InterpolationUtils.hxx"
25 #include "InterpKernelException.hxx"
26 #include "InterpKernelGeo2DEdgeLin.hxx"
27 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
28 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
33 namespace INTERP_KERNEL
35 inline void calculateBarycenterDyn(const double **pts, int nbPts,
36 int dim, double *bary);
38 inline double calculateAreaForPolyg(const double **coords, int nbOfPtsInPolygs,
42 inline double calculateAreaForQPolyg(const double **coords, int nbOfPtsInPolygs,
45 inline double calculateLgthForSeg2(const double *p1, const double *p2, int spaceDim)
52 for(int i=0;i<spaceDim;i++)
53 ret+=(p2[i]-p1[i])*(p2[i]-p1[i]);
58 inline double calculateLgthForSeg3(const double *begin, const double *end, const double *middle, int spaceDim)
62 Edge *ed=Edge::BuildEdgeFrom3Points(begin,middle,end);
63 double ret=ed->getCurveLength(); ed->decrRef();
67 return calculateLgthForSeg2(begin,end,spaceDim);
70 // ===========================
71 // Calculate Area for triangle
72 // ===========================
73 inline double calculateAreaForTria(const double *p1, const double *p2,
74 const double *p3, int spaceDim)
80 area = -((p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1]))/2.0;
84 area = sqrt(((p2[1]-p1[1])*(p3[2]-p1[2]) - (p3[1]-p1[1])*(p2[2]-p1[2]))*
85 ((p2[1]-p1[1])*(p3[2]-p1[2]) - (p3[1]-p1[1])*(p2[2]-p1[2]))
87 ((p3[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p3[2]-p1[2]))*
88 ((p3[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p3[2]-p1[2]))
90 ((p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1]))*
91 ((p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1])))/2.0;
97 // =============================
98 // Calculate Area for quadrangle
99 // =============================
100 inline double calculateAreaForQuad(const double *p1, const double *p2,
101 const double *p3, const double *p4,
108 double a1 = (p2[0]-p1[0])/4.0, a2 = (p2[1]-p1[1])/4.0;
109 double b1 = (p3[0]-p4[0])/4.0, b2 = (p3[1]-p4[1])/4.0;
110 double c1 = (p3[0]-p2[0])/4.0, c2 = (p3[1]-p2[1])/4.0;
111 double d1 = (p4[0]-p1[0])/4.0, d2 = (p4[1]-p1[1])/4.0;
113 area = - 4.0*( b1*c2 - c1*b2 + a1*c2 - c1*a2
114 + b1*d2 - d1*b2 + a1*d2 - d1*a2);
118 area = (sqrt(((p2[1]-p1[1])*(p4[2]-p1[2]) - (p4[1]-p1[1])*(p2[2]-p1[2]))*
119 ((p2[1]-p1[1])*(p4[2]-p1[2]) - (p4[1]-p1[1])*(p2[2]-p1[2]))
120 + ((p4[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p4[2]-p1[2]))*
121 ((p4[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p4[2]-p1[2]))
122 + ((p2[0]-p1[0])*(p4[1]-p1[1]) - (p4[0]-p1[0])*(p2[1]-p1[1]))*
123 ((p2[0]-p1[0])*(p4[1]-p1[1]) - (p4[0]-p1[0])*(p2[1]-p1[1])))
125 sqrt(((p4[1]-p3[1])*(p2[2]-p3[2]) - (p2[1]-p3[1])*(p4[2]-p3[2]))*
126 ((p4[1]-p3[1])*(p2[2]-p3[2]) - (p2[1]-p3[1])*(p4[2]-p3[2]))
127 + ((p2[0]-p3[0])*(p4[2]-p3[2]) - (p4[0]-p3[0])*(p2[2]-p3[2]))*
128 ((p2[0]-p3[0])*(p4[2]-p3[2]) - (p4[0]-p3[0])*(p2[2]-p3[2]))
129 + ((p4[0]-p3[0])*(p2[1]-p3[1]) - (p2[0]-p3[0])*(p4[1]-p3[1]))*
130 ((p4[0]-p3[0])*(p2[1]-p3[1]) - (p2[0]-p3[0])*(p4[1]-p3[1])))
137 // ====================================
138 // Calculate Normal Vector for Triangle
139 // ====================================
140 inline void calculateNormalForTria(const double *p1, const double *p2,
141 const double *p3, double *normal)
143 normal[0] = ((p2[1]-p1[1])*(p3[2]-p1[2]) - (p3[1]-p1[1])*(p2[2]-p1[2]))/2.0;
144 normal[1] = ((p3[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p3[2]-p1[2]))/2.0;
145 normal[2] = ((p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1]))/2.0;
148 // ======================================
149 // Calculate Normal Vector for Quadrangle
150 // ======================================
151 inline void calculateNormalForQuad(const double *p1, const double *p2,
152 const double *p3, const double *p4,
155 double xnormal1 = (p2[1]-p1[1])*(p4[2]-p1[2]) - (p4[1]-p1[1])*(p2[2]-p1[2]);
156 double xnormal2 = (p4[0]-p1[0])*(p2[2]-p1[2]) - (p2[0]-p1[0])*(p4[2]-p1[2]);
157 double xnormal3 = (p2[0]-p1[0])*(p4[1]-p1[1]) - (p4[0]-p1[0])*(p2[1]-p1[1]);
158 double xarea = sqrt(xnormal1*xnormal1 + xnormal2*xnormal2 + xnormal3*xnormal3);
159 xnormal1 = xnormal1/xarea;
160 xnormal2 = xnormal2/xarea;
161 xnormal3 = xnormal3/xarea;
162 xarea = calculateAreaForQuad(p1,p2,p3,p4,3);
163 normal[0] = xnormal1*xarea ;
164 normal[1] = xnormal2*xarea ;
165 normal[2] = xnormal3*xarea ;
168 // ===================================
169 // Calculate Normal Vector for Polygon
170 // ===================================
171 inline void calculateNormalForPolyg(const double **coords, int nbOfPtsInPolygs,
174 double coordOfBary[3];
176 calculateBarycenterDyn(coords,nbOfPtsInPolygs,3,coordOfBary);
177 double xnormal1 = (coords[0][1]-coords[1][1]) * (coordOfBary[2]-coords[1][2])
178 - (coords[0][2]-coords[1][2]) * (coordOfBary[1]-coords[1][1]);
180 double xnormal2 = (coords[0][2]-coords[1][2]) * (coordOfBary[0]-coords[1][0])
181 - (coords[0][0]-coords[1][0]) * (coordOfBary[2]-coords[1][2]);
183 double xnormal3 = (coords[0][0]-coords[1][0]) * (coordOfBary[1]-coords[1][1])
184 - (coords[0][1]-coords[1][1]) * (coordOfBary[0]-coords[1][0]);
186 double xarea=sqrt(xnormal1*xnormal1 + xnormal2*xnormal2 + xnormal3*xnormal3);
190 //std::string diagnosis"Can not calculate normal - polygon is singular";
191 throw std::exception();
194 xnormal1 = -xnormal1/xarea;
195 xnormal2 = -xnormal2/xarea;
196 xnormal3 = -xnormal3/xarea;
197 xarea = calculateAreaForPolyg(coords,nbOfPtsInPolygs,3);
198 normal[0] = xnormal1*xarea ;
199 normal[1] = xnormal2*xarea ;
200 normal[2] = xnormal3*xarea ;
203 // ==========================
204 // Calculate Area for Polygon
205 // ==========================
206 inline double calculateAreaForPolyg(const double **coords, int nbOfPtsInPolygs,
210 double coordOfBary[3];
212 calculateBarycenterDyn(coords,nbOfPtsInPolygs,spaceDim,coordOfBary);
213 for ( int i=0; i<nbOfPtsInPolygs; i++ )
215 double tmp = calculateAreaForTria(coords[i],coords[(i+1)%nbOfPtsInPolygs],
216 coordOfBary,spaceDim);
222 double calculateAreaForQPolyg(const double **coords, int nbOfPtsInPolygs, int spaceDim)
225 if(nbOfPtsInPolygs%2==0)
229 std::vector<Node *> nodes(nbOfPtsInPolygs);
230 for(int i=0;i<nbOfPtsInPolygs;i++)
231 nodes[i]=new Node(coords[i][0],coords[i][1]);
232 QuadraticPolygon *pol=QuadraticPolygon::BuildArcCirclePolygon(nodes);
233 double ret=pol->getArea();
238 return calculateAreaForPolyg(coords,nbOfPtsInPolygs/2,spaceDim);
242 std::ostringstream oss; oss << "INTERP_KERNEL::calculateAreaForQPolyg : nb of points in quadratic polygon is " << nbOfPtsInPolygs << " should be even !";
243 throw INTERP_KERNEL::Exception(oss.str().c_str());
247 // ==========================
248 // Calculate Volume for Tetra
249 // ==========================
250 inline double calculateVolumeForTetra(const double *p1, const double *p2,
251 const double *p3, const double *p4)
253 return ( (p3[0]-p1[0])*( (p2[1]-p1[1])*(p4[2]-p1[2])
254 - (p2[2]-p1[2])*(p4[1]-p1[1]) )
255 - (p2[0]-p1[0])*( (p3[1]-p1[1])*(p4[2]-p1[2])
256 - (p3[2]-p1[2])*(p4[1]-p1[1]) )
257 + (p4[0]-p1[0])*( (p3[1]-p1[1])*(p2[2]-p1[2])
258 - (p3[2]-p1[2])*(p2[1]-p1[1]) )
262 // =========================
263 // Calculate Volume for Pyra
264 // =========================
265 inline double calculateVolumeForPyra(const double *p1, const double *p2,
266 const double *p3, const double *p4,
269 return ( ((p3[0]-p1[0])*( (p2[1]-p1[1])*(p5[2]-p1[2])
270 - (p2[2]-p1[2])*(p5[1]-p1[1]) )
271 -(p2[0]-p1[0])*( (p3[1]-p1[1])*(p5[2]-p1[2])
272 - (p3[2]-p1[2])*(p5[1]-p1[1]) )
273 +(p5[0]-p1[0])*( (p3[1]-p1[1])*(p2[2]-p1[2])
274 - (p3[2]-p1[2])*(p2[1]-p1[1]) ))
276 ((p4[0]-p1[0])*( (p3[1]-p1[1])*(p5[2]-p1[2])
277 - (p3[2]-p1[2])*(p5[1]-p1[1]) )
278 -(p3[0]-p1[0])*( (p4[1]-p1[1])*(p5[2]-p1[2])
279 - (p4[2]-p1[2])*(p5[1]-p1[1]))
280 +(p5[0]-p1[0])*( (p4[1]-p1[1])*(p3[2]-p1[2])
281 - (p4[2]-p1[2])*(p3[1]-p1[1]) ))
285 // ==========================
286 // Calculate Volume for Penta
287 // ==========================
288 inline double calculateVolumeForPenta(const double *p1, const double *p2,
289 const double *p3, const double *p4,
290 const double *p5, const double *p6)
292 double a1 = (p2[0]-p3[0])/2.0, a2 = (p2[1]-p3[1])/2.0, a3 = (p2[2]-p3[2])/2.0;
293 double b1 = (p5[0]-p6[0])/2.0, b2 = (p5[1]-p6[1])/2.0, b3 = (p5[2]-p6[2])/2.0;
294 double c1 = (p4[0]-p1[0])/2.0, c2 = (p4[1]-p1[1])/2.0, c3 = (p4[2]-p1[2])/2.0;
295 double d1 = (p5[0]-p2[0])/2.0, d2 = (p5[1]-p2[1])/2.0, d3 = (p5[2]-p2[2])/2.0;
296 double e1 = (p6[0]-p3[0])/2.0, e2 = (p6[1]-p3[1])/2.0, e3 = (p6[2]-p3[2])/2.0;
297 double f1 = (p1[0]-p3[0])/2.0, f2 = (p1[1]-p3[1])/2.0, f3 = (p1[2]-p3[2])/2.0;
298 double h1 = (p4[0]-p6[0])/2.0, h2 = (p4[1]-p6[1])/2.0, h3 = (p4[2]-p6[2])/2.0;
300 double A = a1*c2*f3 - a1*c3*f2 - a2*c1*f3 + a2*c3*f1 + a3*c1*f2 - a3*c2*f1;
301 double B = b1*c2*h3 - b1*c3*h2 - b2*c1*h3 + b2*c3*h1 + b3*c1*h2 - b3*c2*h1;
302 double C = (a1*c2*h3 + b1*c2*f3) - (a1*c3*h2 + b1*c3*f2)
303 - (a2*c1*h3 + b2*c1*f3) + (a2*c3*h1 + b2*c3*f1)
304 + (a3*c1*h2 + b3*c1*f2) - (a3*c2*h1 + b3*c2*f1);
305 double D = a1*d2*f3 - a1*d3*f2 - a2*d1*f3 + a2*d3*f1 + a3*d1*f2 - a3*d2*f1;
306 double E = b1*d2*h3 - b1*d3*h2 - b2*d1*h3 + b2*d3*h1 + b3*d1*h2 - b3*d2*h1;
307 double F = (a1*d2*h3 + b1*d2*f3) - (a1*d3*h2 + b1*d3*f2)
308 - (a2*d1*h3 + b2*d1*f3) + (a2*d3*h1 + b2*d3*f1)
309 + (a3*d1*h2 + b3*d1*f2) - (a3*d2*h1 + b3*d2*f1);
310 double G = a1*e2*f3 - a1*e3*f2 - a2*e1*f3 + a2*e3*f1 + a3*e1*f2 - a3*e2*f1;
311 double H = b1*e2*h3 - b1*e3*h2 - b2*e1*h3 + b2*e3*h1 + b3*e1*h2 - b3*e2*h1;
312 double P = (a1*e2*h3 + b1*e2*f3) - (a1*e3*h2 + b1*e3*f2)
313 - (a2*e1*h3 + b2*e1*f3) + (a2*e3*h1 + b2*e3*f1)
314 + (a3*e1*h2 + b3*e1*f2) - (a3*e2*h1 + b3*e2*f1);
316 return (-2.0*(2.0*(A + B + D + E + G + H) + C + F + P)/9.0);
319 // =========================
320 // Calculate Volume for Hexa
321 // =========================
322 inline double calculateVolumeForHexa(const double *pt1, const double *pt2,
323 const double *pt3, const double *pt4,
324 const double *pt5, const double *pt6,
325 const double *pt7, const double *pt8)
327 double a1=(pt3[0]-pt4[0])/8.0, a2=(pt3[1]-pt4[1])/8.0, a3=(pt3[2]-pt4[2])/8.0;
328 double b1=(pt2[0]-pt1[0])/8.0, b2=(pt2[1]-pt1[1])/8.0, b3=(pt2[2]-pt1[2])/8.0;
329 double c1=(pt7[0]-pt8[0])/8.0, c2=(pt7[1]-pt8[1])/8.0, c3=(pt7[2]-pt8[2])/8.0;
330 double d1=(pt6[0]-pt5[0])/8.0, d2=(pt6[1]-pt5[1])/8.0, d3=(pt6[2]-pt5[2])/8.0;
331 double e1=(pt3[0]-pt2[0])/8.0, e2=(pt3[1]-pt2[1])/8.0, e3=(pt3[2]-pt2[2])/8.0;
332 double f1=(pt4[0]-pt1[0])/8.0, f2=(pt4[1]-pt1[1])/8.0, f3=(pt4[2]-pt1[2])/8.0;
333 double h1=(pt7[0]-pt6[0])/8.0, h2=(pt7[1]-pt6[1])/8.0, h3=(pt7[2]-pt6[2])/8.0;
334 double p1=(pt8[0]-pt5[0])/8.0, p2=(pt8[1]-pt5[1])/8.0, p3=(pt8[2]-pt5[2])/8.0;
335 double q1=(pt3[0]-pt7[0])/8.0, q2=(pt3[1]-pt7[1])/8.0, q3=(pt3[2]-pt7[2])/8.0;
336 double r1=(pt4[0]-pt8[0])/8.0, r2=(pt4[1]-pt8[1])/8.0, r3=(pt4[2]-pt8[2])/8.0;
337 double s1=(pt2[0]-pt6[0])/8.0, s2=(pt2[1]-pt6[1])/8.0, s3=(pt2[2]-pt6[2])/8.0;
338 double t1=(pt1[0]-pt5[0])/8.0, t2=(pt1[1]-pt5[1])/8.0, t3=(pt1[2]-pt5[2])/8.0;
340 double A = a1*e2*q3 - a1*e3*q2 - a2*e1*q3 + a2*e3*q1 + a3*e1*q2 - a3*e2*q1;
341 double B = c1*h2*q3 - c1*h3*q2 - c2*h1*q3 + c2*h3*q1 + c3*h1*q2 - c3*h2*q1;
342 double C = (a1*h2 + c1*e2)*q3 - (a1*h3 + c1*e3)*q2
343 - (a2*h1 + c2*e1)*q3 + (a2*h3 + c2*e3)*q1
344 + (a3*h1 + c3*e1)*q2 - (a3*h2 + c3*e2)*q1;
345 double D = b1*e2*s3 - b1*e3*s2 - b2*e1*s3 + b2*e3*s1 + b3*e1*s2 - b3*e2*s1;
346 double E = d1*h2*s3 - d1*h3*s2 - d2*h1*s3 + d2*h3*s1 + d3*h1*s2 - d3*h2*s1;
347 double F = (b1*h2 + d1*e2)*s3 - (b1*h3 + d1*e3)*s2
348 - (b2*h1 + d2*e1)*s3 + (b2*h3 + d2*e3)*s1
349 + (b3*h1 + d3*e1)*s2 - (b3*h2 + d3*e2)*s1;
350 double G = (a1*e2*s3 + b1*e2*q3) - (a1*e3*s2 + b1*e3*q2)
351 - (a2*e1*s3 + b2*e1*q3) + (a2*e3*s1 + b2*e3*q1)
352 + (a3*e1*s2 + b3*e1*q2) - (a3*e2*s1 + b3*e2*q1);
353 double H = (c1*h2*s3 + d1*h2*q3) - (c1*h3*s2 + d1*h3*q2)
354 - (c2*h1*s3 + d2*h1*q3) + (c2*h3*s1 + d2*h3*q1)
355 + (c3*h1*s2 + d3*h1*q2) - (c3*h2*s1 + d3*h2*q1);
356 double I = ((a1*h2 + c1*e2)*s3 + (b1*h2 + d1*e2)*q3)
357 - ((a1*h3 + c1*e3)*s2 + (b1*h3 + d1*e3)*q2)
358 - ((a2*h1 + c2*e1)*s3 + (b2*h1 + d2*e1)*q3)
359 + ((a2*h3 + c2*e3)*s1 + (b2*h3 + d2*e3)*q1)
360 + ((a3*h1 + c3*e1)*s2 + (b3*h1 + d3*e1)*q2)
361 - ((a3*h2 + c3*e2)*s1 + (b3*h2 + d3*e2)*q1);
362 double J = a1*f2*r3 - a1*f3*r2 - a2*f1*r3 + a2*f3*r1 + a3*f1*r2 - a3*f2*r1;
363 double K = c1*p2*r3 - c1*p3*r2 - c2*p1*r3 + c2*p3*r1 + c3*p1*r2 - c3*p2*r1;
364 double L = (a1*p2 + c1*f2)*r3 - (a1*p3 + c1*f3)*r2
365 - (a2*p1 + c2*f1)*r3 + (a2*p3 + c2*f3)*r1
366 + (a3*p1 + c3*f1)*r2 - (a3*p2 + c3*f2)*r1;
367 double M = b1*f2*t3 - b1*f3*t2 - b2*f1*t3 + b2*f3*t1 + b3*f1*t2 - b3*f2*t1;
368 double N = d1*p2*t3 - d1*p3*t2 - d2*p1*t3 + d2*p3*t1 + d3*p1*t2 - d3*p2*t1;
369 double O = (b1*p2 + d1*f2)*t3 - (b1*p3 + d1*f3)*t2
370 - (b2*p1 + d2*f1)*t3 + (b2*p3 + d2*f3)*t1
371 + (b3*p1 + d3*f1)*t2 - (b3*p2 + d3*f2)*t1;
372 double P = (a1*f2*t3 + b1*f2*r3) - (a1*f3*t2 + b1*f3*r2)
373 - (a2*f1*t3 + b2*f1*r3) + (a2*f3*t1 + b2*f3*r1)
374 + (a3*f1*t2 + b3*f1*r2) - (a3*f2*t1 + b3*f2*r1);
375 double Q = (c1*p2*t3 + d1*p2*r3) - (c1*p3*t2 + d1*p3*r2)
376 - (c2*p1*t3 + d2*p1*r3) + (c2*p3*t1 + d2*p3*r1)
377 + (c3*p1*t2 + d3*p1*r2) - (c3*p2*t1 + d3*p2*r1);
378 double R = ((a1*p2 + c1*f2)*t3 + (b1*p2 + d1*f2)*r3)
379 - ((a1*p3 + c1*f3)*t2 + (b1*p3 + d1*f3)*r2)
380 - ((a2*p1 + c2*f1)*t3 + (b2*p1 + d2*f1)*r3)
381 + ((a2*p3 + c2*f3)*t1 + (b2*p3 + d2*f3)*r1)
382 + ((a3*p1 + c3*f1)*t2 + (b3*p1 + d3*f1)*r2)
383 - ((a3*p2 + c3*f2)*t1 + (b3*p2 + d3*f2)*r1);
384 double S = (a1*e2*r3 + a1*f2*q3) - (a1*e3*r2 + a1*f3*q2)
385 - (a2*e1*r3 + a2*f1*q3) + (a2*e3*r1 + a2*f3*q1)
386 + (a3*e1*r2 + a3*f1*q2) - (a3*e2*r1 + a3*f2*q1);
387 double T = (c1*h2*r3 + c1*p2*q3) - (c1*h3*r2 + c1*p3*q2)
388 - (c2*h1*r3 + c2*p1*q3) + (c2*h3*r1 + c2*p3*q1)
389 + (c3*h1*r2 + c3*p1*q2) - (c3*h2*r1 + c3*p2*q1);
390 double U = ((a1*h2 + c1*e2)*r3 + (a1*p2 + c1*f2)*q3)
391 - ((a1*h3 + c1*e3)*r2 + (a1*p3 + c1*f3)*q2)
392 - ((a2*h1 + c2*e1)*r3 + (a2*p1 + c2*f1)*q3)
393 + ((a2*h3 + c2*e3)*r1 + (a2*p3 + c2*f3)*q1)
394 + ((a3*h1 + c3*e1)*r2 + (a3*p1 + c3*f1)*q2)
395 - ((a3*h2 + c3*e2)*r1 + (a3*p2 + c3*f2)*q1);
396 double V = (b1*e2*t3 + b1*f2*s3) - (b1*e3*t2 + b1*f3*s2)
397 - (b2*e1*t3 + b2*f1*s3) + (b2*e3*t1 + b2*f3*s1)
398 + (b3*e1*t2 + b3*f1*s2) - (b3*e2*t1 + b3*f2*s1);
399 double W = (d1*h2*t3 + d1*p2*s3) - (d1*h3*t2 + d1*p3*s2)
400 - (d2*h1*t3 + d2*p1*s3) + (d2*h3*t1 + d2*p3*s1)
401 + (d3*h1*t2 + d3*p1*s2) - (d3*h2*t1 + d3*p2*s1);
402 double X = ((b1*h2 + d1*e2)*t3 + (b1*p2 + d1*f2)*s3)
403 - ((b1*h3 + d1*e3)*t2 + (b1*p3 + d1*f3)*s2)
404 - ((b2*h1 + d2*e1)*t3 + (b2*p1 + d2*f1)*s3)
405 + ((b2*h3 + d2*e3)*t1 + (b2*p3 + d2*f3)*s1)
406 + ((b3*h1 + d3*e1)*t2 + (b3*p1 + d3*f1)*s2)
407 - ((b3*h2 + d3*e2)*t1 + (b3*p2 + d3*f2)*s1);
408 double Y = (a1*e2*t3 + a1*f2*s3 + b1*e2*r3 + b1*f2*q3)
409 - (a1*e3*t2 + a1*f3*s2 + b1*e3*r2 + b1*f3*q2)
410 - (a2*e1*t3 + a2*f1*s3 + b2*e1*r3 + b2*f1*q3)
411 + (a2*e3*t1 + a2*f3*s1 + b2*e3*r1 + b2*f3*q1)
412 + (a3*e1*t2 + a3*f1*s2 + b3*e1*r2 + b3*f1*q2)
413 - (a3*e2*t1 + a3*f2*s1 + b3*e2*r1 + b3*f2*q1);
414 double Z = (c1*h2*t3 + c1*p2*s3 + d1*h2*r3 + d1*p2*q3)
415 - (c1*h3*t2 + c1*p3*s2 + d1*h3*r2 + d1*p3*q2)
416 - (c2*h1*t3 + c2*p1*s3 + d2*h1*r3 + d2*p1*q3)
417 + (c2*h3*t1 + c2*p3*s1 + d2*h3*r1 + d2*p3*q1)
418 + (c3*h1*t2 + c3*p1*s2 + d3*h1*r2 + d3*p1*q2)
419 - (c3*h2*t1 + c3*p2*s1 + d3*h2*r1 + d3*p2*q1);
420 double AA = ((a1*h2 + c1*e2)*t3 + (a1*p2 + c1*f2)*s3
421 +(b1*h2 + d1*e2)*r3 + (b1*p2 + d1*f2)*q3)
422 - ((a1*h3 + c1*e3)*t2 + (a1*p3 + c1*f3)*s2
423 +(b1*h3 + d1*e3)*r2 + (b1*p3 + d1*f3)*q2)
424 - ((a2*h1 + c2*e1)*t3 + (a2*p1 + c2*f1)*s3
425 +(b2*h1 + d2*e1)*r3 + (b2*p1 + d2*f1)*q3)
426 + ((a2*h3 + c2*e3)*t1 + (a2*p3 + c2*f3)*s1
427 +(b2*h3 + d2*e3)*r1 + (b2*p3 + d2*f3)*q1)
428 + ((a3*h1 + c3*e1)*t2 + (a3*p1 + c3*f1)*s2
429 +(b3*h1 + d3*e1)*r2 + (b3*p1 + d3*f1)*q2)
430 - ((a3*h2 + c3*e2)*t1 + (a3*p2 + c3*f2)*s1
431 + (b3*h2 + d3*e2)*r1 + (b3*p2 + d3*f2)*q1);
433 return 64.0*( 8.0*(A + B + D + E + J + K + M + N)
434 + 4.0*(C + F + G + H + L + O + P + Q + S + T + V + W)
435 + 2.0*(I + R + U + X + Y + Z) + AA ) / 27.0 ;
438 // =========================================================================================================================
439 // Calculate Volume for Generic Polyedron, even not convex one, WARNING !!! The polyedron's faces must be correctly ordered
440 // =========================================================================================================================
441 inline double calculateVolumeForPolyh(const double ***pts,
442 const int *nbOfNodesPerFaces,
448 for ( int i=0; i<nbOfFaces; i++ )
453 calculateNormalForPolyg(pts[i],nbOfNodesPerFaces[i],normal);
454 vecForAlt[0]=bary[0]-pts[i][0][0];
455 vecForAlt[1]=bary[1]-pts[i][0][1];
456 vecForAlt[2]=bary[2]-pts[i][0][2];
457 volume+=vecForAlt[0]*normal[0]+vecForAlt[1]*normal[1]+vecForAlt[2]*normal[2];
463 * Calculate Volume for Generic Polyedron, even not convex one, WARNING !!! The polyedron's faces must be correctly ordered.
464 * 2nd API avoiding to create double** arrays. The returned value could be negative if polyhedrons faces are not oriented with normal outside of the
467 template<class ConnType, NumberingPolicy numPol>
468 inline double calculateVolumeForPolyh2(const ConnType *connec, int lgth, const double *coords)
470 std::size_t nbOfFaces=std::count(connec,connec+lgth,-1)+1;
472 const int *work=connec;
473 for(std::size_t iFace=0;iFace<nbOfFaces;iFace++)
475 const int *work2=std::find(work+1,connec+lgth,-1);
476 std::size_t nbOfNodesOfCurFace=std::distance(work,work2);
477 double areaVector[3]={0.,0.,0.};
478 for(std::size_t ptId=0;ptId<nbOfNodesOfCurFace;ptId++)
480 const double *pti=coords+3*OTT<ConnType,numPol>::coo2C(work[ptId]);
481 const double *pti1=coords+3*OTT<ConnType,numPol>::coo2C(work[(ptId+1)%nbOfNodesOfCurFace]);
482 areaVector[0]+=pti[1]*pti1[2]-pti[2]*pti1[1];
483 areaVector[1]+=pti[2]*pti1[0]-pti[0]*pti1[2];
484 areaVector[2]+=pti[0]*pti1[1]-pti[1]*pti1[0];
486 const double *pt=coords+3*work[0];
487 volume+=pt[0]*areaVector[0]+pt[1]*areaVector[1]+pt[2]*areaVector[2];
494 * This method returns the area oriented vector of a polygon. This method is useful for normal computation without
495 * any troubles if several edges are colinears.
496 * @param res must be of size at least 3 to store the result.
498 template<class ConnType, NumberingPolicy numPol>
499 inline void areaVectorOfPolygon(const ConnType *connec, int lgth, const double *coords, double *res)
501 res[0]=0.; res[1]=0.; res[2]=0.;
502 for(int ptId=0;ptId<lgth;ptId++)
504 const double *pti=coords+3*OTT<ConnType,numPol>::coo2C(connec[ptId]);
505 const double *pti1=coords+3*OTT<ConnType,numPol>::coo2C(connec[(ptId+1)%lgth]);
506 res[0]+=pti[1]*pti1[2]-pti[2]*pti1[1];
507 res[1]+=pti[2]*pti1[0]-pti[0]*pti1[2];
508 res[2]+=pti[0]*pti1[1]-pti[1]*pti1[0];
512 template<class ConnType, NumberingPolicy numPol>
513 inline void computePolygonBarycenter3D(const ConnType *connec, int lgth, const double *coords, double *res)
516 areaVectorOfPolygon<ConnType,numPol>(connec,lgth,coords,area);
517 double norm=sqrt(area[0]*area[0]+area[1]*area[1]+area[2]*area[2]);
518 if(norm>std::numeric_limits<double>::min())
520 area[0]/=norm; area[1]/=norm; area[2]/=norm;
521 res[0]=0.; res[1]=0.; res[2]=0.;
522 for(int i=1;i<lgth-1;i++)
526 v[0]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[0])]+
527 coords[3*OTT<ConnType,numPol>::coo2C(connec[i])]+
528 coords[3*OTT<ConnType,numPol>::coo2C(connec[i+1])])/3.;
529 v[1]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[0])+1]+
530 coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+1]+
531 coords[3*OTT<ConnType,numPol>::coo2C(connec[i+1])+1])/3.;
532 v[2]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[0])+2]+
533 coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+2]+
534 coords[3*OTT<ConnType,numPol>::coo2C(connec[i+1])+2])/3.;
535 ConnType tmpConn[3]={connec[0],connec[i],connec[i+1]};
536 areaVectorOfPolygon<ConnType,numPol>(tmpConn,3,coords,tmpArea);
537 double norm2=sqrt(tmpArea[0]*tmpArea[0]+tmpArea[1]*tmpArea[1]+tmpArea[2]*tmpArea[2]);
540 tmpArea[0]/=norm2; tmpArea[1]/=norm2; tmpArea[2]/=norm2;
541 double signOfArea=area[0]*tmpArea[0]+area[1]*tmpArea[1]+area[2]*tmpArea[2];
542 res[0]+=signOfArea*norm2*v[0]/norm; res[1]+=signOfArea*norm2*v[1]/norm; res[2]+=signOfArea*norm2*v[2]/norm;
548 res[0]=0.; res[1]=0.; res[2]=0.;
550 throw INTERP_KERNEL::Exception("computePolygonBarycenter3D : lgth of polygon is < 1 !");
553 for(int i=0;i<lgth;i++)
555 v[0]=coords[3*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])]-coords[3*OTT<ConnType,numPol>::coo2C(connec[i])];
556 v[1]=coords[3*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])+1]-coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+1];
557 v[2]=coords[3*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])+2]-coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+2];
558 double norm2=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
559 res[0]+=(coords[3*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])]+coords[3*OTT<ConnType,numPol>::coo2C(connec[i])])/2.*norm2;
560 res[1]+=(coords[3*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])+1]+coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+1])/2.*norm2;
561 res[2]+=(coords[3*OTT<ConnType,numPol>::coo2C(connec[(i+1)%lgth])+2]+coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+2])/2.*norm2;
564 if(norm>std::numeric_limits<double>::min())
566 res[0]/=norm; res[1]/=norm; res[2]/=norm;
571 res[0]=0.; res[1]=0.; res[2]=0.;
572 for(int i=0;i<lgth;i++)
574 res[0]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[i])];
575 res[1]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+1];
576 res[2]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[i])+2];
578 res[0]/=lgth; res[1]/=lgth; res[2]/=lgth;
584 inline double integrationOverA3DLine(double u1, double v1, double u2, double v2, double A, double B, double C)
586 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))+
587 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.;
590 template<class ConnType, NumberingPolicy numPol>
591 inline void barycenterOfPolyhedron(const ConnType *connec, int lgth, const double *coords, double *res)
593 std::size_t nbOfFaces=std::count(connec,connec+lgth,-1)+1;
594 res[0]=0.; res[1]=0.; res[2]=0.;
595 const int *work=connec;
596 for(std::size_t i=0;i<nbOfFaces;i++)
598 const int *work2=std::find(work+1,connec+lgth,-1);
599 int nbOfNodesOfCurFace=(int)std::distance(work,work2);
600 // projection to (u,v) of each faces of polyh to compute integral(x^2/2) on each faces.
602 areaVectorOfPolygon<ConnType,numPol>(work,nbOfNodesOfCurFace,coords,normal);
603 double normOfNormal=sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
604 if(normOfNormal<std::numeric_limits<double>::min())
606 normal[0]/=normOfNormal; normal[1]/=normOfNormal; normal[2]/=normOfNormal;
607 double u[2]={normal[1],-normal[0]};
608 double s=sqrt(u[0]*u[0]+u[1]*u[1]);
612 u[0]/=std::abs(s); u[1]/=std::abs(s);
615 { u[0]=1.; u[1]=0.; }
616 //C : high in plane of polyhedron face : always constant
617 double w=normal[0]*coords[3*OTT<ConnType,numPol>::coo2C(work[0])]+
618 normal[1]*coords[3*OTT<ConnType,numPol>::coo2C(work[0])+1]+
619 normal[2]*coords[3*OTT<ConnType,numPol>::coo2C(work[0])+2];
620 // A,B,D,F,G,H,L,M,N coeffs of rotation matrix defined by (u,c,s)
621 double A=u[0]*u[0]*(1-c)+c;
622 double B=u[0]*u[1]*(1-c);
625 double G=u[1]*u[1]*(1-c)+c;
633 for(int j=0;j<nbOfNodesOfCurFace;j++)
635 const double *p1=coords+3*OTT<ConnType,numPol>::coo2C(work[j]);
636 const double *p2=coords+3*OTT<ConnType,numPol>::coo2C(work[(j+1)%nbOfNodesOfCurFace]);
637 double u1=A*p1[0]+B*p1[1]+D*p1[2];
638 double u2=A*p2[0]+B*p2[1]+D*p2[2];
639 double v1=F*p1[0]+G*p1[1]+H*p1[2];
640 double v2=F*p2[0]+G*p2[1]+H*p2[2];
642 double gx=integrationOverA3DLine(u1,v1,u2,v2,A,B,CX);
643 double gy=integrationOverA3DLine(u1,v1,u2,v2,F,G,CY);
644 double gz=integrationOverA3DLine(u1,v1,u2,v2,L,M,CZ);
645 res[0]+=gx*normal[0];
646 res[1]+=gy*normal[1];
647 res[2]+=gz*normal[2];
651 double vol=calculateVolumeForPolyh2<ConnType,numPol>(connec,lgth,coords);
652 if(fabs(vol)>std::numeric_limits<double>::min())
654 res[0]/=vol; res[1]/=vol; res[2]/=vol;
659 res[0]=0.; res[1]=0.; res[2]=0.;
661 for(std::size_t i=0;i<nbOfFaces;i++)
663 const int *work2=std::find(work+1,connec+lgth,-1);
664 int nbOfNodesOfCurFace=(int)std::distance(work,work2);
666 areaVectorOfPolygon<ConnType,numPol>(work,nbOfNodesOfCurFace,coords,normal);
667 double normOfNormal=sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
668 if(normOfNormal<std::numeric_limits<double>::min())
672 computePolygonBarycenter3D<ConnType,numPol>(work,nbOfNodesOfCurFace,coords,tmpBary);
673 res[0]+=normOfNormal*tmpBary[0]; res[1]+=normOfNormal*tmpBary[1]; res[2]+=normOfNormal*tmpBary[2];
676 res[0]/=sum; res[1]/=sum; res[2]/=sum;
680 // ============================================================================================================================================
681 // Calculate Volume for NON Generic Polyedron. Only polydrons with bary included in pts is supported by this method. Result is always positive.
682 // ============================================================================================================================================
683 inline double calculateVolumeForPolyhAbs(const double ***pts,
684 const int *nbOfNodesPerFaces,
690 for ( int i=0; i<nbOfFaces; i++ )
695 calculateNormalForPolyg(pts[i],nbOfNodesPerFaces[i],normal);
696 vecForAlt[0]=bary[0]-pts[i][0][0];
697 vecForAlt[1]=bary[1]-pts[i][0][1];
698 vecForAlt[2]=bary[2]-pts[i][0][2];
699 volume+=fabs(vecForAlt[0]*normal[0]+vecForAlt[1]*normal[1]+vecForAlt[2]*normal[2]);
705 inline double addComponentsOfVec(const double **pts, int rk)
707 return pts[N-1][rk]+addComponentsOfVec<N-1>(pts,rk);
711 inline double addComponentsOfVec<1>(const double **pts, int rk)
716 template<int N, int DIM>
717 inline void calculateBarycenter(const double **pts, double *bary)
719 bary[DIM-1]=addComponentsOfVec<N>(pts,DIM-1)/N;
720 calculateBarycenter<N,DIM-1>(pts,bary);
724 inline void calculateBarycenter<2,0>(const double ** /*pts*/, double * /*bary*/)
729 inline void calculateBarycenter<3,0>(const double ** /*pts*/, double * /*bary*/)
734 inline void calculateBarycenter<4,0>(const double ** /*pts*/, double * /*bary*/)
739 inline void calculateBarycenter<5,0>(const double ** /*pts*/, double * /*bary*/)
744 inline void calculateBarycenter<6,0>(const double ** /*pts*/, double * /*bary*/)
749 inline void calculateBarycenter<7,0>(const double ** /*pts*/, double * /*bary*/)
754 inline void calculateBarycenter<8,0>(const double ** /*pts*/, double * /*bary*/)
758 inline void calculateBarycenterDyn(const double **pts, int nbPts,
759 int dim, double *bary)
761 for(int i=0;i<dim;i++)
764 for(int j=0;j<nbPts;j++)
772 template<int SPACEDIM>
773 inline void calculateBarycenterDyn2(const double *pts, int nbPts, double *bary)
775 for(int i=0;i<SPACEDIM;i++)
778 for(int j=0;j<nbPts;j++)
780 temp+=pts[j*SPACEDIM+i];
786 inline void computePolygonBarycenter2DEngine(double **coords, int lgth, double *res)
789 res[0]=0.; res[1]=0.;
790 for(int i=0;i<lgth;i++)
792 double cp=coords[i][0]*coords[(i+1)%lgth][1]-coords[i][1]*coords[(i+1)%lgth][0];
794 res[0]+=cp*(coords[i][0]+coords[(i+1)%lgth][0]);
795 res[1]+=cp*(coords[i][1]+coords[(i+1)%lgth][1]);
801 template<class ConnType, NumberingPolicy numPol>
802 inline void computePolygonBarycenter2D(const ConnType *connec, int lgth, const double *coords, double *res)
804 double **coords2=new double *[lgth];
805 for(int i=0;i<lgth;i++)
806 coords2[i]=const_cast<double *>(coords+2*OTT<ConnType,numPol>::coo2C(connec[i]));
807 computePolygonBarycenter2DEngine(coords2,lgth,res);
811 inline void computeQPolygonBarycenter2D(double **coords, int nbOfPtsInPolygs, int spaceDim, double *res)
813 if(nbOfPtsInPolygs%2==0)
817 std::vector<Node *> nodes(nbOfPtsInPolygs);
818 for(int i=0;i<nbOfPtsInPolygs;i++)
819 nodes[i]=new Node(coords[i][0],coords[i][1]);
820 QuadraticPolygon *pol=QuadraticPolygon::BuildArcCirclePolygon(nodes);
821 pol->getBarycenter(res);
825 return computePolygonBarycenter2DEngine(coords,nbOfPtsInPolygs/2,res);
829 std::ostringstream oss; oss << "INTERP_KERNEL::computeQPolygonBarycenter2D : nb of points in quadratic polygon is " << nbOfPtsInPolygs << " should be even !";
830 throw INTERP_KERNEL::Exception(oss.str().c_str());