1 #ifndef MEDMEM_FORMULAE
2 #define MEDMEM_FORMULAE
6 inline void CalculateBarycenterDyn(const double **pts, int nbPts, int dim, double *bary);
8 inline double CalculateAreaForPolyg(const double **coords, int nbOfPtsInPolygs, int spaceDim);
10 inline double CalculateAreaForTria(const double *p1, const double *p2, const double *p3,int spaceDim)
13 return (-((p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1]))/2.0);
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);
23 inline double CalculateAreaForQuad(const double *p1, const double *p2, const double *p3, const double *p4, int spaceDim)
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;
32 return (- 4.0*(b1*c2 - c1*b2 + a1*c2 - c1*a2 + b1*d2 -
33 d1*b2 + a1*d2 - d1*a2));
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);
52 inline void CalculateNormalForTria(const double *p1, const double *p2, const double *p3, double *normal)
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;
59 inline void CalculateNormalForQuad(const double *p1, const double *p2, const double *p3, const double *p4, double *normal)
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 ;
74 inline void CalculateNormalForPolyg(const double **coords, int nbOfPtsInPolygs, double *normal)
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 ;
91 inline double CalculateAreaForPolyg(const double **coords, int nbOfPtsInPolygs, int spaceDim)
93 double coordOfBary[3];
94 CalculateBarycenterDyn(coords,nbOfPtsInPolygs,spaceDim,coordOfBary);
96 for(int i=0;i<nbOfPtsInPolygs;i++)
98 double tmp=CalculateAreaForTria(coords[i],coords[(i+1)%nbOfPtsInPolygs],coordOfBary,spaceDim);
104 inline double CalculateVolumeForTetra(const double *p1, const double *p2, const double *p3, const double *p4)
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;
111 inline double CalculateVolumeForPyra(const double *p1, const double *p2, const double *p3, const double *p4, const double *p5)
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])))
122 inline double CalculateVolumeForPenta(const double *p1, const double *p2, const double *p3, const double *p4, const double *p5, const double *p6)
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;
132 double A = a1*c2*f3 - a1*c3*f2 - a2*c1*f3 + a2*c3*f1 +
134 double B = b1*c2*h3 - b1*c3*h2 - b2*c1*h3 + b2*c3*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 +
141 double E = b1*d2*h3 - b1*d3*h2 - b2*d1*h3 + b2*d3*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 +
148 double H = b1*e2*h3 - b1*e3*h2 - b2*e1*h3 + b2*e3*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);
154 return (-2.0*(2.0*(A + B + D + E + G + H) + C + F + P)/9.0);
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)
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;
172 double A = a1*e2*q3 - a1*e3*q2 - a2*e1*q3 + a2*e3*q1 +
174 double B = c1*h2*q3 - c1*h3*q2 - c2*h1*q3 + c2*h3*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 +
181 double E = d1*h2*s3 - d1*h3*s2 - d2*h1*s3 + d2*h3*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 +
200 double K = c1*p2*r3 - c1*p3*r2 - c2*p1*r3 + c2*p3*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 +
207 double N = d1*p2*t3 - d1*p3*t2 - d2*p1*t3 + d2*p3*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);
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) +
279 inline double CalculateVolumeForPolyh(const double ***pts, const int *nbOfNodesPerFaces, int nbOfFaces, const double *bary)
282 for(int i=0;i<nbOfFaces;i++)
285 CalculateNormalForPolyg(pts[i],nbOfNodesPerFaces[i],normal);
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]);
296 inline double addComponentsOfVec(const double **pts, int rk)
298 return pts[N-1][rk]+addComponentsOfVec<N-1>(pts,rk);
302 inline double addComponentsOfVec<1>(const double **pts, int rk)
307 template<int N, int DIM>
308 inline void CalculateBarycenter(const double **pts, double *bary)
310 bary[DIM-1]=addComponentsOfVec<N>(pts,DIM-1)/N;
311 CalculateBarycenter<N,DIM-1>(pts,bary);
315 inline void CalculateBarycenter<2,0>(const double **pts, double *bary) { }
318 inline void CalculateBarycenter<3,0>(const double **pts, double *bary) { }
321 inline void CalculateBarycenter<4,0>(const double **pts, double *bary) { }
324 inline void CalculateBarycenter<5,0>(const double **pts, double *bary) { }
327 inline void CalculateBarycenter<6,0>(const double **pts, double *bary) { }
330 inline void CalculateBarycenter<7,0>(const double **pts, double *bary) { }
333 inline void CalculateBarycenter<8,0>(const double **pts, double *bary) { }
335 inline void CalculateBarycenterDyn(const double **pts, int nbPts, int dim, double *bary)
337 for(int i=0;i<dim;i++)
340 for(int j=0;j<nbPts;j++)