 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 /*=========================================================================
2
3 Module: $RCSfile: V_QuadMetric.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 * QuadMetric.cpp contains quality calculations for Quads
18 *
19 * This file is part of VERDICT
20 *
21 */
22
23 #define VERDICT_EXPORTS
24
25 #include "moab/verdict.h"
26 #include "VerdictVector.hpp"
27 #include "verdict_defines.hpp"
28 #include "V_GaussIntegration.hpp"
29 #include <memory.h>
30
31 //! the average area of a quad
32 double verdict_quad_size = 0;
33
34 /*!
35 weights based on the average size of a quad
36 */
37 int get_weight( double& m11, double& m21, double& m12, double& m22 )
38 {
39
40 m11 = 1;
41 m21 = 0;
42 m12 = 0;
43 m22 = 1;
44
45 double scale = sqrt( verdict_quad_size / ( m11 * m22 - m21 * m12 ) );
46
47 m11 *= scale;
48 m21 *= scale;
49 m12 *= scale;
50 m22 *= scale;
51
52 return 1;
53 }
54
55 //! return the average area of a quad
56 C_FUNC_DEF void v_set_quad_size( double size )
57 {
58 verdict_quad_size = size;
59 }
60
61 //! returns whether the quad is collapsed or not
62 VerdictBoolean is_collapsed_quad( double coordinates[][3] )
63 {
64 if( coordinates[3][0] == coordinates[2][0] && coordinates[3][1] == coordinates[2][1] &&
65 coordinates[3][2] == coordinates[2][2] )
66 return VERDICT_TRUE;
67
68 else
69 return VERDICT_FALSE;
70 }
71
72 void make_quad_edges( VerdictVector edges[4], double coordinates[][3] )
73 {
74
75 edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
76 coordinates[1][2] - coordinates[0][2] );
77 edges[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
78 coordinates[2][2] - coordinates[1][2] );
79 edges[2].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
80 coordinates[3][2] - coordinates[2][2] );
81 edges[3].set( coordinates[0][0] - coordinates[3][0], coordinates[0][1] - coordinates[3][1],
82 coordinates[0][2] - coordinates[3][2] );
83 }
84
85 void signed_corner_areas( double areas[4], double coordinates[][3] )
86 {
87 VerdictVector edges[4];
88 make_quad_edges( edges, coordinates );
89
90 VerdictVector corner_normals[4];
91 corner_normals[0] = edges[3] * edges[0];
92 corner_normals[1] = edges[0] * edges[1];
93 corner_normals[2] = edges[1] * edges[2];
94 corner_normals[3] = edges[2] * edges[3];
95
96 // principal axes
97 VerdictVector principal_axes[2];
98 principal_axes[0] = edges[0] - edges[2];
99 principal_axes[1] = edges[1] - edges[3];
100
101 // quad center unit normal
102 VerdictVector unit_center_normal;
103 unit_center_normal = principal_axes[0] * principal_axes[1];
104 unit_center_normal.normalize();
105
106 areas[0] = unit_center_normal % corner_normals[0];
107 areas[1] = unit_center_normal % corner_normals[1];
108 areas[2] = unit_center_normal % corner_normals[2];
109 areas[3] = unit_center_normal % corner_normals[3];
110 }
111
112 /*!
113 localize the coordinates of a quad
114
115 localizing puts the centriod of the quad
116 at the orgin and also rotates the quad
117 such that edge (0,1) is aligned with the x axis
118 and the quad normal lines up with the y axis.
119
120 */
121 void localize_quad_coordinates( VerdictVector nodes[4] )
122 {
123 int i;
124 VerdictVector global[4] = { nodes[0], nodes[1], nodes[2], nodes[3] };
125
126 VerdictVector center = ( global[0] + global[1] + global[2] + global[3] ) / 4.0;
127 for( i = 0; i < 4; i++ )
128 global[i] -= center;
129
130 VerdictVector vector_diff;
131 VerdictVector vector_sum;
132 VerdictVector ref_point( 0.0, 0.0, 0.0 );
133 VerdictVector tmp_vector, normal( 0.0, 0.0, 0.0 );
134 VerdictVector vector1, vector2;
135
136 for( i = 0; i < 4; i++ )
137 {
138 vector1 = global[i];
139 vector2 = global[( i + 1 ) % 4];
140 vector_diff = vector2 - vector1;
141 ref_point += vector1;
142 vector_sum = vector1 + vector2;
143
144 tmp_vector.set( vector_diff.y() * vector_sum.z(), vector_diff.z() * vector_sum.x(),
145 vector_diff.x() * vector_sum.y() );
146 normal += tmp_vector;
147 }
148
149 normal.normalize();
150 normal *= -1.0;
151
152 VerdictVector local_x_axis = global[1] - global[0];
153 local_x_axis.normalize();
154
155 VerdictVector local_y_axis = normal * local_x_axis;
156 local_y_axis.normalize();
157
158 for( i = 0; i < 4; i++ )
159 {
160 nodes[i].x( global[i] % local_x_axis );
161 nodes[i].y( global[i] % local_y_axis );
162 nodes[i].z( global[i] % normal );
163 }
164 }
165
166 /*!
167 moves and rotates the quad such that it enables us to
168 use components of ef's
169 */
170 void localize_quad_for_ef( VerdictVector node_pos[4] )
171 {
172
173 VerdictVector centroid( node_pos[0] );
174 centroid += node_pos[1];
175 centroid += node_pos[2];
176 centroid += node_pos[3];
177
178 centroid /= 4.0;
179
180 node_pos[0] -= centroid;
181 node_pos[1] -= centroid;
182 node_pos[2] -= centroid;
183 node_pos[3] -= centroid;
184
185 VerdictVector rotate = node_pos[1] + node_pos[2] - node_pos[3] - node_pos[0];
186 rotate.normalize();
187
188 double cosine = rotate.x();
189 double sine = rotate.y();
190
191 double xnew;
192
193 for( int i = 0; i < 4; i++ )
194 {
195 xnew = cosine * node_pos[i].x() + sine * node_pos[i].y();
196 node_pos[i].y( -sine * node_pos[i].x() + cosine * node_pos[i].y() );
197 node_pos[i].x( xnew );
198 }
199 }
200
201 /*!
202 returns the normal vector of a quad
203 */
204 VerdictVector quad_normal( double coordinates[][3] )
205 {
206 // get normal at node 0
207 VerdictVector edge0, edge1;
208
209 edge0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
210 coordinates[1][2] - coordinates[0][2] );
211
212 edge1.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
213 coordinates[3][2] - coordinates[0][2] );
214
215 VerdictVector norm0 = edge0 * edge1;
216 norm0.normalize();
217
218 // because some faces may have obtuse angles, check another normal at
219 // node 2 for consistent sense
220
221 edge0.set( coordinates[2][0] - coordinates[3][0], coordinates[2][1] - coordinates[3][1],
222 coordinates[2][2] - coordinates[3][2] );
223
224 edge1.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
225 coordinates[2][2] - coordinates[1][2] );
226
227 VerdictVector norm2 = edge0 * edge1;
228 norm2.normalize();
229
230 // if these two agree, we are done, else test a third to decide
231
232 if( ( norm0 % norm2 ) > 0.0 )
233 {
234 norm0 += norm2;
235 norm0 *= 0.5;
236 return norm0;
237 }
238
239 // test normal at node1
240
241 edge0.set( coordinates[1][0] - coordinates[2][0], coordinates[1][1] - coordinates[2][1],
242 coordinates[1][2] - coordinates[2][2] );
243
244 edge1.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
245 coordinates[1][2] - coordinates[0][2] );
246
247 VerdictVector norm1 = edge0 * edge1;
248 norm1.normalize();
249
250 if( ( norm0 % norm1 ) > 0.0 )
251 {
252 norm0 += norm1;
253 norm0 *= 0.5;
254 return norm0;
255 }
256 else
257 {
258 norm2 += norm1;
259 norm2 *= 0.5;
260 return norm2;
261 }
262 }
263
264 /*!
265 the edge ratio of a quad
266
267 NB (P. Pebay 01/19/07):
268 Hmax / Hmin where Hmax and Hmin are respectively the maximum and the
269 minimum edge lengths
270 */
271 C_FUNC_DEF double v_quad_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
272 {
273 VerdictVector edges[4];
274 make_quad_edges( edges, coordinates );
275
276 double a2 = edges[0].length_squared();
277 double b2 = edges[1].length_squared();
278 double c2 = edges[2].length_squared();
279 double d2 = edges[3].length_squared();
280
281 double mab, Mab, mcd, Mcd, m2, M2;
282 if( a2 < b2 )
283 {
284 mab = a2;
285 Mab = b2;
286 }
287 else // b2 <= a2
288 {
289 mab = b2;
290 Mab = a2;
291 }
292 if( c2 < d2 )
293 {
294 mcd = c2;
295 Mcd = d2;
296 }
297 else // d2 <= c2
298 {
299 mcd = d2;
300 Mcd = c2;
301 }
302 m2 = mab < mcd ? mab : mcd;
303 M2 = Mab > Mcd ? Mab : Mcd;
304
305 if( m2 < VERDICT_DBL_MIN )
306 return (double)VERDICT_DBL_MAX;
307 else
308 {
309 double edge_ratio = sqrt( M2 / m2 );
310
311 if( edge_ratio > 0 ) return (double)VERDICT_MIN( edge_ratio, VERDICT_DBL_MAX );
312 return (double)VERDICT_MAX( edge_ratio, -VERDICT_DBL_MAX );
313 }
314 }
315
316 /*!
317 maximum of edge ratio of a quad
318
319 maximum edge length ratio at quad center
320 */
321 C_FUNC_DEF double v_quad_max_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
322 {
323 VerdictVector quad_nodes[4];
324 quad_nodes[0].set( coordinates[0][0], coordinates[0][1], coordinates[0][2] );
325 quad_nodes[1].set( coordinates[1][0], coordinates[1][1], coordinates[1][2] );
326 quad_nodes[2].set( coordinates[2][0], coordinates[2][1], coordinates[2][2] );
327 quad_nodes[3].set( coordinates[3][0], coordinates[3][1], coordinates[3][2] );
328
329 VerdictVector principal_axes[2];
330 principal_axes[0] = quad_nodes[1] + quad_nodes[2] - quad_nodes[0] - quad_nodes[3];
331 principal_axes[1] = quad_nodes[2] + quad_nodes[3] - quad_nodes[0] - quad_nodes[1];
332
333 double len1 = principal_axes[0].length();
334 double len2 = principal_axes[1].length();
335
336 if( len1 < VERDICT_DBL_MIN || len2 < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
337
338 double max_edge_ratio = VERDICT_MAX( len1 / len2, len2 / len1 );
339
340 if( max_edge_ratio > 0 ) return (double)VERDICT_MIN( max_edge_ratio, VERDICT_DBL_MAX );
341 return (double)VERDICT_MAX( max_edge_ratio, -VERDICT_DBL_MAX );
342 }
343
344 /*!
345 the aspect ratio of a quad
346
347 NB (P. Pebay 01/20/07):
348 this is a generalization of the triangle aspect ratio
349 using Heron's formula.
350 */
351 C_FUNC_DEF double v_quad_aspect_ratio( int /*num_nodes*/, double coordinates[][3] )
352 {
353
354 VerdictVector edges[4];
355 make_quad_edges( edges, coordinates );
356
357 double a1 = edges[0].length();
358 double b1 = edges[1].length();
359 double c1 = edges[2].length();
360 double d1 = edges[3].length();
361
362 double ma = a1 > b1 ? a1 : b1;
363 double mb = c1 > d1 ? c1 : d1;
364 double hm = ma > mb ? ma : mb;
365
366 VerdictVector ab = edges[0] * edges[1];
367 VerdictVector cd = edges[2] * edges[3];
368 double denominator = ab.length() + cd.length();
369
370 if( denominator < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
371
372 double aspect_ratio = .5 * hm * ( a1 + b1 + c1 + d1 ) / denominator;
373
374 if( aspect_ratio > 0 ) return (double)VERDICT_MIN( aspect_ratio, VERDICT_DBL_MAX );
375 return (double)VERDICT_MAX( aspect_ratio, -VERDICT_DBL_MAX );
376 }
377
378 /*!
379 the radius ratio of a quad
380
381 NB (P. Pebay 01/19/07):
382 this metric is called "radius ratio" by extension of a concept that does
383 not exist in general with quads -- although a different name should probably
384 be used in the future.
385 */
386 C_FUNC_DEF double v_quad_radius_ratio( int /*num_nodes*/, double coordinates[][3] )
387 {
388 static const double normal_coeff = 1. / ( 2. * sqrt( 2. ) );
389
390 VerdictVector edges[4];
391 make_quad_edges( edges, coordinates );
392
393 double a2 = edges[0].length_squared();
394 double b2 = edges[1].length_squared();
395 double c2 = edges[2].length_squared();
396 double d2 = edges[3].length_squared();
397
398 VerdictVector diag;
399 diag.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
400 coordinates[2][2] - coordinates[0][2] );
401 double m2 = diag.length_squared();
402
403 diag.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
404 coordinates[3][2] - coordinates[1][2] );
405 double n2 = diag.length_squared();
406
407 double t0 = a2 > b2 ? a2 : b2;
408 double t1 = c2 > d2 ? c2 : d2;
409 double t2 = m2 > n2 ? m2 : n2;
410 double h2 = t0 > t1 ? t0 : t1;
411 h2 = h2 > t2 ? h2 : t2;
412
413 VerdictVector ab = edges[0] * edges[1];
414 VerdictVector bc = edges[1] * edges[2];
415 VerdictVector cd = edges[2] * edges[3];
416 VerdictVector da = edges[3] * edges[0];
417
418 t0 = da.length();
419 t1 = ab.length();
420 t2 = bc.length();
421 double t3 = cd.length();
422
423 t0 = t0 < t1 ? t0 : t1;
424 t2 = t2 < t3 ? t2 : t3;
425 t0 = t0 < t2 ? t0 : t2;
426
427 if( t0 < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
428
429 double radius_ratio = normal_coeff * sqrt( ( a2 + b2 + c2 + d2 ) * h2 ) / t0;
430
431 if( radius_ratio > 0 ) return (double)VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
432 return (double)VERDICT_MAX( radius_ratio, -VERDICT_DBL_MAX );
433 }
434
435 /*!
436 the average Frobenius aspect of a quad
437
438 NB (P. Pebay 01/20/07):
439 this metric is calculated by averaging the 4 Frobenius aspects at
440 each corner of the quad, when the reference triangle is right isosceles.
441 */
442 C_FUNC_DEF double v_quad_med_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
443 {
444
445 VerdictVector edges[4];
446 make_quad_edges( edges, coordinates );
447
448 double a2 = edges[0].length_squared();
449 double b2 = edges[1].length_squared();
450 double c2 = edges[2].length_squared();
451 double d2 = edges[3].length_squared();
452
453 VerdictVector ab = edges[0] * edges[1];
454 VerdictVector bc = edges[1] * edges[2];
455 VerdictVector cd = edges[2] * edges[3];
456 VerdictVector da = edges[3] * edges[0];
457
458 double ab1 = ab.length();
459 double bc1 = bc.length();
460 double cd1 = cd.length();
461 double da1 = da.length();
462
463 if( ab1 < VERDICT_DBL_MIN || bc1 < VERDICT_DBL_MIN || cd1 < VERDICT_DBL_MIN || da1 < VERDICT_DBL_MIN )
464 return (double)VERDICT_DBL_MAX;
465
466 double qsum = ( a2 + b2 ) / ab1;
467 qsum += ( b2 + c2 ) / bc1;
468 qsum += ( c2 + d2 ) / cd1;
469 qsum += ( d2 + a2 ) / da1;
470
471 double med_aspect_frobenius = .125 * qsum;
472
473 if( med_aspect_frobenius > 0 ) return (double)VERDICT_MIN( med_aspect_frobenius, VERDICT_DBL_MAX );
474 return (double)VERDICT_MAX( med_aspect_frobenius, -VERDICT_DBL_MAX );
475 }
476
477 /*!
478 the maximum Frobenius aspect of a quad
479
480 NB (P. Pebay 01/20/07):
481 this metric is calculated by taking the maximum of the 4 Frobenius aspects at
482 each corner of the quad, when the reference triangle is right isosceles.
483 */
484 C_FUNC_DEF double v_quad_max_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
485 {
486
487 VerdictVector edges[4];
488 make_quad_edges( edges, coordinates );
489
490 double a2 = edges[0].length_squared();
491 double b2 = edges[1].length_squared();
492 double c2 = edges[2].length_squared();
493 double d2 = edges[3].length_squared();
494
495 VerdictVector ab = edges[0] * edges[1];
496 VerdictVector bc = edges[1] * edges[2];
497 VerdictVector cd = edges[2] * edges[3];
498 VerdictVector da = edges[3] * edges[0];
499
500 double ab1 = ab.length();
501 double bc1 = bc.length();
502 double cd1 = cd.length();
503 double da1 = da.length();
504
505 if( ab1 < VERDICT_DBL_MIN || bc1 < VERDICT_DBL_MIN || cd1 < VERDICT_DBL_MIN || da1 < VERDICT_DBL_MIN )
506 return (double)VERDICT_DBL_MAX;
507
508 double qmax = ( a2 + b2 ) / ab1;
509
510 double qcur = ( b2 + c2 ) / bc1;
511 qmax = qmax > qcur ? qmax : qcur;
512
513 qcur = ( c2 + d2 ) / cd1;
514 qmax = qmax > qcur ? qmax : qcur;
515
516 qcur = ( d2 + a2 ) / da1;
517 qmax = qmax > qcur ? qmax : qcur;
518
519 double max_aspect_frobenius = .5 * qmax;
520
521 if( max_aspect_frobenius > 0 ) return (double)VERDICT_MIN( max_aspect_frobenius, VERDICT_DBL_MAX );
522 return (double)VERDICT_MAX( max_aspect_frobenius, -VERDICT_DBL_MAX );
523 }
524
525 /*!
526 skew of a quad
527
528 maximum ||cos A|| where A is the angle between edges at quad center
529 */
530 C_FUNC_DEF double v_quad_skew( int /*num_nodes*/, double coordinates[][3] )
531 {
532 VerdictVector node_pos[4];
533 for( int i = 0; i < 4; i++ )
534 node_pos[i].set( coordinates[i][0], coordinates[i][1], coordinates[i][2] );
535
536 VerdictVector principle_axes[2];
537 principle_axes[0] = node_pos[1] + node_pos[2] - node_pos[3] - node_pos[0];
538 principle_axes[1] = node_pos[2] + node_pos[3] - node_pos[0] - node_pos[1];
539
540 if( principle_axes[0].normalize() < VERDICT_DBL_MIN ) return 0.0;
541 if( principle_axes[1].normalize() < VERDICT_DBL_MIN ) return 0.0;
542
543 double skew = fabs( principle_axes[0] % principle_axes[1] );
544
545 return (double)VERDICT_MIN( skew, VERDICT_DBL_MAX );
546 }
547
548 /*!
549 taper of a quad
550
551 maximum ratio of lengths derived from opposite edges
552 */
553 C_FUNC_DEF double v_quad_taper( int /*num_nodes*/, double coordinates[][3] )
554 {
555 VerdictVector node_pos[4];
556 for( int i = 0; i < 4; i++ )
557 node_pos[i].set( coordinates[i][0], coordinates[i][1], coordinates[i][2] );
558
559 VerdictVector principle_axes[2];
560 principle_axes[0] = node_pos[1] + node_pos[2] - node_pos[3] - node_pos[0];
561 principle_axes[1] = node_pos[2] + node_pos[3] - node_pos[0] - node_pos[1];
562
563 VerdictVector cross_derivative = node_pos[0] + node_pos[2] - node_pos[1] - node_pos[3];
564
565 double lengths[2];
566 lengths[0] = principle_axes[0].length();
567 lengths[1] = principle_axes[1].length();
568
569 // get min length
570 lengths[0] = VERDICT_MIN( lengths[0], lengths[1] );
571
572 if( lengths[0] < VERDICT_DBL_MIN ) return VERDICT_DBL_MAX;
573
574 double taper = cross_derivative.length() / lengths[0];
575 return (double)VERDICT_MIN( taper, VERDICT_DBL_MAX );
576 }
577
578 /*!
579 warpage of a quad
580
581 deviation of element from planarity
582 */
583 C_FUNC_DEF double v_quad_warpage( int /*num_nodes*/, double coordinates[][3] )
584 {
585
586 VerdictVector edges[4];
587 make_quad_edges( edges, coordinates );
588
589 VerdictVector corner_normals[4];
590 corner_normals[0] = edges[3] * edges[0];
591 corner_normals[1] = edges[0] * edges[1];
592 corner_normals[2] = edges[1] * edges[2];
593 corner_normals[3] = edges[2] * edges[3];
594
595 if( corner_normals[0].normalize() < VERDICT_DBL_MIN || corner_normals[1].normalize() < VERDICT_DBL_MIN ||
596 corner_normals[2].normalize() < VERDICT_DBL_MIN || corner_normals[3].normalize() < VERDICT_DBL_MIN )
597 return (double)VERDICT_DBL_MIN;
598
599 double warpage =
600 pow( VERDICT_MIN( corner_normals[0] % corner_normals[2], corner_normals[1] % corner_normals[3] ), 3 );
601
602 if( warpage > 0 ) return (double)VERDICT_MIN( warpage, VERDICT_DBL_MAX );
603 return (double)VERDICT_MAX( warpage, -VERDICT_DBL_MAX );
604 }
605
606 /*!
607 the area of a quad
608
609 jacobian at quad center
610 */
611 C_FUNC_DEF double v_quad_area( int /*num_nodes*/, double coordinates[][3] )
612 {
613
614 double corner_areas[4];
615 signed_corner_areas( corner_areas, coordinates );
616
617 double area = 0.25 * ( corner_areas[0] + corner_areas[1] + corner_areas[2] + corner_areas[3] );
618
619 if( area > 0 ) return (double)VERDICT_MIN( area, VERDICT_DBL_MAX );
620 return (double)VERDICT_MAX( area, -VERDICT_DBL_MAX );
621 }
622
623 /*!
624 the stretch of a quad
625
626 sqrt(2) * minimum edge length / maximum diagonal length
627 */
628 C_FUNC_DEF double v_quad_stretch( int /*num_nodes*/, double coordinates[][3] )
629 {
630 VerdictVector edges[4], temp;
631 make_quad_edges( edges, coordinates );
632
633 double lengths_squared[4];
634 lengths_squared[0] = edges[0].length_squared();
635 lengths_squared[1] = edges[1].length_squared();
636 lengths_squared[2] = edges[2].length_squared();
637 lengths_squared[3] = edges[3].length_squared();
638
639 temp.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
640 coordinates[2][2] - coordinates[0][2] );
641 double diag02 = temp.length_squared();
642
643 temp.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
644 coordinates[3][2] - coordinates[1][2] );
645 double diag13 = temp.length_squared();
646
647 static const double QUAD_STRETCH_FACTOR = sqrt( 2.0 );
648
649 // 'diag02' is now the max diagonal of the quad
650 diag02 = VERDICT_MAX( diag02, diag13 );
651
652 if( diag02 < VERDICT_DBL_MIN )
653 return (double)VERDICT_DBL_MAX;
654 else
655 {
656 double stretch =
657 (double)( QUAD_STRETCH_FACTOR * sqrt( VERDICT_MIN( VERDICT_MIN( lengths_squared[0], lengths_squared[1] ),
658 VERDICT_MIN( lengths_squared[2], lengths_squared[3] ) ) /
659 diag02 ) );
660
661 return (double)VERDICT_MIN( stretch, VERDICT_DBL_MAX );
662 }
663 }
664
665 /*!
666 the largest angle of a quad
667
668 largest included quad area (degrees)
669 */
670 C_FUNC_DEF double v_quad_maximum_angle( int /*num_nodes*/, double coordinates[][3] )
671 {
672
673 // if this is a collapsed quad, just pass it on to
674 // the tri_largest_angle routine
675 if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_maximum_angle( 3, coordinates );
676
677 double angle;
678 double max_angle = 0.0;
679
680 VerdictVector edges[4];
681 edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
682 coordinates[1][2] - coordinates[0][2] );
683 edges[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
684 coordinates[2][2] - coordinates[1][2] );
685 edges[2].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
686 coordinates[3][2] - coordinates[2][2] );
687 edges[3].set( coordinates[0][0] - coordinates[3][0], coordinates[0][1] - coordinates[3][1],
688 coordinates[0][2] - coordinates[3][2] );
689
690 // go around each node and calculate the angle
691 // at each node
692 double length[4];
693 length[0] = edges[0].length();
694 length[1] = edges[1].length();
695 length[2] = edges[2].length();
696 length[3] = edges[3].length();
697
698 if( length[0] <= VERDICT_DBL_MIN || length[1] <= VERDICT_DBL_MIN || length[2] <= VERDICT_DBL_MIN ||
699 length[3] <= VERDICT_DBL_MIN )
700 return 0.0;
701
702 angle = acos( -( edges[0] % edges[1] ) / ( length[0] * length[1] ) );
703 max_angle = VERDICT_MAX( angle, max_angle );
704
705 angle = acos( -( edges[1] % edges[2] ) / ( length[1] * length[2] ) );
706 max_angle = VERDICT_MAX( angle, max_angle );
707
708 angle = acos( -( edges[2] % edges[3] ) / ( length[2] * length[3] ) );
709 max_angle = VERDICT_MAX( angle, max_angle );
710
711 angle = acos( -( edges[3] % edges[0] ) / ( length[3] * length[0] ) );
712 max_angle = VERDICT_MAX( angle, max_angle );
713
714 max_angle = max_angle * 180.0 / VERDICT_PI;
715
716 // if any signed areas are < 0, then you are getting the wrong angle
717 double areas[4];
718 signed_corner_areas( areas, coordinates );
719
720 if( areas[0] < 0 || areas[1] < 0 || areas[2] < 0 || areas[3] < 0 )
721 {
722 max_angle = 360 - max_angle;
723 }
724
725 if( max_angle > 0 ) return (double)VERDICT_MIN( max_angle, VERDICT_DBL_MAX );
726 return (double)VERDICT_MAX( max_angle, -VERDICT_DBL_MAX );
727 }
728
729 /*!
730 the smallest angle of a quad
731
732 smallest included quad angle (degrees)
733 */
734 C_FUNC_DEF double v_quad_minimum_angle( int /*num_nodes*/, double coordinates[][3] )
735 {
736 // if this quad is a collapsed quad, then just
737 // send it to the tri_smallest_angle routine
738 if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_minimum_angle( 3, coordinates );
739
740 double angle;
741 double min_angle = 360.0;
742
743 VerdictVector edges[4];
744 edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
745 coordinates[1][2] - coordinates[0][2] );
746 edges[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
747 coordinates[2][2] - coordinates[1][2] );
748 edges[2].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
749 coordinates[3][2] - coordinates[2][2] );
750 edges[3].set( coordinates[0][0] - coordinates[3][0], coordinates[0][1] - coordinates[3][1],
751 coordinates[0][2] - coordinates[3][2] );
752
753 // go around each node and calculate the angle
754 // at each node
755 double length[4];
756 length[0] = edges[0].length();
757 length[1] = edges[1].length();
758 length[2] = edges[2].length();
759 length[3] = edges[3].length();
760
761 if( length[0] <= VERDICT_DBL_MIN || length[1] <= VERDICT_DBL_MIN || length[2] <= VERDICT_DBL_MIN ||
762 length[3] <= VERDICT_DBL_MIN )
763 return 360.0;
764
765 angle = acos( -( edges[0] % edges[1] ) / ( length[0] * length[1] ) );
766 min_angle = VERDICT_MIN( angle, min_angle );
767
768 angle = acos( -( edges[1] % edges[2] ) / ( length[1] * length[2] ) );
769 min_angle = VERDICT_MIN( angle, min_angle );
770
771 angle = acos( -( edges[2] % edges[3] ) / ( length[2] * length[3] ) );
772 min_angle = VERDICT_MIN( angle, min_angle );
773
774 angle = acos( -( edges[3] % edges[0] ) / ( length[3] * length[0] ) );
775 min_angle = VERDICT_MIN( angle, min_angle );
776
777 min_angle = min_angle * 180.0 / VERDICT_PI;
778
779 if( min_angle > 0 ) return (double)VERDICT_MIN( min_angle, VERDICT_DBL_MAX );
780 return (double)VERDICT_MAX( min_angle, -VERDICT_DBL_MAX );
781 }
782
783 /*!
784 the oddy of a quad
785
786 general distortion measure based on left Cauchy-Green Tensor
787 */
788 C_FUNC_DEF double v_quad_oddy( int /*num_nodes*/, double coordinates[][3] )
789 {
790
791 double max_oddy = 0.;
792
793 VerdictVector first, second, node_pos[4];
794
795 double g, g11, g12, g22, cur_oddy;
796 int i;
797
798 for( i = 0; i < 4; i++ )
799 node_pos[i].set( coordinates[i][0], coordinates[i][1], coordinates[i][2] );
800
801 for( i = 0; i < 4; i++ )
802 {
803 first = node_pos[i] - node_pos[( i + 1 ) % 4];
804 second = node_pos[i] - node_pos[( i + 3 ) % 4];
805
806 g11 = first % first;
807 g12 = first % second;
808 g22 = second % second;
809 g = g11 * g22 - g12 * g12;
810
811 if( g < VERDICT_DBL_MIN )
812 cur_oddy = VERDICT_DBL_MAX;
813 else
814 cur_oddy = ( ( g11 - g22 ) * ( g11 - g22 ) + 4. * g12 * g12 ) / 2. / g;
815
816 max_oddy = VERDICT_MAX( max_oddy, cur_oddy );
817 }
818
819 if( max_oddy > 0 ) return (double)VERDICT_MIN( max_oddy, VERDICT_DBL_MAX );
820 return (double)VERDICT_MAX( max_oddy, -VERDICT_DBL_MAX );
821 }
822
823 /*!
824 the condition of a quad
825
826 maximum condition number of the Jacobian matrix at 4 corners
827 */
828 C_FUNC_DEF double v_quad_condition( int /*num_nodes*/, double coordinates[][3] )
829 {
830
831 if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_condition( 3, coordinates );
832
833 double areas[4];
834 signed_corner_areas( areas, coordinates );
835
836 double max_condition = 0.;
837
838 VerdictVector xxi, xet;
839
840 double condition;
841
842 for( int i = 0; i < 4; i++ )
843 {
844
845 xxi.set( coordinates[i][0] - coordinates[( i + 1 ) % 4][0], coordinates[i][1] - coordinates[( i + 1 ) % 4][1],
846 coordinates[i][2] - coordinates[( i + 1 ) % 4][2] );
847
848 xet.set( coordinates[i][0] - coordinates[( i + 3 ) % 4][0], coordinates[i][1] - coordinates[( i + 3 ) % 4][1],
849 coordinates[i][2] - coordinates[( i + 3 ) % 4][2] );
850
851 if( areas[i] < VERDICT_DBL_MIN )
852 condition = VERDICT_DBL_MAX;
853 else
854 condition = ( xxi % xxi + xet % xet ) / areas[i];
855
856 max_condition = VERDICT_MAX( max_condition, condition );
857 }
858
859 max_condition /= 2;
860
861 if( max_condition > 0 ) return (double)VERDICT_MIN( max_condition, VERDICT_DBL_MAX );
862 return (double)VERDICT_MAX( max_condition, -VERDICT_DBL_MAX );
863 }
864
865 /*!
866 the jacobian of a quad
867
868 minimum pointwise volume of local map at 4 corners and center of quad
869 */
870 C_FUNC_DEF double v_quad_jacobian( int /*num_nodes*/, double coordinates[][3] )
871 {
872
873 if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return (double)( v_tri_area( 3, coordinates ) * 2.0 );
874
875 double areas[4];
876 signed_corner_areas( areas, coordinates );
877
878 double jacobian = VERDICT_MIN( VERDICT_MIN( areas[0], areas[1] ), VERDICT_MIN( areas[2], areas[3] ) );
879 if( jacobian > 0 ) return (double)VERDICT_MIN( jacobian, VERDICT_DBL_MAX );
880 return (double)VERDICT_MAX( jacobian, -VERDICT_DBL_MAX );
881 }
882
883 /*!
884 scaled jacobian of a quad
885
886 Minimum Jacobian divided by the lengths of the 2 edge vector
887 */
888 C_FUNC_DEF double v_quad_scaled_jacobian( int /*num_nodes*/, double coordinates[][3] )
889 {
890 if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_scaled_jacobian( 3, coordinates );
891
892 double corner_areas[4], min_scaled_jac = VERDICT_DBL_MAX, scaled_jac;
893 signed_corner_areas( corner_areas, coordinates );
894
895 VerdictVector edges[4];
896 make_quad_edges( edges, coordinates );
897
898 double length[4];
899 length[0] = edges[0].length();
900 length[1] = edges[1].length();
901 length[2] = edges[2].length();
902 length[3] = edges[3].length();
903
904 if( length[0] < VERDICT_DBL_MIN || length[1] < VERDICT_DBL_MIN || length[2] < VERDICT_DBL_MIN ||
905 length[3] < VERDICT_DBL_MIN )
906 return 0.0;
907
908 scaled_jac = corner_areas[0] / ( length[0] * length[3] );
909 min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
910
911 scaled_jac = corner_areas[1] / ( length[1] * length[0] );
912 min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
913
914 scaled_jac = corner_areas[2] / ( length[2] * length[1] );
915 min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
916
917 scaled_jac = corner_areas[3] / ( length[3] * length[2] );
918 min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
919
920 if( min_scaled_jac > 0 ) return (double)VERDICT_MIN( min_scaled_jac, VERDICT_DBL_MAX );
921 return (double)VERDICT_MAX( min_scaled_jac, -VERDICT_DBL_MAX );
922 }
923
924 /*!
925 the shear of a quad
926
927 2/Condition number of Jacobian Skew matrix
928 */
929 C_FUNC_DEF double v_quad_shear( int /*num_nodes*/, double coordinates[][3] )
930 {
931 double scaled_jacobian = v_quad_scaled_jacobian( 4, coordinates );
932
933 if( scaled_jacobian <= VERDICT_DBL_MIN )
934 return 0.0;
935 else
936 return (double)VERDICT_MIN( scaled_jacobian, VERDICT_DBL_MAX );
937 }
938
939 /*!
940 the shape of a quad
941
942 2/Condition number of weighted Jacobian matrix
943 */
944 C_FUNC_DEF double v_quad_shape( int /*num_nodes*/, double coordinates[][3] )
945 {
946
947 double corner_areas[4], min_shape = VERDICT_DBL_MAX, shape;
948 signed_corner_areas( corner_areas, coordinates );
949
950 VerdictVector edges[4];
951 make_quad_edges( edges, coordinates );
952
953 double length_squared[4];
954 length_squared[0] = edges[0].length_squared();
955 length_squared[1] = edges[1].length_squared();
956 length_squared[2] = edges[2].length_squared();
957 length_squared[3] = edges[3].length_squared();
958
959 if( length_squared[0] <= VERDICT_DBL_MIN || length_squared[1] <= VERDICT_DBL_MIN ||
960 length_squared[2] <= VERDICT_DBL_MIN || length_squared[3] <= VERDICT_DBL_MIN )
961 return 0.0;
962
963 shape = corner_areas[0] / ( length_squared[0] + length_squared[3] );
964 min_shape = VERDICT_MIN( shape, min_shape );
965
966 shape = corner_areas[1] / ( length_squared[1] + length_squared[0] );
967 min_shape = VERDICT_MIN( shape, min_shape );
968
969 shape = corner_areas[2] / ( length_squared[2] + length_squared[1] );
970 min_shape = VERDICT_MIN( shape, min_shape );
971
972 shape = corner_areas[3] / ( length_squared[3] + length_squared[2] );
973 min_shape = VERDICT_MIN( shape, min_shape );
974
975 min_shape *= 2;
976
977 if( min_shape < VERDICT_DBL_MIN ) min_shape = 0;
978
979 if( min_shape > 0 ) return (double)VERDICT_MIN( min_shape, VERDICT_DBL_MAX );
980 return (double)VERDICT_MAX( min_shape, -VERDICT_DBL_MAX );
981 }
982
983 /*!
984 the relative size of a quad
985
986 Min( J, 1/J ), where J is determinant of weighted Jacobian matrix
987 */
988 C_FUNC_DEF double v_quad_relative_size_squared( int /*num_nodes*/, double coordinates[][3] )
989 {
990
991 double quad_area = v_quad_area( 4, coordinates );
992 double rel_size = 0;
993
994 v_set_quad_size( quad_area );
995 double w11, w21, w12, w22;
996 get_weight( w11, w21, w12, w22 );
997 double avg_area = determinant( w11, w21, w12, w22 );
998
999 if( avg_area > VERDICT_DBL_MIN )
1000 {
1001
1002 w11 = quad_area / avg_area;
1003
1004 if( w11 > VERDICT_DBL_MIN )
1005 {
1006 rel_size = VERDICT_MIN( w11, 1 / w11 );
1007 rel_size *= rel_size;
1008 }
1009 }
1010
1011 if( rel_size > 0 ) return (double)VERDICT_MIN( rel_size, VERDICT_DBL_MAX );
1012 return (double)VERDICT_MAX( rel_size, -VERDICT_DBL_MAX );
1013 }
1014
1015 /*!
1016 the relative shape and size of a quad
1017
1018 Product of Shape and Relative Size
1019 */
1020 C_FUNC_DEF double v_quad_shape_and_size( int num_nodes, double coordinates[][3] )
1021 {
1022 double shape, size;
1023 size = v_quad_relative_size_squared( num_nodes, coordinates );
1024 shape = v_quad_shape( num_nodes, coordinates );
1025
1026 double shape_and_size = shape * size;
1027
1028 if( shape_and_size > 0 ) return (double)VERDICT_MIN( shape_and_size, VERDICT_DBL_MAX );
1029 return (double)VERDICT_MAX( shape_and_size, -VERDICT_DBL_MAX );
1030 }
1031
1032 /*!
1033 the shear and size of a quad
1034
1035 product of shear and relative size
1036 */
1037 C_FUNC_DEF double v_quad_shear_and_size( int num_nodes, double coordinates[][3] )
1038 {
1039 double shear, size;
1040 shear = v_quad_shear( num_nodes, coordinates );
1041 size = v_quad_relative_size_squared( num_nodes, coordinates );
1042
1043 double shear_and_size = shear * size;
1044
1045 if( shear_and_size > 0 ) return (double)VERDICT_MIN( shear_and_size, VERDICT_DBL_MAX );
1046 return (double)VERDICT_MAX( shear_and_size, -VERDICT_DBL_MAX );
1047 }
1048
1049 /*!
1050 the distortion of a quad
1051 */
1052 C_FUNC_DEF double v_quad_distortion( int num_nodes, double coordinates[][3] )
1053 {
1054 // To calculate distortion for linear and 2nd order quads
1055 // distortion = {min(|J|)/actual area}*{parent area}
1056 // parent area = 4 for a quad.
1057 // min |J| is the minimum over nodes and gaussian integration points
1058 // created by Ling Pan, CAT on 4/30/01
1059
1060 double element_area = 0.0, distrt, thickness_gauss;
1061 double cur_jacobian = 0., sign_jacobian, jacobian;
1062 VerdictVector aa, bb, cc, normal_at_point, xin;
1063
1064 // use 2x2 gauss points for linear quads and 3x3 for 2nd order quads
1065 int number_of_gauss_points = 0;
1066 if( num_nodes == 4 )
1067 { // 2x2 quadrature rule
1068 number_of_gauss_points = 2;
1069 }
1070 else if( num_nodes == 8 )
1071 { // 3x3 quadrature rule
1072 number_of_gauss_points = 3;
1073 }
1074
1075 int total_number_of_gauss_points = number_of_gauss_points * number_of_gauss_points;
1076
1077 VerdictVector face_normal = quad_normal( coordinates );
1078
1079 double distortion = VERDICT_DBL_MAX;
1080
1081 VerdictVector first, second;
1082
1083 int i;
1084 // Will work out the case for collapsed quad later
1085 if( is_collapsed_quad( coordinates ) == VERDICT_TRUE )
1086 {
1087 for( i = 0; i < 3; i++ )
1088 {
1089
1090 first.set( coordinates[i][0] - coordinates[( i + 1 ) % 3][0],
1091 coordinates[i][1] - coordinates[( i + 1 ) % 3][1],
1092 coordinates[i][2] - coordinates[( i + 1 ) % 3][2] );
1093
1094 second.set( coordinates[i][0] - coordinates[( i + 2 ) % 3][0],
1095 coordinates[i][1] - coordinates[( i + 2 ) % 3][1],
1096 coordinates[i][2] - coordinates[( i + 2 ) % 3][2] );
1097
1098 sign_jacobian = ( face_normal % ( first * second ) ) > 0 ? 1. : -1.;
1099 cur_jacobian = sign_jacobian * ( first * second ).length();
1100 distortion = VERDICT_MIN( distortion, cur_jacobian );
1101 }
1102 element_area = ( first * second ).length() / 2.0;
1103 distortion /= element_area;
1104 }
1105 else
1106 {
1107 double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
1108 double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
1109 double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
1110 double weight[maxTotalNumberGaussPoints];
1111
1112 // create an object of GaussIntegration
1113 GaussIntegration::initialize( number_of_gauss_points, num_nodes );
1114 GaussIntegration::calculate_shape_function_2d_quad();
1115 GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], weight );
1116
1117 // calculate element area
1118 int ife, ja;
1119 for( ife = 0; ife < total_number_of_gauss_points; ife++ )
1120 {
1121 aa.set( 0.0, 0.0, 0.0 );
1122 bb.set( 0.0, 0.0, 0.0 );
1123
1124 for( ja = 0; ja < num_nodes; ja++ )
1125 {
1126 xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
1127 aa += dndy1[ife][ja] * xin;
1128 bb += dndy2[ife][ja] * xin;
1129 }
1130 normal_at_point = aa * bb;
1131 jacobian = normal_at_point.length();
1132 element_area += weight[ife] * jacobian;
1133 }
1134
1135 double dndy1_at_node[maxNumberNodes][maxNumberNodes];
1136 double dndy2_at_node[maxNumberNodes][maxNumberNodes];
1137
1138 GaussIntegration::calculate_derivative_at_nodes( dndy1_at_node, dndy2_at_node );
1139
1140 VerdictVector normal_at_nodes[9];
1141
1142 // evaluate normal at nodes and distortion values at nodes
1143 int jai;
1144 for( ja = 0; ja < num_nodes; ja++ )
1145 {
1146 aa.set( 0.0, 0.0, 0.0 );
1147 bb.set( 0.0, 0.0, 0.0 );
1148 for( jai = 0; jai < num_nodes; jai++ )
1149 {
1150 xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
1151 aa += dndy1_at_node[ja][jai] * xin;
1152 bb += dndy2_at_node[ja][jai] * xin;
1153 }
1154 normal_at_nodes[ja] = aa * bb;
1155 normal_at_nodes[ja].normalize();
1156 }
1157
1158 // determine if element is flat
1159 bool flat_element = true;
1160 double dot_product;
1161
1162 for( ja = 0; ja < num_nodes; ja++ )
1163 {
1164 dot_product = normal_at_nodes[0] % normal_at_nodes[ja];
1165 if( fabs( dot_product ) < 0.99 )
1166 {
1167 flat_element = false;
1168 break;
1169 }
1170 }
1171
1172 // take into consideration of the thickness of the element
1173 double thickness;
1174 // get_quad_thickness(face, element_area, thickness );
1175 thickness = 0.001 * sqrt( element_area );
1176
1177 // set thickness gauss point location
1178 double zl = 0.5773502691896;
1179 if( flat_element ) zl = 0.0;
1180
1181 int no_gauss_pts_z = ( flat_element ) ? 1 : 2;
1182 double thickness_z;
1183 int igz;
1184 // loop on Gauss points
1185 for( ife = 0; ife < total_number_of_gauss_points; ife++ )
1186 {
1187 // loop on the thickness direction gauss points
1188 for( igz = 0; igz < no_gauss_pts_z; igz++ )
1189 {
1190 zl = -zl;
1191 thickness_z = zl * thickness / 2.0;
1192
1193 aa.set( 0.0, 0.0, 0.0 );
1194 bb.set( 0.0, 0.0, 0.0 );
1195 cc.set( 0.0, 0.0, 0.0 );
1196
1197 for( ja = 0; ja < num_nodes; ja++ )
1198 {
1199 xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
1200 xin += thickness_z * normal_at_nodes[ja];
1201 aa += dndy1[ife][ja] * xin;
1202 bb += dndy2[ife][ja] * xin;
1203 thickness_gauss = shape_function[ife][ja] * thickness / 2.0;
1204 cc += thickness_gauss * normal_at_nodes[ja];
1205 }
1206
1207 normal_at_point = aa * bb;
1208 // jacobian = normal_at_point.length();
1209 distrt = cc % normal_at_point;
1210 if( distrt < distortion ) distortion = distrt;
1211 }
1212 }
1213
1214 // loop through nodal points
1215 for( ja = 0; ja < num_nodes; ja++ )
1216 {
1217 for( igz = 0; igz < no_gauss_pts_z; igz++ )
1218 {
1219 zl = -zl;
1220 thickness_z = zl * thickness / 2.0;
1221
1222 aa.set( 0.0, 0.0, 0.0 );
1223 bb.set( 0.0, 0.0, 0.0 );
1224 cc.set( 0.0, 0.0, 0.0 );
1225
1226 for( jai = 0; jai < num_nodes; jai++ )
1227 {
1228 xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
1229 xin += thickness_z * normal_at_nodes[ja];
1230 aa += dndy1_at_node[ja][jai] * xin;
1231 bb += dndy2_at_node[ja][jai] * xin;
1232 if( jai == ja )
1233 thickness_gauss = thickness / 2.0;
1234 else
1235 thickness_gauss = 0.;
1236 cc += thickness_gauss * normal_at_nodes[jai];
1237 }
1238 }
1239 normal_at_point = aa * bb;
1240 sign_jacobian = ( face_normal % normal_at_point ) > 0 ? 1. : -1.;
1241 distrt = sign_jacobian * ( cc % normal_at_point );
1242
1243 if( distrt < distortion ) distortion = distrt;
1244 }
1245
1246 if( element_area * thickness != 0 )
1247 distortion *= 8. / ( element_area * thickness );
1248 else
1249 distortion *= 8.;
1250 }
1251
1252 return (double)distortion;
1253 }
1254
1255 /*!
1256 multiple quality measures of a quad
1257 */
1258 C_FUNC_DEF void v_quad_quality( int num_nodes,
1259 double coordinates[][3],
1260 unsigned int metrics_request_flag,
1261 QuadMetricVals* metric_vals )
1262 {
1263
1264 memset( metric_vals, 0, sizeof( QuadMetricVals ) );
1265
1266 // for starts, lets set up some basic and common information
1267
1268 /* node numbers and side numbers used below
1269
1270 2
1271 3 +--------- 2
1272 / +
1273 / |
1274 3 / | 1
1275 / |
1276 + |
1277 0 -------------+ 1
1278 0
1279 */
1280
1281 // vectors for each side
1282 VerdictVector edges[4];
1283 make_quad_edges( edges, coordinates );
1284
1285 double areas[4];
1286 signed_corner_areas( areas, coordinates );
1287
1288 double lengths[4];
1289 lengths[0] = edges[0].length();
1290 lengths[1] = edges[1].length();
1291 lengths[2] = edges[2].length();
1292 lengths[3] = edges[3].length();
1293
1294 VerdictBoolean is_collapsed = is_collapsed_quad( coordinates );
1295
1296 // handle collapsed quads metrics here
1297 if( is_collapsed == VERDICT_TRUE && metrics_request_flag & ( V_QUAD_MINIMUM_ANGLE | V_QUAD_MAXIMUM_ANGLE |
1298 V_QUAD_JACOBIAN | V_QUAD_SCALED_JACOBIAN ) )
1299 {
1300 if( metrics_request_flag & V_QUAD_MINIMUM_ANGLE )
1301 metric_vals->minimum_angle = v_tri_minimum_angle( 3, coordinates );
1302 if( metrics_request_flag & V_QUAD_MAXIMUM_ANGLE )
1303 metric_vals->maximum_angle = v_tri_maximum_angle( 3, coordinates );
1304 if( metrics_request_flag & V_QUAD_JACOBIAN )
1305 metric_vals->jacobian = (double)( v_tri_area( 3, coordinates ) * 2.0 );
1306 if( metrics_request_flag & V_QUAD_SCALED_JACOBIAN )
1307 metric_vals->jacobian = (double)( v_tri_scaled_jacobian( 3, coordinates ) * 2.0 );
1308 }
1309
1310 // calculate both largest and smallest angles
1311 if( metrics_request_flag & ( V_QUAD_MINIMUM_ANGLE | V_QUAD_MAXIMUM_ANGLE ) && is_collapsed == VERDICT_FALSE )
1312 {
1313 // gather the angles
1314 double angles[4];
1315 angles[0] = acos( -( edges[0] % edges[1] ) / ( lengths[0] * lengths[1] ) );
1316 angles[1] = acos( -( edges[1] % edges[2] ) / ( lengths[1] * lengths[2] ) );
1317 angles[2] = acos( -( edges[2] % edges[3] ) / ( lengths[2] * lengths[3] ) );
1318 angles[3] = acos( -( edges[3] % edges[0] ) / ( lengths[3] * lengths[0] ) );
1319
1320 if( lengths[0] <= VERDICT_DBL_MIN || lengths[1] <= VERDICT_DBL_MIN || lengths[2] <= VERDICT_DBL_MIN ||
1321 lengths[3] <= VERDICT_DBL_MIN )
1322 {
1323 metric_vals->minimum_angle = 360.0;
1324 metric_vals->maximum_angle = 0.0;
1325 }
1326 else
1327 {
1328 // if smallest angle, find the smallest angle
1329 if( metrics_request_flag & V_QUAD_MINIMUM_ANGLE )
1330 {
1331 metric_vals->minimum_angle = VERDICT_DBL_MAX;
1332 for( int i = 0; i < 4; i++ )
1333 metric_vals->minimum_angle = VERDICT_MIN( angles[i], metric_vals->minimum_angle );
1334 metric_vals->minimum_angle *= 180.0 / VERDICT_PI;
1335 }
1336 // if largest angle, find the largest angle
1337 if( metrics_request_flag & V_QUAD_MAXIMUM_ANGLE )
1338 {
1339 metric_vals->maximum_angle = 0.0;
1340 for( int i = 0; i < 4; i++ )
1341 metric_vals->maximum_angle = VERDICT_MAX( angles[i], metric_vals->maximum_angle );
1342 metric_vals->maximum_angle *= 180.0 / VERDICT_PI;
1343
1344 if( areas[0] < 0 || areas[1] < 0 || areas[2] < 0 || areas[3] < 0 )
1345 metric_vals->maximum_angle = 360 - metric_vals->maximum_angle;
1346 }
1347 }
1348 }
1349
1350 // handle max_edge_ratio, skew, taper, and area together
1351 if( metrics_request_flag & ( V_QUAD_MAX_EDGE_RATIO | V_QUAD_SKEW | V_QUAD_TAPER ) )
1352 {
1353 // get principle axes
1354 VerdictVector principal_axes[2];
1355 principal_axes[0] = edges[0] - edges[2];
1356 principal_axes[1] = edges[1] - edges[3];
1357
1358 if( metrics_request_flag & ( V_QUAD_MAX_EDGE_RATIO | V_QUAD_SKEW | V_QUAD_TAPER ) )
1359 {
1360 double len1 = principal_axes[0].length();
1361 double len2 = principal_axes[1].length();
1362
1363 // calculate the max_edge_ratio ratio
1364 if( metrics_request_flag & V_QUAD_MAX_EDGE_RATIO )
1365 {
1366 if( len1 < VERDICT_DBL_MIN || len2 < VERDICT_DBL_MIN )
1367 metric_vals->max_edge_ratio = VERDICT_DBL_MAX;
1368 else
1369 metric_vals->max_edge_ratio = VERDICT_MAX( len1 / len2, len2 / len1 );
1370 }
1371
1372 // calculate the taper
1373 if( metrics_request_flag & V_QUAD_TAPER )
1374 {
1375 double min_length = VERDICT_MIN( len1, len2 );
1376
1377 VerdictVector cross_derivative = edges[1] + edges[3];
1378
1379 if( min_length < VERDICT_DBL_MIN )
1380 metric_vals->taper = VERDICT_DBL_MAX;
1381 else
1382 metric_vals->taper = cross_derivative.length() / min_length;
1383 }
1384
1385 // calculate the skew
1386 if( metrics_request_flag & V_QUAD_SKEW )
1387 {
1388 if( principal_axes[0].normalize() < VERDICT_DBL_MIN || principal_axes[1].normalize() < VERDICT_DBL_MIN )
1389 metric_vals->skew = 0.0;
1390 else
1391 metric_vals->skew = fabs( principal_axes[0] % principal_axes[1] );
1392 }
1393 }
1394 }
1395
1396 // calculate the area
1397 if( metrics_request_flag & ( V_QUAD_AREA | V_QUAD_RELATIVE_SIZE_SQUARED ) )
1398 {
1399 metric_vals->area = 0.25 * ( areas[0] + areas[1] + areas[2] + areas[3] );
1400 }
1401
1402 // calculate the relative size
1403 if( metrics_request_flag & ( V_QUAD_RELATIVE_SIZE_SQUARED | V_QUAD_SHAPE_AND_SIZE | V_QUAD_SHEAR_AND_SIZE ) )
1404 {
1405 double quad_area = v_quad_area( 4, coordinates );
1406 v_set_quad_size( quad_area );
1407 double w11, w21, w12, w22;
1408 get_weight( w11, w21, w12, w22 );
1409 double avg_area = determinant( w11, w21, w12, w22 );
1410
1411 if( avg_area < VERDICT_DBL_MIN )
1412 metric_vals->relative_size_squared = 0.0;
1413 else
1414 metric_vals->relative_size_squared =
1415 pow( VERDICT_MIN( metric_vals->area / avg_area, avg_area / metric_vals->area ), 2 );
1416 }
1417
1418 // calculate the jacobian
1419 if( metrics_request_flag & V_QUAD_JACOBIAN )
1420 {
1421 metric_vals->jacobian = VERDICT_MIN( VERDICT_MIN( areas[0], areas[1] ), VERDICT_MIN( areas[2], areas[3] ) );
1422 }
1423
1424 if( metrics_request_flag & ( V_QUAD_SCALED_JACOBIAN | V_QUAD_SHEAR | V_QUAD_SHEAR_AND_SIZE ) )
1425 {
1426 double scaled_jac, min_scaled_jac = VERDICT_DBL_MAX;
1427
1428 if( lengths[0] < VERDICT_DBL_MIN || lengths[1] < VERDICT_DBL_MIN || lengths[2] < VERDICT_DBL_MIN ||
1429 lengths[3] < VERDICT_DBL_MIN )
1430 {
1431 metric_vals->scaled_jacobian = 0.0;
1432 metric_vals->shear = 0.0;
1433 }
1434 else
1435 {
1436 scaled_jac = areas[0] / ( lengths[0] * lengths[3] );
1437 min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
1438
1439 scaled_jac = areas[1] / ( lengths[1] * lengths[0] );
1440 min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
1441
1442 scaled_jac = areas[2] / ( lengths[2] * lengths[1] );
1443 min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
1444
1445 scaled_jac = areas[3] / ( lengths[3] * lengths[2] );
1446 min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
1447
1448 metric_vals->scaled_jacobian = min_scaled_jac;
1449
1450 // what the heck...set shear as well
1451 if( min_scaled_jac <= VERDICT_DBL_MIN )
1452 metric_vals->shear = 0.0;
1453 else
1454 metric_vals->shear = min_scaled_jac;
1455 }
1456 }
1457
1458 if( metrics_request_flag & ( V_QUAD_WARPAGE | V_QUAD_ODDY ) )
1459 {
1460 VerdictVector corner_normals[4];
1461 corner_normals[0] = edges[3] * edges[0];
1462 corner_normals[1] = edges[0] * edges[1];
1463 corner_normals[2] = edges[1] * edges[2];
1464 corner_normals[3] = edges[2] * edges[3];
1465
1466 if( metrics_request_flag & V_QUAD_ODDY )
1467 {
1468 double oddy, max_oddy = 0.0;
1469
1470 double diff, dot_prod;
1471
1472 double length_squared[4];
1473 length_squared[0] = corner_normals[0].length_squared();
1474 length_squared[1] = corner_normals[1].length_squared();
1475 length_squared[2] = corner_normals[2].length_squared();
1476 length_squared[3] = corner_normals[3].length_squared();
1477
1478 if( length_squared[0] < VERDICT_DBL_MIN || length_squared[1] < VERDICT_DBL_MIN ||
1479 length_squared[2] < VERDICT_DBL_MIN || length_squared[3] < VERDICT_DBL_MIN )
1480 metric_vals->oddy = VERDICT_DBL_MAX;
1481 else
1482 {
1483 diff = ( lengths[0] * lengths[0] ) - ( lengths[1] * lengths[1] );
1484 dot_prod = edges[0] % edges[1];
1485 oddy = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[1] );
1486 max_oddy = VERDICT_MAX( oddy, max_oddy );
1487
1488 diff = ( lengths[1] * lengths[1] ) - ( lengths[2] * lengths[2] );
1489 dot_prod = edges[1] % edges[2];
1490 oddy = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[2] );
1491 max_oddy = VERDICT_MAX( oddy, max_oddy );
1492
1493 diff = ( lengths[2] * lengths[2] ) - ( lengths[3] * lengths[3] );
1494 dot_prod = edges[2] % edges[3];
1495 oddy = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[3] );
1496 max_oddy = VERDICT_MAX( oddy, max_oddy );
1497
1498 diff = ( lengths[3] * lengths[3] ) - ( lengths[0] * lengths[0] );
1499 dot_prod = edges[3] % edges[0];
1500 oddy = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[0] );
1501 max_oddy = VERDICT_MAX( oddy, max_oddy );
1502
1503 metric_vals->oddy = max_oddy;
1504 }
1505 }
1506
1507 if( metrics_request_flag & V_QUAD_WARPAGE )
1508 {
1509 if( corner_normals[0].normalize() < VERDICT_DBL_MIN || corner_normals[1].normalize() < VERDICT_DBL_MIN ||
1510 corner_normals[2].normalize() < VERDICT_DBL_MIN || corner_normals[3].normalize() < VERDICT_DBL_MIN )
1511 metric_vals->warpage = VERDICT_DBL_MAX;
1512 else
1513 {
1514 metric_vals->warpage =
1515 pow( VERDICT_MIN( corner_normals[0] % corner_normals[2], corner_normals[1] % corner_normals[3] ),
1516 3 );
1517 }
1518 }
1519 }
1520
1521 if( metrics_request_flag & V_QUAD_STRETCH )
1522 {
1523 VerdictVector temp;
1524
1525 temp.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
1526 coordinates[2][2] - coordinates[0][2] );
1527 double diag02 = temp.length_squared();
1528
1529 temp.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
1530 coordinates[3][2] - coordinates[1][2] );
1531 double diag13 = temp.length_squared();
1532
1533 static const double QUAD_STRETCH_FACTOR = sqrt( 2.0 );
1534
1535 // 'diag02' is now the max diagonal of the quad
1536 diag02 = VERDICT_MAX( diag02, diag13 );
1537
1538 if( diag02 < VERDICT_DBL_MIN )
1539 metric_vals->stretch = VERDICT_DBL_MAX;
1540 else
1541 metric_vals->stretch =
1542 QUAD_STRETCH_FACTOR *
1543 VERDICT_MIN( VERDICT_MIN( lengths[0], lengths[1] ), VERDICT_MIN( lengths[2], lengths[3] ) ) /
1544 sqrt( diag02 );
1545 }
1546
1547 if( metrics_request_flag & ( V_QUAD_CONDITION | V_QUAD_SHAPE | V_QUAD_SHAPE_AND_SIZE ) )
1548 {
1549 double lengths_squared[4];
1550 lengths_squared[0] = edges[0].length_squared();
1551 lengths_squared[1] = edges[1].length_squared();
1552 lengths_squared[2] = edges[2].length_squared();
1553 lengths_squared[3] = edges[3].length_squared();
1554
1555 if( areas[0] < VERDICT_DBL_MIN || areas[1] < VERDICT_DBL_MIN || areas[2] < VERDICT_DBL_MIN ||
1556 areas[3] < VERDICT_DBL_MIN )
1557 {
1558 metric_vals->condition = VERDICT_DBL_MAX;
1559 metric_vals->shape = 0.0;
1560 }
1561 else
1562 {
1563 double max_condition = 0.0, condition;
1564
1565 condition = ( lengths_squared[0] + lengths_squared[3] ) / areas[0];
1566 max_condition = VERDICT_MAX( max_condition, condition );
1567
1568 condition = ( lengths_squared[1] + lengths_squared[0] ) / areas[1];
1569 max_condition = VERDICT_MAX( max_condition, condition );
1570
1571 condition = ( lengths_squared[2] + lengths_squared[1] ) / areas[2];
1572 max_condition = VERDICT_MAX( max_condition, condition );
1573
1574 condition = ( lengths_squared[3] + lengths_squared[2] ) / areas[3];
1575 max_condition = VERDICT_MAX( max_condition, condition );
1576
1577 metric_vals->condition = 0.5 * max_condition;
1578 metric_vals->shape = 2 / max_condition;
1579 }
1580 }
1581
1582 if( metrics_request_flag & V_QUAD_AREA )
1583 {
1584 if( metric_vals->area > 0 ) metric_vals->area = (double)VERDICT_MIN( metric_vals->area, VERDICT_DBL_MAX );
1585 metric_vals->area = (double)VERDICT_MAX( metric_vals->area, -VERDICT_DBL_MAX );
1586 }
1587
1588 if( metrics_request_flag & V_QUAD_MAX_EDGE_RATIO )
1589 {
1590 if( metric_vals->max_edge_ratio > 0 )
1591 metric_vals->max_edge_ratio = (double)VERDICT_MIN( metric_vals->max_edge_ratio, VERDICT_DBL_MAX );
1592 metric_vals->max_edge_ratio = (double)VERDICT_MAX( metric_vals->max_edge_ratio, -VERDICT_DBL_MAX );
1593 }
1594
1595 if( metrics_request_flag & V_QUAD_CONDITION )
1596 {
1597 if( metric_vals->condition > 0 )
1598 metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
1599 metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
1600 }
1601
1602 // calculate distortion
1603 if( metrics_request_flag & V_QUAD_DISTORTION )
1604 {
1605 metric_vals->distortion = v_quad_distortion( num_nodes, coordinates );
1606
1607 if( metric_vals->distortion > 0 )
1608 metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
1609 metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
1610 }
1611
1612 if( metrics_request_flag & V_QUAD_JACOBIAN )
1613 {
1614 if( metric_vals->jacobian > 0 )
1615 metric_vals->jacobian = (double)VERDICT_MIN( metric_vals->jacobian, VERDICT_DBL_MAX );
1616 metric_vals->jacobian = (double)VERDICT_MAX( metric_vals->jacobian, -VERDICT_DBL_MAX );
1617 }
1618
1619 if( metrics_request_flag & V_QUAD_MAXIMUM_ANGLE )
1620 {
1621 if( metric_vals->maximum_angle > 0 )
1622 metric_vals->maximum_angle = (double)VERDICT_MIN( metric_vals->maximum_angle, VERDICT_DBL_MAX );
1623 metric_vals->maximum_angle = (double)VERDICT_MAX( metric_vals->maximum_angle, -VERDICT_DBL_MAX );
1624 }
1625
1626 if( metrics_request_flag & V_QUAD_MINIMUM_ANGLE )
1627 {
1628 if( metric_vals->minimum_angle > 0 )
1629 metric_vals->minimum_angle = (double)VERDICT_MIN( metric_vals->minimum_angle, VERDICT_DBL_MAX );
1630 metric_vals->minimum_angle = (double)VERDICT_MAX( metric_vals->minimum_angle, -VERDICT_DBL_MAX );
1631 }
1632
1633 if( metrics_request_flag & V_QUAD_ODDY )
1634 {
1635 if( metric_vals->oddy > 0 ) metric_vals->oddy = (double)VERDICT_MIN( metric_vals->oddy, VERDICT_DBL_MAX );
1636 metric_vals->oddy = (double)VERDICT_MAX( metric_vals->oddy, -VERDICT_DBL_MAX );
1637 }
1638
1639 if( metrics_request_flag & V_QUAD_RELATIVE_SIZE_SQUARED )
1640 {
1641 if( metric_vals->relative_size_squared > 0 )
1642 metric_vals->relative_size_squared =
1643 (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
1644 metric_vals->relative_size_squared =
1645 (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
1646 }
1647
1648 if( metrics_request_flag & V_QUAD_SCALED_JACOBIAN )
1649 {
1650 if( metric_vals->scaled_jacobian > 0 )
1651 metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
1652 metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
1653 }
1654
1655 if( metrics_request_flag & V_QUAD_SHEAR )
1656 {
1657 if( metric_vals->shear > 0 ) metric_vals->shear = (double)VERDICT_MIN( metric_vals->shear, VERDICT_DBL_MAX );
1658 metric_vals->shear = (double)VERDICT_MAX( metric_vals->shear, -VERDICT_DBL_MAX );
1659 }
1660
1661 // calculate shear and size
1662 // reuse values from above
1663 if( metrics_request_flag & V_QUAD_SHEAR_AND_SIZE )
1664 {
1665 metric_vals->shear_and_size = metric_vals->shear * metric_vals->relative_size_squared;
1666
1667 if( metric_vals->shear_and_size > 0 )
1668 metric_vals->shear_and_size = (double)VERDICT_MIN( metric_vals->shear_and_size, VERDICT_DBL_MAX );
1669 metric_vals->shear_and_size = (double)VERDICT_MAX( metric_vals->shear_and_size, -VERDICT_DBL_MAX );
1670 }
1671
1672 if( metrics_request_flag & V_QUAD_SHAPE )
1673 {
1674 if( metric_vals->shape > 0 ) metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
1675 metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
1676 }
1677
1678 // calculate shape and size
1679 // reuse values from above
1680 if( metrics_request_flag & V_QUAD_SHAPE_AND_SIZE )
1681 {
1682 metric_vals->shape_and_size = metric_vals->shape * metric_vals->relative_size_squared;
1683
1684 if( metric_vals->shape_and_size > 0 )
1685 metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
1686 metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
1687 }
1688
1689 if( metrics_request_flag & V_QUAD_SKEW )
1690 {
1691 if( metric_vals->skew > 0 ) metric_vals->skew = (double)VERDICT_MIN( metric_vals->skew, VERDICT_DBL_MAX );
1692 metric_vals->skew = (double)VERDICT_MAX( metric_vals->skew, -VERDICT_DBL_MAX );
1693 }
1694
1695 if( metrics_request_flag & V_QUAD_STRETCH )
1696 {
1697 if( metric_vals->stretch > 0 )
1698 metric_vals->stretch = (double)VERDICT_MIN( metric_vals->stretch, VERDICT_DBL_MAX );
1699 metric_vals->stretch = (double)VERDICT_MAX( metric_vals->stretch, -VERDICT_DBL_MAX );
1700 }
1701
1702 if( metrics_request_flag & V_QUAD_TAPER )
1703 {
1704 if( metric_vals->taper > 0 ) metric_vals->taper = (double)VERDICT_MIN( metric_vals->taper, VERDICT_DBL_MAX );
1705 metric_vals->taper = (double)VERDICT_MAX( metric_vals->taper, -VERDICT_DBL_MAX );
1706 }
1707
1708 if( metrics_request_flag & V_QUAD_WARPAGE )
1709 {
1710 if( metric_vals->warpage > 0 )
1711 metric_vals->warpage = (double)VERDICT_MIN( metric_vals->warpage, VERDICT_DBL_MAX );
1712 metric_vals->warpage = (double)VERDICT_MAX( metric_vals->warpage, -VERDICT_DBL_MAX );
1713 }
1714
1715 if( metrics_request_flag & V_QUAD_MED_ASPECT_FROBENIUS )
1716 metric_vals->med_aspect_frobenius = v_quad_med_aspect_frobenius( 4, coordinates );
1717 if( metrics_request_flag & V_QUAD_MAX_ASPECT_FROBENIUS )
1718 metric_vals->max_aspect_frobenius = v_quad_max_aspect_frobenius( 4, coordinates );
1719 if( metrics_request_flag & V_QUAD_EDGE_RATIO ) metric_vals->edge_ratio = v_quad_edge_ratio( 4, coordinates );
1720 if( metrics_request_flag & V_QUAD_ASPECT_RATIO ) metric_vals->aspect_ratio = v_quad_aspect_ratio( 4, coordinates );
1721 if( metrics_request_flag & V_QUAD_RADIUS_RATIO ) metric_vals->radius_ratio = v_quad_radius_ratio( 4, coordinates );
1722 }
1.9.1.