 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 /*=========================================================================
2
3 Module: $RCSfile: V_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 }
1.9.1.