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