Branch data Line data Source code
1 : : /*=========================================================================
2 : :
3 : : Module: $RCSfile: V_TetMetric.cpp,v $
4 : :
5 : : Copyright (c) 2006 Sandia Corporation.
6 : : All rights reserved.
7 : : See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
8 : :
9 : : This software is distributed WITHOUT ANY WARRANTY; without even
10 : : the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 : : PURPOSE. See the above copyright notice for more information.
12 : :
13 : : =========================================================================*/
14 : :
15 : : /*
16 : : *
17 : : * TetMetric.cpp contains quality calculations for Tets
18 : : *
19 : : * This file is part of VERDICT
20 : : *
21 : : */
22 : :
23 : : #define VERDICT_EXPORTS
24 : :
25 : : #include "moab/verdict.h"
26 : : #include "verdict_defines.hpp"
27 : : #include "VerdictVector.hpp"
28 : : #include "V_GaussIntegration.hpp"
29 : : #include <memory.h>
30 : :
31 : : //! the average volume of a tet
32 : : double verdict_tet_size = 0;
33 : :
34 : : /*!
35 : : set the average volume of a tet
36 : : */
37 : 1 : C_FUNC_DEF void v_set_tet_size( double size )
38 : : {
39 : 1 : verdict_tet_size = size;
40 : 1 : }
41 : :
42 : : /*!
43 : : get the weights based on the average size
44 : : of a tet
45 : : */
46 : 46 : int get_weight( VerdictVector& w1, VerdictVector& w2, VerdictVector& w3 )
47 : : {
48 : : static const double rt3 = sqrt( 3.0 );
49 : : static const double root_of_2 = sqrt( 2.0 );
50 : :
51 : 46 : w1.set( 1, 0, 0 );
52 : 46 : w2.set( 0.5, 0.5 * rt3, 0 );
53 : 46 : w3.set( 0.5, rt3 / 6.0, root_of_2 / rt3 );
54 : :
55 [ + - ][ + - ]: 46 : double scale = pow( 6. * verdict_tet_size / determinant( w1, w2, w3 ), 0.3333333333333 );
[ + - ]
56 : :
57 : 46 : w1 *= scale;
58 : 46 : w2 *= scale;
59 : 46 : w3 *= scale;
60 : :
61 : 46 : return 1;
62 : : }
63 : :
64 : : /*!
65 : : the edge ratio of a tet
66 : :
67 : : NB (P. Pebay 01/22/07):
68 : : Hmax / Hmin where Hmax and Hmin are respectively the maximum and the
69 : : minimum edge lengths
70 : : */
71 : 22 : C_FUNC_DEF double v_tet_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
72 : : {
73 [ + - ][ + - ]: 22 : VerdictVector a, b, c, d, e, f;
[ + - ][ + - ]
[ + - ][ + - ]
74 : :
75 : 44 : a.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
76 [ + - ]: 22 : coordinates[1][2] - coordinates[0][2] );
77 : :
78 : 44 : b.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
79 [ + - ]: 22 : coordinates[2][2] - coordinates[1][2] );
80 : :
81 : 44 : c.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
82 [ + - ]: 22 : coordinates[0][2] - coordinates[2][2] );
83 : :
84 : 44 : d.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
85 [ + - ]: 22 : coordinates[3][2] - coordinates[0][2] );
86 : :
87 : 44 : e.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
88 [ + - ]: 22 : coordinates[3][2] - coordinates[1][2] );
89 : :
90 : 44 : f.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
91 [ + - ]: 22 : coordinates[3][2] - coordinates[2][2] );
92 : :
93 [ + - ]: 22 : double a2 = a.length_squared();
94 [ + - ]: 22 : double b2 = b.length_squared();
95 [ + - ]: 22 : double c2 = c.length_squared();
96 [ + - ]: 22 : double d2 = d.length_squared();
97 [ + - ]: 22 : double e2 = e.length_squared();
98 [ + - ]: 22 : double f2 = f.length_squared();
99 : :
100 : : double m2, M2, mab, mcd, mef, Mab, Mcd, Mef;
101 : :
102 [ + + ]: 22 : if( a2 < b2 )
103 : : {
104 : 6 : mab = a2;
105 : 6 : Mab = b2;
106 : : }
107 : : else // b2 <= a2
108 : : {
109 : 16 : mab = b2;
110 : 16 : Mab = a2;
111 : : }
112 [ + + ]: 22 : if( c2 < d2 )
113 : : {
114 : 6 : mcd = c2;
115 : 6 : Mcd = d2;
116 : : }
117 : : else // d2 <= c2
118 : : {
119 : 16 : mcd = d2;
120 : 16 : Mcd = c2;
121 : : }
122 [ + + ]: 22 : if( e2 < f2 )
123 : : {
124 : 8 : mef = e2;
125 : 8 : Mef = f2;
126 : : }
127 : : else // f2 <= e2
128 : : {
129 : 14 : mef = f2;
130 : 14 : Mef = e2;
131 : : }
132 : :
133 [ + + ]: 22 : m2 = mab < mcd ? mab : mcd;
134 [ + + ]: 22 : m2 = m2 < mef ? m2 : mef;
135 : :
136 [ - + ]: 22 : if( m2 < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
137 : :
138 [ + + ]: 22 : M2 = Mab > Mcd ? Mab : Mcd;
139 [ + + ]: 22 : M2 = M2 > Mef ? M2 : Mef;
140 : :
141 : 22 : double edge_ratio = sqrt( M2 / m2 );
142 : :
143 [ + - ][ + - ]: 22 : if( edge_ratio > 0 ) return (double)VERDICT_MIN( edge_ratio, VERDICT_DBL_MAX );
144 [ # # ]: 22 : return (double)VERDICT_MAX( edge_ratio, -VERDICT_DBL_MAX );
145 : : }
146 : :
147 : : /*!
148 : : the scaled jacobian of a tet
149 : :
150 : : minimum of the jacobian divided by the lengths of 3 edge vectors
151 : :
152 : : */
153 : 22 : C_FUNC_DEF double v_tet_scaled_jacobian( int /*num_nodes*/, double coordinates[][3] )
154 : : {
155 : :
156 [ + - ][ + - ]: 22 : VerdictVector side0, side1, side2, side3, side4, side5;
[ + - ][ + - ]
[ + - ][ + - ]
157 : :
158 : 44 : side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
159 [ + - ]: 22 : coordinates[1][2] - coordinates[0][2] );
160 : :
161 : 44 : side1.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
162 [ + - ]: 22 : coordinates[2][2] - coordinates[1][2] );
163 : :
164 : 44 : side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
165 [ + - ]: 22 : coordinates[0][2] - coordinates[2][2] );
166 : :
167 : 44 : side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
168 [ + - ]: 22 : coordinates[3][2] - coordinates[0][2] );
169 : :
170 : 44 : side4.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
171 [ + - ]: 22 : coordinates[3][2] - coordinates[1][2] );
172 : :
173 : 44 : side5.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
174 [ + - ]: 22 : coordinates[3][2] - coordinates[2][2] );
175 : :
176 : : double jacobi;
177 : :
178 [ + - ][ + - ]: 22 : jacobi = side3 % ( side2 * side0 );
179 : :
180 : : // products of lengths squared of each edge attached to a node.
181 [ + - ][ + - ]: 22 : double length_squared[4] = { side0.length_squared() * side2.length_squared() * side3.length_squared(),
[ + - ]
182 [ + - ][ + - ]: 22 : side0.length_squared() * side1.length_squared() * side4.length_squared(),
[ + - ]
183 [ + - ][ + - ]: 22 : side1.length_squared() * side2.length_squared() * side5.length_squared(),
[ + - ]
184 [ + - ][ + - ]: 22 : side3.length_squared() * side4.length_squared() * side5.length_squared() };
[ + - ]
185 : 22 : int which_node = 0;
186 [ + + ]: 22 : if( length_squared[1] > length_squared[which_node] ) which_node = 1;
187 [ + + ]: 22 : if( length_squared[2] > length_squared[which_node] ) which_node = 2;
188 [ + + ]: 22 : if( length_squared[3] > length_squared[which_node] ) which_node = 3;
189 : :
190 : 22 : double length_product = sqrt( length_squared[which_node] );
191 [ - + ]: 22 : if( length_product < fabs( jacobi ) ) length_product = fabs( jacobi );
192 : :
193 [ - + ]: 22 : if( length_product < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
194 : :
195 : : static const double root_of_2 = sqrt( 2.0 );
196 : :
197 : 22 : return (double)( root_of_2 * jacobi / length_product );
198 : : }
199 : :
200 : : /*!
201 : : The radius ratio of a tet
202 : :
203 : : NB (P. Pebay 04/16/07):
204 : : CR / (3.0 * IR) where CR is the circumsphere radius and IR is the inscribed
205 : : sphere radius.
206 : : Note that this metric is similar to the aspect beta of a tet, except that
207 : : it does not return VERDICT_DBL_MAX if the element has negative orientation.
208 : : */
209 : 44 : C_FUNC_DEF double v_tet_radius_ratio( int /*num_nodes*/, double coordinates[][3] )
210 : : {
211 : :
212 : : // Determine side vectors
213 [ + - ][ + + ]: 308 : VerdictVector side[6];
214 : :
215 : 88 : side[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
216 [ + - ]: 44 : coordinates[1][2] - coordinates[0][2] );
217 : :
218 : 88 : side[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
219 [ + - ]: 44 : coordinates[2][2] - coordinates[1][2] );
220 : :
221 : 88 : side[2].set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
222 [ + - ]: 44 : coordinates[0][2] - coordinates[2][2] );
223 : :
224 : 88 : side[3].set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
225 [ + - ]: 44 : coordinates[3][2] - coordinates[0][2] );
226 : :
227 : 88 : side[4].set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
228 [ + - ]: 44 : coordinates[3][2] - coordinates[1][2] );
229 : :
230 : 88 : side[5].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
231 [ + - ]: 44 : coordinates[3][2] - coordinates[2][2] );
232 : :
233 [ + - ][ + - ]: 44 : VerdictVector numerator = side[3].length_squared() * ( side[2] * side[0] ) +
[ + - ][ + - ]
234 [ + - ][ + - ]: 44 : side[2].length_squared() * ( side[3] * side[0] ) +
[ + - ]
235 [ + - ][ + - ]: 88 : side[0].length_squared() * ( side[3] * side[2] );
[ + - ][ + - ]
236 : :
237 : : double area_sum;
238 [ + - ][ + - ]: 44 : area_sum = ( ( side[2] * side[0] ).length() + ( side[3] * side[0] ).length() + ( side[4] * side[1] ).length() +
[ + - ][ + - ]
[ + - ][ + - ]
239 [ + - ][ + - ]: 44 : ( side[3] * side[2] ).length() ) *
240 : 44 : 0.5;
241 : :
242 [ + - ]: 44 : double volume = v_tet_volume( 4, coordinates );
243 : :
244 [ - + ]: 44 : if( fabs( volume ) < VERDICT_DBL_MIN )
245 : 0 : return (double)VERDICT_DBL_MAX;
246 : : else
247 : : {
248 : : double radius_ratio;
249 [ + - ]: 44 : radius_ratio = numerator.length() * area_sum / ( 108 * volume * volume );
250 : :
251 [ + - ]: 44 : return (double)VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
252 : : }
253 : : }
254 : :
255 : : /*!
256 : : The radius ratio of a positively-oriented tet, a.k.a. "aspect beta"
257 : :
258 : : NB (P. Pebay 04/16/07):
259 : : CR / (3.0 * IR) where CR is the circumsphere radius and IR is the inscribed
260 : : sphere radius if the element has positive orientation.
261 : : Note that this metric is similar to the radius ratio of a tet, except that
262 : : it returns VERDICT_DBL_MAX if the element has negative orientation.
263 : :
264 : : */
265 : 22 : C_FUNC_DEF double v_tet_aspect_beta( int /*num_nodes*/, double coordinates[][3] )
266 : : {
267 : :
268 : : // Determine side vectors
269 [ + - ][ + + ]: 154 : VerdictVector side[6];
270 : :
271 : 44 : side[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
272 [ + - ]: 22 : coordinates[1][2] - coordinates[0][2] );
273 : :
274 : 44 : side[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
275 [ + - ]: 22 : coordinates[2][2] - coordinates[1][2] );
276 : :
277 : 44 : side[2].set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
278 [ + - ]: 22 : coordinates[0][2] - coordinates[2][2] );
279 : :
280 : 44 : side[3].set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
281 [ + - ]: 22 : coordinates[3][2] - coordinates[0][2] );
282 : :
283 : 44 : side[4].set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
284 [ + - ]: 22 : coordinates[3][2] - coordinates[1][2] );
285 : :
286 : 44 : side[5].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
287 [ + - ]: 22 : coordinates[3][2] - coordinates[2][2] );
288 : :
289 [ + - ][ + - ]: 22 : VerdictVector numerator = side[3].length_squared() * ( side[2] * side[0] ) +
[ + - ][ + - ]
290 [ + - ][ + - ]: 22 : side[2].length_squared() * ( side[3] * side[0] ) +
[ + - ]
291 [ + - ][ + - ]: 44 : side[0].length_squared() * ( side[3] * side[2] );
[ + - ][ + - ]
292 : :
293 : : double area_sum;
294 [ + - ][ + - ]: 22 : area_sum = ( ( side[2] * side[0] ).length() + ( side[3] * side[0] ).length() + ( side[4] * side[1] ).length() +
[ + - ][ + - ]
[ + - ][ + - ]
295 [ + - ][ + - ]: 22 : ( side[3] * side[2] ).length() ) *
296 : 22 : 0.5;
297 : :
298 [ + - ]: 22 : double volume = v_tet_volume( 4, coordinates );
299 : :
300 [ - + ]: 22 : if( volume < VERDICT_DBL_MIN )
301 : 0 : return (double)VERDICT_DBL_MAX;
302 : : else
303 : : {
304 : : double radius_ratio;
305 [ + - ]: 22 : radius_ratio = numerator.length() * area_sum / ( 108 * volume * volume );
306 : :
307 [ + - ]: 22 : return (double)VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
308 : : }
309 : : }
310 : :
311 : : /*!
312 : : The aspect ratio of a tet
313 : :
314 : : NB (P. Pebay 01/22/07):
315 : : Hmax / (2 sqrt(6) r) where Hmax and r respectively denote the greatest edge
316 : : length and the inradius of the tetrahedron
317 : : */
318 : 44 : C_FUNC_DEF double v_tet_aspect_ratio( int /*num_nodes*/, double coordinates[][3] )
319 : : {
320 : : static const double normal_coeff = sqrt( 6. ) / 12.;
321 : :
322 : : // Determine side vectors
323 [ + - ][ + - ]: 44 : VerdictVector ab, bc, ac, ad, bd, cd;
[ + - ][ + - ]
[ + - ][ + - ]
324 : :
325 : 88 : ab.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
326 [ + - ]: 44 : coordinates[1][2] - coordinates[0][2] );
327 : :
328 : 88 : ac.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
329 [ + - ]: 44 : coordinates[2][2] - coordinates[0][2] );
330 : :
331 : 88 : ad.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
332 [ + - ]: 44 : coordinates[3][2] - coordinates[0][2] );
333 : :
334 [ + - ][ + - ]: 44 : double detTet = ab % ( ac * ad );
335 : :
336 [ - + ]: 44 : if( detTet < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
337 : :
338 : 88 : bc.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
339 [ + - ]: 44 : coordinates[2][2] - coordinates[1][2] );
340 : :
341 : 88 : bd.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
342 [ + - ]: 44 : coordinates[3][2] - coordinates[1][2] );
343 : :
344 : 88 : cd.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
345 [ + - ]: 44 : coordinates[3][2] - coordinates[2][2] );
346 : :
347 [ + - ]: 44 : double ab2 = ab.length_squared();
348 [ + - ]: 44 : double bc2 = bc.length_squared();
349 [ + - ]: 44 : double ac2 = ac.length_squared();
350 [ + - ]: 44 : double ad2 = ad.length_squared();
351 [ + - ]: 44 : double bd2 = bd.length_squared();
352 [ + - ]: 44 : double cd2 = cd.length_squared();
353 : :
354 [ + + ]: 44 : double A = ab2 > bc2 ? ab2 : bc2;
355 [ + + ]: 44 : double B = ac2 > ad2 ? ac2 : ad2;
356 [ + + ]: 44 : double C = bd2 > cd2 ? bd2 : cd2;
357 [ + + ]: 44 : double D = A > B ? A : B;
358 [ + + ]: 44 : double hm = D > C ? sqrt( D ) : sqrt( C );
359 : :
360 [ + - ][ + - ]: 44 : bd = ab * bc;
361 [ + - ]: 44 : A = bd.length();
362 [ + - ][ + - ]: 44 : bd = ab * ad;
363 [ + - ]: 44 : B = bd.length();
364 [ + - ][ + - ]: 44 : bd = ac * ad;
365 [ + - ]: 44 : C = bd.length();
366 [ + - ][ + - ]: 44 : bd = bc * cd;
367 [ + - ]: 44 : D = bd.length();
368 : :
369 : : double aspect_ratio;
370 : 44 : aspect_ratio = normal_coeff * hm * ( A + B + C + D ) / fabs( detTet );
371 : :
372 [ + - ][ + - ]: 44 : if( aspect_ratio > 0 ) return (double)VERDICT_MIN( aspect_ratio, VERDICT_DBL_MAX );
373 [ # # ]: 44 : return (double)VERDICT_MAX( aspect_ratio, -VERDICT_DBL_MAX );
374 : : }
375 : :
376 : : /*!
377 : : the aspect gamma of a tet
378 : :
379 : : srms^3 / (8.48528137423857*V) where srms = sqrt(sum(Si^2)/6), where Si is the edge length
380 : : */
381 : 22 : C_FUNC_DEF double v_tet_aspect_gamma( int /*num_nodes*/, double coordinates[][3] )
382 : : {
383 : :
384 : : // Determine side vectors
385 [ + - ][ + - ]: 22 : VerdictVector side0, side1, side2, side3, side4, side5;
[ + - ][ + - ]
[ + - ][ + - ]
386 : :
387 : 44 : side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
388 [ + - ]: 22 : coordinates[1][2] - coordinates[0][2] );
389 : :
390 : 44 : side1.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
391 [ + - ]: 22 : coordinates[2][2] - coordinates[1][2] );
392 : :
393 : 44 : side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
394 [ + - ]: 22 : coordinates[0][2] - coordinates[2][2] );
395 : :
396 : 44 : side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
397 [ + - ]: 22 : coordinates[3][2] - coordinates[0][2] );
398 : :
399 : 44 : side4.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
400 [ + - ]: 22 : coordinates[3][2] - coordinates[1][2] );
401 : :
402 : 44 : side5.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
403 [ + - ]: 22 : coordinates[3][2] - coordinates[2][2] );
404 : :
405 [ + - ]: 22 : double volume = fabs( v_tet_volume( 4, coordinates ) );
406 : :
407 [ - + ]: 22 : if( volume < VERDICT_DBL_MIN )
408 : 0 : return (double)VERDICT_DBL_MAX;
409 : : else
410 : : {
411 [ + - ][ + - ]: 22 : double srms = sqrt( ( side0.length_squared() + side1.length_squared() + side2.length_squared() +
[ + - ]
412 [ + - ][ + - ]: 22 : side3.length_squared() + side4.length_squared() + side5.length_squared() ) /
[ + - ]
413 : 22 : 6.0 );
414 : :
415 : 22 : double aspect_ratio_gamma = pow( srms, 3 ) / ( 8.48528137423857 * volume );
416 : 22 : return (double)aspect_ratio_gamma;
417 : : }
418 : : }
419 : :
420 : : /*!
421 : : The aspect frobenius of a tet
422 : :
423 : : NB (P. Pebay 01/22/07):
424 : : Frobenius condition number when the reference element is regular
425 : : */
426 : 44 : C_FUNC_DEF double v_tet_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
427 : : {
428 : : static const double normal_exp = 1. / 3.;
429 : :
430 [ + - ][ + - ]: 44 : VerdictVector ab, ac, ad;
[ + - ]
431 : :
432 : 88 : ab.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
433 [ + - ]: 44 : coordinates[1][2] - coordinates[0][2] );
434 : :
435 : 88 : ac.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
436 [ + - ]: 44 : coordinates[2][2] - coordinates[0][2] );
437 : :
438 : 88 : ad.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
439 [ + - ]: 44 : coordinates[3][2] - coordinates[0][2] );
440 : :
441 [ + - ][ + - ]: 44 : double denominator = ab % ( ac * ad );
442 : 44 : denominator *= denominator;
443 : 44 : denominator *= 2.;
444 : 44 : denominator = 3. * pow( denominator, normal_exp );
445 : :
446 [ - + ]: 44 : if( denominator < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
447 : :
448 : : double u[3];
449 [ + - ]: 44 : ab.get_xyz( u );
450 : : double v[3];
451 [ + - ]: 44 : ac.get_xyz( v );
452 : : double w[3];
453 [ + - ]: 44 : ad.get_xyz( w );
454 : :
455 : 44 : double numerator = u[0] * u[0] + u[1] * u[1] + u[2] * u[2];
456 : 44 : numerator += v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
457 : 44 : numerator += w[0] * w[0] + w[1] * w[1] + w[2] * w[2];
458 : 44 : numerator *= 1.5;
459 : 44 : numerator -= v[0] * u[0] + v[1] * u[1] + v[2] * u[2];
460 : 44 : numerator -= w[0] * u[0] + w[1] * u[1] + w[2] * u[2];
461 : 44 : numerator -= w[0] * v[0] + w[1] * v[1] + w[2] * v[2];
462 : :
463 : 44 : double aspect_frobenius = numerator / denominator;
464 : :
465 [ + - ][ + - ]: 44 : if( aspect_frobenius > 0 ) return (double)VERDICT_MIN( aspect_frobenius, VERDICT_DBL_MAX );
466 [ # # ]: 44 : return (double)VERDICT_MAX( aspect_frobenius, -VERDICT_DBL_MAX );
467 : : }
468 : :
469 : : /*!
470 : : The minimum angle of a tet
471 : :
472 : : NB (P. Pebay 01/22/07):
473 : : minimum nonoriented dihedral angle
474 : : */
475 : 44 : C_FUNC_DEF double v_tet_minimum_angle( int /*num_nodes*/, double coordinates[][3] )
476 : : {
477 : : static const double normal_coeff = 180. * .3183098861837906715377675267450287;
478 : :
479 : : // Determine side vectors
480 [ + - ][ + - ]: 44 : VerdictVector ab, bc, ad, cd;
[ + - ][ + - ]
481 : :
482 : 88 : ab.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
483 [ + - ]: 44 : coordinates[1][2] - coordinates[0][2] );
484 : :
485 : 88 : ad.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
486 [ + - ]: 44 : coordinates[3][2] - coordinates[0][2] );
487 : :
488 : 88 : bc.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
489 [ + - ]: 44 : coordinates[2][2] - coordinates[1][2] );
490 : :
491 : 88 : cd.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
492 [ + - ]: 44 : coordinates[3][2] - coordinates[2][2] );
493 : :
494 [ + - ]: 44 : VerdictVector abc = ab * bc;
495 [ + - ]: 44 : double nabc = abc.length();
496 [ + - ]: 44 : VerdictVector abd = ab * ad;
497 [ + - ]: 44 : double nabd = abd.length();
498 [ + - ]: 44 : VerdictVector acd = ad * cd;
499 [ + - ]: 44 : double nacd = acd.length();
500 [ + - ]: 44 : VerdictVector bcd = bc * cd;
501 [ + - ]: 44 : double nbcd = bcd.length();
502 : :
503 [ + - ]: 44 : double alpha = acos( ( abc % abd ) / ( nabc * nabd ) );
504 [ + - ]: 44 : double beta = acos( ( abc % acd ) / ( nabc * nacd ) );
505 [ + - ]: 44 : double gamma = acos( ( abc % bcd ) / ( nabc * nbcd ) );
506 [ + - ]: 44 : double delta = acos( ( abd % acd ) / ( nabd * nacd ) );
507 [ + - ]: 44 : double epsilon = acos( ( abd % bcd ) / ( nabd * nbcd ) );
508 [ + - ]: 44 : double zeta = acos( ( acd % bcd ) / ( nacd * nbcd ) );
509 : :
510 [ + + ]: 44 : alpha = alpha < beta ? alpha : beta;
511 [ + + ]: 44 : alpha = alpha < gamma ? alpha : gamma;
512 [ + + ]: 44 : alpha = alpha < delta ? alpha : delta;
513 [ + - ]: 44 : alpha = alpha < epsilon ? alpha : epsilon;
514 [ + + ]: 44 : alpha = alpha < zeta ? alpha : zeta;
515 : 44 : alpha *= normal_coeff;
516 : :
517 [ - + ]: 44 : if( alpha < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
518 : :
519 [ + - ][ + - ]: 44 : if( alpha > 0 ) return (double)VERDICT_MIN( alpha, VERDICT_DBL_MAX );
520 [ # # ]: 44 : return (double)VERDICT_MAX( alpha, -VERDICT_DBL_MAX );
521 : : }
522 : :
523 : : /*!
524 : : The collapse ratio of a tet
525 : : */
526 : 44 : C_FUNC_DEF double v_tet_collapse_ratio( int /*num_nodes*/, double coordinates[][3] )
527 : : {
528 : : // Determine side vectors
529 [ + - ][ + - ]: 44 : VerdictVector e01, e02, e03, e12, e13, e23;
[ + - ][ + - ]
[ + - ][ + - ]
530 : :
531 : 88 : e01.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
532 [ + - ]: 44 : coordinates[1][2] - coordinates[0][2] );
533 : :
534 : 88 : e02.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
535 [ + - ]: 44 : coordinates[2][2] - coordinates[0][2] );
536 : :
537 : 88 : e03.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
538 [ + - ]: 44 : coordinates[3][2] - coordinates[0][2] );
539 : :
540 : 88 : e12.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
541 [ + - ]: 44 : coordinates[2][2] - coordinates[1][2] );
542 : :
543 : 88 : e13.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
544 [ + - ]: 44 : coordinates[3][2] - coordinates[1][2] );
545 : :
546 : 88 : e23.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
547 [ + - ]: 44 : coordinates[3][2] - coordinates[2][2] );
548 : :
549 : : double l[6];
550 [ + - ]: 44 : l[0] = e01.length();
551 [ + - ]: 44 : l[1] = e02.length();
552 [ + - ]: 44 : l[2] = e03.length();
553 [ + - ]: 44 : l[3] = e12.length();
554 [ + - ]: 44 : l[4] = e13.length();
555 [ + - ]: 44 : l[5] = e23.length();
556 : :
557 : : // Find longest edge for each bounding triangle of tetrahedron
558 [ + + ]: 44 : double l012 = l[4] > l[0] ? l[4] : l[0];
559 [ + + ]: 44 : l012 = l[1] > l012 ? l[1] : l012;
560 [ + + ]: 44 : double l031 = l[0] > l[2] ? l[0] : l[2];
561 [ + + ]: 44 : l031 = l[3] > l031 ? l[3] : l031;
562 [ + + ]: 44 : double l023 = l[2] > l[1] ? l[2] : l[1];
563 [ + + ]: 44 : l023 = l[5] > l023 ? l[5] : l023;
564 [ + + ]: 44 : double l132 = l[4] > l[3] ? l[4] : l[3];
565 [ + + ]: 44 : l132 = l[5] > l132 ? l[5] : l132;
566 : :
567 : : // Compute collapse ratio for each vertex/triangle pair
568 [ + - ]: 44 : VerdictVector N;
569 : : double h, magN;
570 : : double cr;
571 : : double crMin;
572 : :
573 [ + - ][ + - ]: 44 : N = e01 * e02;
574 [ + - ]: 44 : magN = N.length();
575 [ + - ]: 44 : h = ( e03 % N ) / magN; // height of vertex 3 above 0-1-2
576 : 44 : crMin = h / l012; // ratio of height to longest edge of 0-1-2
577 : :
578 [ + - ][ + - ]: 44 : N = e03 * e01;
579 [ + - ]: 44 : magN = N.length();
580 [ + - ]: 44 : h = ( e02 % N ) / magN; // height of vertex 2 above 0-3-1
581 : 44 : cr = h / l031; // ratio of height to longest edge of 0-3-1
582 [ + + ]: 44 : if( cr < crMin ) crMin = cr;
583 : :
584 [ + - ][ + - ]: 44 : N = e02 * e03;
585 [ + - ]: 44 : magN = N.length();
586 [ + - ]: 44 : h = ( e01 % N ) / magN; // height of vertex 1 above 0-2-3
587 : 44 : cr = h / l023; // ratio of height to longest edge of 0-2-3
588 [ + + ]: 44 : if( cr < crMin ) crMin = cr;
589 : :
590 [ + - ][ + - ]: 44 : N = e12 * e13;
591 [ + - ]: 44 : magN = N.length();
592 [ + - ]: 44 : h = ( e01 % N ) / magN; // height of vertex 0 above 1-3-2
593 : 44 : cr = h / l132; // ratio of height to longest edge of 1-3-2
594 [ + + ]: 44 : if( cr < crMin ) crMin = cr;
595 : :
596 [ - + ]: 44 : if( crMin < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
597 [ + - ][ + - ]: 44 : if( crMin > 0 ) return (double)VERDICT_MIN( crMin, VERDICT_DBL_MAX );
598 [ # # ]: 44 : return (double)VERDICT_MAX( crMin, -VERDICT_DBL_MAX );
599 : : }
600 : :
601 : : /*!
602 : : the volume of a tet
603 : :
604 : : 1/6 * jacobian at a corner node
605 : : */
606 : 135 : C_FUNC_DEF double v_tet_volume( int /*num_nodes*/, double coordinates[][3] )
607 : : {
608 : :
609 : : // Determine side vectors
610 [ + - ][ + - ]: 135 : VerdictVector side0, side2, side3;
[ + - ]
611 : :
612 : 270 : side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
613 [ + - ]: 135 : coordinates[1][2] - coordinates[0][2] );
614 : :
615 : 270 : side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
616 [ + - ]: 135 : coordinates[0][2] - coordinates[2][2] );
617 : :
618 : 270 : side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
619 [ + - ]: 135 : coordinates[3][2] - coordinates[0][2] );
620 : :
621 [ + - ][ + - ]: 135 : return (double)( ( side3 % ( side2 * side0 ) ) / 6.0 );
622 : : }
623 : :
624 : : /*!
625 : : the condition of a tet
626 : :
627 : : condition number of the jacobian matrix at any corner
628 : : */
629 : 23 : C_FUNC_DEF double v_tet_condition( int /*num_nodes*/, double coordinates[][3] )
630 : : {
631 : :
632 : : double condition, term1, term2, det;
633 : 23 : double rt3 = sqrt( 3.0 );
634 : 23 : double rt6 = sqrt( 6.0 );
635 : :
636 [ + - ][ + - ]: 23 : VerdictVector side0, side2, side3;
[ + - ]
637 : :
638 : 46 : side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
639 [ + - ]: 23 : coordinates[1][2] - coordinates[0][2] );
640 : :
641 : 46 : side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
642 [ + - ]: 23 : coordinates[0][2] - coordinates[2][2] );
643 : :
644 : 46 : side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
645 [ + - ]: 23 : coordinates[3][2] - coordinates[0][2] );
646 : :
647 [ + - ][ + - ]: 23 : VerdictVector c_1, c_2, c_3;
[ + - ]
648 : :
649 [ + - ]: 23 : c_1 = side0;
650 [ + - ][ + - ]: 23 : c_2 = ( -2 * side2 - side0 ) / rt3;
[ + - ][ + - ]
651 [ + - ][ + - ]: 23 : c_3 = ( 3 * side3 + side2 - side0 ) / rt6;
[ + - ][ + - ]
[ + - ]
652 : :
653 [ + - ][ + - ]: 23 : term1 = c_1 % c_1 + c_2 % c_2 + c_3 % c_3;
[ + - ]
654 [ + - ][ + - ]: 23 : term2 = ( c_1 * c_2 ) % ( c_1 * c_2 ) + ( c_2 * c_3 ) % ( c_2 * c_3 ) + ( c_1 * c_3 ) % ( c_1 * c_3 );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
655 [ + - ][ + - ]: 23 : det = c_1 % ( c_2 * c_3 );
656 : :
657 [ - + ]: 23 : if( det <= VERDICT_DBL_MIN )
658 : 0 : return VERDICT_DBL_MAX;
659 : : else
660 : 23 : condition = sqrt( term1 * term2 ) / ( 3.0 * det );
661 : :
662 : 23 : return (double)condition;
663 : : }
664 : :
665 : : /*!
666 : : the jacobian of a tet
667 : :
668 : : TODO
669 : : */
670 : 23 : C_FUNC_DEF double v_tet_jacobian( int /*num_nodes*/, double coordinates[][3] )
671 : : {
672 [ + - ][ + - ]: 23 : VerdictVector side0, side2, side3;
[ + - ]
673 : :
674 : 46 : side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
675 [ + - ]: 23 : coordinates[1][2] - coordinates[0][2] );
676 : :
677 : 46 : side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
678 [ + - ]: 23 : coordinates[0][2] - coordinates[2][2] );
679 : :
680 : 46 : side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
681 [ + - ]: 23 : coordinates[3][2] - coordinates[0][2] );
682 : :
683 [ + - ][ + - ]: 23 : return (double)( side3 % ( side2 * side0 ) );
684 : : }
685 : :
686 : : /*!
687 : : the shape of a tet
688 : :
689 : : 3/ condition number of weighted jacobian matrix
690 : : */
691 : 24 : C_FUNC_DEF double v_tet_shape( int /*num_nodes*/, double coordinates[][3] )
692 : : {
693 : :
694 : : static const double two_thirds = 2.0 / 3.0;
695 : : static const double root_of_2 = sqrt( 2.0 );
696 : :
697 [ + - ][ + - ]: 24 : VerdictVector edge0, edge2, edge3;
[ + - ]
698 : :
699 : 48 : edge0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
700 [ + - ]: 24 : coordinates[1][2] - coordinates[0][2] );
701 : :
702 : 48 : edge2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
703 [ + - ]: 24 : coordinates[0][2] - coordinates[2][2] );
704 : :
705 : 48 : edge3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
706 [ + - ]: 24 : coordinates[3][2] - coordinates[0][2] );
707 : :
708 [ + - ][ + - ]: 24 : double jacobian = edge3 % ( edge2 * edge0 );
709 [ - + ]: 24 : if( jacobian < VERDICT_DBL_MIN ) { return (double)0.0; }
710 : 24 : double num = 3 * pow( root_of_2 * jacobian, two_thirds );
711 : : double den =
712 [ + - ][ + - ]: 24 : 1.5 * ( edge0 % edge0 + edge2 % edge2 + edge3 % edge3 ) - ( edge0 % -edge2 + -edge2 % edge3 + edge3 % edge0 );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
713 : :
714 [ - + ]: 24 : if( den < VERDICT_DBL_MIN ) return (double)0.0;
715 : :
716 [ + - ]: 24 : return (double)VERDICT_MAX( num / den, 0 );
717 : : }
718 : :
719 : : /*!
720 : : the relative size of a tet
721 : :
722 : : Min(J,1/J), where J is the determinant of the weighted Jacobian matrix
723 : : */
724 : 24 : C_FUNC_DEF double v_tet_relative_size_squared( int /*num_nodes*/, double coordinates[][3] )
725 : : {
726 : : double size;
727 [ + - ][ + - ]: 24 : VerdictVector w1, w2, w3;
[ + - ]
728 [ + - ]: 24 : get_weight( w1, w2, w3 );
729 [ + - ][ + - ]: 24 : double avg_volume = ( w1 % ( w2 * w3 ) ) / 6.0;
730 : :
731 [ + - ]: 24 : double volume = v_tet_volume( 4, coordinates );
732 : :
733 [ - + ]: 24 : if( avg_volume < VERDICT_DBL_MIN )
734 : 0 : return 0.0;
735 : : else
736 : : {
737 : 24 : size = volume / avg_volume;
738 [ - + ]: 24 : if( size <= VERDICT_DBL_MIN ) return 0.0;
739 [ + + ]: 24 : if( size > 1 ) size = (double)( 1 ) / size;
740 : : }
741 : 24 : return (double)( size * size );
742 : : }
743 : :
744 : : /*!
745 : : the shape and size of a tet
746 : :
747 : : Product of the shape and relative size
748 : : */
749 : 1 : C_FUNC_DEF double v_tet_shape_and_size( int num_nodes, double coordinates[][3] )
750 : : {
751 : :
752 : : double shape, size;
753 : 1 : shape = v_tet_shape( num_nodes, coordinates );
754 : 1 : size = v_tet_relative_size_squared( num_nodes, coordinates );
755 : :
756 : 1 : return (double)( shape * size );
757 : : }
758 : :
759 : : /*!
760 : : the distortion of a tet
761 : : */
762 : 23 : C_FUNC_DEF double v_tet_distortion( int num_nodes, double coordinates[][3] )
763 : : {
764 : :
765 : 23 : double distortion = VERDICT_DBL_MAX;
766 : 23 : int number_of_gauss_points = 0;
767 [ + - ]: 23 : if( num_nodes == 4 )
768 : : // for linear tet, the distortion is always 1 because
769 : : // straight edge tets are the target shape for tet
770 : 23 : return 1.0;
771 : :
772 [ # # ]: 0 : else if( num_nodes == 10 )
773 : : // use four integration points for quadratic tet
774 : 0 : number_of_gauss_points = 4;
775 : :
776 : 0 : int number_dims = 3;
777 : 0 : int total_number_of_gauss_points = number_of_gauss_points;
778 : : // use is_tri=1 to indicate this is for tet in 3D
779 : 0 : int is_tri = 1;
780 : :
781 : : double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
782 : : double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
783 : : double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
784 : : double dndy3[maxTotalNumberGaussPoints][maxNumberNodes];
785 : : double weight[maxTotalNumberGaussPoints];
786 : :
787 : : // create an object of GaussIntegration for tet
788 [ # # ]: 0 : GaussIntegration::initialize( number_of_gauss_points, num_nodes, number_dims, is_tri );
789 [ # # ]: 0 : GaussIntegration::calculate_shape_function_3d_tet();
790 [ # # ]: 0 : GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], dndy3[0], weight );
791 : :
792 : : // vector xxi is the derivative vector of coordinates w.r.t local xi coordinate in the
793 : : // computation space
794 : : // vector xet is the derivative vector of coordinates w.r.t local et coordinate in the
795 : : // computation space
796 : : // vector xze is the derivative vector of coordinates w.r.t local ze coordinate in the
797 : : // computation space
798 [ # # ][ # # ]: 0 : VerdictVector xxi, xet, xze, xin;
[ # # ][ # # ]
799 : :
800 : : double jacobian, minimum_jacobian;
801 : 0 : double element_volume = 0.0;
802 : 0 : minimum_jacobian = VERDICT_DBL_MAX;
803 : :
804 : : // calculate element volume
805 : : int ife, ja;
806 [ # # ]: 0 : for( ife = 0; ife < total_number_of_gauss_points; ife++ )
807 : : {
808 [ # # ]: 0 : xxi.set( 0.0, 0.0, 0.0 );
809 [ # # ]: 0 : xet.set( 0.0, 0.0, 0.0 );
810 [ # # ]: 0 : xze.set( 0.0, 0.0, 0.0 );
811 : :
812 [ # # ]: 0 : for( ja = 0; ja < num_nodes; ja++ )
813 : : {
814 [ # # ]: 0 : xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
815 [ # # ][ # # ]: 0 : xxi += dndy1[ife][ja] * xin;
816 [ # # ][ # # ]: 0 : xet += dndy2[ife][ja] * xin;
817 [ # # ][ # # ]: 0 : xze += dndy3[ife][ja] * xin;
818 : : }
819 : :
820 : : // determinant
821 [ # # ][ # # ]: 0 : jacobian = xxi % ( xet * xze );
822 [ # # ]: 0 : if( minimum_jacobian > jacobian ) minimum_jacobian = jacobian;
823 : :
824 : 0 : element_volume += weight[ife] * jacobian;
825 : : } // element_volume is 6 times the actual volume
826 : :
827 : : // loop through all nodes
828 : : double dndy1_at_node[maxNumberNodes][maxNumberNodes];
829 : : double dndy2_at_node[maxNumberNodes][maxNumberNodes];
830 : : double dndy3_at_node[maxNumberNodes][maxNumberNodes];
831 : :
832 [ # # ]: 0 : GaussIntegration::calculate_derivative_at_nodes_3d_tet( dndy1_at_node, dndy2_at_node, dndy3_at_node );
833 : : int node_id;
834 [ # # ]: 0 : for( node_id = 0; node_id < num_nodes; node_id++ )
835 : : {
836 [ # # ]: 0 : xxi.set( 0.0, 0.0, 0.0 );
837 [ # # ]: 0 : xet.set( 0.0, 0.0, 0.0 );
838 [ # # ]: 0 : xze.set( 0.0, 0.0, 0.0 );
839 : :
840 [ # # ]: 0 : for( ja = 0; ja < num_nodes; ja++ )
841 : : {
842 [ # # ]: 0 : xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
843 [ # # ][ # # ]: 0 : xxi += dndy1_at_node[node_id][ja] * xin;
844 [ # # ][ # # ]: 0 : xet += dndy2_at_node[node_id][ja] * xin;
845 [ # # ][ # # ]: 0 : xze += dndy3_at_node[node_id][ja] * xin;
846 : : }
847 : :
848 [ # # ][ # # ]: 0 : jacobian = xxi % ( xet * xze );
849 [ # # ]: 0 : if( minimum_jacobian > jacobian ) minimum_jacobian = jacobian;
850 : : }
851 : 0 : distortion = minimum_jacobian / element_volume;
852 : :
853 : 23 : return (double)distortion;
854 : : }
855 : :
856 : : /*!
857 : : the quality metrics of a tet
858 : : */
859 : 22 : C_FUNC_DEF void v_tet_quality( int num_nodes, double coordinates[][3], unsigned int metrics_request_flag,
860 : : TetMetricVals* metric_vals )
861 : : {
862 : :
863 : 22 : memset( metric_vals, 0, sizeof( TetMetricVals ) );
864 : :
865 : : /*
866 : :
867 : : node numbers and edge numbers below
868 : :
869 : :
870 : :
871 : : 3
872 : : + edge 0 is node 0 to 1
873 : : +|+ edge 1 is node 1 to 2
874 : : 3/ | \5 edge 2 is node 0 to 2
875 : : / 4| \ edge 3 is node 0 to 3
876 : : 0 - -|- + 2 edge 4 is node 1 to 3
877 : : \ | + edge 5 is node 2 to 3
878 : : 0\ | /1
879 : : +|/ edge 2 is behind edge 4
880 : : 1
881 : :
882 : :
883 : : */
884 : :
885 : : // lets start with making the vectors
886 [ + - ][ + + ]: 154 : VerdictVector edges[6];
887 : 44 : edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
888 [ + - ]: 22 : coordinates[1][2] - coordinates[0][2] );
889 : :
890 : 44 : edges[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
891 [ + - ]: 22 : coordinates[2][2] - coordinates[1][2] );
892 : :
893 : 44 : edges[2].set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
894 [ + - ]: 22 : coordinates[0][2] - coordinates[2][2] );
895 : :
896 : 44 : edges[3].set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
897 [ + - ]: 22 : coordinates[3][2] - coordinates[0][2] );
898 : :
899 : 44 : edges[4].set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
900 [ + - ]: 22 : coordinates[3][2] - coordinates[1][2] );
901 : :
902 : 44 : edges[5].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
903 [ + - ]: 22 : coordinates[3][2] - coordinates[2][2] );
904 : :
905 : : // common numbers
906 : : static const double root_of_2 = sqrt( 2.0 );
907 : :
908 : : // calculate the jacobian
909 : : static const int do_jacobian = V_TET_JACOBIAN | V_TET_VOLUME | V_TET_ASPECT_BETA | V_TET_ASPECT_GAMMA |
910 : : V_TET_SHAPE | V_TET_RELATIVE_SIZE_SQUARED | V_TET_SHAPE_AND_SIZE |
911 : : V_TET_SCALED_JACOBIAN | V_TET_CONDITION;
912 [ + - ][ + - ]: 22 : if( metrics_request_flag & do_jacobian ) { metric_vals->jacobian = (double)( edges[3] % ( edges[2] * edges[0] ) ); }
[ + - ]
913 : :
914 : : // calculate the volume
915 [ + - ]: 22 : if( metrics_request_flag & V_TET_VOLUME ) { metric_vals->volume = (double)( metric_vals->jacobian / 6.0 ); }
916 : :
917 : : // calculate aspect ratio
918 [ + - ]: 22 : if( metrics_request_flag & V_TET_ASPECT_BETA )
919 : : {
920 [ + - ][ + - ]: 22 : double surface_area = ( ( edges[2] * edges[0] ).length() + ( edges[3] * edges[0] ).length() +
[ + - ][ + - ]
921 [ + - ][ + - ]: 22 : ( edges[4] * edges[1] ).length() + ( edges[3] * edges[2] ).length() ) *
[ + - ][ + - ]
922 : 22 : 0.5;
923 : :
924 [ + - ][ + - ]: 22 : VerdictVector numerator = edges[3].length_squared() * ( edges[2] * edges[0] ) +
[ + - ][ + - ]
925 [ + - ][ + - ]: 22 : edges[2].length_squared() * ( edges[3] * edges[0] ) +
[ + - ]
926 [ + - ][ + - ]: 44 : edges[0].length_squared() * ( edges[3] * edges[2] );
[ + - ][ + - ]
927 : :
928 : 22 : double volume = metric_vals->jacobian / 6.0;
929 : :
930 [ - + ]: 22 : if( volume < VERDICT_DBL_MIN )
931 : 0 : metric_vals->aspect_beta = (double)( VERDICT_DBL_MAX );
932 : : else
933 [ + - ]: 22 : metric_vals->aspect_beta = (double)( numerator.length() * surface_area / ( 108 * volume * volume ) );
934 : : }
935 : :
936 : : // calculate the aspect gamma
937 [ + - ]: 22 : if( metrics_request_flag & V_TET_ASPECT_GAMMA )
938 : : {
939 : 22 : double volume = fabs( metric_vals->jacobian / 6.0 );
940 [ - + ]: 22 : if( fabs( volume ) < VERDICT_DBL_MIN )
941 : 0 : metric_vals->aspect_gamma = VERDICT_DBL_MAX;
942 : : else
943 : : {
944 [ + - ][ + - ]: 22 : double srms = sqrt( ( edges[0].length_squared() + edges[1].length_squared() + edges[2].length_squared() +
[ + - ]
945 [ + - ][ + - ]: 22 : edges[3].length_squared() + edges[4].length_squared() + edges[5].length_squared() ) /
[ + - ]
946 : 22 : 6.0 );
947 : :
948 : : // cube the srms
949 : 22 : srms *= ( srms * srms );
950 : 22 : metric_vals->aspect_gamma = (double)( srms / ( 8.48528137423857 * volume ) );
951 : : }
952 : : }
953 : :
954 : : // calculate the shape of the tet
955 [ + - ]: 22 : if( metrics_request_flag & ( V_TET_SHAPE | V_TET_SHAPE_AND_SIZE ) )
956 : : {
957 : : // if the jacobian is non-positive, the shape is 0
958 [ - + ]: 22 : if( metric_vals->jacobian < VERDICT_DBL_MIN ) { metric_vals->shape = (double)0.0; }
959 : : else
960 : : {
961 : : static const double two_thirds = 2.0 / 3.0;
962 : 22 : double num = 3.0 * pow( root_of_2 * metric_vals->jacobian, two_thirds );
963 [ + - ][ + - ]: 22 : double den = 1.5 * ( edges[0] % edges[0] + edges[2] % edges[2] + edges[3] % edges[3] ) -
[ + - ]
964 [ + - ][ + - ]: 22 : ( edges[0] % -edges[2] + -edges[2] % edges[3] + edges[3] % edges[0] );
[ + - ][ + - ]
[ + - ]
965 : :
966 [ - + ]: 22 : if( den < VERDICT_DBL_MIN )
967 : 0 : metric_vals->shape = (double)0.0;
968 : : else
969 [ + - ]: 22 : metric_vals->shape = (double)VERDICT_MAX( num / den, 0 );
970 : : }
971 : : }
972 : :
973 : : // calculate the relative size of the tet
974 [ + - ]: 22 : if( metrics_request_flag & ( V_TET_RELATIVE_SIZE_SQUARED | V_TET_SHAPE_AND_SIZE ) )
975 : : {
976 [ + - ][ + - ]: 22 : VerdictVector w1, w2, w3;
[ + - ]
977 [ + - ]: 22 : get_weight( w1, w2, w3 );
978 [ + - ][ + - ]: 22 : double avg_vol = ( w1 % ( w2 * w3 ) ) / 6;
979 : :
980 [ - + ]: 22 : if( avg_vol < VERDICT_DBL_MIN )
981 : 0 : metric_vals->relative_size_squared = 0.0;
982 : : else
983 : : {
984 : 22 : double tmp = metric_vals->jacobian / ( 6 * avg_vol );
985 [ - + ]: 22 : if( tmp < VERDICT_DBL_MIN )
986 : 0 : metric_vals->relative_size_squared = 0.0;
987 : : else
988 : : {
989 : 22 : tmp *= tmp;
990 [ + - ]: 22 : metric_vals->relative_size_squared = (double)VERDICT_MIN( tmp, 1 / tmp );
991 : : }
992 : : }
993 : : }
994 : :
995 : : // calculate the shape and size
996 [ + - ]: 22 : if( metrics_request_flag & V_TET_SHAPE_AND_SIZE )
997 : 22 : { metric_vals->shape_and_size = (double)( metric_vals->shape * metric_vals->relative_size_squared ); }
998 : :
999 : : // calculate the scaled jacobian
1000 [ + - ]: 22 : if( metrics_request_flag & V_TET_SCALED_JACOBIAN )
1001 : : {
1002 : : // find out which node the normalized jacobian can be calculated at
1003 : : // and it will be the smaller than at other nodes
1004 [ + - ][ + - ]: 22 : double length_squared[4] = { edges[0].length_squared() * edges[2].length_squared() * edges[3].length_squared(),
[ + - ]
1005 [ + - ][ + - ]: 22 : edges[0].length_squared() * edges[1].length_squared() * edges[4].length_squared(),
[ + - ]
1006 [ + - ][ + - ]: 22 : edges[1].length_squared() * edges[2].length_squared() * edges[5].length_squared(),
[ + - ]
1007 [ + - ][ + - ]: 22 : edges[3].length_squared() * edges[4].length_squared() *
1008 [ + - ]: 22 : edges[5].length_squared() };
1009 : :
1010 : 22 : int which_node = 0;
1011 [ + + ]: 22 : if( length_squared[1] > length_squared[which_node] ) which_node = 1;
1012 [ + + ]: 22 : if( length_squared[2] > length_squared[which_node] ) which_node = 2;
1013 [ + + ]: 22 : if( length_squared[3] > length_squared[which_node] ) which_node = 3;
1014 : :
1015 : : // find the scaled jacobian at this node
1016 : 22 : double length_product = sqrt( length_squared[which_node] );
1017 [ - + ]: 22 : if( length_product < fabs( metric_vals->jacobian ) ) length_product = fabs( metric_vals->jacobian );
1018 : :
1019 [ - + ]: 22 : if( length_product < VERDICT_DBL_MIN )
1020 : 0 : metric_vals->scaled_jacobian = (double)VERDICT_DBL_MAX;
1021 : : else
1022 : 22 : metric_vals->scaled_jacobian = (double)( root_of_2 * metric_vals->jacobian / length_product );
1023 : : }
1024 : :
1025 : : // calculate the condition number
1026 [ + - ]: 22 : if( metrics_request_flag & V_TET_CONDITION )
1027 : : {
1028 : : static const double root_of_3 = sqrt( 3.0 );
1029 : : static const double root_of_6 = sqrt( 6.0 );
1030 : :
1031 [ + - ][ + - ]: 22 : VerdictVector c_1, c_2, c_3;
[ + - ]
1032 [ + - ]: 22 : c_1 = edges[0];
1033 [ + - ][ + - ]: 22 : c_2 = ( -2 * edges[2] - edges[0] ) / root_of_3;
[ + - ][ + - ]
1034 [ + - ][ + - ]: 22 : c_3 = ( 3 * edges[3] + edges[2] - edges[0] ) / root_of_6;
[ + - ][ + - ]
[ + - ]
1035 : :
1036 [ + - ][ + - ]: 22 : double term1 = c_1 % c_1 + c_2 % c_2 + c_3 % c_3;
[ + - ]
1037 [ + - ][ + - ]: 22 : double term2 = ( c_1 * c_2 ) % ( c_1 * c_2 ) + ( c_2 * c_3 ) % ( c_2 * c_3 ) + ( c_3 * c_1 ) % ( c_3 * c_1 );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
1038 : :
1039 [ + - ][ + - ]: 22 : double det = c_1 % ( c_2 * c_3 );
1040 : :
1041 [ - + ]: 22 : if( det <= VERDICT_DBL_MIN )
1042 : 0 : metric_vals->condition = (double)VERDICT_DBL_MAX;
1043 : : else
1044 : 22 : metric_vals->condition = (double)( sqrt( term1 * term2 ) / ( 3.0 * det ) );
1045 : : }
1046 : :
1047 : : // calculate the distortion
1048 [ + - ]: 22 : if( metrics_request_flag & V_TET_DISTORTION )
1049 [ + - ]: 22 : { metric_vals->distortion = v_tet_distortion( num_nodes, coordinates ); }
1050 : :
1051 : : // check for overflow
1052 [ + - ]: 22 : if( metrics_request_flag & V_TET_ASPECT_BETA )
1053 : : {
1054 [ + - ]: 22 : if( metric_vals->aspect_beta > 0 )
1055 [ + - ]: 22 : metric_vals->aspect_beta = (double)VERDICT_MIN( metric_vals->aspect_beta, VERDICT_DBL_MAX );
1056 [ + - ]: 22 : metric_vals->aspect_beta = (double)VERDICT_MAX( metric_vals->aspect_beta, -VERDICT_DBL_MAX );
1057 : : }
1058 : :
1059 [ + - ]: 22 : if( metrics_request_flag & V_TET_ASPECT_GAMMA )
1060 : : {
1061 [ + - ]: 22 : if( metric_vals->aspect_gamma > 0 )
1062 [ + - ]: 22 : metric_vals->aspect_gamma = (double)VERDICT_MIN( metric_vals->aspect_gamma, VERDICT_DBL_MAX );
1063 [ + - ]: 22 : metric_vals->aspect_gamma = (double)VERDICT_MAX( metric_vals->aspect_gamma, -VERDICT_DBL_MAX );
1064 : : }
1065 : :
1066 [ + - ]: 22 : if( metrics_request_flag & V_TET_VOLUME )
1067 : : {
1068 [ + - ][ + - ]: 22 : if( metric_vals->volume > 0 ) metric_vals->volume = (double)VERDICT_MIN( metric_vals->volume, VERDICT_DBL_MAX );
1069 [ + - ]: 22 : metric_vals->volume = (double)VERDICT_MAX( metric_vals->volume, -VERDICT_DBL_MAX );
1070 : : }
1071 : :
1072 [ + - ]: 22 : if( metrics_request_flag & V_TET_CONDITION )
1073 : : {
1074 [ + - ]: 22 : if( metric_vals->condition > 0 )
1075 [ + - ]: 22 : metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
1076 [ + - ]: 22 : metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
1077 : : }
1078 : :
1079 [ + - ]: 22 : if( metrics_request_flag & V_TET_JACOBIAN )
1080 : : {
1081 [ + - ]: 22 : if( metric_vals->jacobian > 0 )
1082 [ + - ]: 22 : metric_vals->jacobian = (double)VERDICT_MIN( metric_vals->jacobian, VERDICT_DBL_MAX );
1083 [ + - ]: 22 : metric_vals->jacobian = (double)VERDICT_MAX( metric_vals->jacobian, -VERDICT_DBL_MAX );
1084 : : }
1085 : :
1086 [ + - ]: 22 : if( metrics_request_flag & V_TET_SCALED_JACOBIAN )
1087 : : {
1088 [ + - ]: 22 : if( metric_vals->scaled_jacobian > 0 )
1089 [ + - ]: 22 : metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
1090 [ + - ]: 22 : metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
1091 : : }
1092 : :
1093 [ + - ]: 22 : if( metrics_request_flag & V_TET_SHAPE )
1094 : : {
1095 [ + - ][ + - ]: 22 : if( metric_vals->shape > 0 ) metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
1096 [ + - ]: 22 : metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
1097 : : }
1098 : :
1099 [ + - ]: 22 : if( metrics_request_flag & V_TET_RELATIVE_SIZE_SQUARED )
1100 : : {
1101 [ + - ]: 22 : if( metric_vals->relative_size_squared > 0 )
1102 : : metric_vals->relative_size_squared =
1103 [ + - ]: 22 : (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
1104 : : metric_vals->relative_size_squared =
1105 [ + - ]: 22 : (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
1106 : : }
1107 : :
1108 [ + - ]: 22 : if( metrics_request_flag & V_TET_SHAPE_AND_SIZE )
1109 : : {
1110 [ + - ]: 22 : if( metric_vals->shape_and_size > 0 )
1111 [ + - ]: 22 : metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
1112 [ + - ]: 22 : metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
1113 : : }
1114 : :
1115 [ + - ]: 22 : if( metrics_request_flag & V_TET_DISTORTION )
1116 : : {
1117 [ + - ]: 22 : if( metric_vals->distortion > 0 )
1118 [ + - ]: 22 : metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
1119 [ + - ]: 22 : metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
1120 : : }
1121 : :
1122 [ + - ][ + - ]: 22 : if( metrics_request_flag & V_TET_ASPECT_RATIO ) metric_vals->aspect_ratio = v_tet_aspect_ratio( 4, coordinates );
1123 : :
1124 [ + - ]: 22 : if( metrics_request_flag & V_TET_ASPECT_FROBENIUS )
1125 [ + - ]: 22 : metric_vals->aspect_frobenius = v_tet_aspect_frobenius( 4, coordinates );
1126 : :
1127 [ + - ][ + - ]: 22 : if( metrics_request_flag & V_TET_MINIMUM_ANGLE ) metric_vals->minimum_angle = v_tet_minimum_angle( 4, coordinates );
1128 : :
1129 [ + - ]: 22 : if( metrics_request_flag & V_TET_COLLAPSE_RATIO )
1130 [ + - ]: 22 : metric_vals->collapse_ratio = v_tet_collapse_ratio( 4, coordinates );
1131 : :
1132 [ + - ][ + - ]: 22 : if( metrics_request_flag & V_TET_RADIUS_RATIO ) metric_vals->radius_ratio = v_tet_radius_ratio( 4, coordinates );
1133 : 22 : }
|
