Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
V_QuadMetric.cpp
Go to the documentation of this file.
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 }