Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
V_HexMetric.cpp
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Module: $RCSfile: V_HexMetric.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  * HexMetric.cpp contains quality calculations for hexes
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 "V_GaussIntegration.hpp"
28 #include "verdict_defines.hpp"
29 #include <memory.h>
30 
31 #if defined( __BORLANDC__ )
32 #pragma warn - 8004 /* "assigned a value that is never used" */
33 #endif
34 
35 //! the average volume of a hex
36 double verdict_hex_size = 0;
37 
38 //! weights based on the average size of a hex
39 int v_hex_get_weight( VerdictVector& v1, VerdictVector& v2, VerdictVector& v3 )
40 {
41 
42  if( verdict_hex_size == 0 ) return 0;
43 
44  v1.set( 1, 0, 0 );
45  v2.set( 0, 1, 0 );
46  v3.set( 0, 0, 1 );
47 
48  double scale = pow( verdict_hex_size / ( v1 % ( v2 * v3 ) ), 0.33333333333 );
49  v1 *= scale;
50  v2 *= scale;
51  v3 *= scale;
52 
53  return 1;
54 }
55 
56 //! returns the average volume of a hex
57 C_FUNC_DEF void v_set_hex_size( double size )
58 {
59  verdict_hex_size = size;
60 }
61 
62 #define make_hex_nodes( coord, pos ) \
63  for( int mhcii = 0; mhcii < 8; mhcii++ ) \
64  { \
65  ( pos )[mhcii].set( ( coord )[mhcii][0], ( coord )[mhcii][1], ( coord )[mhcii][2] ); \
66  }
67 
68 #define make_edge_length_squares( edges, lengths ) \
69  { \
70  for( int melii = 0; melii < 12; melii++ ) \
71  ( lengths )[melii] = ( edges )[melii].length_squared(); \
72  }
73 
74 //! make VerdictVectors from coordinates
75 void make_hex_edges( double coordinates[][3], VerdictVector edges[12] )
76 {
77  edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
78  coordinates[1][2] - coordinates[0][2] );
79  edges[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
80  coordinates[2][2] - coordinates[1][2] );
81  edges[2].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
82  coordinates[3][2] - coordinates[2][2] );
83  edges[3].set( coordinates[0][0] - coordinates[3][0], coordinates[0][1] - coordinates[3][1],
84  coordinates[0][2] - coordinates[3][2] );
85  edges[4].set( coordinates[5][0] - coordinates[4][0], coordinates[5][1] - coordinates[4][1],
86  coordinates[5][2] - coordinates[4][2] );
87  edges[5].set( coordinates[6][0] - coordinates[5][0], coordinates[6][1] - coordinates[5][1],
88  coordinates[6][2] - coordinates[5][2] );
89  edges[6].set( coordinates[7][0] - coordinates[6][0], coordinates[7][1] - coordinates[6][1],
90  coordinates[7][2] - coordinates[6][2] );
91  edges[7].set( coordinates[4][0] - coordinates[7][0], coordinates[4][1] - coordinates[7][1],
92  coordinates[4][2] - coordinates[7][2] );
93  edges[8].set( coordinates[4][0] - coordinates[0][0], coordinates[4][1] - coordinates[0][1],
94  coordinates[4][2] - coordinates[0][2] );
95  edges[9].set( coordinates[5][0] - coordinates[1][0], coordinates[5][1] - coordinates[1][1],
96  coordinates[5][2] - coordinates[1][2] );
97  edges[10].set( coordinates[6][0] - coordinates[2][0], coordinates[6][1] - coordinates[2][1],
98  coordinates[6][2] - coordinates[2][2] );
99  edges[11].set( coordinates[7][0] - coordinates[3][0], coordinates[7][1] - coordinates[3][1],
100  coordinates[7][2] - coordinates[3][2] );
101 }
102 
103 /*!
104  localizes hex coordinates
105 */
106 void localize_hex_coordinates( double coordinates[][3], VerdictVector position[8] )
107 {
108 
109  int ii;
110  for( ii = 0; ii < 8; ii++ )
111  {
112  position[ii].set( coordinates[ii][0], coordinates[ii][1], coordinates[ii][2] );
113  }
114 
115  // ... Make centroid of element the center of coordinate system
116  VerdictVector point_1256 = position[1];
117  point_1256 += position[2];
118  point_1256 += position[5];
119  point_1256 += position[6];
120 
121  VerdictVector point_0374 = position[0];
122  point_0374 += position[3];
123  point_0374 += position[7];
124  point_0374 += position[4];
125 
126  VerdictVector centroid = point_1256;
127  centroid += point_0374;
128  centroid /= 8.0;
129 
130  int i;
131  for( i = 0; i < 8; i++ )
132  position[i] -= centroid;
133 
134  // ... Rotate element such that center of side 1-2 and 0-3 define X axis
135  double DX = point_1256.x() - point_0374.x();
136  double DY = point_1256.y() - point_0374.y();
137  double DZ = point_1256.z() - point_0374.z();
138 
139  double AMAGX = sqrt( DX * DX + DZ * DZ );
140  double AMAGY = sqrt( DX * DX + DY * DY + DZ * DZ );
141  double FMAGX = AMAGX == 0.0 ? 1.0 : 0.0;
142  double FMAGY = AMAGY == 0.0 ? 1.0 : 0.0;
143 
144  double CZ = DX / ( AMAGX + FMAGX ) + FMAGX;
145  double SZ = DZ / ( AMAGX + FMAGX );
146  double CY = AMAGX / ( AMAGY + FMAGY ) + FMAGY;
147  double SY = DY / ( AMAGY + FMAGY );
148 
149  double temp;
150 
151  for( i = 0; i < 8; i++ )
152  {
153  temp = CY * CZ * position[i].x() + CY * SZ * position[i].z() + SY * position[i].y();
154  position[i].y( -SY * CZ * position[i].x() - SY * SZ * position[i].z() + CY * position[i].y() );
155  position[i].z( -SZ * position[i].x() + CZ * position[i].z() );
156  position[i].x( temp );
157  }
158 
159  // ... Now, rotate about Y
160  VerdictVector delta = -position[0];
161  delta -= position[1];
162  delta += position[2];
163  delta += position[3];
164  delta -= position[4];
165  delta -= position[5];
166  delta += position[6];
167  delta += position[7];
168 
169  DY = delta.y();
170  DZ = delta.z();
171 
172  AMAGY = sqrt( DY * DY + DZ * DZ );
173  FMAGY = AMAGY == 0.0 ? 1.0 : 0.0;
174 
175  double CX = DY / ( AMAGY + FMAGY ) + FMAGY;
176  double SX = DZ / ( AMAGY + FMAGY );
177 
178  for( i = 0; i < 8; i++ )
179  {
180  temp = CX * position[i].y() + SX * position[i].z();
181  position[i].z( -SX * position[i].y() + CX * position[i].z() );
182  position[i].y( temp );
183  }
184 }
185 
186 double safe_ratio3( const double numerator, const double denominator, const double max_ratio )
187 {
188  // this filter is essential for good running time in practice
189  double return_value;
190 
191  const double filter_n = max_ratio * 1.0e-16;
192  const double filter_d = 1.0e-16;
193  if( fabs( numerator ) <= filter_n && fabs( denominator ) >= filter_d )
194  {
195  return_value = numerator / denominator;
196  }
197  else
198  {
199  return_value = fabs( numerator ) / max_ratio >= fabs( denominator )
200  ? ( ( numerator >= 0.0 && denominator >= 0.0 ) || ( numerator < 0.0 && denominator < 0.0 )
201  ? max_ratio
202  : -max_ratio )
203  : numerator / denominator;
204  }
205 
206  return return_value;
207 }
208 
209 double safe_ratio( const double numerator, const double denominator )
210 {
211 
212  double return_value;
213  const double filter_n = VERDICT_DBL_MAX;
214  const double filter_d = VERDICT_DBL_MIN;
215  if( fabs( numerator ) <= filter_n && fabs( denominator ) >= filter_d )
216  {
217  return_value = numerator / denominator;
218  }
219  else
220  {
221  return_value = VERDICT_DBL_MAX;
222  }
223 
224  return return_value;
225 }
226 
227 double condition_comp( const VerdictVector& xxi, const VerdictVector& xet, const VerdictVector& xze )
228 {
229  double det = xxi % ( xet * xze );
230 
231  if( det <= VERDICT_DBL_MIN )
232  {
233  return VERDICT_DBL_MAX;
234  }
235 
236  double term1 = xxi % xxi + xet % xet + xze % xze;
237  double term2 =
238  ( ( xxi * xet ) % ( xxi * xet ) ) + ( ( xet * xze ) % ( xet * xze ) ) + ( ( xze * xxi ) % ( xze * xxi ) );
239 
240  return sqrt( term1 * term2 ) / det;
241 }
242 
243 double oddy_comp( const VerdictVector& xxi, const VerdictVector& xet, const VerdictVector& xze )
244 {
245  static const double third = 1.0 / 3.0;
246 
247  double g11, g12, g13, g22, g23, g33, rt_g;
248 
249  g11 = xxi % xxi;
250  g12 = xxi % xet;
251  g13 = xxi % xze;
252  g22 = xet % xet;
253  g23 = xet % xze;
254  g33 = xze % xze;
255  rt_g = xxi % ( xet * xze );
256 
257  double oddy_metric;
258  if( rt_g > VERDICT_DBL_MIN )
259  {
260  double norm_G_squared = g11 * g11 + 2.0 * g12 * g12 + 2.0 * g13 * g13 + g22 * g22 + 2.0 * g23 * g23 + g33 * g33;
261 
262  double norm_J_squared = g11 + g22 + g33;
263 
264  oddy_metric = ( norm_G_squared - third * norm_J_squared * norm_J_squared ) / pow( rt_g, 4. * third );
265  }
266 
267  else
268  oddy_metric = VERDICT_DBL_MAX;
269 
270  return oddy_metric;
271 }
272 
273 //! calcualates edge lengths of a hex
274 double hex_edge_length( int max_min, double coordinates[][3] )
275 {
276  double temp[3], edge[12];
277  int i;
278 
279  // lengths^2 of edges
280  for( i = 0; i < 3; i++ )
281  {
282  temp[i] = coordinates[1][i] - coordinates[0][i];
283  temp[i] = temp[i] * temp[i];
284  }
285  edge[0] = sqrt( temp[0] + temp[1] + temp[2] );
286 
287  for( i = 0; i < 3; i++ )
288  {
289  temp[i] = coordinates[2][i] - coordinates[1][i];
290  temp[i] = temp[i] * temp[i];
291  }
292  edge[1] = sqrt( temp[0] + temp[1] + temp[2] );
293 
294  for( i = 0; i < 3; i++ )
295  {
296  temp[i] = coordinates[3][i] - coordinates[2][i];
297  temp[i] = temp[i] * temp[i];
298  }
299  edge[2] = sqrt( temp[0] + temp[1] + temp[2] );
300 
301  for( i = 0; i < 3; i++ )
302  {
303  temp[i] = coordinates[0][i] - coordinates[3][i];
304  temp[i] = temp[i] * temp[i];
305  }
306  edge[3] = sqrt( temp[0] + temp[1] + temp[2] );
307 
308  for( i = 0; i < 3; i++ )
309  {
310  temp[i] = coordinates[5][i] - coordinates[4][i];
311  temp[i] = temp[i] * temp[i];
312  }
313  edge[4] = sqrt( temp[0] + temp[1] + temp[2] );
314 
315  for( i = 0; i < 3; i++ )
316  {
317  temp[i] = coordinates[6][i] - coordinates[5][i];
318  temp[i] = temp[i] * temp[i];
319  }
320  edge[5] = sqrt( temp[0] + temp[1] + temp[2] );
321 
322  for( i = 0; i < 3; i++ )
323  {
324  temp[i] = coordinates[7][i] - coordinates[6][i];
325  temp[i] = temp[i] * temp[i];
326  }
327  edge[6] = sqrt( temp[0] + temp[1] + temp[2] );
328 
329  for( i = 0; i < 3; i++ )
330  {
331  temp[i] = coordinates[4][i] - coordinates[7][i];
332  temp[i] = temp[i] * temp[i];
333  }
334  edge[7] = sqrt( temp[0] + temp[1] + temp[2] );
335 
336  for( i = 0; i < 3; i++ )
337  {
338  temp[i] = coordinates[4][i] - coordinates[0][i];
339  temp[i] = temp[i] * temp[i];
340  }
341  edge[8] = sqrt( temp[0] + temp[1] + temp[2] );
342 
343  for( i = 0; i < 3; i++ )
344  {
345  temp[i] = coordinates[5][i] - coordinates[1][i];
346  temp[i] = temp[i] * temp[i];
347  }
348  edge[9] = sqrt( temp[0] + temp[1] + temp[2] );
349 
350  for( i = 0; i < 3; i++ )
351  {
352  temp[i] = coordinates[6][i] - coordinates[2][i];
353  temp[i] = temp[i] * temp[i];
354  }
355  edge[10] = sqrt( temp[0] + temp[1] + temp[2] );
356 
357  for( i = 0; i < 3; i++ )
358  {
359  temp[i] = coordinates[7][i] - coordinates[3][i];
360  temp[i] = temp[i] * temp[i];
361  }
362  edge[11] = sqrt( temp[0] + temp[1] + temp[2] );
363 
364  double _edge = edge[0];
365 
366  if( max_min == 0 )
367  {
368  for( i = 1; i < 12; i++ )
369  _edge = VERDICT_MIN( _edge, edge[i] );
370  return (double)_edge;
371  }
372  else
373  {
374  for( i = 1; i < 12; i++ )
375  _edge = VERDICT_MAX( _edge, edge[i] );
376  return (double)_edge;
377  }
378 }
379 
380 double diag_length( int max_min, double coordinates[][3] )
381 {
382  double temp[3], diag[4];
383  int i;
384 
385  // lengths^2 f diag nals
386  for( i = 0; i < 3; i++ )
387  {
388  temp[i] = coordinates[6][i] - coordinates[0][i];
389  temp[i] = temp[i] * temp[i];
390  }
391  diag[0] = sqrt( temp[0] + temp[1] + temp[2] );
392 
393  for( i = 0; i < 3; i++ )
394  {
395  temp[i] = coordinates[4][i] - coordinates[2][i];
396  temp[i] = temp[i] * temp[i];
397  }
398  diag[1] = sqrt( temp[0] + temp[1] + temp[2] );
399 
400  for( i = 0; i < 3; i++ )
401  {
402  temp[i] = coordinates[7][i] - coordinates[1][i];
403  temp[i] = temp[i] * temp[i];
404  }
405  diag[2] = sqrt( temp[0] + temp[1] + temp[2] );
406 
407  for( i = 0; i < 3; i++ )
408  {
409  temp[i] = coordinates[5][i] - coordinates[3][i];
410  temp[i] = temp[i] * temp[i];
411  }
412  diag[3] = sqrt( temp[0] + temp[1] + temp[2] );
413 
414  double diagonal = diag[0];
415  if( max_min == 0 ) // Return min diagonal
416  {
417  for( i = 1; i < 4; i++ )
418  diagonal = VERDICT_MIN( diagonal, diag[i] );
419  return (double)diagonal;
420  }
421  else // Return max diagonal
422  {
423  for( i = 1; i < 4; i++ )
424  diagonal = VERDICT_MAX( diagonal, diag[i] );
425  return (double)diagonal;
426  }
427 }
428 
429 //! calculates efg values
430 VerdictVector calc_hex_efg( int efg_index, VerdictVector coordinates[8] )
431 {
432 
433  VerdictVector efg;
434 
435  switch( efg_index )
436  {
437 
438  case 1:
439  efg = coordinates[1];
440  efg += coordinates[2];
441  efg += coordinates[5];
442  efg += coordinates[6];
443  efg -= coordinates[0];
444  efg -= coordinates[3];
445  efg -= coordinates[4];
446  efg -= coordinates[7];
447  break;
448 
449  case 2:
450  efg = coordinates[2];
451  efg += coordinates[3];
452  efg += coordinates[6];
453  efg += coordinates[7];
454  efg -= coordinates[0];
455  efg -= coordinates[1];
456  efg -= coordinates[4];
457  efg -= coordinates[5];
458  break;
459 
460  case 3:
461  efg = coordinates[4];
462  efg += coordinates[5];
463  efg += coordinates[6];
464  efg += coordinates[7];
465  efg -= coordinates[0];
466  efg -= coordinates[1];
467  efg -= coordinates[2];
468  efg -= coordinates[3];
469  break;
470 
471  case 12:
472  efg = coordinates[0];
473  efg += coordinates[2];
474  efg += coordinates[4];
475  efg += coordinates[6];
476  efg -= coordinates[1];
477  efg -= coordinates[3];
478  efg -= coordinates[5];
479  efg -= coordinates[7];
480  break;
481 
482  case 13:
483  efg = coordinates[0];
484  efg += coordinates[3];
485  efg += coordinates[5];
486  efg += coordinates[6];
487  efg -= coordinates[1];
488  efg -= coordinates[2];
489  efg -= coordinates[4];
490  efg -= coordinates[7];
491  break;
492 
493  case 23:
494  efg = coordinates[0];
495  efg += coordinates[1];
496  efg += coordinates[6];
497  efg += coordinates[7];
498  efg -= coordinates[2];
499  efg -= coordinates[3];
500  efg -= coordinates[4];
501  efg -= coordinates[5];
502  break;
503 
504  case 123:
505  efg = coordinates[0];
506  efg += coordinates[2];
507  efg += coordinates[5];
508  efg += coordinates[7];
509  efg -= coordinates[1];
510  efg -= coordinates[5];
511  efg -= coordinates[6];
512  efg -= coordinates[2];
513  break;
514 
515  default:
516  efg.set( 0, 0, 0 );
517  }
518 
519  return efg;
520 }
521 
522 /*!
523  the edge ratio of a hex
524 
525  NB (P. Pebay 01/23/07):
526  Hmax / Hmin where Hmax and Hmin are respectively the maximum and the
527  minimum edge lengths
528 */
529 C_FUNC_DEF double v_hex_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
530 {
531 
532  VerdictVector edges[12];
533  make_hex_edges( coordinates, edges );
534 
535  double a2 = edges[0].length_squared();
536  double b2 = edges[1].length_squared();
537  double c2 = edges[2].length_squared();
538  double d2 = edges[3].length_squared();
539  double e2 = edges[4].length_squared();
540  double f2 = edges[5].length_squared();
541  double g2 = edges[6].length_squared();
542  double h2 = edges[7].length_squared();
543  double i2 = edges[8].length_squared();
544  double j2 = edges[9].length_squared();
545  double k2 = edges[10].length_squared();
546  double l2 = edges[11].length_squared();
547 
548  double mab, mcd, mef, Mab, Mcd, Mef;
549  double mgh, mij, mkl, Mgh, Mij, Mkl;
550 
551  if( a2 < b2 )
552  {
553  mab = a2;
554  Mab = b2;
555  }
556  else // b2 <= a2
557  {
558  mab = b2;
559  Mab = a2;
560  }
561  if( c2 < d2 )
562  {
563  mcd = c2;
564  Mcd = d2;
565  }
566  else // d2 <= c2
567  {
568  mcd = d2;
569  Mcd = c2;
570  }
571  if( e2 < f2 )
572  {
573  mef = e2;
574  Mef = f2;
575  }
576  else // f2 <= e2
577  {
578  mef = f2;
579  Mef = e2;
580  }
581  if( g2 < h2 )
582  {
583  mgh = g2;
584  Mgh = h2;
585  }
586  else // h2 <= g2
587  {
588  mgh = h2;
589  Mgh = g2;
590  }
591  if( i2 < j2 )
592  {
593  mij = i2;
594  Mij = j2;
595  }
596  else // j2 <= i2
597  {
598  mij = j2;
599  Mij = i2;
600  }
601  if( k2 < l2 )
602  {
603  mkl = k2;
604  Mkl = l2;
605  }
606  else // l2 <= k2
607  {
608  mkl = l2;
609  Mkl = k2;
610  }
611 
612  double m2;
613  m2 = mab < mcd ? mab : mcd;
614  m2 = m2 < mef ? m2 : mef;
615  m2 = m2 < mgh ? m2 : mgh;
616  m2 = m2 < mij ? m2 : mij;
617  m2 = m2 < mkl ? m2 : mkl;
618 
619  if( m2 < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
620 
621  double M2;
622  M2 = Mab > Mcd ? Mab : Mcd;
623  M2 = M2 > Mef ? M2 : Mef;
624  M2 = M2 > Mgh ? M2 : Mgh;
625  M2 = M2 > Mij ? M2 : Mij;
626  M2 = M2 > Mkl ? M2 : Mkl;
627  m2 = m2 < mef ? m2 : mef;
628 
629  M2 = Mab > Mcd ? Mab : Mcd;
630  M2 = M2 > Mef ? M2 : Mef;
631 
632  double edge_ratio = sqrt( M2 / m2 );
633 
634  if( edge_ratio > 0 ) return (double)VERDICT_MIN( edge_ratio, VERDICT_DBL_MAX );
635  return (double)VERDICT_MAX( edge_ratio, -VERDICT_DBL_MAX );
636 }
637 
638 /*!
639  max edge ratio of a hex
640 
641  Maximum edge length ratio at hex center
642 */
643 C_FUNC_DEF double v_hex_max_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
644 {
645  double aspect;
646  VerdictVector node_pos[8];
647  make_hex_nodes( coordinates, node_pos );
648 
649  double aspect_12, aspect_13, aspect_23;
650 
651  VerdictVector efg1 = calc_hex_efg( 1, node_pos );
652  VerdictVector efg2 = calc_hex_efg( 2, node_pos );
653  VerdictVector efg3 = calc_hex_efg( 3, node_pos );
654 
655  double mag_efg1 = efg1.length();
656  double mag_efg2 = efg2.length();
657  double mag_efg3 = efg3.length();
658 
659  aspect_12 = safe_ratio( VERDICT_MAX( mag_efg1, mag_efg2 ), VERDICT_MIN( mag_efg1, mag_efg2 ) );
660  aspect_13 = safe_ratio( VERDICT_MAX( mag_efg1, mag_efg3 ), VERDICT_MIN( mag_efg1, mag_efg3 ) );
661  aspect_23 = safe_ratio( VERDICT_MAX( mag_efg2, mag_efg3 ), VERDICT_MIN( mag_efg2, mag_efg3 ) );
662 
663  aspect = VERDICT_MAX( aspect_12, VERDICT_MAX( aspect_13, aspect_23 ) );
664 
665  if( aspect > 0 ) return (double)VERDICT_MIN( aspect, VERDICT_DBL_MAX );
666  return (double)VERDICT_MAX( aspect, -VERDICT_DBL_MAX );
667 }
668 
669 /*!
670  skew of a hex
671 
672  Maximum ||cosA|| where A is the angle between edges at hex center.
673 */
674 C_FUNC_DEF double v_hex_skew( int /*num_nodes*/, double coordinates[][3] )
675 {
676  VerdictVector node_pos[8];
677  make_hex_nodes( coordinates, node_pos );
678 
679  double skew_1, skew_2, skew_3;
680 
681  VerdictVector efg1 = calc_hex_efg( 1, node_pos );
682  VerdictVector efg2 = calc_hex_efg( 2, node_pos );
683  VerdictVector efg3 = calc_hex_efg( 3, node_pos );
684 
685  if( efg1.normalize() <= VERDICT_DBL_MIN ) return VERDICT_DBL_MAX;
686  if( efg2.normalize() <= VERDICT_DBL_MIN ) return VERDICT_DBL_MAX;
687  if( efg3.normalize() <= VERDICT_DBL_MIN ) return VERDICT_DBL_MAX;
688 
689  skew_1 = fabs( efg1 % efg2 );
690  skew_2 = fabs( efg1 % efg3 );
691  skew_3 = fabs( efg2 % efg3 );
692 
693  double skew = ( VERDICT_MAX( skew_1, VERDICT_MAX( skew_2, skew_3 ) ) );
694 
695  if( skew > 0 ) return (double)VERDICT_MIN( skew, VERDICT_DBL_MAX );
696  return (double)VERDICT_MAX( skew, -VERDICT_DBL_MAX );
697 }
698 
699 /*!
700  taper of a hex
701 
702  Maximum ratio of lengths derived from opposite edges.
703 */
704 C_FUNC_DEF double v_hex_taper( int /*num_nodes*/, double coordinates[][3] )
705 {
706  VerdictVector node_pos[8];
707  make_hex_nodes( coordinates, node_pos );
708 
709  VerdictVector efg1 = calc_hex_efg( 1, node_pos );
710  VerdictVector efg2 = calc_hex_efg( 2, node_pos );
711  VerdictVector efg3 = calc_hex_efg( 3, node_pos );
712 
713  VerdictVector efg12 = calc_hex_efg( 12, node_pos );
714  VerdictVector efg13 = calc_hex_efg( 13, node_pos );
715  VerdictVector efg23 = calc_hex_efg( 23, node_pos );
716 
717  double taper_1 = fabs( safe_ratio( efg12.length(), VERDICT_MIN( efg1.length(), efg2.length() ) ) );
718  double taper_2 = fabs( safe_ratio( efg13.length(), VERDICT_MIN( efg1.length(), efg3.length() ) ) );
719  double taper_3 = fabs( safe_ratio( efg23.length(), VERDICT_MIN( efg2.length(), efg3.length() ) ) );
720 
721  double taper = (double)VERDICT_MAX( taper_1, VERDICT_MAX( taper_2, taper_3 ) );
722 
723  if( taper > 0 ) return (double)VERDICT_MIN( taper, VERDICT_DBL_MAX );
724  return (double)VERDICT_MAX( taper, -VERDICT_DBL_MAX );
725 }
726 
727 /*!
728  volume of a hex
729 
730  Jacobian at hex center
731 */
732 C_FUNC_DEF double v_hex_volume( int /*num_nodes*/, double coordinates[][3] )
733 {
734  VerdictVector node_pos[8];
735  make_hex_nodes( coordinates, node_pos );
736 
737  VerdictVector efg1 = calc_hex_efg( 1, node_pos );
738  VerdictVector efg2 = calc_hex_efg( 2, node_pos );
739  VerdictVector efg3 = calc_hex_efg( 3, node_pos );
740 
741  double volume;
742  volume = (double)( efg1 % ( efg2 * efg3 ) ) / 64.0;
743 
744  if( volume > 0 ) return (double)VERDICT_MIN( volume, VERDICT_DBL_MAX );
745  return (double)VERDICT_MAX( volume, -VERDICT_DBL_MAX );
746 }
747 
748 /*!
749  stretch of a hex
750 
751  sqrt(3) * minimum edge length / maximum diagonal length
752 */
753 C_FUNC_DEF double v_hex_stretch( int /*num_nodes*/, double coordinates[][3] )
754 {
755  static const double HEX_STRETCH_SCALE_FACTOR = sqrt( 3.0 );
756 
757  double min_edge = hex_edge_length( 0, coordinates );
758  double max_diag = diag_length( 1, coordinates );
759 
760  double stretch = HEX_STRETCH_SCALE_FACTOR * safe_ratio( min_edge, max_diag );
761 
762  if( stretch > 0 ) return (double)VERDICT_MIN( stretch, VERDICT_DBL_MAX );
763  return (double)VERDICT_MAX( stretch, -VERDICT_DBL_MAX );
764 }
765 
766 /*!
767  diagonal ratio of a hex
768 
769  Minimum diagonal length / maximum diagonal length
770 */
771 C_FUNC_DEF double v_hex_diagonal( int /*num_nodes*/, double coordinates[][3] )
772 {
773 
774  double min_diag = diag_length( 0, coordinates );
775  double max_diag = diag_length( 1, coordinates );
776 
777  double diagonal = safe_ratio( min_diag, max_diag );
778 
779  if( diagonal > 0 ) return (double)VERDICT_MIN( diagonal, VERDICT_DBL_MAX );
780  return (double)VERDICT_MAX( diagonal, -VERDICT_DBL_MAX );
781 }
782 
783 #define SQR( x ) ( ( x ) * ( x ) )
784 
785 /*!
786  dimension of a hex
787 
788  Pronto-specific characteristic length for stable time step calculation.
789  Char_length = Volume / 2 grad Volume
790 */
791 C_FUNC_DEF double v_hex_dimension( int /*num_nodes*/, double coordinates[][3] )
792 {
793 
794  double gradop[9][4];
795 
796  double x1 = coordinates[0][0];
797  double x2 = coordinates[1][0];
798  double x3 = coordinates[2][0];
799  double x4 = coordinates[3][0];
800  double x5 = coordinates[4][0];
801  double x6 = coordinates[5][0];
802  double x7 = coordinates[6][0];
803  double x8 = coordinates[7][0];
804 
805  double y1 = coordinates[0][1];
806  double y2 = coordinates[1][1];
807  double y3 = coordinates[2][1];
808  double y4 = coordinates[3][1];
809  double y5 = coordinates[4][1];
810  double y6 = coordinates[5][1];
811  double y7 = coordinates[6][1];
812  double y8 = coordinates[7][1];
813 
814  double z1 = coordinates[0][2];
815  double z2 = coordinates[1][2];
816  double z3 = coordinates[2][2];
817  double z4 = coordinates[3][2];
818  double z5 = coordinates[4][2];
819  double z6 = coordinates[5][2];
820  double z7 = coordinates[6][2];
821  double z8 = coordinates[7][2];
822 
823  double z24 = z2 - z4;
824  double z52 = z5 - z2;
825  double z45 = z4 - z5;
826  gradop[1][1] =
827  ( y2 * ( z6 - z3 - z45 ) + y3 * z24 + y4 * ( z3 - z8 - z52 ) + y5 * ( z8 - z6 - z24 ) + y6 * z52 + y8 * z45 ) /
828  12.0;
829 
830  double z31 = z3 - z1;
831  double z63 = z6 - z3;
832  double z16 = z1 - z6;
833  gradop[2][1] =
834  ( y3 * ( z7 - z4 - z16 ) + y4 * z31 + y1 * ( z4 - z5 - z63 ) + y6 * ( z5 - z7 - z31 ) + y7 * z63 + y5 * z16 ) /
835  12.0;
836 
837  double z42 = z4 - z2;
838  double z74 = z7 - z4;
839  double z27 = z2 - z7;
840  gradop[3][1] =
841  ( y4 * ( z8 - z1 - z27 ) + y1 * z42 + y2 * ( z1 - z6 - z74 ) + y7 * ( z6 - z8 - z42 ) + y8 * z74 + y6 * z27 ) /
842  12.0;
843 
844  double z13 = z1 - z3;
845  double z81 = z8 - z1;
846  double z38 = z3 - z8;
847  gradop[4][1] =
848  ( y1 * ( z5 - z2 - z38 ) + y2 * z13 + y3 * ( z2 - z7 - z81 ) + y8 * ( z7 - z5 - z13 ) + y5 * z81 + y7 * z38 ) /
849  12.0;
850 
851  double z86 = z8 - z6;
852  double z18 = z1 - z8;
853  double z61 = z6 - z1;
854  gradop[5][1] =
855  ( y8 * ( z4 - z7 - z61 ) + y7 * z86 + y6 * ( z7 - z2 - z18 ) + y1 * ( z2 - z4 - z86 ) + y4 * z18 + y2 * z61 ) /
856  12.0;
857 
858  double z57 = z5 - z7;
859  double z25 = z2 - z5;
860  double z72 = z7 - z2;
861  gradop[6][1] =
862  ( y5 * ( z1 - z8 - z72 ) + y8 * z57 + y7 * ( z8 - z3 - z25 ) + y2 * ( z3 - z1 - z57 ) + y1 * z25 + y3 * z72 ) /
863  12.0;
864 
865  double z68 = z6 - z8;
866  double z36 = z3 - z6;
867  double z83 = z8 - z3;
868  gradop[7][1] =
869  ( y6 * ( z2 - z5 - z83 ) + y5 * z68 + y8 * ( z5 - z4 - z36 ) + y3 * ( z4 - z2 - z68 ) + y2 * z36 + y4 * z83 ) /
870  12.0;
871 
872  double z75 = z7 - z5;
873  double z47 = z4 - z7;
874  double z54 = z5 - z4;
875  gradop[8][1] =
876  ( y7 * ( z3 - z6 - z54 ) + y6 * z75 + y5 * ( z6 - z1 - z47 ) + y4 * ( z1 - z3 - z75 ) + y3 * z47 + y1 * z54 ) /
877  12.0;
878 
879  double x24 = x2 - x4;
880  double x52 = x5 - x2;
881  double x45 = x4 - x5;
882  gradop[1][2] =
883  ( z2 * ( x6 - x3 - x45 ) + z3 * x24 + z4 * ( x3 - x8 - x52 ) + z5 * ( x8 - x6 - x24 ) + z6 * x52 + z8 * x45 ) /
884  12.0;
885 
886  double x31 = x3 - x1;
887  double x63 = x6 - x3;
888  double x16 = x1 - x6;
889  gradop[2][2] =
890  ( z3 * ( x7 - x4 - x16 ) + z4 * x31 + z1 * ( x4 - x5 - x63 ) + z6 * ( x5 - x7 - x31 ) + z7 * x63 + z5 * x16 ) /
891  12.0;
892 
893  double x42 = x4 - x2;
894  double x74 = x7 - x4;
895  double x27 = x2 - x7;
896  gradop[3][2] =
897  ( z4 * ( x8 - x1 - x27 ) + z1 * x42 + z2 * ( x1 - x6 - x74 ) + z7 * ( x6 - x8 - x42 ) + z8 * x74 + z6 * x27 ) /
898  12.0;
899 
900  double x13 = x1 - x3;
901  double x81 = x8 - x1;
902  double x38 = x3 - x8;
903  gradop[4][2] =
904  ( z1 * ( x5 - x2 - x38 ) + z2 * x13 + z3 * ( x2 - x7 - x81 ) + z8 * ( x7 - x5 - x13 ) + z5 * x81 + z7 * x38 ) /
905  12.0;
906 
907  double x86 = x8 - x6;
908  double x18 = x1 - x8;
909  double x61 = x6 - x1;
910  gradop[5][2] =
911  ( z8 * ( x4 - x7 - x61 ) + z7 * x86 + z6 * ( x7 - x2 - x18 ) + z1 * ( x2 - x4 - x86 ) + z4 * x18 + z2 * x61 ) /
912  12.0;
913 
914  double x57 = x5 - x7;
915  double x25 = x2 - x5;
916  double x72 = x7 - x2;
917  gradop[6][2] =
918  ( z5 * ( x1 - x8 - x72 ) + z8 * x57 + z7 * ( x8 - x3 - x25 ) + z2 * ( x3 - x1 - x57 ) + z1 * x25 + z3 * x72 ) /
919  12.0;
920 
921  double x68 = x6 - x8;
922  double x36 = x3 - x6;
923  double x83 = x8 - x3;
924  gradop[7][2] =
925  ( z6 * ( x2 - x5 - x83 ) + z5 * x68 + z8 * ( x5 - x4 - x36 ) + z3 * ( x4 - x2 - x68 ) + z2 * x36 + z4 * x83 ) /
926  12.0;
927 
928  double x75 = x7 - x5;
929  double x47 = x4 - x7;
930  double x54 = x5 - x4;
931  gradop[8][2] =
932  ( z7 * ( x3 - x6 - x54 ) + z6 * x75 + z5 * ( x6 - x1 - x47 ) + z4 * ( x1 - x3 - x75 ) + z3 * x47 + z1 * x54 ) /
933  12.0;
934 
935  double y24 = y2 - y4;
936  double y52 = y5 - y2;
937  double y45 = y4 - y5;
938  gradop[1][3] =
939  ( x2 * ( y6 - y3 - y45 ) + x3 * y24 + x4 * ( y3 - y8 - y52 ) + x5 * ( y8 - y6 - y24 ) + x6 * y52 + x8 * y45 ) /
940  12.0;
941 
942  double y31 = y3 - y1;
943  double y63 = y6 - y3;
944  double y16 = y1 - y6;
945  gradop[2][3] =
946  ( x3 * ( y7 - y4 - y16 ) + x4 * y31 + x1 * ( y4 - y5 - y63 ) + x6 * ( y5 - y7 - y31 ) + x7 * y63 + x5 * y16 ) /
947  12.0;
948 
949  double y42 = y4 - y2;
950  double y74 = y7 - y4;
951  double y27 = y2 - y7;
952  gradop[3][3] =
953  ( x4 * ( y8 - y1 - y27 ) + x1 * y42 + x2 * ( y1 - y6 - y74 ) + x7 * ( y6 - y8 - y42 ) + x8 * y74 + x6 * y27 ) /
954  12.0;
955 
956  double y13 = y1 - y3;
957  double y81 = y8 - y1;
958  double y38 = y3 - y8;
959  gradop[4][3] =
960  ( x1 * ( y5 - y2 - y38 ) + x2 * y13 + x3 * ( y2 - y7 - y81 ) + x8 * ( y7 - y5 - y13 ) + x5 * y81 + x7 * y38 ) /
961  12.0;
962 
963  double y86 = y8 - y6;
964  double y18 = y1 - y8;
965  double y61 = y6 - y1;
966  gradop[5][3] =
967  ( x8 * ( y4 - y7 - y61 ) + x7 * y86 + x6 * ( y7 - y2 - y18 ) + x1 * ( y2 - y4 - y86 ) + x4 * y18 + x2 * y61 ) /
968  12.0;
969 
970  double y57 = y5 - y7;
971  double y25 = y2 - y5;
972  double y72 = y7 - y2;
973  gradop[6][3] =
974  ( x5 * ( y1 - y8 - y72 ) + x8 * y57 + x7 * ( y8 - y3 - y25 ) + x2 * ( y3 - y1 - y57 ) + x1 * y25 + x3 * y72 ) /
975  12.0;
976 
977  double y68 = y6 - y8;
978  double y36 = y3 - y6;
979  double y83 = y8 - y3;
980  gradop[7][3] =
981  ( x6 * ( y2 - y5 - y83 ) + x5 * y68 + x8 * ( y5 - y4 - y36 ) + x3 * ( y4 - y2 - y68 ) + x2 * y36 + x4 * y83 ) /
982  12.0;
983 
984  double y75 = y7 - y5;
985  double y47 = y4 - y7;
986  double y54 = y5 - y4;
987  gradop[8][3] =
988  ( x7 * ( y3 - y6 - y54 ) + x6 * y75 + x5 * ( y6 - y1 - y47 ) + x4 * ( y1 - y3 - y75 ) + x3 * y47 + x1 * y54 ) /
989  12.0;
990 
991  // calculate element volume and characteristic element aspect ratio
992  // (used in time step and hourglass control) -
993 
994  double volume = coordinates[0][0] * gradop[1][1] + coordinates[1][0] * gradop[2][1] +
995  coordinates[2][0] * gradop[3][1] + coordinates[3][0] * gradop[4][1] +
996  coordinates[4][0] * gradop[5][1] + coordinates[5][0] * gradop[6][1] +
997  coordinates[6][0] * gradop[7][1] + coordinates[7][0] * gradop[8][1];
998  double aspect =
999  .5 * SQR( volume ) /
1000  ( SQR( gradop[1][1] ) + SQR( gradop[2][1] ) + SQR( gradop[3][1] ) + SQR( gradop[4][1] ) + SQR( gradop[5][1] ) +
1001  SQR( gradop[6][1] ) + SQR( gradop[7][1] ) + SQR( gradop[8][1] ) + SQR( gradop[1][2] ) + SQR( gradop[2][2] ) +
1002  SQR( gradop[3][2] ) + SQR( gradop[4][2] ) + SQR( gradop[5][2] ) + SQR( gradop[6][2] ) + SQR( gradop[7][2] ) +
1003  SQR( gradop[8][2] ) + SQR( gradop[1][3] ) + SQR( gradop[2][3] ) + SQR( gradop[3][3] ) + SQR( gradop[4][3] ) +
1004  SQR( gradop[5][3] ) + SQR( gradop[6][3] ) + SQR( gradop[7][3] ) + SQR( gradop[8][3] ) );
1005 
1006  return (double)sqrt( aspect );
1007 }
1008 
1009 /*!
1010  oddy of a hex
1011 
1012  General distortion measure based on left Cauchy-Green Tensor
1013 */
1014 C_FUNC_DEF double v_hex_oddy( int /*num_nodes*/, double coordinates[][3] )
1015 {
1016 
1017  double oddy = 0.0, current_oddy;
1018  VerdictVector xxi, xet, xze;
1019 
1020  VerdictVector node_pos[8];
1021  make_hex_nodes( coordinates, node_pos );
1022 
1023  xxi = calc_hex_efg( 1, node_pos );
1024  xet = calc_hex_efg( 2, node_pos );
1025  xze = calc_hex_efg( 3, node_pos );
1026 
1027  current_oddy = oddy_comp( xxi, xet, xze );
1028  if( current_oddy > oddy )
1029  {
1030  oddy = current_oddy;
1031  }
1032 
1033  xxi.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
1034  coordinates[1][2] - coordinates[0][2] );
1035 
1036  xet.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
1037  coordinates[3][2] - coordinates[0][2] );
1038 
1039  xze.set( coordinates[4][0] - coordinates[0][0], coordinates[4][1] - coordinates[0][1],
1040  coordinates[4][2] - coordinates[0][2] );
1041 
1042  current_oddy = oddy_comp( xxi, xet, xze );
1043  if( current_oddy > oddy )
1044  {
1045  oddy = current_oddy;
1046  }
1047 
1048  xxi.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
1049  coordinates[2][2] - coordinates[1][2] );
1050 
1051  xet.set( coordinates[0][0] - coordinates[1][0], coordinates[0][1] - coordinates[1][1],
1052  coordinates[0][2] - coordinates[1][2] );
1053 
1054  xze.set( coordinates[5][0] - coordinates[1][0], coordinates[5][1] - coordinates[1][1],
1055  coordinates[5][2] - coordinates[1][2] );
1056 
1057  current_oddy = oddy_comp( xxi, xet, xze );
1058  if( current_oddy > oddy )
1059  {
1060  oddy = current_oddy;
1061  }
1062 
1063  xxi.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
1064  coordinates[3][2] - coordinates[2][2] );
1065 
1066  xet.set( coordinates[1][0] - coordinates[2][0], coordinates[1][1] - coordinates[2][1],
1067  coordinates[1][2] - coordinates[2][2] );
1068 
1069  xze.set( coordinates[6][0] - coordinates[2][0], coordinates[6][1] - coordinates[2][1],
1070  coordinates[6][2] - coordinates[2][2] );
1071 
1072  current_oddy = oddy_comp( xxi, xet, xze );
1073  if( current_oddy > oddy )
1074  {
1075  oddy = current_oddy;
1076  }
1077 
1078  xxi.set( coordinates[0][0] - coordinates[3][0], coordinates[0][1] - coordinates[3][1],
1079  coordinates[0][2] - coordinates[3][2] );
1080 
1081  xet.set( coordinates[2][0] - coordinates[3][0], coordinates[2][1] - coordinates[3][1],
1082  coordinates[2][2] - coordinates[3][2] );
1083 
1084  xze.set( coordinates[7][0] - coordinates[3][0], coordinates[7][1] - coordinates[3][1],
1085  coordinates[7][2] - coordinates[3][2] );
1086 
1087  current_oddy = oddy_comp( xxi, xet, xze );
1088  if( current_oddy > oddy )
1089  {
1090  oddy = current_oddy;
1091  }
1092 
1093  xxi.set( coordinates[7][0] - coordinates[4][0], coordinates[7][1] - coordinates[4][1],
1094  coordinates[7][2] - coordinates[4][2] );
1095 
1096  xet.set( coordinates[5][0] - coordinates[4][0], coordinates[5][1] - coordinates[4][1],
1097  coordinates[5][2] - coordinates[4][2] );
1098 
1099  xze.set( coordinates[0][0] - coordinates[4][0], coordinates[0][1] - coordinates[4][1],
1100  coordinates[0][2] - coordinates[4][2] );
1101 
1102  current_oddy = oddy_comp( xxi, xet, xze );
1103  if( current_oddy > oddy )
1104  {
1105  oddy = current_oddy;
1106  }
1107 
1108  xxi.set( coordinates[4][0] - coordinates[5][0], coordinates[4][1] - coordinates[5][1],
1109  coordinates[4][2] - coordinates[5][2] );
1110 
1111  xet.set( coordinates[6][0] - coordinates[5][0], coordinates[6][1] - coordinates[5][1],
1112  coordinates[6][2] - coordinates[5][2] );
1113 
1114  xze.set( coordinates[1][0] - coordinates[5][0], coordinates[1][1] - coordinates[5][1],
1115  coordinates[1][2] - coordinates[5][2] );
1116 
1117  current_oddy = oddy_comp( xxi, xet, xze );
1118  if( current_oddy > oddy )
1119  {
1120  oddy = current_oddy;
1121  }
1122 
1123  xxi.set( coordinates[5][0] - coordinates[6][0], coordinates[5][1] - coordinates[6][1],
1124  coordinates[5][2] - coordinates[6][2] );
1125 
1126  xet.set( coordinates[7][0] - coordinates[6][0], coordinates[7][1] - coordinates[6][1],
1127  coordinates[7][2] - coordinates[6][2] );
1128 
1129  xze.set( coordinates[2][0] - coordinates[6][0], coordinates[2][1] - coordinates[6][1],
1130  coordinates[2][2] - coordinates[6][2] );
1131 
1132  current_oddy = oddy_comp( xxi, xet, xze );
1133  if( current_oddy > oddy )
1134  {
1135  oddy = current_oddy;
1136  }
1137 
1138  xxi.set( coordinates[6][0] - coordinates[7][0], coordinates[6][1] - coordinates[7][1],
1139  coordinates[6][2] - coordinates[7][2] );
1140 
1141  xet.set( coordinates[4][0] - coordinates[7][0], coordinates[4][1] - coordinates[7][1],
1142  coordinates[4][2] - coordinates[7][2] );
1143 
1144  xze.set( coordinates[3][0] - coordinates[7][0], coordinates[3][1] - coordinates[7][1],
1145  coordinates[3][2] - coordinates[7][2] );
1146 
1147  current_oddy = oddy_comp( xxi, xet, xze );
1148  if( current_oddy > oddy )
1149  {
1150  oddy = current_oddy;
1151  }
1152 
1153  if( oddy > 0 ) return (double)VERDICT_MIN( oddy, VERDICT_DBL_MAX );
1154  return (double)VERDICT_MAX( oddy, -VERDICT_DBL_MAX );
1155 }
1156 
1157 /*!
1158  the average Frobenius aspect of a hex
1159 
1160  NB (P. Pebay 01/20/07):
1161  this metric is calculated by averaging the 8 Frobenius aspects at
1162  each corner of the hex, when the reference corner is right isosceles.
1163 */
1164 C_FUNC_DEF double v_hex_med_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
1165 {
1166 
1167  VerdictVector node_pos[8];
1168  make_hex_nodes( coordinates, node_pos );
1169 
1170  double med_aspect_frobenius = 0.;
1171  VerdictVector xxi, xet, xze;
1172 
1173  // J(0,0,0):
1174 
1175  xxi = node_pos[1] - node_pos[0];
1176  xet = node_pos[3] - node_pos[0];
1177  xze = node_pos[4] - node_pos[0];
1178 
1179  med_aspect_frobenius += condition_comp( xxi, xet, xze );
1180  // J(1,0,0):
1181 
1182  xxi = node_pos[2] - node_pos[1];
1183  xet = node_pos[0] - node_pos[1];
1184  xze = node_pos[5] - node_pos[1];
1185 
1186  med_aspect_frobenius += condition_comp( xxi, xet, xze );
1187  // J(1,1,0):
1188 
1189  xxi = node_pos[3] - node_pos[2];
1190  xet = node_pos[1] - node_pos[2];
1191  xze = node_pos[6] - node_pos[2];
1192 
1193  med_aspect_frobenius += condition_comp( xxi, xet, xze );
1194  // J(0,1,0):
1195 
1196  xxi = node_pos[0] - node_pos[3];
1197  xet = node_pos[2] - node_pos[3];
1198  xze = node_pos[7] - node_pos[3];
1199 
1200  med_aspect_frobenius += condition_comp( xxi, xet, xze );
1201  // J(0,0,1):
1202 
1203  xxi = node_pos[7] - node_pos[4];
1204  xet = node_pos[5] - node_pos[4];
1205  xze = node_pos[0] - node_pos[4];
1206 
1207  med_aspect_frobenius += condition_comp( xxi, xet, xze );
1208  // J(1,0,1):
1209 
1210  xxi = node_pos[4] - node_pos[5];
1211  xet = node_pos[6] - node_pos[5];
1212  xze = node_pos[1] - node_pos[5];
1213 
1214  med_aspect_frobenius += condition_comp( xxi, xet, xze );
1215  // J(1,1,1):
1216 
1217  xxi = node_pos[5] - node_pos[6];
1218  xet = node_pos[7] - node_pos[6];
1219  xze = node_pos[2] - node_pos[6];
1220 
1221  med_aspect_frobenius += condition_comp( xxi, xet, xze );
1222  // J(1,1,1):
1223 
1224  xxi = node_pos[6] - node_pos[7];
1225  xet = node_pos[4] - node_pos[7];
1226  xze = node_pos[3] - node_pos[7];
1227 
1228  med_aspect_frobenius += condition_comp( xxi, xet, xze );
1229  med_aspect_frobenius /= 24.;
1230 
1231  if( med_aspect_frobenius > 0 ) return (double)VERDICT_MIN( med_aspect_frobenius, VERDICT_DBL_MAX );
1232  return (double)VERDICT_MAX( med_aspect_frobenius, -VERDICT_DBL_MAX );
1233 }
1234 
1235 /*!
1236  maximum Frobenius condition number of a hex
1237 
1238  Maximum Frobenius condition number of the Jacobian matrix at 8 corners
1239  NB (P. Pebay 01/25/07):
1240  this metric is calculated by taking the maximum of the 8 Frobenius aspects at
1241  each corner of the hex, when the reference corner is right isosceles.
1242 */
1243 C_FUNC_DEF double v_hex_max_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
1244 {
1245 
1246  VerdictVector node_pos[8];
1247  make_hex_nodes( coordinates, node_pos );
1248 
1249  double condition = 0.0, current_condition;
1250  VerdictVector xxi, xet, xze;
1251 
1252  xxi = calc_hex_efg( 1, node_pos );
1253  xet = calc_hex_efg( 2, node_pos );
1254  xze = calc_hex_efg( 3, node_pos );
1255 
1256  current_condition = condition_comp( xxi, xet, xze );
1257  if( current_condition > condition )
1258  {
1259  condition = current_condition;
1260  }
1261 
1262  // J(0,0,0):
1263 
1264  xxi = node_pos[1] - node_pos[0];
1265  xet = node_pos[3] - node_pos[0];
1266  xze = node_pos[4] - node_pos[0];
1267 
1268  current_condition = condition_comp( xxi, xet, xze );
1269  if( current_condition > condition )
1270  {
1271  condition = current_condition;
1272  }
1273 
1274  // J(1,0,0):
1275 
1276  xxi = node_pos[2] - node_pos[1];
1277  xet = node_pos[0] - node_pos[1];
1278  xze = node_pos[5] - node_pos[1];
1279 
1280  current_condition = condition_comp( xxi, xet, xze );
1281  if( current_condition > condition )
1282  {
1283  condition = current_condition;
1284  }
1285 
1286  // J(1,1,0):
1287 
1288  xxi = node_pos[3] - node_pos[2];
1289  xet = node_pos[1] - node_pos[2];
1290  xze = node_pos[6] - node_pos[2];
1291 
1292  current_condition = condition_comp( xxi, xet, xze );
1293  if( current_condition > condition )
1294  {
1295  condition = current_condition;
1296  }
1297 
1298  // J(0,1,0):
1299 
1300  xxi = node_pos[0] - node_pos[3];
1301  xet = node_pos[2] - node_pos[3];
1302  xze = node_pos[7] - node_pos[3];
1303 
1304  current_condition = condition_comp( xxi, xet, xze );
1305  if( current_condition > condition )
1306  {
1307  condition = current_condition;
1308  }
1309 
1310  // J(0,0,1):
1311 
1312  xxi = node_pos[7] - node_pos[4];
1313  xet = node_pos[5] - node_pos[4];
1314  xze = node_pos[0] - node_pos[4];
1315 
1316  current_condition = condition_comp( xxi, xet, xze );
1317  if( current_condition > condition )
1318  {
1319  condition = current_condition;
1320  }
1321 
1322  // J(1,0,1):
1323 
1324  xxi = node_pos[4] - node_pos[5];
1325  xet = node_pos[6] - node_pos[5];
1326  xze = node_pos[1] - node_pos[5];
1327 
1328  current_condition = condition_comp( xxi, xet, xze );
1329  if( current_condition > condition )
1330  {
1331  condition = current_condition;
1332  }
1333 
1334  // J(1,1,1):
1335 
1336  xxi = node_pos[5] - node_pos[6];
1337  xet = node_pos[7] - node_pos[6];
1338  xze = node_pos[2] - node_pos[6];
1339 
1340  current_condition = condition_comp( xxi, xet, xze );
1341  if( current_condition > condition )
1342  {
1343  condition = current_condition;
1344  }
1345 
1346  // J(1,1,1):
1347 
1348  xxi = node_pos[6] - node_pos[7];
1349  xet = node_pos[4] - node_pos[7];
1350  xze = node_pos[3] - node_pos[7];
1351 
1352  current_condition = condition_comp( xxi, xet, xze );
1353  if( current_condition > condition )
1354  {
1355  condition = current_condition;
1356  }
1357 
1358  condition /= 3.0;
1359 
1360  if( condition > 0 ) return (double)VERDICT_MIN( condition, VERDICT_DBL_MAX );
1361  return (double)VERDICT_MAX( condition, -VERDICT_DBL_MAX );
1362 }
1363 
1364 /*!
1365  The maximum Frobenius condition of a hex, a.k.a. condition
1366  NB (P. Pebay 01/25/07):
1367  this method is maintained for backwards compatibility only.
1368  It will become deprecated at some point.
1369 
1370 */
1371 C_FUNC_DEF double v_hex_condition( int /*num_nodes*/, double coordinates[][3] )
1372 {
1373 
1374  return v_hex_max_aspect_frobenius( 8, coordinates );
1375 }
1376 
1377 /*!
1378  jacobian of a hex
1379 
1380  Minimum pointwise volume of local map at 8 corners & center of hex
1381 */
1382 C_FUNC_DEF double v_hex_jacobian( int /*num_nodes*/, double coordinates[][3] )
1383 {
1384 
1385  VerdictVector node_pos[8];
1386  make_hex_nodes( coordinates, node_pos );
1387 
1388  double jacobian = VERDICT_DBL_MAX;
1389  double current_jacobian;
1390  VerdictVector xxi, xet, xze;
1391 
1392  xxi = calc_hex_efg( 1, node_pos );
1393  xet = calc_hex_efg( 2, node_pos );
1394  xze = calc_hex_efg( 3, node_pos );
1395 
1396  current_jacobian = xxi % ( xet * xze ) / 64.0;
1397  if( current_jacobian < jacobian )
1398  {
1399  jacobian = current_jacobian;
1400  }
1401 
1402  // J(0,0,0):
1403 
1404  xxi = node_pos[1] - node_pos[0];
1405  xet = node_pos[3] - node_pos[0];
1406  xze = node_pos[4] - node_pos[0];
1407 
1408  current_jacobian = xxi % ( xet * xze );
1409  if( current_jacobian < jacobian )
1410  {
1411  jacobian = current_jacobian;
1412  }
1413 
1414  // J(1,0,0):
1415 
1416  xxi = node_pos[2] - node_pos[1];
1417  xet = node_pos[0] - node_pos[1];
1418  xze = node_pos[5] - node_pos[1];
1419 
1420  current_jacobian = xxi % ( xet * xze );
1421  if( current_jacobian < jacobian )
1422  {
1423  jacobian = current_jacobian;
1424  }
1425 
1426  // J(1,1,0):
1427 
1428  xxi = node_pos[3] - node_pos[2];
1429  xet = node_pos[1] - node_pos[2];
1430  xze = node_pos[6] - node_pos[2];
1431 
1432  current_jacobian = xxi % ( xet * xze );
1433  if( current_jacobian < jacobian )
1434  {
1435  jacobian = current_jacobian;
1436  }
1437 
1438  // J(0,1,0):
1439 
1440  xxi = node_pos[0] - node_pos[3];
1441  xet = node_pos[2] - node_pos[3];
1442  xze = node_pos[7] - node_pos[3];
1443 
1444  current_jacobian = xxi % ( xet * xze );
1445  if( current_jacobian < jacobian )
1446  {
1447  jacobian = current_jacobian;
1448  }
1449 
1450  // J(0,0,1):
1451 
1452  xxi = node_pos[7] - node_pos[4];
1453  xet = node_pos[5] - node_pos[4];
1454  xze = node_pos[0] - node_pos[4];
1455 
1456  current_jacobian = xxi % ( xet * xze );
1457  if( current_jacobian < jacobian )
1458  {
1459  jacobian = current_jacobian;
1460  }
1461 
1462  // J(1,0,1):
1463 
1464  xxi = node_pos[4] - node_pos[5];
1465  xet = node_pos[6] - node_pos[5];
1466  xze = node_pos[1] - node_pos[5];
1467 
1468  current_jacobian = xxi % ( xet * xze );
1469  if( current_jacobian < jacobian )
1470  {
1471  jacobian = current_jacobian;
1472  }
1473 
1474  // J(1,1,1):
1475 
1476  xxi = node_pos[5] - node_pos[6];
1477  xet = node_pos[7] - node_pos[6];
1478  xze = node_pos[2] - node_pos[6];
1479 
1480  current_jacobian = xxi % ( xet * xze );
1481  if( current_jacobian < jacobian )
1482  {
1483  jacobian = current_jacobian;
1484  }
1485 
1486  // J(0,1,1):
1487 
1488  xxi = node_pos[6] - node_pos[7];
1489  xet = node_pos[4] - node_pos[7];
1490  xze = node_pos[3] - node_pos[7];
1491 
1492  current_jacobian = xxi % ( xet * xze );
1493  if( current_jacobian < jacobian )
1494  {
1495  jacobian = current_jacobian;
1496  }
1497 
1498  if( jacobian > 0 ) return (double)VERDICT_MIN( jacobian, VERDICT_DBL_MAX );
1499  return (double)VERDICT_MAX( jacobian, -VERDICT_DBL_MAX );
1500 }
1501 
1502 /*!
1503  scaled jacobian of a hex
1504 
1505  Minimum Jacobian divided by the lengths of the 3 edge vectors
1506 */
1507 C_FUNC_DEF double v_hex_scaled_jacobian( int /*num_nodes*/, double coordinates[][3] )
1508 {
1509 
1510  double jacobi, min_norm_jac = VERDICT_DBL_MAX;
1511  // double min_jacobi = VERDICT_DBL_MAX;
1512  double temp_norm_jac, lengths;
1513  double len1_sq, len2_sq, len3_sq;
1514  VerdictVector xxi, xet, xze;
1515 
1516  VerdictVector node_pos[8];
1517  make_hex_nodes( coordinates, node_pos );
1518 
1519  xxi = calc_hex_efg( 1, node_pos );
1520  xet = calc_hex_efg( 2, node_pos );
1521  xze = calc_hex_efg( 3, node_pos );
1522 
1523  jacobi = xxi % ( xet * xze );
1524  // if( jacobi < min_jacobi) { min_jacobi = jacobi; }
1525 
1526  len1_sq = xxi.length_squared();
1527  len2_sq = xet.length_squared();
1528  len3_sq = xze.length_squared();
1529 
1530  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
1531  return (double)VERDICT_DBL_MAX;
1532 
1533  lengths = sqrt( len1_sq * len2_sq * len3_sq );
1534  temp_norm_jac = jacobi / lengths;
1535 
1536  if( temp_norm_jac < min_norm_jac )
1537  min_norm_jac = temp_norm_jac;
1538  else
1539  temp_norm_jac = jacobi;
1540 
1541  // J(0,0,0):
1542 
1543  xxi = node_pos[1] - node_pos[0];
1544  xet = node_pos[3] - node_pos[0];
1545  xze = node_pos[4] - node_pos[0];
1546 
1547  jacobi = xxi % ( xet * xze );
1548  // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
1549 
1550  len1_sq = xxi.length_squared();
1551  len2_sq = xet.length_squared();
1552  len3_sq = xze.length_squared();
1553 
1554  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
1555  return (double)VERDICT_DBL_MAX;
1556 
1557  lengths = sqrt( len1_sq * len2_sq * len3_sq );
1558  temp_norm_jac = jacobi / lengths;
1559  if( temp_norm_jac < min_norm_jac )
1560  min_norm_jac = temp_norm_jac;
1561  else
1562  temp_norm_jac = jacobi;
1563 
1564  // J(1,0,0):
1565 
1566  xxi = node_pos[2] - node_pos[1];
1567  xet = node_pos[0] - node_pos[1];
1568  xze = node_pos[5] - node_pos[1];
1569 
1570  jacobi = xxi % ( xet * xze );
1571  // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
1572 
1573  len1_sq = xxi.length_squared();
1574  len2_sq = xet.length_squared();
1575  len3_sq = xze.length_squared();
1576 
1577  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
1578  return (double)VERDICT_DBL_MAX;
1579 
1580  lengths = sqrt( len1_sq * len2_sq * len3_sq );
1581  temp_norm_jac = jacobi / lengths;
1582  if( temp_norm_jac < min_norm_jac )
1583  min_norm_jac = temp_norm_jac;
1584  else
1585  temp_norm_jac = jacobi;
1586 
1587  // J(1,1,0):
1588 
1589  xxi = node_pos[3] - node_pos[2];
1590  xet = node_pos[1] - node_pos[2];
1591  xze = node_pos[6] - node_pos[2];
1592 
1593  jacobi = xxi % ( xet * xze );
1594  // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
1595 
1596  len1_sq = xxi.length_squared();
1597  len2_sq = xet.length_squared();
1598  len3_sq = xze.length_squared();
1599 
1600  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
1601  return (double)VERDICT_DBL_MAX;
1602 
1603  lengths = sqrt( len1_sq * len2_sq * len3_sq );
1604  temp_norm_jac = jacobi / lengths;
1605  if( temp_norm_jac < min_norm_jac )
1606  min_norm_jac = temp_norm_jac;
1607  else
1608  temp_norm_jac = jacobi;
1609 
1610  // J(0,1,0):
1611 
1612  xxi = node_pos[0] - node_pos[3];
1613  xet = node_pos[2] - node_pos[3];
1614  xze = node_pos[7] - node_pos[3];
1615 
1616  jacobi = xxi % ( xet * xze );
1617  // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
1618 
1619  len1_sq = xxi.length_squared();
1620  len2_sq = xet.length_squared();
1621  len3_sq = xze.length_squared();
1622 
1623  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
1624  return (double)VERDICT_DBL_MAX;
1625 
1626  lengths = sqrt( len1_sq * len2_sq * len3_sq );
1627  temp_norm_jac = jacobi / lengths;
1628  if( temp_norm_jac < min_norm_jac )
1629  min_norm_jac = temp_norm_jac;
1630  else
1631  temp_norm_jac = jacobi;
1632 
1633  // J(0,0,1):
1634 
1635  xxi = node_pos[7] - node_pos[4];
1636  xet = node_pos[5] - node_pos[4];
1637  xze = node_pos[0] - node_pos[4];
1638 
1639  jacobi = xxi % ( xet * xze );
1640  // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
1641 
1642  len1_sq = xxi.length_squared();
1643  len2_sq = xet.length_squared();
1644  len3_sq = xze.length_squared();
1645 
1646  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
1647  return (double)VERDICT_DBL_MAX;
1648 
1649  lengths = sqrt( len1_sq * len2_sq * len3_sq );
1650  temp_norm_jac = jacobi / lengths;
1651  if( temp_norm_jac < min_norm_jac )
1652  min_norm_jac = temp_norm_jac;
1653  else
1654  temp_norm_jac = jacobi;
1655 
1656  // J(1,0,1):
1657 
1658  xxi = node_pos[4] - node_pos[5];
1659  xet = node_pos[6] - node_pos[5];
1660  xze = node_pos[1] - node_pos[5];
1661 
1662  jacobi = xxi % ( xet * xze );
1663  // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
1664 
1665  len1_sq = xxi.length_squared();
1666  len2_sq = xet.length_squared();
1667  len3_sq = xze.length_squared();
1668 
1669  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
1670  return (double)VERDICT_DBL_MAX;
1671 
1672  lengths = sqrt( len1_sq * len2_sq * len3_sq );
1673  temp_norm_jac = jacobi / lengths;
1674  if( temp_norm_jac < min_norm_jac )
1675  min_norm_jac = temp_norm_jac;
1676  else
1677  temp_norm_jac = jacobi;
1678 
1679  // J(1,1,1):
1680 
1681  xxi = node_pos[5] - node_pos[6];
1682  xet = node_pos[7] - node_pos[6];
1683  xze = node_pos[2] - node_pos[6];
1684 
1685  jacobi = xxi % ( xet * xze );
1686  // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
1687 
1688  len1_sq = xxi.length_squared();
1689  len2_sq = xet.length_squared();
1690  len3_sq = xze.length_squared();
1691 
1692  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
1693  return (double)VERDICT_DBL_MAX;
1694 
1695  lengths = sqrt( len1_sq * len2_sq * len3_sq );
1696  temp_norm_jac = jacobi / lengths;
1697  if( temp_norm_jac < min_norm_jac )
1698  min_norm_jac = temp_norm_jac;
1699  else
1700  temp_norm_jac = jacobi;
1701 
1702  // J(0,1,1):
1703 
1704  xxi = node_pos[6] - node_pos[7];
1705  xet = node_pos[4] - node_pos[7];
1706  xze = node_pos[3] - node_pos[7];
1707 
1708  jacobi = xxi % ( xet * xze );
1709  // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
1710 
1711  len1_sq = xxi.length_squared();
1712  len2_sq = xet.length_squared();
1713  len3_sq = xze.length_squared();
1714 
1715  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
1716  return (double)VERDICT_DBL_MAX;
1717 
1718  lengths = sqrt( len1_sq * len2_sq * len3_sq );
1719  temp_norm_jac = jacobi / lengths;
1720  if( temp_norm_jac < min_norm_jac ) min_norm_jac = temp_norm_jac;
1721  // else
1722  // temp_norm_jac = jacobi;
1723 
1724  if( min_norm_jac > 0 ) return (double)VERDICT_MIN( min_norm_jac, VERDICT_DBL_MAX );
1725  return (double)VERDICT_MAX( min_norm_jac, -VERDICT_DBL_MAX );
1726 }
1727 
1728 /*!
1729  shear of a hex
1730 
1731  3/Condition number of Jacobian Skew matrix
1732 */
1733 C_FUNC_DEF double v_hex_shear( int /*num_nodes*/, double coordinates[][3] )
1734 {
1735 
1736  double shear;
1737  double min_shear = 1.0;
1738  VerdictVector xxi, xet, xze;
1739  double det, len1_sq, len2_sq, len3_sq, lengths;
1740 
1741  VerdictVector node_pos[8];
1742  make_hex_nodes( coordinates, node_pos );
1743 
1744  // J(0,0,0):
1745 
1746  xxi = node_pos[1] - node_pos[0];
1747  xet = node_pos[3] - node_pos[0];
1748  xze = node_pos[4] - node_pos[0];
1749 
1750  len1_sq = xxi.length_squared();
1751  len2_sq = xet.length_squared();
1752  len3_sq = xze.length_squared();
1753 
1754  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;
1755 
1756  lengths = sqrt( len1_sq * len2_sq * len3_sq );
1757  det = xxi % ( xet * xze );
1758  if( det < VERDICT_DBL_MIN )
1759  {
1760  return 0;
1761  }
1762 
1763  shear = det / lengths;
1764  min_shear = VERDICT_MIN( shear, min_shear );
1765 
1766  // J(1,0,0):
1767 
1768  xxi = node_pos[2] - node_pos[1];
1769  xet = node_pos[0] - node_pos[1];
1770  xze = node_pos[5] - node_pos[1];
1771 
1772  len1_sq = xxi.length_squared();
1773  len2_sq = xet.length_squared();
1774  len3_sq = xze.length_squared();
1775 
1776  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;
1777 
1778  lengths = sqrt( len1_sq * len2_sq * len3_sq );
1779  det = xxi % ( xet * xze );
1780  if( det < VERDICT_DBL_MIN )
1781  {
1782  return 0;
1783  }
1784 
1785  shear = det / lengths;
1786  min_shear = VERDICT_MIN( shear, min_shear );
1787 
1788  // J(1,1,0):
1789 
1790  xxi = node_pos[3] - node_pos[2];
1791  xet = node_pos[1] - node_pos[2];
1792  xze = node_pos[6] - node_pos[2];
1793 
1794  len1_sq = xxi.length_squared();
1795  len2_sq = xet.length_squared();
1796  len3_sq = xze.length_squared();
1797 
1798  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;
1799 
1800  lengths = sqrt( len1_sq * len2_sq * len3_sq );
1801  det = xxi % ( xet * xze );
1802  if( det < VERDICT_DBL_MIN )
1803  {
1804  return 0;
1805  }
1806 
1807  shear = det / lengths;
1808  min_shear = VERDICT_MIN( shear, min_shear );
1809 
1810  // J(0,1,0):
1811 
1812  xxi = node_pos[0] - node_pos[3];
1813  xet = node_pos[2] - node_pos[3];
1814  xze = node_pos[7] - node_pos[3];
1815 
1816  len1_sq = xxi.length_squared();
1817  len2_sq = xet.length_squared();
1818  len3_sq = xze.length_squared();
1819 
1820  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;
1821 
1822  lengths = sqrt( len1_sq * len2_sq * len3_sq );
1823  det = xxi % ( xet * xze );
1824  if( det < VERDICT_DBL_MIN )
1825  {
1826  return 0;
1827  }
1828 
1829  shear = det / lengths;
1830  min_shear = VERDICT_MIN( shear, min_shear );
1831 
1832  // J(0,0,1):
1833 
1834  xxi = node_pos[7] - node_pos[4];
1835  xet = node_pos[5] - node_pos[4];
1836  xze = node_pos[0] - node_pos[4];
1837 
1838  len1_sq = xxi.length_squared();
1839  len2_sq = xet.length_squared();
1840  len3_sq = xze.length_squared();
1841 
1842  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;
1843 
1844  lengths = sqrt( len1_sq * len2_sq * len3_sq );
1845  det = xxi % ( xet * xze );
1846  if( det < VERDICT_DBL_MIN )
1847  {
1848  return 0;
1849  }
1850 
1851  shear = det / lengths;
1852  min_shear = VERDICT_MIN( shear, min_shear );
1853 
1854  // J(1,0,1):
1855 
1856  xxi = node_pos[4] - node_pos[5];
1857  xet = node_pos[6] - node_pos[5];
1858  xze = node_pos[1] - node_pos[5];
1859 
1860  len1_sq = xxi.length_squared();
1861  len2_sq = xet.length_squared();
1862  len3_sq = xze.length_squared();
1863 
1864  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;
1865 
1866  lengths = sqrt( len1_sq * len2_sq * len3_sq );
1867  det = xxi % ( xet * xze );
1868  if( det < VERDICT_DBL_MIN )
1869  {
1870  return 0;
1871  }
1872 
1873  shear = det / lengths;
1874  min_shear = VERDICT_MIN( shear, min_shear );
1875 
1876  // J(1,1,1):
1877 
1878  xxi = node_pos[5] - node_pos[6];
1879  xet = node_pos[7] - node_pos[6];
1880  xze = node_pos[2] - node_pos[6];
1881 
1882  len1_sq = xxi.length_squared();
1883  len2_sq = xet.length_squared();
1884  len3_sq = xze.length_squared();
1885 
1886  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;
1887 
1888  lengths = sqrt( len1_sq * len2_sq * len3_sq );
1889  det = xxi % ( xet * xze );
1890  if( det < VERDICT_DBL_MIN )
1891  {
1892  return 0;
1893  }
1894 
1895  shear = det / lengths;
1896  min_shear = VERDICT_MIN( shear, min_shear );
1897 
1898  // J(0,1,1):
1899 
1900  xxi = node_pos[6] - node_pos[7];
1901  xet = node_pos[4] - node_pos[7];
1902  xze = node_pos[3] - node_pos[7];
1903 
1904  len1_sq = xxi.length_squared();
1905  len2_sq = xet.length_squared();
1906  len3_sq = xze.length_squared();
1907 
1908  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;
1909 
1910  lengths = sqrt( len1_sq * len2_sq * len3_sq );
1911  det = xxi % ( xet * xze );
1912  if( det < VERDICT_DBL_MIN )
1913  {
1914  return 0;
1915  }
1916 
1917  shear = det / lengths;
1918  min_shear = VERDICT_MIN( shear, min_shear );
1919 
1920  if( min_shear <= VERDICT_DBL_MIN ) min_shear = 0;
1921 
1922  if( min_shear > 0 ) return (double)VERDICT_MIN( min_shear, VERDICT_DBL_MAX );
1923  return (double)VERDICT_MAX( min_shear, -VERDICT_DBL_MAX );
1924 }
1925 
1926 /*!
1927  shape of a hex
1928 
1929  3/Condition number of weighted Jacobian matrix
1930 */
1931 C_FUNC_DEF double v_hex_shape( int /*num_nodes*/, double coordinates[][3] )
1932 {
1933 
1934  double det, shape;
1935  double min_shape = 1.0;
1936  static const double two_thirds = 2.0 / 3.0;
1937 
1938  VerdictVector xxi, xet, xze;
1939 
1940  VerdictVector node_pos[8];
1941  make_hex_nodes( coordinates, node_pos );
1942 
1943  // J(0,0,0):
1944 
1945  xxi = node_pos[1] - node_pos[0];
1946  xet = node_pos[3] - node_pos[0];
1947  xze = node_pos[4] - node_pos[0];
1948 
1949  det = xxi % ( xet * xze );
1950  if( det > VERDICT_DBL_MIN )
1951  shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
1952  else
1953  return 0;
1954 
1955  if( shape < min_shape )
1956  {
1957  min_shape = shape;
1958  }
1959 
1960  // J(1,0,0):
1961 
1962  xxi = node_pos[2] - node_pos[1];
1963  xet = node_pos[0] - node_pos[1];
1964  xze = node_pos[5] - node_pos[1];
1965 
1966  det = xxi % ( xet * xze );
1967  if( det > VERDICT_DBL_MIN )
1968  shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
1969  else
1970  return 0;
1971 
1972  if( shape < min_shape )
1973  {
1974  min_shape = shape;
1975  }
1976 
1977  // J(1,1,0):
1978 
1979  xxi = node_pos[3] - node_pos[2];
1980  xet = node_pos[1] - node_pos[2];
1981  xze = node_pos[6] - node_pos[2];
1982 
1983  det = xxi % ( xet * xze );
1984  if( det > VERDICT_DBL_MIN )
1985  shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
1986  else
1987  return 0;
1988 
1989  if( shape < min_shape )
1990  {
1991  min_shape = shape;
1992  }
1993 
1994  // J(0,1,0):
1995 
1996  xxi = node_pos[0] - node_pos[3];
1997  xet = node_pos[2] - node_pos[3];
1998  xze = node_pos[7] - node_pos[3];
1999 
2000  det = xxi % ( xet * xze );
2001  if( det > VERDICT_DBL_MIN )
2002  shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
2003  else
2004  return 0;
2005 
2006  if( shape < min_shape )
2007  {
2008  min_shape = shape;
2009  }
2010 
2011  // J(0,0,1):
2012 
2013  xxi = node_pos[7] - node_pos[4];
2014  xet = node_pos[5] - node_pos[4];
2015  xze = node_pos[0] - node_pos[4];
2016 
2017  det = xxi % ( xet * xze );
2018  if( det > VERDICT_DBL_MIN )
2019  shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
2020  else
2021  return 0;
2022 
2023  if( shape < min_shape )
2024  {
2025  min_shape = shape;
2026  }
2027 
2028  // J(1,0,1):
2029 
2030  xxi = node_pos[4] - node_pos[5];
2031  xet = node_pos[6] - node_pos[5];
2032  xze = node_pos[1] - node_pos[5];
2033 
2034  det = xxi % ( xet * xze );
2035  if( det > VERDICT_DBL_MIN )
2036  shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
2037  else
2038  return 0;
2039 
2040  if( shape < min_shape )
2041  {
2042  min_shape = shape;
2043  }
2044 
2045  // J(1,1,1):
2046 
2047  xxi = node_pos[5] - node_pos[6];
2048  xet = node_pos[7] - node_pos[6];
2049  xze = node_pos[2] - node_pos[6];
2050 
2051  det = xxi % ( xet * xze );
2052  if( det > VERDICT_DBL_MIN )
2053  shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
2054  else
2055  return 0;
2056 
2057  if( shape < min_shape )
2058  {
2059  min_shape = shape;
2060  }
2061 
2062  // J(1,1,1):
2063 
2064  xxi = node_pos[6] - node_pos[7];
2065  xet = node_pos[4] - node_pos[7];
2066  xze = node_pos[3] - node_pos[7];
2067 
2068  det = xxi % ( xet * xze );
2069  if( det > VERDICT_DBL_MIN )
2070  shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
2071  else
2072  return 0;
2073 
2074  if( shape < min_shape )
2075  {
2076  min_shape = shape;
2077  }
2078 
2079  if( min_shape <= VERDICT_DBL_MIN ) min_shape = 0;
2080 
2081  if( min_shape > 0 ) return (double)VERDICT_MIN( min_shape, VERDICT_DBL_MAX );
2082  return (double)VERDICT_MAX( min_shape, -VERDICT_DBL_MAX );
2083 }
2084 
2085 /*!
2086  relative size of a hex
2087 
2088  Min( J, 1/J ), where J is determinant of weighted Jacobian matrix
2089 */
2090 C_FUNC_DEF double v_hex_relative_size_squared( int /*num_nodes*/, double coordinates[][3] )
2091 {
2092  double size = 0;
2093  double tau;
2094 
2095  VerdictVector xxi, xet, xze;
2096  double det, det_sum = 0;
2097 
2098  v_hex_get_weight( xxi, xet, xze );
2099 
2100  // This is the average relative size
2101  double detw = xxi % ( xet * xze );
2102 
2103  if( detw < VERDICT_DBL_MIN ) return 0;
2104 
2105  VerdictVector node_pos[8];
2106  make_hex_nodes( coordinates, node_pos );
2107 
2108  // J(0,0,0):
2109 
2110  xxi = node_pos[1] - node_pos[0];
2111  xet = node_pos[3] - node_pos[0];
2112  xze = node_pos[4] - node_pos[0];
2113 
2114  det = xxi % ( xet * xze );
2115  det_sum += det;
2116 
2117  // J(1,0,0):
2118 
2119  xxi = node_pos[2] - node_pos[1];
2120  xet = node_pos[0] - node_pos[1];
2121  xze = node_pos[5] - node_pos[1];
2122 
2123  det = xxi % ( xet * xze );
2124  det_sum += det;
2125 
2126  // J(0,1,0):
2127 
2128  xxi = node_pos[3] - node_pos[2];
2129  xet = node_pos[1] - node_pos[2];
2130  xze = node_pos[6] - node_pos[2];
2131 
2132  det = xxi % ( xet * xze );
2133  det_sum += det;
2134 
2135  // J(1,1,0):
2136 
2137  xxi = node_pos[0] - node_pos[3];
2138  xet = node_pos[2] - node_pos[3];
2139  xze = node_pos[7] - node_pos[3];
2140 
2141  det = xxi % ( xet * xze );
2142  det_sum += det;
2143 
2144  // J(0,1,0):
2145 
2146  xxi = node_pos[7] - node_pos[4];
2147  xet = node_pos[5] - node_pos[4];
2148  xze = node_pos[0] - node_pos[4];
2149 
2150  det = xxi % ( xet * xze );
2151  det_sum += det;
2152 
2153  // J(1,0,1):
2154 
2155  xxi = node_pos[4] - node_pos[5];
2156  xet = node_pos[6] - node_pos[5];
2157  xze = node_pos[1] - node_pos[5];
2158 
2159  det = xxi % ( xet * xze );
2160  det_sum += det;
2161 
2162  // J(1,1,1):
2163 
2164  xxi = node_pos[5] - node_pos[6];
2165  xet = node_pos[7] - node_pos[6];
2166  xze = node_pos[2] - node_pos[6];
2167 
2168  det = xxi % ( xet * xze );
2169  det_sum += det;
2170 
2171  // J(1,1,1):
2172 
2173  xxi = node_pos[6] - node_pos[7];
2174  xet = node_pos[4] - node_pos[7];
2175  xze = node_pos[3] - node_pos[7];
2176 
2177  det = xxi % ( xet * xze );
2178  det_sum += det;
2179 
2180  if( det_sum > VERDICT_DBL_MIN )
2181  {
2182  tau = det_sum / ( 8 * detw );
2183 
2184  tau = VERDICT_MIN( tau, 1.0 / tau );
2185 
2186  size = tau * tau;
2187  }
2188 
2189  if( size > 0 ) return (double)VERDICT_MIN( size, VERDICT_DBL_MAX );
2190  return (double)VERDICT_MAX( size, -VERDICT_DBL_MAX );
2191 }
2192 
2193 /*!
2194  shape and size of a hex
2195 
2196  Product of Shape and Relative Size
2197 */
2198 C_FUNC_DEF double v_hex_shape_and_size( int num_nodes, double coordinates[][3] )
2199 {
2200  double size = v_hex_relative_size_squared( num_nodes, coordinates );
2201  double shape = v_hex_shape( num_nodes, coordinates );
2202 
2203  double shape_size = size * shape;
2204 
2205  if( shape_size > 0 ) return (double)VERDICT_MIN( shape_size, VERDICT_DBL_MAX );
2206  return (double)VERDICT_MAX( shape_size, -VERDICT_DBL_MAX );
2207 }
2208 
2209 /*!
2210  shear and size of a hex
2211 
2212  Product of Shear and Relative Size
2213 */
2214 C_FUNC_DEF double v_hex_shear_and_size( int num_nodes, double coordinates[][3] )
2215 {
2216  double size = v_hex_relative_size_squared( num_nodes, coordinates );
2217  double shear = v_hex_shear( num_nodes, coordinates );
2218 
2219  double shear_size = shear * size;
2220 
2221  if( shear_size > 0 ) return (double)VERDICT_MIN( shear_size, VERDICT_DBL_MAX );
2222  return (double)VERDICT_MAX( shear_size, -VERDICT_DBL_MAX );
2223 }
2224 
2225 /*!
2226  distortion of a hex
2227 */
2228 C_FUNC_DEF double v_hex_distortion( int num_nodes, double coordinates[][3] )
2229 {
2230 
2231  // use 2x2 gauss points for linear hex and 3x3 for 2nd order hex
2232  int number_of_gauss_points = 0;
2233  if( num_nodes == 8 )
2234  // 2x2 quadrature rule
2235  number_of_gauss_points = 2;
2236  else if( num_nodes == 20 )
2237  // 3x3 quadrature rule
2238  number_of_gauss_points = 3;
2239 
2240  int number_dimension = 3;
2241  int total_number_of_gauss_points = number_of_gauss_points * number_of_gauss_points * number_of_gauss_points;
2242  double distortion = VERDICT_DBL_MAX;
2243 
2244  // maxTotalNumberGaussPoints =27, maxNumberNodes = 20
2245  // they are defined in GaussIntegration.hpp
2246  // This is used to make these arrays static.
2247  // I tried dynamically allocated arrays but the new and delete
2248  // was very expensive
2249 
2250  double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
2251  double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
2252  double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
2253  double dndy3[maxTotalNumberGaussPoints][maxNumberNodes];
2254  double weight[maxTotalNumberGaussPoints];
2255 
2256  // create an object of GaussIntegration
2257  GaussIntegration::initialize( number_of_gauss_points, num_nodes, number_dimension );
2258  GaussIntegration::calculate_shape_function_3d_hex();
2259  GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], dndy3[0], weight );
2260 
2261  VerdictVector xxi, xet, xze, xin;
2262 
2263  double jacobian, minimum_jacobian;
2264  double element_volume = 0.0;
2265  minimum_jacobian = VERDICT_DBL_MAX;
2266  // calculate element volume
2267  int ife, ja;
2268  for( ife = 0; ife < total_number_of_gauss_points; ife++ )
2269  {
2270 
2271  xxi.set( 0.0, 0.0, 0.0 );
2272  xet.set( 0.0, 0.0, 0.0 );
2273  xze.set( 0.0, 0.0, 0.0 );
2274 
2275  for( ja = 0; ja < num_nodes; ja++ )
2276  {
2277  xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
2278  xxi += dndy1[ife][ja] * xin;
2279  xet += dndy2[ife][ja] * xin;
2280  xze += dndy3[ife][ja] * xin;
2281  }
2282 
2283  jacobian = xxi % ( xet * xze );
2284  if( minimum_jacobian > jacobian ) minimum_jacobian = jacobian;
2285 
2286  element_volume += weight[ife] * jacobian;
2287  }
2288 
2289  // loop through all nodes
2290  double dndy1_at_node[maxNumberNodes][maxNumberNodes];
2291  double dndy2_at_node[maxNumberNodes][maxNumberNodes];
2292  double dndy3_at_node[maxNumberNodes][maxNumberNodes];
2293 
2294  GaussIntegration::calculate_derivative_at_nodes_3d( dndy1_at_node, dndy2_at_node, dndy3_at_node );
2295  int node_id;
2296  for( node_id = 0; node_id < num_nodes; node_id++ )
2297  {
2298 
2299  xxi.set( 0.0, 0.0, 0.0 );
2300  xet.set( 0.0, 0.0, 0.0 );
2301  xze.set( 0.0, 0.0, 0.0 );
2302 
2303  for( ja = 0; ja < num_nodes; ja++ )
2304  {
2305  xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
2306  xxi += dndy1_at_node[node_id][ja] * xin;
2307  xet += dndy2_at_node[node_id][ja] * xin;
2308  xze += dndy3_at_node[node_id][ja] * xin;
2309  }
2310 
2311  jacobian = xxi % ( xet * xze );
2312  if( minimum_jacobian > jacobian ) minimum_jacobian = jacobian;
2313  }
2314  distortion = minimum_jacobian / element_volume * 8.;
2315  return (double)distortion;
2316 }
2317 
2318 /*
2319 C_FUNC_DEF double hex_jac_normjac_oddy_cond( int choices[],
2320  double coordinates[][3],
2321  double answers[4] )
2322 {
2323 
2324  //Define variables
2325  int i;
2326 
2327  double xxi[3], xet[3], xze[3];
2328  double norm_jacobian = 0.0, current_norm_jac = 0.0;
2329  double jacobian = 0.0, current_jacobian = 0.0;
2330  double oddy = 0.0, current_oddy = 0.0;
2331  double condition = 0.0, current_condition = 0.0;
2332 
2333 
2334  for( i=0; i<3; i++)
2335  xxi[i] = calc_hex_efg( 2, i, coordinates );
2336  for( i=0; i<3; i++)
2337  xet[i] = calc_hex_efg( 3, i, coordinates );
2338  for( i=0; i<3; i++)
2339  xze[i] = calc_hex_efg( 6, i, coordinates );
2340 
2341  // norm jacobian and jacobian
2342  if( choices[0] || choices[1] )
2343  {
2344  current_jacobian = determinant( xxi, xet, xze );
2345  current_norm_jac = normalize_jacobian( current_jacobian,
2346  xxi, xet, xze );
2347 
2348  if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; }
2349 
2350  current_jacobian /= 64.0;
2351  if (current_jacobian < jacobian) { jacobian = current_jacobian; }
2352  }
2353 
2354  // oddy
2355  if( choices[2] )
2356  {
2357  current_oddy = oddy_comp( xxi, xet, xze );
2358  if ( current_oddy > oddy ) { oddy = current_oddy; }
2359  }
2360 
2361  // condition
2362  if( choices[3] )
2363  {
2364  current_condition = condition_comp( xxi, xet, xze );
2365  if ( current_condition > condition ) { condition = current_condition; }
2366  }
2367 
2368 
2369  for( i=0; i<3; i++)
2370  {
2371  xxi[i] = coordinates[1][i] - coordinates[0][i];
2372  xet[i] = coordinates[3][i] - coordinates[0][i];
2373  xze[i] = coordinates[4][i] - coordinates[0][i];
2374  }
2375 
2376  // norm jacobian and jacobian
2377  if( choices[0] || choices[1] )
2378  {
2379  current_jacobian = determinant( xxi, xet, xze );
2380  current_norm_jac = normalize_jacobian( current_jacobian,
2381  xxi, xet, xze );
2382 
2383  if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; }
2384  if (current_jacobian < jacobian) { jacobian = current_jacobian; }
2385  }
2386 
2387  // oddy
2388  if( choices[2] )
2389  {
2390  current_oddy = oddy_comp( xxi, xet, xze );
2391  if ( current_oddy > oddy ) { oddy = current_oddy; }
2392  }
2393 
2394  // condition
2395  if( choices[3] )
2396  {
2397  current_condition = condition_comp( xxi, xet, xze );
2398  if ( current_condition > condition ) { condition = current_condition; }
2399  }
2400 
2401 
2402  for( i=0; i<3; i++)
2403  {
2404  xxi[i] = coordinates[1][i] - coordinates[0][i];
2405  xet[i] = coordinates[2][i] - coordinates[1][i];
2406  xze[i] = coordinates[5][i] - coordinates[1][i];
2407  }
2408 
2409  // norm jacobian and jacobian
2410  if( choices[0] || choices[1] )
2411  {
2412  current_jacobian = determinant( xxi, xet, xze );
2413  current_norm_jac = normalize_jacobian( current_jacobian,
2414  xxi, xet, xze );
2415 
2416  if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; }
2417  if (current_jacobian < jacobian) { jacobian = current_jacobian; }
2418  }
2419 
2420  // oddy
2421  if( choices[2] )
2422  {
2423  current_oddy = oddy_comp( xxi, xet, xze );
2424  if ( current_oddy > oddy ) { oddy = current_oddy; }
2425  }
2426 
2427  // condition
2428  if( choices[3] )
2429  {
2430  current_condition = condition_comp( xxi, xet, xze );
2431  if ( current_condition > condition ) { condition = current_condition; }
2432  }
2433 
2434 
2435  for( i=0; i<3; i++)
2436  {
2437  xxi[i] = coordinates[2][i] - coordinates[3][i];
2438  xet[i] = coordinates[3][i] - coordinates[0][i];
2439  xze[i] = coordinates[7][i] - coordinates[3][i];
2440  }
2441 
2442  // norm jacobian and jacobian
2443  if( choices[0] || choices[1] )
2444  {
2445  current_jacobian = determinant( xxi, xet, xze );
2446  current_norm_jac = normalize_jacobian( current_jacobian,
2447  xxi, xet, xze );
2448 
2449  if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; }
2450  if (current_jacobian < jacobian) { jacobian = current_jacobian; }
2451  }
2452 
2453  // oddy
2454  if( choices[2] )
2455  {
2456  current_oddy = oddy_comp( xxi, xet, xze );
2457  if ( current_oddy > oddy ) { oddy = current_oddy; }
2458  }
2459 
2460  // condition
2461  if( choices[3] )
2462  {
2463  current_condition = condition_comp( xxi, xet, xze );
2464  if ( current_condition > condition ) { condition = current_condition; }
2465  }
2466 
2467 
2468  for( i=0; i<3; i++)
2469  {
2470  xxi[i] = coordinates[5][i] - coordinates[4][i];
2471  xet[i] = coordinates[7][i] - coordinates[4][i];
2472  xze[i] = coordinates[4][i] - coordinates[0][i];
2473  }
2474 
2475 
2476  // norm jacobian and jacobian
2477  if( choices[0] || choices[1] )
2478  {
2479  current_jacobian = determinant( xxi, xet, xze );
2480  current_norm_jac = normalize_jacobian( current_jacobian,
2481  xxi, xet, xze );
2482 
2483  if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; }
2484  if (current_jacobian < jacobian) { jacobian = current_jacobian; }
2485  }
2486 
2487  // oddy
2488  if( choices[2] )
2489  {
2490  current_oddy = oddy_comp( xxi, xet, xze );
2491  if ( current_oddy > oddy ) { oddy = current_oddy; }
2492  }
2493 
2494  // condition
2495  if( choices[3] )
2496  {
2497  current_condition = condition_comp( xxi, xet, xze );
2498  if ( current_condition > condition ) { condition = current_condition; }
2499  }
2500 
2501 
2502  for( i=0; i<3; i++)
2503  {
2504  xxi[i] = coordinates[2][i] - coordinates[3][i];
2505  xet[i] = coordinates[2][i] - coordinates[1][i];
2506  xze[i] = coordinates[6][i] - coordinates[2][i];
2507  }
2508 
2509  // norm jacobian and jacobian
2510  if( choices[0] || choices[1] )
2511  {
2512  current_jacobian = determinant( xxi, xet, xze );
2513  current_norm_jac = normalize_jacobian( current_jacobian,
2514  xxi, xet, xze );
2515 
2516  if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; }
2517  if (current_jacobian < jacobian) { jacobian = current_jacobian; }
2518  }
2519 
2520  // oddy
2521  if( choices[2] )
2522  {
2523  current_oddy = oddy_comp( xxi, xet, xze );
2524  if ( current_oddy > oddy ) { oddy = current_oddy; }
2525  }
2526 
2527  // condition
2528  if( choices[3] )
2529  {
2530  current_condition = condition_comp( xxi, xet, xze );
2531  if ( current_condition > condition ) { condition = current_condition; }
2532  }
2533 
2534 
2535  for( i=0; i<3; i++)
2536  {
2537  xxi[i] = coordinates[5][i] - coordinates[4][i];
2538  xet[i] = coordinates[6][i] - coordinates[5][i];
2539  xze[i] = coordinates[5][i] - coordinates[1][i];
2540  }
2541 
2542 
2543  // norm jacobian and jacobian
2544  if( choices[0] || choices[1] )
2545  {
2546  current_jacobian = determinant( xxi, xet, xze );
2547  current_norm_jac = normalize_jacobian( current_jacobian,
2548  xxi, xet, xze );
2549 
2550  if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; }
2551  if (current_jacobian < jacobian) { jacobian = current_jacobian; }
2552  }
2553 
2554  // oddy
2555  if( choices[2] )
2556  {
2557  current_oddy = oddy_comp( xxi, xet, xze );
2558  if ( current_oddy > oddy ) { oddy = current_oddy; }
2559  }
2560 
2561  // condition
2562  if( choices[3] )
2563  {
2564  current_condition = condition_comp( xxi, xet, xze );
2565  if ( current_condition > condition ) { condition = current_condition; }
2566  }
2567 
2568 
2569  for( i=0; i<3; i++)
2570  {
2571  xxi[i] = coordinates[6][i] - coordinates[7][i];
2572  xet[i] = coordinates[7][i] - coordinates[4][i];
2573  xze[i] = coordinates[7][i] - coordinates[3][i];
2574  }
2575 
2576 
2577  // norm jacobian and jacobian
2578  if( choices[0] || choices[1] )
2579  {
2580  current_jacobian = determinant( xxi, xet, xze );
2581  current_norm_jac = normalize_jacobian( current_jacobian,
2582  xxi, xet, xze );
2583 
2584  if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; }
2585  if (current_jacobian < jacobian) { jacobian = current_jacobian; }
2586  }
2587 
2588  // oddy
2589  if( choices[2] )
2590  {
2591  current_oddy = oddy_comp( xxi, xet, xze );
2592  if ( current_oddy > oddy ) { oddy = current_oddy; }
2593  }
2594 
2595  // condition
2596  if( choices[3] )
2597  {
2598  current_condition = condition_comp( xxi, xet, xze );
2599  if ( current_condition > condition ) { condition = current_condition; }
2600  }
2601 
2602 
2603  for( i=0; i<3; i++)
2604  {
2605  xxi[i] = coordinates[6][i] - coordinates[7][i];
2606  xet[i] = coordinates[6][i] - coordinates[5][i];
2607  xze[i] = coordinates[6][i] - coordinates[2][i];
2608  }
2609 
2610  // norm jacobian and jacobian
2611  if( choices[0] || choices[1] )
2612  {
2613  current_jacobian = determinant( xxi, xet, xze );
2614  current_norm_jac = normalize_jacobian( current_jacobian,
2615  xxi, xet, xze );
2616 
2617  if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; }
2618  if (current_jacobian < jacobian) { jacobian = current_jacobian; }
2619  }
2620 
2621  // oddy
2622  if( choices[2] )
2623  {
2624  current_oddy = oddy_comp( xxi, xet, xze );
2625  if ( current_oddy > oddy ) { oddy = current_oddy; }
2626  }
2627 
2628  // condition
2629  if( choices[3] )
2630  {
2631  current_condition = condition_comp( xxi, xet, xze );
2632  if ( current_condition > condition ) { condition = current_condition; }
2633 
2634  condition /= 3.0;
2635  }
2636 
2637 
2638  answers[0] = jacobian;
2639  answers[1] = norm_jacobian;
2640  answers[2] = oddy;
2641  answers[3] = condition;
2642 
2643  return 1.0;
2644 
2645 }
2646 */
2647 
2648 /*!
2649  multiple quality metrics of a hex
2650 */
2651 C_FUNC_DEF void v_hex_quality( int num_nodes,
2652  double coordinates[][3],
2653  unsigned int metrics_request_flag,
2654  HexMetricVals* metric_vals )
2655 {
2656  memset( metric_vals, 0, sizeof( HexMetricVals ) );
2657 
2658  // max edge ratio, skew, taper, and volume
2659  if( metrics_request_flag & ( V_HEX_MAX_EDGE_RATIO | V_HEX_SKEW | V_HEX_TAPER ) )
2660  {
2661  VerdictVector node_pos[8];
2662  make_hex_nodes( coordinates, node_pos );
2663 
2664  VerdictVector efg1, efg2, efg3;
2665  efg1 = calc_hex_efg( 1, node_pos );
2666  efg2 = calc_hex_efg( 2, node_pos );
2667  efg3 = calc_hex_efg( 3, node_pos );
2668 
2669  if( metrics_request_flag & V_HEX_MAX_EDGE_RATIO )
2670  {
2671  double max_edge_ratio_12, max_edge_ratio_13, max_edge_ratio_23;
2672 
2673  double mag_efg1 = efg1.length();
2674  double mag_efg2 = efg2.length();
2675  double mag_efg3 = efg3.length();
2676 
2677  max_edge_ratio_12 = safe_ratio( VERDICT_MAX( mag_efg1, mag_efg2 ), VERDICT_MIN( mag_efg1, mag_efg2 ) );
2678  max_edge_ratio_13 = safe_ratio( VERDICT_MAX( mag_efg1, mag_efg3 ), VERDICT_MIN( mag_efg1, mag_efg3 ) );
2679  max_edge_ratio_23 = safe_ratio( VERDICT_MAX( mag_efg2, mag_efg3 ), VERDICT_MIN( mag_efg2, mag_efg3 ) );
2680 
2681  metric_vals->max_edge_ratio =
2682  (double)VERDICT_MAX( max_edge_ratio_12, VERDICT_MAX( max_edge_ratio_13, max_edge_ratio_23 ) );
2683  }
2684 
2685  if( metrics_request_flag & V_HEX_SKEW )
2686  {
2687 
2688  VerdictVector vec1 = efg1;
2689  VerdictVector vec2 = efg2;
2690  VerdictVector vec3 = efg3;
2691 
2692  if( vec1.normalize() <= VERDICT_DBL_MIN || vec2.normalize() <= VERDICT_DBL_MIN ||
2693  vec3.normalize() <= VERDICT_DBL_MIN )
2694  {
2695  metric_vals->skew = (double)VERDICT_DBL_MAX;
2696  }
2697  else
2698  {
2699  double skewx = fabs( vec1 % vec2 );
2700  double skewy = fabs( vec1 % vec3 );
2701  double skewz = fabs( vec2 % vec3 );
2702 
2703  metric_vals->skew = (double)( VERDICT_MAX( skewx, VERDICT_MAX( skewy, skewz ) ) );
2704  }
2705  }
2706 
2707  if( metrics_request_flag & V_HEX_TAPER )
2708  {
2709  VerdictVector efg12 = calc_hex_efg( 12, node_pos );
2710  VerdictVector efg13 = calc_hex_efg( 13, node_pos );
2711  VerdictVector efg23 = calc_hex_efg( 23, node_pos );
2712 
2713  double taperx = fabs( safe_ratio( efg12.length(), VERDICT_MIN( efg1.length(), efg2.length() ) ) );
2714  double tapery = fabs( safe_ratio( efg13.length(), VERDICT_MIN( efg1.length(), efg3.length() ) ) );
2715  double taperz = fabs( safe_ratio( efg23.length(), VERDICT_MIN( efg2.length(), efg3.length() ) ) );
2716 
2717  metric_vals->taper = (double)VERDICT_MAX( taperx, VERDICT_MAX( tapery, taperz ) );
2718  }
2719  }
2720 
2721  if( metrics_request_flag & V_HEX_VOLUME )
2722  {
2723  metric_vals->volume = v_hex_volume( 8, coordinates );
2724  }
2725 
2726  if( metrics_request_flag &
2727  ( V_HEX_JACOBIAN | V_HEX_SCALED_JACOBIAN | V_HEX_CONDITION | V_HEX_SHEAR | V_HEX_SHAPE |
2728  V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE | V_HEX_STRETCH ) )
2729  {
2730 
2731  static const double two_thirds = 2.0 / 3.0;
2732  VerdictVector edges[12];
2733  // the length squares
2734  double length_squared[12];
2735  // make vectors from coordinates
2736  make_hex_edges( coordinates, edges );
2737 
2738  // calculate the length squares if we need to
2739  if( metrics_request_flag &
2740  ( V_HEX_JACOBIAN | V_HEX_SHEAR | V_HEX_SCALED_JACOBIAN | V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE |
2741  V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHEAR_AND_SIZE | V_HEX_STRETCH ) )
2742  make_edge_length_squares( edges, length_squared );
2743 
2744  double jacobian = VERDICT_DBL_MAX, scaled_jacobian = VERDICT_DBL_MAX, condition = 0.0, shear = 1.0, shape = 1.0,
2745  oddy = 0.0;
2746  double current_jacobian, current_scaled_jacobian, current_condition, current_shape, detw = 0, det_sum = 0,
2747  current_oddy;
2748  VerdictBoolean rel_size_error = VERDICT_FALSE;
2749 
2750  VerdictVector xxi, xet, xze;
2751 
2752  // get weights if we need based on average size of a hex
2753  if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
2754  {
2755  v_hex_get_weight( xxi, xet, xze );
2756  detw = xxi % ( xet * xze );
2757  if( detw < VERDICT_DBL_MIN ) rel_size_error = VERDICT_TRUE;
2758  }
2759 
2760  xxi = edges[0] - edges[2] + edges[4] - edges[6];
2761  xet = edges[1] - edges[3] + edges[5] - edges[7];
2762  xze = edges[8] + edges[9] + edges[10] + edges[11];
2763 
2764  current_jacobian = xxi % ( xet * xze ) / 64.0;
2765  if( current_jacobian < jacobian ) jacobian = current_jacobian;
2766 
2767  if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
2768  {
2769  current_jacobian *= 64.0;
2770  current_scaled_jacobian =
2771  current_jacobian / sqrt( xxi.length_squared() * xet.length_squared() * xze.length_squared() );
2772  if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
2773  }
2774 
2775  if( metrics_request_flag & V_HEX_CONDITION )
2776  {
2777  current_condition = condition_comp( xxi, xet, xze );
2778  if( current_condition > condition )
2779  {
2780  condition = current_condition;
2781  }
2782  }
2783 
2784  if( metrics_request_flag & V_HEX_ODDY )
2785  {
2786  current_oddy = oddy_comp( xxi, xet, xze );
2787  if( current_oddy > oddy )
2788  {
2789  oddy = current_oddy;
2790  }
2791  }
2792 
2793  // J(0,0,0)
2794  current_jacobian = edges[0] % ( -edges[3] * edges[8] );
2795  if( current_jacobian < jacobian ) jacobian = current_jacobian;
2796 
2797  if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
2798  {
2799  det_sum += current_jacobian;
2800  }
2801 
2802  if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
2803  {
2804  if( length_squared[0] <= VERDICT_DBL_MIN || length_squared[3] <= VERDICT_DBL_MIN ||
2805  length_squared[8] <= VERDICT_DBL_MIN )
2806  {
2807  current_scaled_jacobian = VERDICT_DBL_MAX;
2808  }
2809  else
2810  {
2811  current_scaled_jacobian =
2812  current_jacobian / sqrt( length_squared[0] * length_squared[3] * length_squared[8] );
2813  }
2814  if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
2815  }
2816 
2817  if( metrics_request_flag & V_HEX_CONDITION )
2818  {
2819  current_condition = condition_comp( edges[0], -edges[3], edges[8] );
2820  if( current_condition > condition )
2821  {
2822  condition = current_condition;
2823  }
2824  }
2825 
2826  if( metrics_request_flag & V_HEX_ODDY )
2827  {
2828  current_oddy = oddy_comp( edges[0], -edges[3], edges[8] );
2829  if( current_oddy > oddy )
2830  {
2831  oddy = current_oddy;
2832  }
2833  }
2834 
2835  if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
2836  {
2837  if( current_jacobian > VERDICT_DBL_MIN )
2838  current_shape = 3 * pow( current_jacobian, two_thirds ) /
2839  ( length_squared[0] + length_squared[3] + length_squared[8] );
2840  else
2841  current_shape = 0;
2842 
2843  if( current_shape < shape )
2844  {
2845  shape = current_shape;
2846  }
2847  }
2848 
2849  // J(1,0,0)
2850  current_jacobian = edges[1] % ( -edges[0] * edges[9] );
2851  if( current_jacobian < jacobian ) jacobian = current_jacobian;
2852 
2853  if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
2854  {
2855  det_sum += current_jacobian;
2856  }
2857 
2858  if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
2859  {
2860  if( length_squared[1] <= VERDICT_DBL_MIN || length_squared[0] <= VERDICT_DBL_MIN ||
2861  length_squared[9] <= VERDICT_DBL_MIN )
2862  {
2863  current_scaled_jacobian = VERDICT_DBL_MAX;
2864  }
2865  else
2866  {
2867  current_scaled_jacobian =
2868  current_jacobian / sqrt( length_squared[1] * length_squared[0] * length_squared[9] );
2869  }
2870  if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
2871  }
2872 
2873  if( metrics_request_flag & V_HEX_CONDITION )
2874  {
2875  current_condition = condition_comp( edges[1], -edges[0], edges[9] );
2876  if( current_condition > condition )
2877  {
2878  condition = current_condition;
2879  }
2880  }
2881 
2882  if( metrics_request_flag & V_HEX_ODDY )
2883  {
2884  current_oddy = oddy_comp( edges[1], -edges[0], edges[9] );
2885  if( current_oddy > oddy )
2886  {
2887  oddy = current_oddy;
2888  }
2889  }
2890 
2891  if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
2892  {
2893  if( current_jacobian > VERDICT_DBL_MIN )
2894  current_shape = 3 * pow( current_jacobian, two_thirds ) /
2895  ( length_squared[1] + length_squared[0] + length_squared[9] );
2896  else
2897  current_shape = 0;
2898 
2899  if( current_shape < shape )
2900  {
2901  shape = current_shape;
2902  }
2903  }
2904 
2905  // J(1,1,0)
2906  current_jacobian = edges[2] % ( -edges[1] * edges[10] );
2907  if( current_jacobian < jacobian ) jacobian = current_jacobian;
2908 
2909  if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
2910  {
2911  det_sum += current_jacobian;
2912  }
2913 
2914  if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
2915  {
2916  if( length_squared[2] <= VERDICT_DBL_MIN || length_squared[1] <= VERDICT_DBL_MIN ||
2917  length_squared[10] <= VERDICT_DBL_MIN )
2918  {
2919  current_scaled_jacobian = VERDICT_DBL_MAX;
2920  }
2921  else
2922  {
2923  current_scaled_jacobian =
2924  current_jacobian / sqrt( length_squared[2] * length_squared[1] * length_squared[10] );
2925  }
2926  if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
2927  }
2928 
2929  if( metrics_request_flag & V_HEX_CONDITION )
2930  {
2931  current_condition = condition_comp( edges[2], -edges[1], edges[10] );
2932  if( current_condition > condition )
2933  {
2934  condition = current_condition;
2935  }
2936  }
2937 
2938  if( metrics_request_flag & V_HEX_ODDY )
2939  {
2940  current_oddy = oddy_comp( edges[2], -edges[1], edges[10] );
2941  if( current_oddy > oddy )
2942  {
2943  oddy = current_oddy;
2944  }
2945  }
2946 
2947  if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
2948  {
2949  if( current_jacobian > VERDICT_DBL_MIN )
2950  current_shape = 3 * pow( current_jacobian, two_thirds ) /
2951  ( length_squared[2] + length_squared[1] + length_squared[10] );
2952  else
2953  current_shape = 0;
2954 
2955  if( current_shape < shape )
2956  {
2957  shape = current_shape;
2958  }
2959  }
2960 
2961  // J(0,1,0)
2962  current_jacobian = edges[3] % ( -edges[2] * edges[11] );
2963  if( current_jacobian < jacobian ) jacobian = current_jacobian;
2964 
2965  if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
2966  {
2967  det_sum += current_jacobian;
2968  }
2969 
2970  if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
2971  {
2972  if( length_squared[3] <= VERDICT_DBL_MIN || length_squared[2] <= VERDICT_DBL_MIN ||
2973  length_squared[11] <= VERDICT_DBL_MIN )
2974  {
2975  current_scaled_jacobian = VERDICT_DBL_MAX;
2976  }
2977  else
2978  {
2979  current_scaled_jacobian =
2980  current_jacobian / sqrt( length_squared[3] * length_squared[2] * length_squared[11] );
2981  }
2982  if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
2983  }
2984 
2985  if( metrics_request_flag & V_HEX_CONDITION )
2986  {
2987  current_condition = condition_comp( edges[3], -edges[2], edges[11] );
2988  if( current_condition > condition )
2989  {
2990  condition = current_condition;
2991  }
2992  }
2993 
2994  if( metrics_request_flag & V_HEX_ODDY )
2995  {
2996  current_oddy = oddy_comp( edges[3], -edges[2], edges[11] );
2997  if( current_oddy > oddy )
2998  {
2999  oddy = current_oddy;
3000  }
3001  }
3002 
3003  if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
3004  {
3005  if( current_jacobian > VERDICT_DBL_MIN )
3006  current_shape = 3 * pow( current_jacobian, two_thirds ) /
3007  ( length_squared[3] + length_squared[2] + length_squared[11] );
3008  else
3009  current_shape = 0;
3010 
3011  if( current_shape < shape )
3012  {
3013  shape = current_shape;
3014  }
3015  }
3016 
3017  // J(0,0,1)
3018  current_jacobian = edges[4] % ( -edges[8] * -edges[7] );
3019  if( current_jacobian < jacobian ) jacobian = current_jacobian;
3020 
3021  if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
3022  {
3023  det_sum += current_jacobian;
3024  }
3025 
3026  if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
3027  {
3028  if( length_squared[4] <= VERDICT_DBL_MIN || length_squared[8] <= VERDICT_DBL_MIN ||
3029  length_squared[7] <= VERDICT_DBL_MIN )
3030  {
3031  current_scaled_jacobian = VERDICT_DBL_MAX;
3032  }
3033  else
3034  {
3035  current_scaled_jacobian =
3036  current_jacobian / sqrt( length_squared[4] * length_squared[8] * length_squared[7] );
3037  }
3038  if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
3039  }
3040 
3041  if( metrics_request_flag & V_HEX_CONDITION )
3042  {
3043  current_condition = condition_comp( edges[4], -edges[8], -edges[7] );
3044  if( current_condition > condition )
3045  {
3046  condition = current_condition;
3047  }
3048  }
3049 
3050  if( metrics_request_flag & V_HEX_ODDY )
3051  {
3052  current_oddy = oddy_comp( edges[4], -edges[8], -edges[7] );
3053  if( current_oddy > oddy )
3054  {
3055  oddy = current_oddy;
3056  }
3057  }
3058 
3059  if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
3060  {
3061  if( current_jacobian > VERDICT_DBL_MIN )
3062  current_shape = 3 * pow( current_jacobian, two_thirds ) /
3063  ( length_squared[4] + length_squared[8] + length_squared[7] );
3064  else
3065  current_shape = 0;
3066 
3067  if( current_shape < shape )
3068  {
3069  shape = current_shape;
3070  }
3071  }
3072 
3073  // J(1,0,1)
3074  current_jacobian = -edges[4] % ( edges[5] * -edges[9] );
3075  if( current_jacobian < jacobian ) jacobian = current_jacobian;
3076 
3077  if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
3078  {
3079  det_sum += current_jacobian;
3080  }
3081 
3082  if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
3083  {
3084  if( length_squared[4] <= VERDICT_DBL_MIN || length_squared[5] <= VERDICT_DBL_MIN ||
3085  length_squared[9] <= VERDICT_DBL_MIN )
3086  {
3087  current_scaled_jacobian = VERDICT_DBL_MAX;
3088  }
3089  else
3090  {
3091  current_scaled_jacobian =
3092  current_jacobian / sqrt( length_squared[4] * length_squared[5] * length_squared[9] );
3093  }
3094  if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
3095  }
3096 
3097  if( metrics_request_flag & V_HEX_CONDITION )
3098  {
3099  current_condition = condition_comp( -edges[4], edges[5], -edges[9] );
3100  if( current_condition > condition )
3101  {
3102  condition = current_condition;
3103  }
3104  }
3105 
3106  if( metrics_request_flag & V_HEX_ODDY )
3107  {
3108  current_oddy = oddy_comp( -edges[4], edges[5], -edges[9] );
3109  if( current_oddy > oddy )
3110  {
3111  oddy = current_oddy;
3112  }
3113  }
3114 
3115  if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
3116  {
3117  if( current_jacobian > VERDICT_DBL_MIN )
3118  current_shape = 3 * pow( current_jacobian, two_thirds ) /
3119  ( length_squared[4] + length_squared[5] + length_squared[9] );
3120  else
3121  current_shape = 0;
3122 
3123  if( current_shape < shape )
3124  {
3125  shape = current_shape;
3126  }
3127  }
3128 
3129  // J(1,1,1)
3130  current_jacobian = -edges[5] % ( edges[6] * -edges[10] );
3131  if( current_jacobian < jacobian ) jacobian = current_jacobian;
3132 
3133  if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
3134  {
3135  det_sum += current_jacobian;
3136  }
3137 
3138  if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
3139  {
3140  if( length_squared[5] <= VERDICT_DBL_MIN || length_squared[6] <= VERDICT_DBL_MIN ||
3141  length_squared[10] <= VERDICT_DBL_MIN )
3142  {
3143  current_scaled_jacobian = VERDICT_DBL_MAX;
3144  }
3145  else
3146  {
3147  current_scaled_jacobian =
3148  current_jacobian / sqrt( length_squared[5] * length_squared[6] * length_squared[10] );
3149  }
3150  if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
3151  }
3152 
3153  if( metrics_request_flag & V_HEX_CONDITION )
3154  {
3155  current_condition = condition_comp( -edges[5], edges[6], -edges[10] );
3156  if( current_condition > condition )
3157  {
3158  condition = current_condition;
3159  }
3160  }
3161 
3162  if( metrics_request_flag & V_HEX_ODDY )
3163  {
3164  current_oddy = oddy_comp( -edges[5], edges[6], -edges[10] );
3165  if( current_oddy > oddy )
3166  {
3167  oddy = current_oddy;
3168  }
3169  }
3170 
3171  if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
3172  {
3173  if( current_jacobian > VERDICT_DBL_MIN )
3174  current_shape = 3 * pow( current_jacobian, two_thirds ) /
3175  ( length_squared[5] + length_squared[6] + length_squared[10] );
3176  else
3177  current_shape = 0;
3178 
3179  if( current_shape < shape )
3180  {
3181  shape = current_shape;
3182  }
3183  }
3184 
3185  // J(0,1,1)
3186  current_jacobian = -edges[6] % ( edges[7] * -edges[11] );
3187  if( current_jacobian < jacobian ) jacobian = current_jacobian;
3188 
3189  if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
3190  {
3191  det_sum += current_jacobian;
3192  }
3193 
3194  if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
3195  {
3196  if( length_squared[6] <= VERDICT_DBL_MIN || length_squared[7] <= VERDICT_DBL_MIN ||
3197  length_squared[11] <= VERDICT_DBL_MIN )
3198  {
3199  current_scaled_jacobian = VERDICT_DBL_MAX;
3200  }
3201  else
3202  {
3203  current_scaled_jacobian =
3204  current_jacobian / sqrt( length_squared[6] * length_squared[7] * length_squared[11] );
3205  }
3206  if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
3207  }
3208 
3209  if( metrics_request_flag & V_HEX_CONDITION )
3210  {
3211  current_condition = condition_comp( -edges[6], edges[7], -edges[11] );
3212  if( current_condition > condition )
3213  {
3214  condition = current_condition;
3215  }
3216  }
3217 
3218  if( metrics_request_flag & V_HEX_ODDY )
3219  {
3220  current_oddy = oddy_comp( -edges[6], edges[7], -edges[11] );
3221  if( current_oddy > oddy )
3222  {
3223  oddy = current_oddy;
3224  }
3225  }
3226 
3227  if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
3228  {
3229  if( current_jacobian > VERDICT_DBL_MIN )
3230  current_shape = 3 * pow( current_jacobian, two_thirds ) /
3231  ( length_squared[6] + length_squared[7] + length_squared[11] );
3232  else
3233  current_shape = 0;
3234 
3235  if( current_shape < shape )
3236  {
3237  shape = current_shape;
3238  }
3239  }
3240 
3241  if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
3242  {
3243  if( det_sum > VERDICT_DBL_MIN && rel_size_error != VERDICT_TRUE )
3244  {
3245  double tau = det_sum / ( 8 * detw );
3246  metric_vals->relative_size_squared = (double)VERDICT_MIN( tau * tau, 1.0 / tau / tau );
3247  }
3248  else
3249  metric_vals->relative_size_squared = 0.0;
3250  }
3251 
3252  // set values from above calculations
3253  if( metrics_request_flag & V_HEX_JACOBIAN ) metric_vals->jacobian = (double)jacobian;
3254 
3255  if( metrics_request_flag & V_HEX_SCALED_JACOBIAN ) metric_vals->scaled_jacobian = (double)scaled_jacobian;
3256 
3257  if( metrics_request_flag & V_HEX_CONDITION ) metric_vals->condition = (double)( condition / 3.0 );
3258 
3259  if( metrics_request_flag & V_HEX_SHEAR )
3260  {
3261  if( shear < VERDICT_DBL_MIN ) // shear has range 0 to +1
3262  shear = 0;
3263  metric_vals->shear = (double)shear;
3264  }
3265 
3266  if( metrics_request_flag & V_HEX_SHAPE ) metric_vals->shape = (double)shape;
3267 
3268  if( metrics_request_flag & V_HEX_SHAPE_AND_SIZE )
3269  metric_vals->shape_and_size = (double)( shape * metric_vals->relative_size_squared );
3270 
3271  if( metrics_request_flag & V_HEX_SHEAR_AND_SIZE )
3272  metric_vals->shear_and_size = (double)( shear * metric_vals->relative_size_squared );
3273 
3274  if( metrics_request_flag & V_HEX_ODDY ) metric_vals->oddy = (double)oddy;
3275 
3276  if( metrics_request_flag & V_HEX_STRETCH )
3277  {
3278  static const double HEX_STRETCH_SCALE_FACTOR = sqrt( 3.0 );
3279  double min_edge = length_squared[0];
3280  for( int j = 1; j < 12; j++ )
3281  min_edge = VERDICT_MIN( min_edge, length_squared[j] );
3282 
3283  double max_diag = diag_length( 1, coordinates );
3284 
3285  metric_vals->stretch = (double)( HEX_STRETCH_SCALE_FACTOR * ( safe_ratio( sqrt( min_edge ), max_diag ) ) );
3286  }
3287  }
3288 
3289  if( metrics_request_flag & V_HEX_DIAGONAL ) metric_vals->diagonal = v_hex_diagonal( num_nodes, coordinates );
3290  if( metrics_request_flag & V_HEX_DIMENSION ) metric_vals->dimension = v_hex_dimension( num_nodes, coordinates );
3291  if( metrics_request_flag & V_HEX_DISTORTION ) metric_vals->distortion = v_hex_distortion( num_nodes, coordinates );
3292 
3293  // take care of any overflow problems
3294  // max_edge_ratio
3295  if( metric_vals->max_edge_ratio > 0 )
3296  metric_vals->max_edge_ratio = (double)VERDICT_MIN( metric_vals->max_edge_ratio, VERDICT_DBL_MAX );
3297  else
3298  metric_vals->max_edge_ratio = (double)VERDICT_MAX( metric_vals->max_edge_ratio, -VERDICT_DBL_MAX );
3299 
3300  // skew
3301  if( metric_vals->skew > 0 )
3302  metric_vals->skew = (double)VERDICT_MIN( metric_vals->skew, VERDICT_DBL_MAX );
3303  else
3304  metric_vals->skew = (double)VERDICT_MAX( metric_vals->skew, -VERDICT_DBL_MAX );
3305 
3306  // taper
3307  if( metric_vals->taper > 0 )
3308  metric_vals->taper = (double)VERDICT_MIN( metric_vals->taper, VERDICT_DBL_MAX );
3309  else
3310  metric_vals->taper = (double)VERDICT_MAX( metric_vals->taper, -VERDICT_DBL_MAX );
3311 
3312  // volume
3313  if( metric_vals->volume > 0 )
3314  metric_vals->volume = (double)VERDICT_MIN( metric_vals->volume, VERDICT_DBL_MAX );
3315  else
3316  metric_vals->volume = (double)VERDICT_MAX( metric_vals->volume, -VERDICT_DBL_MAX );
3317 
3318  // stretch
3319  if( metric_vals->stretch > 0 )
3320  metric_vals->stretch = (double)VERDICT_MIN( metric_vals->stretch, VERDICT_DBL_MAX );
3321  else
3322  metric_vals->stretch = (double)VERDICT_MAX( metric_vals->stretch, -VERDICT_DBL_MAX );
3323 
3324  // diagonal
3325  if( metric_vals->diagonal > 0 )
3326  metric_vals->diagonal = (double)VERDICT_MIN( metric_vals->diagonal, VERDICT_DBL_MAX );
3327  else
3328  metric_vals->diagonal = (double)VERDICT_MAX( metric_vals->diagonal, -VERDICT_DBL_MAX );
3329 
3330  // dimension
3331  if( metric_vals->dimension > 0 )
3332  metric_vals->dimension = (double)VERDICT_MIN( metric_vals->dimension, VERDICT_DBL_MAX );
3333  else
3334  metric_vals->dimension = (double)VERDICT_MAX( metric_vals->dimension, -VERDICT_DBL_MAX );
3335 
3336  // oddy
3337  if( metric_vals->oddy > 0 )
3338  metric_vals->oddy = (double)VERDICT_MIN( metric_vals->oddy, VERDICT_DBL_MAX );
3339  else
3340  metric_vals->oddy = (double)VERDICT_MAX( metric_vals->oddy, -VERDICT_DBL_MAX );
3341 
3342  // condition
3343  if( metric_vals->condition > 0 )
3344  metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
3345  else
3346  metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
3347 
3348  // jacobian
3349  if( metric_vals->jacobian > 0 )
3350  metric_vals->jacobian = (double)VERDICT_MIN( metric_vals->jacobian, VERDICT_DBL_MAX );
3351  else
3352  metric_vals->jacobian = (double)VERDICT_MAX( metric_vals->jacobian, -VERDICT_DBL_MAX );
3353 
3354  // scaled_jacobian
3355  if( metric_vals->scaled_jacobian > 0 )
3356  metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
3357  else
3358  metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
3359 
3360  // shear
3361  if( metric_vals->shear > 0 )
3362  metric_vals->shear = (double)VERDICT_MIN( metric_vals->shear, VERDICT_DBL_MAX );
3363  else
3364  metric_vals->shear = (double)VERDICT_MAX( metric_vals->shear, -VERDICT_DBL_MAX );
3365 
3366  // shape
3367  if( metric_vals->shape > 0 )
3368  metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
3369  else
3370  metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
3371 
3372  // relative_size_squared
3373  if( metric_vals->relative_size_squared > 0 )
3374  metric_vals->relative_size_squared = (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
3375  else
3376  metric_vals->relative_size_squared =
3377  (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
3378 
3379  // shape_and_size
3380  if( metric_vals->shape_and_size > 0 )
3381  metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
3382  else
3383  metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
3384 
3385  // shear_and_size
3386  if( metric_vals->shear_and_size > 0 )
3387  metric_vals->shear_and_size = (double)VERDICT_MIN( metric_vals->shear_and_size, VERDICT_DBL_MAX );
3388  else
3389  metric_vals->shear_and_size = (double)VERDICT_MAX( metric_vals->shear_and_size, -VERDICT_DBL_MAX );
3390 
3391  // distortion
3392  if( metric_vals->distortion > 0 )
3393  metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
3394  else
3395  metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
3396 
3397  if( metrics_request_flag & V_HEX_MED_ASPECT_FROBENIUS )
3398  metric_vals->med_aspect_frobenius = v_hex_med_aspect_frobenius( 8, coordinates );
3399 
3400  if( metrics_request_flag & V_HEX_EDGE_RATIO ) metric_vals->edge_ratio = v_hex_edge_ratio( 8, coordinates );
3401 }