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