Branch data Line data Source code
1 : : /*=========================================================================
2 : :
3 : : Module: $RCSfile: V_TriMetric.cpp,v $
4 : :
5 : : Copyright (c) 2007 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 : : * TriMetric.cpp contains quality calculations for Tris
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 "V_GaussIntegration.hpp"
28 : : #include "VerdictVector.hpp"
29 : : #include <memory.h>
30 : : #include <stddef.h>
31 : :
32 : : // the average area of a tri
33 : : static double verdict_tri_size = 0;
34 : : static ComputeNormal compute_normal = NULL;
35 : :
36 : : /*!
37 : : get weights based on the average area of a set of
38 : : tris
39 : : */
40 : 56 : static int v_tri_get_weight( double& m11, double& m21, double& m12, double& m22 )
41 : : {
42 : : static const double rootOf3 = sqrt( 3.0 );
43 : 56 : m11 = 1;
44 : 56 : m21 = 0;
45 : 56 : m12 = 0.5;
46 : 56 : m22 = 0.5 * rootOf3;
47 : 56 : double scale = sqrt( 2.0 * verdict_tri_size / ( m11 * m22 - m21 * m12 ) );
48 : :
49 : 56 : m11 *= scale;
50 : 56 : m21 *= scale;
51 : 56 : m12 *= scale;
52 : 56 : m22 *= scale;
53 : :
54 : 56 : return 1;
55 : : }
56 : :
57 : : /*! sets the average area of a tri */
58 : 1 : C_FUNC_DEF void v_set_tri_size( double size )
59 : : {
60 : 1 : verdict_tri_size = size;
61 : 1 : }
62 : :
63 : 0 : C_FUNC_DEF void v_set_tri_normal_func( ComputeNormal func )
64 : : {
65 : 0 : compute_normal = func;
66 : 0 : }
67 : :
68 : : /*!
69 : : the edge ratio of a triangle
70 : :
71 : : NB (P. Pebay 01/14/07):
72 : : Hmax / Hmin where Hmax and Hmin are respectively the maximum and the
73 : : minimum edge lengths
74 : :
75 : : */
76 : 36 : C_FUNC_DEF double v_tri_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
77 : : {
78 : :
79 : : // three vectors for each side
80 : 72 : VerdictVector a( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
81 [ + - ]: 36 : coordinates[1][2] - coordinates[0][2] );
82 : :
83 : 72 : VerdictVector b( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
84 [ + - ]: 36 : coordinates[2][2] - coordinates[1][2] );
85 : :
86 : 72 : VerdictVector c( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
87 [ + - ]: 36 : coordinates[0][2] - coordinates[2][2] );
88 : :
89 [ + - ]: 36 : double a2 = a.length_squared();
90 [ + - ]: 36 : double b2 = b.length_squared();
91 [ + - ]: 36 : double c2 = c.length_squared();
92 : :
93 : : double m2, M2;
94 [ + + ]: 36 : if( a2 < b2 )
95 : : {
96 [ - + ]: 8 : if( b2 < c2 )
97 : : {
98 : 0 : m2 = a2;
99 : 0 : M2 = c2;
100 : : }
101 : : else // b2 <= a2
102 : : {
103 [ - + ]: 8 : if( a2 < c2 )
104 : : {
105 : 0 : m2 = a2;
106 : 0 : M2 = b2;
107 : : }
108 : : else // c2 <= a2
109 : : {
110 : 8 : m2 = c2;
111 : 8 : M2 = b2;
112 : : }
113 : : }
114 : : }
115 : : else // b2 <= a2
116 : : {
117 [ + + ]: 28 : if( a2 < c2 )
118 : : {
119 : 18 : m2 = b2;
120 : 18 : M2 = c2;
121 : : }
122 : : else // c2 <= a2
123 : : {
124 [ - + ]: 10 : if( b2 < c2 )
125 : : {
126 : 0 : m2 = b2;
127 : 0 : M2 = a2;
128 : : }
129 : : else // c2 <= b2
130 : : {
131 : 10 : m2 = c2;
132 : 10 : M2 = a2;
133 : : }
134 : : }
135 : : }
136 : :
137 [ - + ]: 36 : if( m2 < VERDICT_DBL_MIN )
138 : 0 : return (double)VERDICT_DBL_MAX;
139 : : else
140 : : {
141 : : double edge_ratio;
142 : 36 : edge_ratio = sqrt( M2 / m2 );
143 : :
144 [ + - ][ + - ]: 36 : if( edge_ratio > 0 ) return (double)VERDICT_MIN( edge_ratio, VERDICT_DBL_MAX );
145 [ # # ]: 36 : return (double)VERDICT_MAX( edge_ratio, -VERDICT_DBL_MAX );
146 : : }
147 : : }
148 : :
149 : : /*!
150 : : the aspect ratio of a triangle
151 : :
152 : : NB (P. Pebay 01/14/07):
153 : : Hmax / ( 2.0 * sqrt(3.0) * IR) where Hmax is the maximum edge length
154 : : and IR is the inradius
155 : :
156 : : note that previous incarnations of verdict used "v_tri_aspect_ratio" to denote
157 : : what is now called "v_tri_aspect_frobenius"
158 : :
159 : : */
160 : 18 : C_FUNC_DEF double v_tri_aspect_ratio( int /*num_nodes*/, double coordinates[][3] )
161 : : {
162 : : static const double normal_coeff = sqrt( 3. ) / 6.;
163 : :
164 : : // three vectors for each side
165 : 36 : VerdictVector a( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
166 [ + - ]: 18 : coordinates[1][2] - coordinates[0][2] );
167 : :
168 : 36 : VerdictVector b( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
169 [ + - ]: 18 : coordinates[2][2] - coordinates[1][2] );
170 : :
171 : 36 : VerdictVector c( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
172 [ + - ]: 18 : coordinates[0][2] - coordinates[2][2] );
173 : :
174 [ + - ]: 18 : double a1 = a.length();
175 [ + - ]: 18 : double b1 = b.length();
176 [ + - ]: 18 : double c1 = c.length();
177 : :
178 [ + + ]: 18 : double hm = a1 > b1 ? a1 : b1;
179 [ + + ]: 18 : hm = hm > c1 ? hm : c1;
180 : :
181 [ + - ]: 18 : VerdictVector ab = a * b;
182 [ + - ]: 18 : double denominator = ab.length();
183 : :
184 [ - + ]: 18 : if( denominator < VERDICT_DBL_MIN )
185 : 0 : return (double)VERDICT_DBL_MAX;
186 : : else
187 : : {
188 : : double aspect_ratio;
189 : 18 : aspect_ratio = normal_coeff * hm * ( a1 + b1 + c1 ) / denominator;
190 : :
191 [ + - ][ + - ]: 18 : if( aspect_ratio > 0 ) return (double)VERDICT_MIN( aspect_ratio, VERDICT_DBL_MAX );
192 [ # # ]: 18 : return (double)VERDICT_MAX( aspect_ratio, -VERDICT_DBL_MAX );
193 : : }
194 : : }
195 : :
196 : : /*!
197 : : the radius ratio of a triangle
198 : :
199 : : NB (P. Pebay 01/13/07):
200 : : CR / (3.0*IR) where CR is the circumradius and IR is the inradius
201 : :
202 : : this quality metric is also known to VERDICT, for tetrahedral elements only,
203 : : a the "aspect beta"
204 : :
205 : : */
206 : 54 : C_FUNC_DEF double v_tri_radius_ratio( int /*num_nodes*/, double coordinates[][3] )
207 : : {
208 : :
209 : : // three vectors for each side
210 : 108 : VerdictVector a( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
211 [ + - ]: 54 : coordinates[1][2] - coordinates[0][2] );
212 : :
213 : 108 : VerdictVector b( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
214 [ + - ]: 54 : coordinates[2][2] - coordinates[1][2] );
215 : :
216 : 108 : VerdictVector c( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
217 [ + - ]: 54 : coordinates[0][2] - coordinates[2][2] );
218 : :
219 [ + - ]: 54 : double a2 = a.length_squared();
220 [ + - ]: 54 : double b2 = b.length_squared();
221 [ + - ]: 54 : double c2 = c.length_squared();
222 : :
223 [ + - ]: 54 : VerdictVector ab = a * b;
224 [ + - ]: 54 : double denominator = ab.length_squared();
225 : :
226 [ - + ]: 54 : if( denominator < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
227 : :
228 : : double radius_ratio;
229 : 54 : radius_ratio = .25 * a2 * b2 * c2 * ( a2 + b2 + c2 ) / denominator;
230 : :
231 [ + - ][ + - ]: 54 : if( radius_ratio > 0 ) return (double)VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
232 [ # # ]: 54 : return (double)VERDICT_MAX( radius_ratio, -VERDICT_DBL_MAX );
233 : : }
234 : :
235 : : /*!
236 : : the Frobenius aspect of a tri
237 : :
238 : : srms^2/(2 * sqrt(3.0) * area)
239 : : where srms^2 is sum of the lengths squared
240 : :
241 : : NB (P. Pebay 01/14/07):
242 : : this method was called "aspect ratio" in earlier incarnations of VERDICT
243 : :
244 : : */
245 : :
246 : 18 : C_FUNC_DEF double v_tri_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
247 : : {
248 : : static const double two_times_root_of_3 = 2 * sqrt( 3.0 );
249 : :
250 : : // three vectors for each side
251 : 36 : VerdictVector side1( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
252 [ + - ]: 18 : coordinates[1][2] - coordinates[0][2] );
253 : :
254 : 36 : VerdictVector side2( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
255 [ + - ]: 18 : coordinates[2][2] - coordinates[1][2] );
256 : :
257 : 36 : VerdictVector side3( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
258 [ + - ]: 18 : coordinates[0][2] - coordinates[2][2] );
259 : :
260 : : // sum the lengths squared of each side
261 [ + - ][ + - ]: 18 : double srms = ( side1.length_squared() + side2.length_squared() + side3.length_squared() );
[ + - ]
262 : :
263 : : // find two times the area of the triangle by cross product
264 [ + - ][ + - ]: 18 : double areaX2 = ( ( side1 * ( -side3 ) ).length() );
[ + - ]
265 : :
266 [ - + ]: 18 : if( areaX2 == 0.0 ) return (double)VERDICT_DBL_MAX;
267 : :
268 : 18 : double aspect = (double)( srms / ( two_times_root_of_3 * ( areaX2 ) ) );
269 [ + - ][ + - ]: 18 : if( aspect > 0 ) return (double)VERDICT_MIN( aspect, VERDICT_DBL_MAX );
270 [ # # ]: 18 : return (double)VERDICT_MAX( aspect, -VERDICT_DBL_MAX );
271 : : }
272 : :
273 : : /*!
274 : : The area of a tri
275 : :
276 : : 0.5 * jacobian at a node
277 : : */
278 : 19 : C_FUNC_DEF double v_tri_area( int /*num_nodes*/, double coordinates[][3] )
279 : : {
280 : : // two vectors for two sides
281 : 38 : VerdictVector side1( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
282 [ + - ]: 19 : coordinates[1][2] - coordinates[0][2] );
283 : :
284 : 38 : VerdictVector side3( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
285 [ + - ]: 19 : coordinates[2][2] - coordinates[0][2] );
286 : :
287 : : // the cross product of the two vectors representing two sides of the
288 : : // triangle
289 [ + - ]: 19 : VerdictVector tmp = side1 * side3;
290 : :
291 : : // return the magnitude of the vector divided by two
292 [ + - ]: 19 : double area = 0.5 * tmp.length();
293 [ + - ][ + - ]: 19 : if( area > 0 ) return (double)VERDICT_MIN( area, VERDICT_DBL_MAX );
294 [ # # ]: 19 : return (double)VERDICT_MAX( area, -VERDICT_DBL_MAX );
295 : : }
296 : :
297 : : /*!
298 : : The minimum angle of a tri
299 : :
300 : : The minimum angle of a tri is the minimum angle between
301 : : two adjacents sides out of all three corners of the triangle.
302 : : */
303 : 19 : C_FUNC_DEF double v_tri_minimum_angle( int /*num_nodes*/, double coordinates[][3] )
304 : : {
305 : :
306 : : // vectors for all the sides
307 [ + - ][ + + ]: 95 : VerdictVector sides[4];
308 : 38 : sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
309 [ + - ]: 19 : coordinates[1][2] - coordinates[0][2] );
310 : 38 : sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
311 [ + - ]: 19 : coordinates[2][2] - coordinates[1][2] );
312 : 38 : sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
313 [ + - ]: 19 : coordinates[2][2] - coordinates[0][2] );
314 : :
315 : : // in case we need to find the interior angle
316 : : // between sides 0 and 1
317 [ + - ][ + - ]: 19 : sides[3] = -sides[1];
318 : :
319 : : // calculate the lengths squared of the sides
320 : : double sides_lengths[3];
321 [ + - ]: 19 : sides_lengths[0] = sides[0].length_squared();
322 [ + - ]: 19 : sides_lengths[1] = sides[1].length_squared();
323 [ + - ]: 19 : sides_lengths[2] = sides[2].length_squared();
324 : :
325 [ + - ][ + - ]: 19 : if( sides_lengths[0] == 0.0 || sides_lengths[1] == 0.0 || sides_lengths[2] == 0.0 ) return 0.0;
[ - + ]
326 : :
327 : : // using the law of sines, we know that the minimum
328 : : // angle is opposite of the shortest side
329 : :
330 : : // find the shortest side
331 : 19 : int short_side = 0;
332 [ + + ]: 19 : if( sides_lengths[1] < sides_lengths[0] ) short_side = 1;
333 [ - + ]: 19 : if( sides_lengths[2] < sides_lengths[short_side] ) short_side = 2;
334 : :
335 : : // from the shortest side, calculate the angle of the
336 : : // opposite angle
337 : : double min_angle;
338 [ + + ][ + - ]: 19 : if( short_side == 0 ) { min_angle = sides[2].interior_angle( sides[1] ); }
339 [ + - ]: 6 : else if( short_side == 1 )
340 : : {
341 [ + - ]: 6 : min_angle = sides[0].interior_angle( sides[2] );
342 : : }
343 : : else
344 : : {
345 [ # # ]: 0 : min_angle = sides[0].interior_angle( sides[3] );
346 : : }
347 : :
348 [ + - ][ + - ]: 19 : if( min_angle > 0 ) return (double)VERDICT_MIN( min_angle, VERDICT_DBL_MAX );
349 [ # # ]: 19 : return (double)VERDICT_MAX( min_angle, -VERDICT_DBL_MAX );
350 : : }
351 : :
352 : : /*!
353 : : The maximum angle of a tri
354 : :
355 : : The maximum angle of a tri is the maximum angle between
356 : : two adjacents sides out of all three corners of the triangle.
357 : : */
358 : 19 : C_FUNC_DEF double v_tri_maximum_angle( int /*num_nodes*/, double coordinates[][3] )
359 : : {
360 : :
361 : : // vectors for all the sides
362 [ + - ][ + + ]: 95 : VerdictVector sides[4];
363 : 38 : sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
364 [ + - ]: 19 : coordinates[1][2] - coordinates[0][2] );
365 : 38 : sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
366 [ + - ]: 19 : coordinates[2][2] - coordinates[1][2] );
367 : 38 : sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
368 [ + - ]: 19 : coordinates[2][2] - coordinates[0][2] );
369 : :
370 : : // in case we need to find the interior angle
371 : : // between sides 0 and 1
372 [ + - ][ + - ]: 19 : sides[3] = -sides[1];
373 : :
374 : : // calculate the lengths squared of the sides
375 : : double sides_lengths[3];
376 [ + - ]: 19 : sides_lengths[0] = sides[0].length_squared();
377 [ + - ]: 19 : sides_lengths[1] = sides[1].length_squared();
378 [ + - ]: 19 : sides_lengths[2] = sides[2].length_squared();
379 : :
380 [ + - ][ + - ]: 19 : if( sides_lengths[0] == 0.0 || sides_lengths[1] == 0.0 || sides_lengths[2] == 0.0 ) { return 0.0; }
[ - + ]
381 : :
382 : : // using the law of sines, we know that the maximum
383 : : // angle is opposite of the longest side
384 : :
385 : : // find the longest side
386 : 19 : int short_side = 0;
387 [ + + ]: 19 : if( sides_lengths[1] > sides_lengths[0] ) short_side = 1;
388 [ + + ]: 19 : if( sides_lengths[2] > sides_lengths[short_side] ) short_side = 2;
389 : :
390 : : // from the longest side, calculate the angle of the
391 : : // opposite angle
392 : : double max_angle;
393 [ + + ][ + - ]: 19 : if( short_side == 0 ) { max_angle = sides[2].interior_angle( sides[1] ); }
394 [ + + ]: 13 : else if( short_side == 1 )
395 : : {
396 [ + - ]: 4 : max_angle = sides[0].interior_angle( sides[2] );
397 : : }
398 : : else
399 : : {
400 [ + - ]: 9 : max_angle = sides[0].interior_angle( sides[3] );
401 : : }
402 : :
403 [ + - ][ + - ]: 19 : if( max_angle > 0 ) return (double)VERDICT_MIN( max_angle, VERDICT_DBL_MAX );
404 [ # # ]: 19 : return (double)VERDICT_MAX( max_angle, -VERDICT_DBL_MAX );
405 : : }
406 : :
407 : : /*!
408 : : The condition of a tri
409 : :
410 : : Condition number of the jacobian matrix at any corner
411 : : */
412 : 57 : C_FUNC_DEF double v_tri_condition( int /*num_nodes*/, double coordinates[][3] )
413 : : {
414 : : static const double rt3 = sqrt( 3.0 );
415 : :
416 : 114 : VerdictVector v1( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
417 [ + - ]: 57 : coordinates[1][2] - coordinates[0][2] );
418 : :
419 : 114 : VerdictVector v2( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
420 [ + - ]: 57 : coordinates[2][2] - coordinates[0][2] );
421 : :
422 [ + - ]: 57 : VerdictVector tri_normal = v1 * v2;
423 [ + - ]: 57 : double areax2 = tri_normal.length();
424 : :
425 [ - + ]: 57 : if( areax2 == 0.0 ) return (double)VERDICT_DBL_MAX;
426 : :
427 [ + - ][ + - ]: 57 : double condition = (double)( ( ( v1 % v1 ) + ( v2 % v2 ) - ( v1 % v2 ) ) / ( areax2 * rt3 ) );
[ + - ]
428 : :
429 : : // check for inverted if we have access to the normal
430 [ - + ]: 57 : if( compute_normal )
431 : : {
432 : : // center of tri
433 : : double point[3], surf_normal[3];
434 : 0 : point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
435 : 0 : point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
436 : 0 : point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
437 : :
438 : : // dot product
439 [ # # ]: 0 : compute_normal( point, surf_normal );
440 [ # # ][ # # ]: 0 : if( ( tri_normal.x() * surf_normal[0] + tri_normal.y() * surf_normal[1] + tri_normal.z() * surf_normal[2] ) <
[ # # ][ # # ]
441 : : 0 )
442 : 0 : return (double)VERDICT_DBL_MAX;
443 : : }
444 [ + - ]: 57 : return (double)VERDICT_MIN( condition, VERDICT_DBL_MAX );
445 : : }
446 : :
447 : : /*!
448 : : The scaled jacobian of a tri
449 : :
450 : : minimum of the jacobian divided by the lengths of 2 edge vectors
451 : : */
452 : 19 : C_FUNC_DEF double v_tri_scaled_jacobian( int /*num_nodes*/, double coordinates[][3] )
453 : : {
454 : : static const double detw = 2. / sqrt( 3.0 );
455 [ + - ][ + - ]: 19 : VerdictVector first, second;
456 : : double jacobian;
457 : :
458 [ + - ][ + + ]: 76 : VerdictVector edge[3];
459 : 38 : edge[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
460 [ + - ]: 19 : coordinates[1][2] - coordinates[0][2] );
461 : :
462 : 38 : edge[1].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
463 [ + - ]: 19 : coordinates[2][2] - coordinates[0][2] );
464 : :
465 : 38 : edge[2].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
466 [ + - ]: 19 : coordinates[2][2] - coordinates[1][2] );
467 [ + - ][ + - ]: 19 : first = edge[1] - edge[0];
468 [ + - ][ + - ]: 19 : second = edge[2] - edge[0];
469 : :
470 [ + - ]: 19 : VerdictVector cross = first * second;
471 [ + - ]: 19 : jacobian = cross.length();
472 : :
473 : : double max_edge_length_product;
474 : : max_edge_length_product =
475 [ + - ][ + - ]: 57 : VERDICT_MAX( edge[0].length() * edge[1].length(),
[ + - ][ + - ]
[ + - ][ + - ]
[ + + ][ + - ]
[ + - ][ + - ]
[ + - ][ # # ]
[ # # ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
476 [ - + ][ + + ]: 57 : VERDICT_MAX( edge[1].length() * edge[2].length(), edge[0].length() * edge[2].length() ) );
477 : :
478 [ - + ]: 19 : if( max_edge_length_product < VERDICT_DBL_MIN ) return (double)0.0;
479 : :
480 : 19 : jacobian *= detw;
481 : 19 : jacobian /= max_edge_length_product;
482 : :
483 [ - + ]: 19 : if( compute_normal )
484 : : {
485 : : // center of tri
486 : : double point[3], surf_normal[3];
487 : 0 : point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
488 : 0 : point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
489 : 0 : point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
490 : :
491 : : // dot product
492 [ # # ]: 0 : compute_normal( point, surf_normal );
493 [ # # ][ # # ]: 0 : if( ( cross.x() * surf_normal[0] + cross.y() * surf_normal[1] + cross.z() * surf_normal[2] ) < 0 )
[ # # ][ # # ]
494 : 0 : jacobian *= -1;
495 : : }
496 : :
497 [ + - ][ + - ]: 19 : if( jacobian > 0 ) return (double)VERDICT_MIN( jacobian, VERDICT_DBL_MAX );
498 [ # # ]: 19 : return (double)VERDICT_MAX( jacobian, -VERDICT_DBL_MAX );
499 : : }
500 : :
501 : : /*!
502 : : The shape of a tri
503 : :
504 : : 2 / condition number of weighted jacobian matrix
505 : : */
506 : 38 : C_FUNC_DEF double v_tri_shape( int num_nodes, double coordinates[][3] )
507 : : {
508 : 38 : double condition = v_tri_condition( num_nodes, coordinates );
509 : :
510 : : double shape;
511 [ - + ]: 38 : if( condition <= VERDICT_DBL_MIN )
512 : 0 : shape = VERDICT_DBL_MAX;
513 : : else
514 : 38 : shape = ( 1 / condition );
515 : :
516 [ + - ][ + - ]: 38 : if( shape > 0 ) return (double)VERDICT_MIN( shape, VERDICT_DBL_MAX );
517 [ # # ]: 0 : return (double)VERDICT_MAX( shape, -VERDICT_DBL_MAX );
518 : : }
519 : :
520 : : /*!
521 : : The relative size of a tri
522 : :
523 : : Min(J,1/J) where J is the determinant of the weighted jacobian matrix.
524 : : */
525 : 38 : C_FUNC_DEF double v_tri_relative_size_squared( int /*num_nodes*/, double coordinates[][3] )
526 : : {
527 : : double w11, w21, w12, w22;
528 : :
529 [ + - ][ + - ]: 38 : VerdictVector xxi, xet, tri_normal;
[ + - ]
530 : :
531 : 38 : v_tri_get_weight( w11, w21, w12, w22 );
532 : :
533 [ + - ]: 38 : double detw = determinant( w11, w21, w12, w22 );
534 : :
535 [ - + ]: 38 : if( detw == 0.0 ) return 0.0;
536 : :
537 : 76 : xxi.set( coordinates[0][0] - coordinates[1][0], coordinates[0][1] - coordinates[1][1],
538 [ + - ]: 38 : coordinates[0][2] - coordinates[1][2] );
539 : :
540 : 76 : xet.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
541 [ + - ]: 38 : coordinates[0][2] - coordinates[2][2] );
542 : :
543 [ + - ][ + - ]: 38 : tri_normal = xxi * xet;
544 : :
545 [ + - ]: 38 : double deta = tri_normal.length();
546 [ + - ][ - + ]: 38 : if( deta == 0.0 || detw == 0.0 ) return 0.0;
547 : :
548 : 38 : double size = pow( deta / detw, 2 );
549 : :
550 [ + + ]: 38 : double rel_size = VERDICT_MIN( size, 1.0 / size );
551 : :
552 [ + - ][ + - ]: 38 : if( rel_size > 0 ) return (double)VERDICT_MIN( rel_size, VERDICT_DBL_MAX );
553 [ # # ]: 38 : return (double)VERDICT_MAX( rel_size, -VERDICT_DBL_MAX );
554 : : }
555 : :
556 : : /*!
557 : : The shape and size of a tri
558 : :
559 : : Product of the Shape and Relative Size
560 : : */
561 : 19 : C_FUNC_DEF double v_tri_shape_and_size( int num_nodes, double coordinates[][3] )
562 : : {
563 : : double size, shape;
564 : :
565 : 19 : size = v_tri_relative_size_squared( num_nodes, coordinates );
566 : 19 : shape = v_tri_shape( num_nodes, coordinates );
567 : :
568 : 19 : double shape_and_size = size * shape;
569 : :
570 [ + - ][ + - ]: 19 : if( shape_and_size > 0 ) return (double)VERDICT_MIN( shape_and_size, VERDICT_DBL_MAX );
571 [ # # ]: 0 : return (double)VERDICT_MAX( shape_and_size, -VERDICT_DBL_MAX );
572 : : }
573 : :
574 : : /*!
575 : : The distortion of a tri
576 : :
577 : : TODO: make a short definition of the distortion and comment below
578 : : */
579 : 37 : C_FUNC_DEF double v_tri_distortion( int num_nodes, double coordinates[][3] )
580 : : {
581 : :
582 : : double distortion;
583 : 37 : int total_number_of_gauss_points = 0;
584 [ + - ][ + - ]: 37 : VerdictVector aa, bb, cc, normal_at_point, xin;
[ + - ][ + - ]
[ + - ]
585 : 37 : double element_area = 0.;
586 : :
587 : 74 : aa.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
588 [ + - ]: 37 : coordinates[1][2] - coordinates[0][2] );
589 : :
590 : 74 : bb.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
591 [ + - ]: 37 : coordinates[2][2] - coordinates[0][2] );
592 : :
593 [ + - ]: 37 : VerdictVector tri_normal = aa * bb;
594 : :
595 : 37 : int number_of_gauss_points = 0;
596 [ + - ]: 37 : if( num_nodes == 3 )
597 : : {
598 : 37 : distortion = 1.0;
599 : 37 : return (double)distortion;
600 : : }
601 : :
602 [ # # ]: 0 : else if( num_nodes == 6 )
603 : : {
604 : 0 : total_number_of_gauss_points = 6;
605 : 0 : number_of_gauss_points = 6;
606 : : }
607 : :
608 : 0 : distortion = VERDICT_DBL_MAX;
609 : : double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
610 : : double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
611 : : double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
612 : : double weight[maxTotalNumberGaussPoints];
613 : :
614 : : // create an object of GaussIntegration
615 : 0 : int number_dims = 2;
616 : 0 : int is_tri = 1;
617 [ # # ]: 0 : GaussIntegration::initialize( number_of_gauss_points, num_nodes, number_dims, is_tri );
618 [ # # ]: 0 : GaussIntegration::calculate_shape_function_2d_tri();
619 [ # # ]: 0 : GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], weight );
620 : :
621 : : // calculate element area
622 : : int ife, ja;
623 [ # # ]: 0 : for( ife = 0; ife < total_number_of_gauss_points; ife++ )
624 : : {
625 [ # # ]: 0 : aa.set( 0.0, 0.0, 0.0 );
626 [ # # ]: 0 : bb.set( 0.0, 0.0, 0.0 );
627 : :
628 [ # # ]: 0 : for( ja = 0; ja < num_nodes; ja++ )
629 : : {
630 [ # # ]: 0 : xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
631 [ # # ][ # # ]: 0 : aa += dndy1[ife][ja] * xin;
632 [ # # ][ # # ]: 0 : bb += dndy2[ife][ja] * xin;
633 : : }
634 [ # # ][ # # ]: 0 : normal_at_point = aa * bb;
635 [ # # ]: 0 : double jacobian = normal_at_point.length();
636 : 0 : element_area += weight[ife] * jacobian;
637 : : }
638 : :
639 : 0 : element_area *= 0.8660254;
640 : : double dndy1_at_node[maxNumberNodes][maxNumberNodes];
641 : : double dndy2_at_node[maxNumberNodes][maxNumberNodes];
642 : :
643 [ # # ]: 0 : GaussIntegration::calculate_derivative_at_nodes_2d_tri( dndy1_at_node, dndy2_at_node );
644 : :
645 [ # # ][ # # ]: 0 : VerdictVector normal_at_nodes[7];
646 : :
647 : : // evaluate normal at nodes and distortion values at nodes
648 : 0 : int jai = 0;
649 [ # # ]: 0 : for( ja = 0; ja < num_nodes; ja++ )
650 : : {
651 [ # # ]: 0 : aa.set( 0.0, 0.0, 0.0 );
652 [ # # ]: 0 : bb.set( 0.0, 0.0, 0.0 );
653 [ # # ]: 0 : for( jai = 0; jai < num_nodes; jai++ )
654 : : {
655 [ # # ]: 0 : xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
656 [ # # ][ # # ]: 0 : aa += dndy1_at_node[ja][jai] * xin;
657 [ # # ][ # # ]: 0 : bb += dndy2_at_node[ja][jai] * xin;
658 : : }
659 [ # # ][ # # ]: 0 : normal_at_nodes[ja] = aa * bb;
660 [ # # ]: 0 : normal_at_nodes[ja].normalize();
661 : : }
662 : :
663 : : // determine if element is flat
664 : 0 : bool flat_element = true;
665 : : double dot_product;
666 : :
667 [ # # ]: 0 : for( ja = 0; ja < num_nodes; ja++ )
668 : : {
669 [ # # ]: 0 : dot_product = normal_at_nodes[0] % normal_at_nodes[ja];
670 [ # # ]: 0 : if( fabs( dot_product ) < 0.99 )
671 : : {
672 : 0 : flat_element = false;
673 : 0 : break;
674 : : }
675 : : }
676 : :
677 : : // take into consideration of the thickness of the element
678 : : double thickness, thickness_gauss;
679 : : double distrt;
680 : : // get_tri_thickness(tri, element_area, thickness );
681 : 0 : thickness = 0.001 * sqrt( element_area );
682 : :
683 : : // set thickness gauss point location
684 : 0 : double zl = 0.5773502691896;
685 [ # # ]: 0 : if( flat_element ) zl = 0.0;
686 : :
687 [ # # ]: 0 : int no_gauss_pts_z = ( flat_element ) ? 1 : 2;
688 : : double thickness_z;
689 : :
690 : : // loop on integration points
691 : : int igz;
692 [ # # ]: 0 : for( ife = 0; ife < total_number_of_gauss_points; ife++ )
693 : : {
694 : : // loop on the thickness direction gauss points
695 [ # # ]: 0 : for( igz = 0; igz < no_gauss_pts_z; igz++ )
696 : : {
697 : 0 : zl = -zl;
698 : 0 : thickness_z = zl * thickness / 2.0;
699 : :
700 [ # # ]: 0 : aa.set( 0.0, 0.0, 0.0 );
701 [ # # ]: 0 : bb.set( 0.0, 0.0, 0.0 );
702 [ # # ]: 0 : cc.set( 0.0, 0.0, 0.0 );
703 : :
704 [ # # ]: 0 : for( ja = 0; ja < num_nodes; ja++ )
705 : : {
706 [ # # ]: 0 : xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
707 [ # # ][ # # ]: 0 : xin += thickness_z * normal_at_nodes[ja];
708 [ # # ][ # # ]: 0 : aa += dndy1[ife][ja] * xin;
709 [ # # ][ # # ]: 0 : bb += dndy2[ife][ja] * xin;
710 : 0 : thickness_gauss = shape_function[ife][ja] * thickness / 2.0;
711 [ # # ][ # # ]: 0 : cc += thickness_gauss * normal_at_nodes[ja];
712 : : }
713 : :
714 [ # # ][ # # ]: 0 : normal_at_point = aa * bb;
715 [ # # ]: 0 : distrt = cc % normal_at_point;
716 [ # # ]: 0 : if( distrt < distortion ) distortion = distrt;
717 : : }
718 : : }
719 : :
720 : : // loop through nodal points
721 [ # # ]: 0 : for( ja = 0; ja < num_nodes; ja++ )
722 : : {
723 [ # # ]: 0 : for( igz = 0; igz < no_gauss_pts_z; igz++ )
724 : : {
725 : 0 : zl = -zl;
726 : 0 : thickness_z = zl * thickness / 2.0;
727 : :
728 [ # # ]: 0 : aa.set( 0.0, 0.0, 0.0 );
729 [ # # ]: 0 : bb.set( 0.0, 0.0, 0.0 );
730 [ # # ]: 0 : cc.set( 0.0, 0.0, 0.0 );
731 : :
732 [ # # ]: 0 : for( jai = 0; jai < num_nodes; jai++ )
733 : : {
734 [ # # ]: 0 : xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
735 [ # # ][ # # ]: 0 : xin += thickness_z * normal_at_nodes[ja];
736 [ # # ][ # # ]: 0 : aa += dndy1_at_node[ja][jai] * xin;
737 [ # # ][ # # ]: 0 : bb += dndy2_at_node[ja][jai] * xin;
738 [ # # ]: 0 : if( jai == ja )
739 : 0 : thickness_gauss = thickness / 2.0;
740 : : else
741 : 0 : thickness_gauss = 0.;
742 [ # # ][ # # ]: 0 : cc += thickness_gauss * normal_at_nodes[jai];
743 : : }
744 : : }
745 : :
746 [ # # ][ # # ]: 0 : normal_at_point = aa * bb;
747 [ # # ][ # # ]: 0 : double sign_jacobian = ( tri_normal % normal_at_point ) > 0 ? 1. : -1.;
748 [ # # ]: 0 : distrt = sign_jacobian * ( cc % normal_at_point );
749 : :
750 [ # # ]: 0 : if( distrt < distortion ) distortion = distrt;
751 : : }
752 [ # # ]: 0 : if( element_area * thickness != 0 )
753 : 0 : distortion *= 1. / ( element_area * thickness );
754 : : else
755 : 0 : distortion *= 1.;
756 : :
757 [ # # ][ # # ]: 0 : if( distortion > 0 ) return (double)VERDICT_MIN( distortion, VERDICT_DBL_MAX );
758 [ # # ]: 37 : return (double)VERDICT_MAX( distortion, -VERDICT_DBL_MAX );
759 : : }
760 : :
761 : : /*!
762 : : tri_quality for calculating multiple tri metrics at once
763 : :
764 : : using this method is generally faster than using the individual
765 : : method multiple times.
766 : :
767 : : */
768 : 18 : C_FUNC_DEF void v_tri_quality( int num_nodes, double coordinates[][3], unsigned int metrics_request_flag,
769 : : TriMetricVals* metric_vals )
770 : : {
771 : :
772 : 18 : memset( metric_vals, 0, sizeof( TriMetricVals ) );
773 : :
774 : : // for starts, lets set up some basic and common information
775 : :
776 : : /* node numbers and side numbers used below
777 : :
778 : : 2
779 : : ++
780 : : / \
781 : : 2 / \ 1
782 : : / \
783 : : / \
784 : : 0 ---------+ 1
785 : : 0
786 : : */
787 : :
788 : : // vectors for each side
789 [ + - ][ + + ]: 72 : VerdictVector sides[3];
790 : 36 : sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
791 [ + - ]: 18 : coordinates[1][2] - coordinates[0][2] );
792 : 36 : sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
793 [ + - ]: 18 : coordinates[2][2] - coordinates[1][2] );
794 : 36 : sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
795 [ + - ]: 18 : coordinates[2][2] - coordinates[0][2] );
796 [ + - ]: 18 : VerdictVector tri_normal = sides[0] * sides[2];
797 : : // if we have access to normal information, check to see if the
798 : : // element is inverted. If we don't have the normal information
799 : : // that we need for this, assume the element is not inverted.
800 : : // This flag will be used for condition number, jacobian, shape,
801 : : // and size and shape.
802 : 18 : bool is_inverted = false;
803 [ - + ]: 18 : if( compute_normal )
804 : : {
805 : : // center of tri
806 : : double point[3], surf_normal[3];
807 : 0 : point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
808 : 0 : point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
809 : 0 : point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
810 : : // dot product
811 [ # # ]: 0 : compute_normal( point, surf_normal );
812 [ # # ][ # # ]: 0 : if( ( tri_normal.x() * surf_normal[0] + tri_normal.y() * surf_normal[1] + tri_normal.z() * surf_normal[2] ) <
[ # # ][ # # ]
813 : : 0 )
814 : 0 : is_inverted = true;
815 : : }
816 : :
817 : : // lengths squared of each side
818 : : double sides_lengths_squared[3];
819 [ + - ]: 18 : sides_lengths_squared[0] = sides[0].length_squared();
820 [ + - ]: 18 : sides_lengths_squared[1] = sides[1].length_squared();
821 [ + - ]: 18 : sides_lengths_squared[2] = sides[2].length_squared();
822 : :
823 : : // if we are doing angle calcuations
824 [ + - ]: 18 : if( metrics_request_flag & ( V_TRI_MINIMUM_ANGLE | V_TRI_MAXIMUM_ANGLE ) )
825 : : {
826 : : // which is short and long side
827 : 18 : int short_side = 0, long_side = 0;
828 : :
829 [ + + ]: 18 : if( sides_lengths_squared[1] < sides_lengths_squared[0] ) short_side = 1;
830 [ - + ]: 18 : if( sides_lengths_squared[2] < sides_lengths_squared[short_side] ) short_side = 2;
831 : :
832 [ + + ]: 18 : if( sides_lengths_squared[1] > sides_lengths_squared[0] ) long_side = 1;
833 [ + + ]: 18 : if( sides_lengths_squared[2] > sides_lengths_squared[long_side] ) long_side = 2;
834 : :
835 : : // calculate the minimum angle of the tri
836 [ + - ]: 18 : if( metrics_request_flag & V_TRI_MINIMUM_ANGLE )
837 : : {
838 [ + - ][ + - ]: 18 : if( sides_lengths_squared[0] == 0.0 || sides_lengths_squared[1] == 0.0 || sides_lengths_squared[2] == 0.0 )
[ - + ]
839 : 0 : { metric_vals->minimum_angle = 0.0; }
840 [ + + ]: 18 : else if( short_side == 0 )
841 [ + - ]: 13 : metric_vals->minimum_angle = (double)sides[2].interior_angle( sides[1] );
842 [ + - ]: 5 : else if( short_side == 1 )
843 [ + - ]: 5 : metric_vals->minimum_angle = (double)sides[0].interior_angle( sides[2] );
844 : : else
845 [ # # ][ # # ]: 18 : metric_vals->minimum_angle = (double)sides[0].interior_angle( -sides[1] );
846 : : }
847 : :
848 : : // calculate the maximum angle of the tri
849 [ + - ]: 18 : if( metrics_request_flag & V_TRI_MAXIMUM_ANGLE )
850 : : {
851 [ + - ][ + - ]: 18 : if( sides_lengths_squared[0] == 0.0 || sides_lengths_squared[1] == 0.0 || sides_lengths_squared[2] == 0.0 )
[ - + ]
852 : 0 : { metric_vals->maximum_angle = 0.0; }
853 [ + + ]: 18 : else if( long_side == 0 )
854 [ + - ]: 5 : metric_vals->maximum_angle = (double)sides[2].interior_angle( sides[1] );
855 [ + + ]: 13 : else if( long_side == 1 )
856 [ + - ]: 4 : metric_vals->maximum_angle = (double)sides[0].interior_angle( sides[2] );
857 : : else
858 [ + - ][ + - ]: 18 : metric_vals->maximum_angle = (double)sides[0].interior_angle( -sides[1] );
859 : : }
860 : : }
861 : :
862 : : // calculate the area of the tri
863 : : // the following metrics depend on area
864 [ + - ]: 18 : if( metrics_request_flag &
865 : : ( V_TRI_AREA | V_TRI_SCALED_JACOBIAN | V_TRI_SHAPE | V_TRI_RELATIVE_SIZE_SQUARED | V_TRI_SHAPE_AND_SIZE ) )
866 [ + - ][ + - ]: 18 : { metric_vals->area = (double)( ( sides[0] * sides[2] ).length() * 0.5 ); }
867 : :
868 : : // calculate the aspect ratio
869 [ + - ]: 18 : if( metrics_request_flag & V_TRI_ASPECT_FROBENIUS )
870 : : {
871 : : // sum the lengths squared
872 : 18 : double srms = sides_lengths_squared[0] + sides_lengths_squared[1] + sides_lengths_squared[2];
873 : :
874 : : // calculate once and reuse
875 : : static const double twoTimesRootOf3 = 2 * sqrt( 3.0 );
876 : :
877 [ + - ][ + - ]: 18 : double div = ( twoTimesRootOf3 * ( ( sides[0] * sides[2] ).length() ) );
878 : :
879 [ - + ]: 18 : if( div == 0.0 )
880 : 0 : metric_vals->aspect_frobenius = (double)VERDICT_DBL_MAX;
881 : : else
882 : 18 : metric_vals->aspect_frobenius = (double)( srms / div );
883 : : }
884 : :
885 : : // calculate the scaled jacobian
886 [ + - ]: 18 : if( metrics_request_flag & V_TRI_SCALED_JACOBIAN )
887 : : {
888 : : // calculate once and reuse
889 : : static const double twoOverRootOf3 = 2 / sqrt( 3.0 );
890 : : // use the area from above
891 : :
892 [ + - ]: 18 : double tmp = tri_normal.length() * twoOverRootOf3;
893 : :
894 : : // now scale it by the lengths of the sides
895 : 18 : double min_scaled_jac = VERDICT_DBL_MAX;
896 : : double temp_scaled_jac;
897 [ + + ]: 72 : for( int i = 0; i < 3; i++ )
898 : : {
899 [ + - ][ - + ]: 54 : if( sides_lengths_squared[i % 3] == 0.0 || sides_lengths_squared[( i + 2 ) % 3] == 0.0 )
900 : 0 : temp_scaled_jac = 0.0;
901 : : else
902 : : temp_scaled_jac =
903 : 54 : tmp / sqrt( sides_lengths_squared[i % 3] ) / sqrt( sides_lengths_squared[( i + 2 ) % 3] );
904 [ + + ]: 54 : if( temp_scaled_jac < min_scaled_jac ) min_scaled_jac = temp_scaled_jac;
905 : : }
906 : : // multiply by -1 if the normals are in opposite directions
907 [ - + ]: 18 : if( is_inverted ) { min_scaled_jac *= -1; }
908 : 18 : metric_vals->scaled_jacobian = (double)min_scaled_jac;
909 : : }
910 : :
911 : : // calculate the condition number
912 [ + - ]: 18 : if( metrics_request_flag & V_TRI_CONDITION )
913 : : {
914 : : // calculate once and reuse
915 : : static const double rootOf3 = sqrt( 3.0 );
916 : : // if it is inverted, the condition number is considered to be infinity.
917 [ - + ]: 18 : if( is_inverted ) { metric_vals->condition = VERDICT_DBL_MAX; }
918 : : else
919 : : {
920 [ + - ][ + - ]: 18 : double area2x = ( sides[0] * sides[2] ).length();
921 [ - + ]: 18 : if( area2x == 0.0 )
922 : 0 : metric_vals->condition = (double)( VERDICT_DBL_MAX );
923 : : else
924 [ + - ][ + - ]: 18 : metric_vals->condition = (double)( ( sides[0] % sides[0] + sides[2] % sides[2] - sides[0] % sides[2] ) /
[ + - ]
925 : 18 : ( area2x * rootOf3 ) );
926 : : }
927 : : }
928 : :
929 : : // calculate the shape
930 [ - + ][ # # ]: 18 : if( metrics_request_flag & V_TRI_SHAPE || metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
931 : : {
932 : : // if element is inverted, shape is zero. We don't need to
933 : : // calculate anything.
934 [ - + ]: 18 : if( is_inverted ) { metric_vals->shape = 0.0; }
935 : : else
936 : : { // otherwise, we calculate the shape
937 : : // calculate once and reuse
938 : : static const double rootOf3 = sqrt( 3.0 );
939 : : // reuse area from before
940 : 18 : double area2x = metric_vals->area * 2;
941 : : // dot products
942 [ + - ][ + - ]: 18 : double dots[3] = { sides[0] % sides[0], sides[2] % sides[2], sides[0] % sides[2] };
[ + - ]
943 : :
944 : : // add the dots
945 : 18 : double sum_dots = dots[0] + dots[1] - dots[2];
946 : :
947 : : // then the finale
948 [ - + ]: 18 : if( sum_dots == 0.0 )
949 : 0 : metric_vals->shape = 0.0;
950 : : else
951 : 18 : metric_vals->shape = (double)( rootOf3 * area2x / sum_dots );
952 : : }
953 : : }
954 : :
955 : : // calculate relative size squared
956 [ - + ][ # # ]: 18 : if( metrics_request_flag & V_TRI_RELATIVE_SIZE_SQUARED || metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
957 : : {
958 : : // get weights
959 : : double w11, w21, w12, w22;
960 : 18 : v_tri_get_weight( w11, w21, w12, w22 );
961 : : // get the determinant
962 [ + - ]: 18 : double detw = determinant( w11, w21, w12, w22 );
963 : : // use the area from above and divide with the determinant
964 [ + - ][ - + ]: 18 : if( metric_vals->area == 0.0 || detw == 0.0 )
965 : 0 : metric_vals->relative_size_squared = 0.0;
966 : : else
967 : : {
968 : 18 : double size = metric_vals->area * 2.0 / detw;
969 : : // square the size
970 : 18 : size *= size;
971 : : // value ranges between 0 to 1
972 [ + - ]: 18 : metric_vals->relative_size_squared = (double)VERDICT_MIN( size, 1.0 / size );
973 : : }
974 : : }
975 : :
976 : : // calculate shape and size
977 [ + - ]: 18 : if( metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
978 : 18 : { metric_vals->shape_and_size = metric_vals->relative_size_squared * metric_vals->shape; }
979 : :
980 : : // calculate distortion
981 [ + - ][ + - ]: 18 : if( metrics_request_flag & V_TRI_DISTORTION ) metric_vals->distortion = v_tri_distortion( num_nodes, coordinates );
982 : :
983 : : // take care of any over-flow problems
984 [ + - ]: 18 : if( metric_vals->aspect_frobenius > 0 )
985 [ + - ]: 18 : metric_vals->aspect_frobenius = (double)VERDICT_MIN( metric_vals->aspect_frobenius, VERDICT_DBL_MAX );
986 : : else
987 [ # # ]: 0 : metric_vals->aspect_frobenius = (double)VERDICT_MAX( metric_vals->aspect_frobenius, -VERDICT_DBL_MAX );
988 : :
989 [ + - ]: 18 : if( metric_vals->area > 0 )
990 [ + - ]: 18 : metric_vals->area = (double)VERDICT_MIN( metric_vals->area, VERDICT_DBL_MAX );
991 : : else
992 [ # # ]: 0 : metric_vals->area = (double)VERDICT_MAX( metric_vals->area, -VERDICT_DBL_MAX );
993 : :
994 [ + - ]: 18 : if( metric_vals->minimum_angle > 0 )
995 [ + - ]: 18 : metric_vals->minimum_angle = (double)VERDICT_MIN( metric_vals->minimum_angle, VERDICT_DBL_MAX );
996 : : else
997 [ # # ]: 0 : metric_vals->minimum_angle = (double)VERDICT_MAX( metric_vals->minimum_angle, -VERDICT_DBL_MAX );
998 : :
999 [ + - ]: 18 : if( metric_vals->maximum_angle > 0 )
1000 [ + - ]: 18 : metric_vals->maximum_angle = (double)VERDICT_MIN( metric_vals->maximum_angle, VERDICT_DBL_MAX );
1001 : : else
1002 [ # # ]: 0 : metric_vals->maximum_angle = (double)VERDICT_MAX( metric_vals->maximum_angle, -VERDICT_DBL_MAX );
1003 : :
1004 [ + - ]: 18 : if( metric_vals->condition > 0 )
1005 [ + - ]: 18 : metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
1006 : : else
1007 [ # # ]: 0 : metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
1008 : :
1009 [ + - ]: 18 : if( metric_vals->shape > 0 )
1010 [ + - ]: 18 : metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
1011 : : else
1012 [ # # ]: 0 : metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
1013 : :
1014 [ + - ]: 18 : if( metric_vals->scaled_jacobian > 0 )
1015 [ + - ]: 18 : metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
1016 : : else
1017 [ # # ]: 0 : metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
1018 : :
1019 [ + - ]: 18 : if( metric_vals->relative_size_squared > 0 )
1020 [ + - ]: 18 : metric_vals->relative_size_squared = (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
1021 : : else
1022 : : metric_vals->relative_size_squared =
1023 [ # # ]: 0 : (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
1024 : :
1025 [ + - ]: 18 : if( metric_vals->shape_and_size > 0 )
1026 [ + - ]: 18 : metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
1027 : : else
1028 [ # # ]: 0 : metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
1029 : :
1030 [ + - ]: 18 : if( metric_vals->distortion > 0 )
1031 [ + - ]: 18 : metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
1032 : : else
1033 [ # # ]: 0 : metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
1034 : :
1035 [ + - ][ + - ]: 18 : if( metrics_request_flag & V_TRI_EDGE_RATIO ) { metric_vals->edge_ratio = v_tri_edge_ratio( 3, coordinates ); }
1036 [ + - ]: 18 : if( metrics_request_flag & V_TRI_RADIUS_RATIO )
1037 [ + - ]: 18 : { metric_vals->radius_ratio = v_tri_radius_ratio( 3, coordinates ); }
1038 [ + - ]: 18 : if( metrics_request_flag & V_TRI_ASPECT_FROBENIUS ) // there is no V_TRI_ASPECT_RATIO !
1039 [ + - ]: 18 : { metric_vals->aspect_ratio = v_tri_radius_ratio( 3, coordinates ); }
1040 : 18 : }
|
