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