MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /*========================================================================= 00002 00003 Module: $RCSfile: V_HexMetric.cpp,v $ 00004 00005 Copyright (c) 2006 Sandia Corporation. 00006 All rights reserved. 00007 See Copyright.txt or http://www.kitware.com/Copyright.htm for details. 00008 00009 This software is distributed WITHOUT ANY WARRANTY; without even 00010 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 00011 PURPOSE. See the above copyright notice for more information. 00012 00013 =========================================================================*/ 00014 00015 /* 00016 * 00017 * HexMetric.cpp contains quality calculations for hexes 00018 * 00019 * This file is part of VERDICT 00020 * 00021 */ 00022 00023 #define VERDICT_EXPORTS 00024 00025 #include "moab/verdict.h" 00026 #include "VerdictVector.hpp" 00027 #include "V_GaussIntegration.hpp" 00028 #include "verdict_defines.hpp" 00029 #include <memory.h> 00030 00031 #if defined( __BORLANDC__ ) 00032 #pragma warn - 8004 /* "assigned a value that is never used" */ 00033 #endif 00034 00035 //! the average volume of a hex 00036 double verdict_hex_size = 0; 00037 00038 //! weights based on the average size of a hex 00039 int v_hex_get_weight( VerdictVector& v1, VerdictVector& v2, VerdictVector& v3 ) 00040 { 00041 00042 if( verdict_hex_size == 0 ) return 0; 00043 00044 v1.set( 1, 0, 0 ); 00045 v2.set( 0, 1, 0 ); 00046 v3.set( 0, 0, 1 ); 00047 00048 double scale = pow( verdict_hex_size / ( v1 % ( v2 * v3 ) ), 0.33333333333 ); 00049 v1 *= scale; 00050 v2 *= scale; 00051 v3 *= scale; 00052 00053 return 1; 00054 } 00055 00056 //! returns the average volume of a hex 00057 C_FUNC_DEF void v_set_hex_size( double size ) 00058 { 00059 verdict_hex_size = size; 00060 } 00061 00062 #define make_hex_nodes( coord, pos ) \ 00063 for( int mhcii = 0; mhcii < 8; mhcii++ ) \ 00064 { \ 00065 pos[mhcii].set( coord[mhcii][0], coord[mhcii][1], coord[mhcii][2] ); \ 00066 } 00067 00068 #define make_edge_length_squares( edges, lengths ) \ 00069 { \ 00070 for( int melii = 0; melii < 12; melii++ ) \ 00071 lengths[melii] = edges[melii].length_squared(); \ 00072 } 00073 00074 //! make VerdictVectors from coordinates 00075 void make_hex_edges( double coordinates[][3], VerdictVector edges[12] ) 00076 { 00077 edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1], 00078 coordinates[1][2] - coordinates[0][2] ); 00079 edges[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1], 00080 coordinates[2][2] - coordinates[1][2] ); 00081 edges[2].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1], 00082 coordinates[3][2] - coordinates[2][2] ); 00083 edges[3].set( coordinates[0][0] - coordinates[3][0], coordinates[0][1] - coordinates[3][1], 00084 coordinates[0][2] - coordinates[3][2] ); 00085 edges[4].set( coordinates[5][0] - coordinates[4][0], coordinates[5][1] - coordinates[4][1], 00086 coordinates[5][2] - coordinates[4][2] ); 00087 edges[5].set( coordinates[6][0] - coordinates[5][0], coordinates[6][1] - coordinates[5][1], 00088 coordinates[6][2] - coordinates[5][2] ); 00089 edges[6].set( coordinates[7][0] - coordinates[6][0], coordinates[7][1] - coordinates[6][1], 00090 coordinates[7][2] - coordinates[6][2] ); 00091 edges[7].set( coordinates[4][0] - coordinates[7][0], coordinates[4][1] - coordinates[7][1], 00092 coordinates[4][2] - coordinates[7][2] ); 00093 edges[8].set( coordinates[4][0] - coordinates[0][0], coordinates[4][1] - coordinates[0][1], 00094 coordinates[4][2] - coordinates[0][2] ); 00095 edges[9].set( coordinates[5][0] - coordinates[1][0], coordinates[5][1] - coordinates[1][1], 00096 coordinates[5][2] - coordinates[1][2] ); 00097 edges[10].set( coordinates[6][0] - coordinates[2][0], coordinates[6][1] - coordinates[2][1], 00098 coordinates[6][2] - coordinates[2][2] ); 00099 edges[11].set( coordinates[7][0] - coordinates[3][0], coordinates[7][1] - coordinates[3][1], 00100 coordinates[7][2] - coordinates[3][2] ); 00101 } 00102 00103 /*! 00104 localizes hex coordinates 00105 */ 00106 void localize_hex_coordinates( double coordinates[][3], VerdictVector position[8] ) 00107 { 00108 00109 int ii; 00110 for( ii = 0; ii < 8; ii++ ) 00111 { 00112 position[ii].set( coordinates[ii][0], coordinates[ii][1], coordinates[ii][2] ); 00113 } 00114 00115 // ... Make centroid of element the center of coordinate system 00116 VerdictVector point_1256 = position[1]; 00117 point_1256 += position[2]; 00118 point_1256 += position[5]; 00119 point_1256 += position[6]; 00120 00121 VerdictVector point_0374 = position[0]; 00122 point_0374 += position[3]; 00123 point_0374 += position[7]; 00124 point_0374 += position[4]; 00125 00126 VerdictVector centroid = point_1256; 00127 centroid += point_0374; 00128 centroid /= 8.0; 00129 00130 int i; 00131 for( i = 0; i < 8; i++ ) 00132 position[i] -= centroid; 00133 00134 // ... Rotate element such that center of side 1-2 and 0-3 define X axis 00135 double DX = point_1256.x() - point_0374.x(); 00136 double DY = point_1256.y() - point_0374.y(); 00137 double DZ = point_1256.z() - point_0374.z(); 00138 00139 double AMAGX = sqrt( DX * DX + DZ * DZ ); 00140 double AMAGY = sqrt( DX * DX + DY * DY + DZ * DZ ); 00141 double FMAGX = AMAGX == 0.0 ? 1.0 : 0.0; 00142 double FMAGY = AMAGY == 0.0 ? 1.0 : 0.0; 00143 00144 double CZ = DX / ( AMAGX + FMAGX ) + FMAGX; 00145 double SZ = DZ / ( AMAGX + FMAGX ); 00146 double CY = AMAGX / ( AMAGY + FMAGY ) + FMAGY; 00147 double SY = DY / ( AMAGY + FMAGY ); 00148 00149 double temp; 00150 00151 for( i = 0; i < 8; i++ ) 00152 { 00153 temp = CY * CZ * position[i].x() + CY * SZ * position[i].z() + SY * position[i].y(); 00154 position[i].y( -SY * CZ * position[i].x() - SY * SZ * position[i].z() + CY * position[i].y() ); 00155 position[i].z( -SZ * position[i].x() + CZ * position[i].z() ); 00156 position[i].x( temp ); 00157 } 00158 00159 // ... Now, rotate about Y 00160 VerdictVector delta = -position[0]; 00161 delta -= position[1]; 00162 delta += position[2]; 00163 delta += position[3]; 00164 delta -= position[4]; 00165 delta -= position[5]; 00166 delta += position[6]; 00167 delta += position[7]; 00168 00169 DY = delta.y(); 00170 DZ = delta.z(); 00171 00172 AMAGY = sqrt( DY * DY + DZ * DZ ); 00173 FMAGY = AMAGY == 0.0 ? 1.0 : 0.0; 00174 00175 double CX = DY / ( AMAGY + FMAGY ) + FMAGY; 00176 double SX = DZ / ( AMAGY + FMAGY ); 00177 00178 for( i = 0; i < 8; i++ ) 00179 { 00180 temp = CX * position[i].y() + SX * position[i].z(); 00181 position[i].z( -SX * position[i].y() + CX * position[i].z() ); 00182 position[i].y( temp ); 00183 } 00184 } 00185 00186 double safe_ratio3( const double numerator, const double denominator, const double max_ratio ) 00187 { 00188 // this filter is essential for good running time in practice 00189 double return_value; 00190 00191 const double filter_n = max_ratio * 1.0e-16; 00192 const double filter_d = 1.0e-16; 00193 if( fabs( numerator ) <= filter_n && fabs( denominator ) >= filter_d ) { return_value = numerator / denominator; } 00194 else 00195 { 00196 return_value = fabs( numerator ) / max_ratio >= fabs( denominator ) 00197 ? ( ( numerator >= 0.0 && denominator >= 0.0 ) || ( numerator < 0.0 && denominator < 0.0 ) 00198 ? max_ratio 00199 : -max_ratio ) 00200 : numerator / denominator; 00201 } 00202 00203 return return_value; 00204 } 00205 00206 double safe_ratio( const double numerator, const double denominator ) 00207 { 00208 00209 double return_value; 00210 const double filter_n = VERDICT_DBL_MAX; 00211 const double filter_d = VERDICT_DBL_MIN; 00212 if( fabs( numerator ) <= filter_n && fabs( denominator ) >= filter_d ) { return_value = numerator / denominator; } 00213 else 00214 { 00215 return_value = VERDICT_DBL_MAX; 00216 } 00217 00218 return return_value; 00219 } 00220 00221 double condition_comp( const VerdictVector& xxi, const VerdictVector& xet, const VerdictVector& xze ) 00222 { 00223 double det = xxi % ( xet * xze ); 00224 00225 if( det <= VERDICT_DBL_MIN ) { return VERDICT_DBL_MAX; } 00226 00227 double term1 = xxi % xxi + xet % xet + xze % xze; 00228 double term2 = 00229 ( ( xxi * xet ) % ( xxi * xet ) ) + ( ( xet * xze ) % ( xet * xze ) ) + ( ( xze * xxi ) % ( xze * xxi ) ); 00230 00231 return sqrt( term1 * term2 ) / det; 00232 } 00233 00234 double oddy_comp( const VerdictVector& xxi, const VerdictVector& xet, const VerdictVector& xze ) 00235 { 00236 static const double third = 1.0 / 3.0; 00237 00238 double g11, g12, g13, g22, g23, g33, rt_g; 00239 00240 g11 = xxi % xxi; 00241 g12 = xxi % xet; 00242 g13 = xxi % xze; 00243 g22 = xet % xet; 00244 g23 = xet % xze; 00245 g33 = xze % xze; 00246 rt_g = xxi % ( xet * xze ); 00247 00248 double oddy_metric; 00249 if( rt_g > VERDICT_DBL_MIN ) 00250 { 00251 double norm_G_squared = g11 * g11 + 2.0 * g12 * g12 + 2.0 * g13 * g13 + g22 * g22 + 2.0 * g23 * g23 + g33 * g33; 00252 00253 double norm_J_squared = g11 + g22 + g33; 00254 00255 oddy_metric = ( norm_G_squared - third * norm_J_squared * norm_J_squared ) / pow( rt_g, 4. * third ); 00256 } 00257 00258 else 00259 oddy_metric = VERDICT_DBL_MAX; 00260 00261 return oddy_metric; 00262 } 00263 00264 //! calcualates edge lengths of a hex 00265 double hex_edge_length( int max_min, double coordinates[][3] ) 00266 { 00267 double temp[3], edge[12]; 00268 int i; 00269 00270 // lengths^2 of edges 00271 for( i = 0; i < 3; i++ ) 00272 { 00273 temp[i] = coordinates[1][i] - coordinates[0][i]; 00274 temp[i] = temp[i] * temp[i]; 00275 } 00276 edge[0] = sqrt( temp[0] + temp[1] + temp[2] ); 00277 00278 for( i = 0; i < 3; i++ ) 00279 { 00280 temp[i] = coordinates[2][i] - coordinates[1][i]; 00281 temp[i] = temp[i] * temp[i]; 00282 } 00283 edge[1] = sqrt( temp[0] + temp[1] + temp[2] ); 00284 00285 for( i = 0; i < 3; i++ ) 00286 { 00287 temp[i] = coordinates[3][i] - coordinates[2][i]; 00288 temp[i] = temp[i] * temp[i]; 00289 } 00290 edge[2] = sqrt( temp[0] + temp[1] + temp[2] ); 00291 00292 for( i = 0; i < 3; i++ ) 00293 { 00294 temp[i] = coordinates[0][i] - coordinates[3][i]; 00295 temp[i] = temp[i] * temp[i]; 00296 } 00297 edge[3] = sqrt( temp[0] + temp[1] + temp[2] ); 00298 00299 for( i = 0; i < 3; i++ ) 00300 { 00301 temp[i] = coordinates[5][i] - coordinates[4][i]; 00302 temp[i] = temp[i] * temp[i]; 00303 } 00304 edge[4] = sqrt( temp[0] + temp[1] + temp[2] ); 00305 00306 for( i = 0; i < 3; i++ ) 00307 { 00308 temp[i] = coordinates[6][i] - coordinates[5][i]; 00309 temp[i] = temp[i] * temp[i]; 00310 } 00311 edge[5] = sqrt( temp[0] + temp[1] + temp[2] ); 00312 00313 for( i = 0; i < 3; i++ ) 00314 { 00315 temp[i] = coordinates[7][i] - coordinates[6][i]; 00316 temp[i] = temp[i] * temp[i]; 00317 } 00318 edge[6] = sqrt( temp[0] + temp[1] + temp[2] ); 00319 00320 for( i = 0; i < 3; i++ ) 00321 { 00322 temp[i] = coordinates[4][i] - coordinates[7][i]; 00323 temp[i] = temp[i] * temp[i]; 00324 } 00325 edge[7] = sqrt( temp[0] + temp[1] + temp[2] ); 00326 00327 for( i = 0; i < 3; i++ ) 00328 { 00329 temp[i] = coordinates[4][i] - coordinates[0][i]; 00330 temp[i] = temp[i] * temp[i]; 00331 } 00332 edge[8] = sqrt( temp[0] + temp[1] + temp[2] ); 00333 00334 for( i = 0; i < 3; i++ ) 00335 { 00336 temp[i] = coordinates[5][i] - coordinates[1][i]; 00337 temp[i] = temp[i] * temp[i]; 00338 } 00339 edge[9] = sqrt( temp[0] + temp[1] + temp[2] ); 00340 00341 for( i = 0; i < 3; i++ ) 00342 { 00343 temp[i] = coordinates[6][i] - coordinates[2][i]; 00344 temp[i] = temp[i] * temp[i]; 00345 } 00346 edge[10] = sqrt( temp[0] + temp[1] + temp[2] ); 00347 00348 for( i = 0; i < 3; i++ ) 00349 { 00350 temp[i] = coordinates[7][i] - coordinates[3][i]; 00351 temp[i] = temp[i] * temp[i]; 00352 } 00353 edge[11] = sqrt( temp[0] + temp[1] + temp[2] ); 00354 00355 double _edge = edge[0]; 00356 00357 if( max_min == 0 ) 00358 { 00359 for( i = 1; i < 12; i++ ) 00360 _edge = VERDICT_MIN( _edge, edge[i] ); 00361 return (double)_edge; 00362 } 00363 else 00364 { 00365 for( i = 1; i < 12; i++ ) 00366 _edge = VERDICT_MAX( _edge, edge[i] ); 00367 return (double)_edge; 00368 } 00369 } 00370 00371 double diag_length( int max_min, double coordinates[][3] ) 00372 { 00373 double temp[3], diag[4]; 00374 int i; 00375 00376 // lengths^2 f diag nals 00377 for( i = 0; i < 3; i++ ) 00378 { 00379 temp[i] = coordinates[6][i] - coordinates[0][i]; 00380 temp[i] = temp[i] * temp[i]; 00381 } 00382 diag[0] = sqrt( temp[0] + temp[1] + temp[2] ); 00383 00384 for( i = 0; i < 3; i++ ) 00385 { 00386 temp[i] = coordinates[4][i] - coordinates[2][i]; 00387 temp[i] = temp[i] * temp[i]; 00388 } 00389 diag[1] = sqrt( temp[0] + temp[1] + temp[2] ); 00390 00391 for( i = 0; i < 3; i++ ) 00392 { 00393 temp[i] = coordinates[7][i] - coordinates[1][i]; 00394 temp[i] = temp[i] * temp[i]; 00395 } 00396 diag[2] = sqrt( temp[0] + temp[1] + temp[2] ); 00397 00398 for( i = 0; i < 3; i++ ) 00399 { 00400 temp[i] = coordinates[5][i] - coordinates[3][i]; 00401 temp[i] = temp[i] * temp[i]; 00402 } 00403 diag[3] = sqrt( temp[0] + temp[1] + temp[2] ); 00404 00405 double diagonal = diag[0]; 00406 if( max_min == 0 ) // Return min diagonal 00407 { 00408 for( i = 1; i < 4; i++ ) 00409 diagonal = VERDICT_MIN( diagonal, diag[i] ); 00410 return (double)diagonal; 00411 } 00412 else // Return max diagonal 00413 { 00414 for( i = 1; i < 4; i++ ) 00415 diagonal = VERDICT_MAX( diagonal, diag[i] ); 00416 return (double)diagonal; 00417 } 00418 } 00419 00420 //! calculates efg values 00421 VerdictVector calc_hex_efg( int efg_index, VerdictVector coordinates[8] ) 00422 { 00423 00424 VerdictVector efg; 00425 00426 switch( efg_index ) 00427 { 00428 00429 case 1: 00430 efg = coordinates[1]; 00431 efg += coordinates[2]; 00432 efg += coordinates[5]; 00433 efg += coordinates[6]; 00434 efg -= coordinates[0]; 00435 efg -= coordinates[3]; 00436 efg -= coordinates[4]; 00437 efg -= coordinates[7]; 00438 break; 00439 00440 case 2: 00441 efg = coordinates[2]; 00442 efg += coordinates[3]; 00443 efg += coordinates[6]; 00444 efg += coordinates[7]; 00445 efg -= coordinates[0]; 00446 efg -= coordinates[1]; 00447 efg -= coordinates[4]; 00448 efg -= coordinates[5]; 00449 break; 00450 00451 case 3: 00452 efg = coordinates[4]; 00453 efg += coordinates[5]; 00454 efg += coordinates[6]; 00455 efg += coordinates[7]; 00456 efg -= coordinates[0]; 00457 efg -= coordinates[1]; 00458 efg -= coordinates[2]; 00459 efg -= coordinates[3]; 00460 break; 00461 00462 case 12: 00463 efg = coordinates[0]; 00464 efg += coordinates[2]; 00465 efg += coordinates[4]; 00466 efg += coordinates[6]; 00467 efg -= coordinates[1]; 00468 efg -= coordinates[3]; 00469 efg -= coordinates[5]; 00470 efg -= coordinates[7]; 00471 break; 00472 00473 case 13: 00474 efg = coordinates[0]; 00475 efg += coordinates[3]; 00476 efg += coordinates[5]; 00477 efg += coordinates[6]; 00478 efg -= coordinates[1]; 00479 efg -= coordinates[2]; 00480 efg -= coordinates[4]; 00481 efg -= coordinates[7]; 00482 break; 00483 00484 case 23: 00485 efg = coordinates[0]; 00486 efg += coordinates[1]; 00487 efg += coordinates[6]; 00488 efg += coordinates[7]; 00489 efg -= coordinates[2]; 00490 efg -= coordinates[3]; 00491 efg -= coordinates[4]; 00492 efg -= coordinates[5]; 00493 break; 00494 00495 case 123: 00496 efg = coordinates[0]; 00497 efg += coordinates[2]; 00498 efg += coordinates[5]; 00499 efg += coordinates[7]; 00500 efg -= coordinates[1]; 00501 efg -= coordinates[5]; 00502 efg -= coordinates[6]; 00503 efg -= coordinates[2]; 00504 break; 00505 00506 default: 00507 efg.set( 0, 0, 0 ); 00508 } 00509 00510 return efg; 00511 } 00512 00513 /*! 00514 the edge ratio of a hex 00515 00516 NB (P. Pebay 01/23/07): 00517 Hmax / Hmin where Hmax and Hmin are respectively the maximum and the 00518 minimum edge lengths 00519 */ 00520 C_FUNC_DEF double v_hex_edge_ratio( int /*num_nodes*/, double coordinates[][3] ) 00521 { 00522 00523 VerdictVector edges[12]; 00524 make_hex_edges( coordinates, edges ); 00525 00526 double a2 = edges[0].length_squared(); 00527 double b2 = edges[1].length_squared(); 00528 double c2 = edges[2].length_squared(); 00529 double d2 = edges[3].length_squared(); 00530 double e2 = edges[4].length_squared(); 00531 double f2 = edges[5].length_squared(); 00532 double g2 = edges[6].length_squared(); 00533 double h2 = edges[7].length_squared(); 00534 double i2 = edges[8].length_squared(); 00535 double j2 = edges[9].length_squared(); 00536 double k2 = edges[10].length_squared(); 00537 double l2 = edges[11].length_squared(); 00538 00539 double mab, mcd, mef, Mab, Mcd, Mef; 00540 double mgh, mij, mkl, Mgh, Mij, Mkl; 00541 00542 if( a2 < b2 ) 00543 { 00544 mab = a2; 00545 Mab = b2; 00546 } 00547 else // b2 <= a2 00548 { 00549 mab = b2; 00550 Mab = a2; 00551 } 00552 if( c2 < d2 ) 00553 { 00554 mcd = c2; 00555 Mcd = d2; 00556 } 00557 else // d2 <= c2 00558 { 00559 mcd = d2; 00560 Mcd = c2; 00561 } 00562 if( e2 < f2 ) 00563 { 00564 mef = e2; 00565 Mef = f2; 00566 } 00567 else // f2 <= e2 00568 { 00569 mef = f2; 00570 Mef = e2; 00571 } 00572 if( g2 < h2 ) 00573 { 00574 mgh = g2; 00575 Mgh = h2; 00576 } 00577 else // h2 <= g2 00578 { 00579 mgh = h2; 00580 Mgh = g2; 00581 } 00582 if( i2 < j2 ) 00583 { 00584 mij = i2; 00585 Mij = j2; 00586 } 00587 else // j2 <= i2 00588 { 00589 mij = j2; 00590 Mij = i2; 00591 } 00592 if( k2 < l2 ) 00593 { 00594 mkl = k2; 00595 Mkl = l2; 00596 } 00597 else // l2 <= k2 00598 { 00599 mkl = l2; 00600 Mkl = k2; 00601 } 00602 00603 double m2; 00604 m2 = mab < mcd ? mab : mcd; 00605 m2 = m2 < mef ? m2 : mef; 00606 m2 = m2 < mgh ? m2 : mgh; 00607 m2 = m2 < mij ? m2 : mij; 00608 m2 = m2 < mkl ? m2 : mkl; 00609 00610 if( m2 < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX; 00611 00612 double M2; 00613 M2 = Mab > Mcd ? Mab : Mcd; 00614 M2 = M2 > Mef ? M2 : Mef; 00615 M2 = M2 > Mgh ? M2 : Mgh; 00616 M2 = M2 > Mij ? M2 : Mij; 00617 M2 = M2 > Mkl ? M2 : Mkl; 00618 m2 = m2 < mef ? m2 : mef; 00619 00620 M2 = Mab > Mcd ? Mab : Mcd; 00621 M2 = M2 > Mef ? M2 : Mef; 00622 00623 double edge_ratio = sqrt( M2 / m2 ); 00624 00625 if( edge_ratio > 0 ) return (double)VERDICT_MIN( edge_ratio, VERDICT_DBL_MAX ); 00626 return (double)VERDICT_MAX( edge_ratio, -VERDICT_DBL_MAX ); 00627 } 00628 00629 /*! 00630 max edge ratio of a hex 00631 00632 Maximum edge length ratio at hex center 00633 */ 00634 C_FUNC_DEF double v_hex_max_edge_ratio( int /*num_nodes*/, double coordinates[][3] ) 00635 { 00636 double aspect; 00637 VerdictVector node_pos[8]; 00638 make_hex_nodes( coordinates, node_pos ); 00639 00640 double aspect_12, aspect_13, aspect_23; 00641 00642 VerdictVector efg1 = calc_hex_efg( 1, node_pos ); 00643 VerdictVector efg2 = calc_hex_efg( 2, node_pos ); 00644 VerdictVector efg3 = calc_hex_efg( 3, node_pos ); 00645 00646 double mag_efg1 = efg1.length(); 00647 double mag_efg2 = efg2.length(); 00648 double mag_efg3 = efg3.length(); 00649 00650 aspect_12 = safe_ratio( VERDICT_MAX( mag_efg1, mag_efg2 ), VERDICT_MIN( mag_efg1, mag_efg2 ) ); 00651 aspect_13 = safe_ratio( VERDICT_MAX( mag_efg1, mag_efg3 ), VERDICT_MIN( mag_efg1, mag_efg3 ) ); 00652 aspect_23 = safe_ratio( VERDICT_MAX( mag_efg2, mag_efg3 ), VERDICT_MIN( mag_efg2, mag_efg3 ) ); 00653 00654 aspect = VERDICT_MAX( aspect_12, VERDICT_MAX( aspect_13, aspect_23 ) ); 00655 00656 if( aspect > 0 ) return (double)VERDICT_MIN( aspect, VERDICT_DBL_MAX ); 00657 return (double)VERDICT_MAX( aspect, -VERDICT_DBL_MAX ); 00658 } 00659 00660 /*! 00661 skew of a hex 00662 00663 Maximum ||cosA|| where A is the angle between edges at hex center. 00664 */ 00665 C_FUNC_DEF double v_hex_skew( int /*num_nodes*/, double coordinates[][3] ) 00666 { 00667 VerdictVector node_pos[8]; 00668 make_hex_nodes( coordinates, node_pos ); 00669 00670 double skew_1, skew_2, skew_3; 00671 00672 VerdictVector efg1 = calc_hex_efg( 1, node_pos ); 00673 VerdictVector efg2 = calc_hex_efg( 2, node_pos ); 00674 VerdictVector efg3 = calc_hex_efg( 3, node_pos ); 00675 00676 if( efg1.normalize() <= VERDICT_DBL_MIN ) return VERDICT_DBL_MAX; 00677 if( efg2.normalize() <= VERDICT_DBL_MIN ) return VERDICT_DBL_MAX; 00678 if( efg3.normalize() <= VERDICT_DBL_MIN ) return VERDICT_DBL_MAX; 00679 00680 skew_1 = fabs( efg1 % efg2 ); 00681 skew_2 = fabs( efg1 % efg3 ); 00682 skew_3 = fabs( efg2 % efg3 ); 00683 00684 double skew = ( VERDICT_MAX( skew_1, VERDICT_MAX( skew_2, skew_3 ) ) ); 00685 00686 if( skew > 0 ) return (double)VERDICT_MIN( skew, VERDICT_DBL_MAX ); 00687 return (double)VERDICT_MAX( skew, -VERDICT_DBL_MAX ); 00688 } 00689 00690 /*! 00691 taper of a hex 00692 00693 Maximum ratio of lengths derived from opposite edges. 00694 */ 00695 C_FUNC_DEF double v_hex_taper( int /*num_nodes*/, double coordinates[][3] ) 00696 { 00697 VerdictVector node_pos[8]; 00698 make_hex_nodes( coordinates, node_pos ); 00699 00700 VerdictVector efg1 = calc_hex_efg( 1, node_pos ); 00701 VerdictVector efg2 = calc_hex_efg( 2, node_pos ); 00702 VerdictVector efg3 = calc_hex_efg( 3, node_pos ); 00703 00704 VerdictVector efg12 = calc_hex_efg( 12, node_pos ); 00705 VerdictVector efg13 = calc_hex_efg( 13, node_pos ); 00706 VerdictVector efg23 = calc_hex_efg( 23, node_pos ); 00707 00708 double taper_1 = fabs( safe_ratio( efg12.length(), VERDICT_MIN( efg1.length(), efg2.length() ) ) ); 00709 double taper_2 = fabs( safe_ratio( efg13.length(), VERDICT_MIN( efg1.length(), efg3.length() ) ) ); 00710 double taper_3 = fabs( safe_ratio( efg23.length(), VERDICT_MIN( efg2.length(), efg3.length() ) ) ); 00711 00712 double taper = (double)VERDICT_MAX( taper_1, VERDICT_MAX( taper_2, taper_3 ) ); 00713 00714 if( taper > 0 ) return (double)VERDICT_MIN( taper, VERDICT_DBL_MAX ); 00715 return (double)VERDICT_MAX( taper, -VERDICT_DBL_MAX ); 00716 } 00717 00718 /*! 00719 volume of a hex 00720 00721 Jacobian at hex center 00722 */ 00723 C_FUNC_DEF double v_hex_volume( int /*num_nodes*/, double coordinates[][3] ) 00724 { 00725 VerdictVector node_pos[8]; 00726 make_hex_nodes( coordinates, node_pos ); 00727 00728 VerdictVector efg1 = calc_hex_efg( 1, node_pos ); 00729 VerdictVector efg2 = calc_hex_efg( 2, node_pos ); 00730 VerdictVector efg3 = calc_hex_efg( 3, node_pos ); 00731 00732 double volume; 00733 volume = (double)( efg1 % ( efg2 * efg3 ) ) / 64.0; 00734 00735 if( volume > 0 ) return (double)VERDICT_MIN( volume, VERDICT_DBL_MAX ); 00736 return (double)VERDICT_MAX( volume, -VERDICT_DBL_MAX ); 00737 } 00738 00739 /*! 00740 stretch of a hex 00741 00742 sqrt(3) * minimum edge length / maximum diagonal length 00743 */ 00744 C_FUNC_DEF double v_hex_stretch( int /*num_nodes*/, double coordinates[][3] ) 00745 { 00746 static const double HEX_STRETCH_SCALE_FACTOR = sqrt( 3.0 ); 00747 00748 double min_edge = hex_edge_length( 0, coordinates ); 00749 double max_diag = diag_length( 1, coordinates ); 00750 00751 double stretch = HEX_STRETCH_SCALE_FACTOR * safe_ratio( min_edge, max_diag ); 00752 00753 if( stretch > 0 ) return (double)VERDICT_MIN( stretch, VERDICT_DBL_MAX ); 00754 return (double)VERDICT_MAX( stretch, -VERDICT_DBL_MAX ); 00755 } 00756 00757 /*! 00758 diagonal ratio of a hex 00759 00760 Minimum diagonal length / maximum diagonal length 00761 */ 00762 C_FUNC_DEF double v_hex_diagonal( int /*num_nodes*/, double coordinates[][3] ) 00763 { 00764 00765 double min_diag = diag_length( 0, coordinates ); 00766 double max_diag = diag_length( 1, coordinates ); 00767 00768 double diagonal = safe_ratio( min_diag, max_diag ); 00769 00770 if( diagonal > 0 ) return (double)VERDICT_MIN( diagonal, VERDICT_DBL_MAX ); 00771 return (double)VERDICT_MAX( diagonal, -VERDICT_DBL_MAX ); 00772 } 00773 00774 #define SQR( x ) ( ( x ) * ( x ) ) 00775 00776 /*! 00777 dimension of a hex 00778 00779 Pronto-specific characteristic length for stable time step calculation. 00780 Char_length = Volume / 2 grad Volume 00781 */ 00782 C_FUNC_DEF double v_hex_dimension( int /*num_nodes*/, double coordinates[][3] ) 00783 { 00784 00785 double gradop[9][4]; 00786 00787 double x1 = coordinates[0][0]; 00788 double x2 = coordinates[1][0]; 00789 double x3 = coordinates[2][0]; 00790 double x4 = coordinates[3][0]; 00791 double x5 = coordinates[4][0]; 00792 double x6 = coordinates[5][0]; 00793 double x7 = coordinates[6][0]; 00794 double x8 = coordinates[7][0]; 00795 00796 double y1 = coordinates[0][1]; 00797 double y2 = coordinates[1][1]; 00798 double y3 = coordinates[2][1]; 00799 double y4 = coordinates[3][1]; 00800 double y5 = coordinates[4][1]; 00801 double y6 = coordinates[5][1]; 00802 double y7 = coordinates[6][1]; 00803 double y8 = coordinates[7][1]; 00804 00805 double z1 = coordinates[0][2]; 00806 double z2 = coordinates[1][2]; 00807 double z3 = coordinates[2][2]; 00808 double z4 = coordinates[3][2]; 00809 double z5 = coordinates[4][2]; 00810 double z6 = coordinates[5][2]; 00811 double z7 = coordinates[6][2]; 00812 double z8 = coordinates[7][2]; 00813 00814 double z24 = z2 - z4; 00815 double z52 = z5 - z2; 00816 double z45 = z4 - z5; 00817 gradop[1][1] = 00818 ( y2 * ( z6 - z3 - z45 ) + y3 * z24 + y4 * ( z3 - z8 - z52 ) + y5 * ( z8 - z6 - z24 ) + y6 * z52 + y8 * z45 ) / 00819 12.0; 00820 00821 double z31 = z3 - z1; 00822 double z63 = z6 - z3; 00823 double z16 = z1 - z6; 00824 gradop[2][1] = 00825 ( y3 * ( z7 - z4 - z16 ) + y4 * z31 + y1 * ( z4 - z5 - z63 ) + y6 * ( z5 - z7 - z31 ) + y7 * z63 + y5 * z16 ) / 00826 12.0; 00827 00828 double z42 = z4 - z2; 00829 double z74 = z7 - z4; 00830 double z27 = z2 - z7; 00831 gradop[3][1] = 00832 ( y4 * ( z8 - z1 - z27 ) + y1 * z42 + y2 * ( z1 - z6 - z74 ) + y7 * ( z6 - z8 - z42 ) + y8 * z74 + y6 * z27 ) / 00833 12.0; 00834 00835 double z13 = z1 - z3; 00836 double z81 = z8 - z1; 00837 double z38 = z3 - z8; 00838 gradop[4][1] = 00839 ( y1 * ( z5 - z2 - z38 ) + y2 * z13 + y3 * ( z2 - z7 - z81 ) + y8 * ( z7 - z5 - z13 ) + y5 * z81 + y7 * z38 ) / 00840 12.0; 00841 00842 double z86 = z8 - z6; 00843 double z18 = z1 - z8; 00844 double z61 = z6 - z1; 00845 gradop[5][1] = 00846 ( y8 * ( z4 - z7 - z61 ) + y7 * z86 + y6 * ( z7 - z2 - z18 ) + y1 * ( z2 - z4 - z86 ) + y4 * z18 + y2 * z61 ) / 00847 12.0; 00848 00849 double z57 = z5 - z7; 00850 double z25 = z2 - z5; 00851 double z72 = z7 - z2; 00852 gradop[6][1] = 00853 ( y5 * ( z1 - z8 - z72 ) + y8 * z57 + y7 * ( z8 - z3 - z25 ) + y2 * ( z3 - z1 - z57 ) + y1 * z25 + y3 * z72 ) / 00854 12.0; 00855 00856 double z68 = z6 - z8; 00857 double z36 = z3 - z6; 00858 double z83 = z8 - z3; 00859 gradop[7][1] = 00860 ( y6 * ( z2 - z5 - z83 ) + y5 * z68 + y8 * ( z5 - z4 - z36 ) + y3 * ( z4 - z2 - z68 ) + y2 * z36 + y4 * z83 ) / 00861 12.0; 00862 00863 double z75 = z7 - z5; 00864 double z47 = z4 - z7; 00865 double z54 = z5 - z4; 00866 gradop[8][1] = 00867 ( y7 * ( z3 - z6 - z54 ) + y6 * z75 + y5 * ( z6 - z1 - z47 ) + y4 * ( z1 - z3 - z75 ) + y3 * z47 + y1 * z54 ) / 00868 12.0; 00869 00870 double x24 = x2 - x4; 00871 double x52 = x5 - x2; 00872 double x45 = x4 - x5; 00873 gradop[1][2] = 00874 ( z2 * ( x6 - x3 - x45 ) + z3 * x24 + z4 * ( x3 - x8 - x52 ) + z5 * ( x8 - x6 - x24 ) + z6 * x52 + z8 * x45 ) / 00875 12.0; 00876 00877 double x31 = x3 - x1; 00878 double x63 = x6 - x3; 00879 double x16 = x1 - x6; 00880 gradop[2][2] = 00881 ( z3 * ( x7 - x4 - x16 ) + z4 * x31 + z1 * ( x4 - x5 - x63 ) + z6 * ( x5 - x7 - x31 ) + z7 * x63 + z5 * x16 ) / 00882 12.0; 00883 00884 double x42 = x4 - x2; 00885 double x74 = x7 - x4; 00886 double x27 = x2 - x7; 00887 gradop[3][2] = 00888 ( z4 * ( x8 - x1 - x27 ) + z1 * x42 + z2 * ( x1 - x6 - x74 ) + z7 * ( x6 - x8 - x42 ) + z8 * x74 + z6 * x27 ) / 00889 12.0; 00890 00891 double x13 = x1 - x3; 00892 double x81 = x8 - x1; 00893 double x38 = x3 - x8; 00894 gradop[4][2] = 00895 ( z1 * ( x5 - x2 - x38 ) + z2 * x13 + z3 * ( x2 - x7 - x81 ) + z8 * ( x7 - x5 - x13 ) + z5 * x81 + z7 * x38 ) / 00896 12.0; 00897 00898 double x86 = x8 - x6; 00899 double x18 = x1 - x8; 00900 double x61 = x6 - x1; 00901 gradop[5][2] = 00902 ( z8 * ( x4 - x7 - x61 ) + z7 * x86 + z6 * ( x7 - x2 - x18 ) + z1 * ( x2 - x4 - x86 ) + z4 * x18 + z2 * x61 ) / 00903 12.0; 00904 00905 double x57 = x5 - x7; 00906 double x25 = x2 - x5; 00907 double x72 = x7 - x2; 00908 gradop[6][2] = 00909 ( z5 * ( x1 - x8 - x72 ) + z8 * x57 + z7 * ( x8 - x3 - x25 ) + z2 * ( x3 - x1 - x57 ) + z1 * x25 + z3 * x72 ) / 00910 12.0; 00911 00912 double x68 = x6 - x8; 00913 double x36 = x3 - x6; 00914 double x83 = x8 - x3; 00915 gradop[7][2] = 00916 ( z6 * ( x2 - x5 - x83 ) + z5 * x68 + z8 * ( x5 - x4 - x36 ) + z3 * ( x4 - x2 - x68 ) + z2 * x36 + z4 * x83 ) / 00917 12.0; 00918 00919 double x75 = x7 - x5; 00920 double x47 = x4 - x7; 00921 double x54 = x5 - x4; 00922 gradop[8][2] = 00923 ( z7 * ( x3 - x6 - x54 ) + z6 * x75 + z5 * ( x6 - x1 - x47 ) + z4 * ( x1 - x3 - x75 ) + z3 * x47 + z1 * x54 ) / 00924 12.0; 00925 00926 double y24 = y2 - y4; 00927 double y52 = y5 - y2; 00928 double y45 = y4 - y5; 00929 gradop[1][3] = 00930 ( x2 * ( y6 - y3 - y45 ) + x3 * y24 + x4 * ( y3 - y8 - y52 ) + x5 * ( y8 - y6 - y24 ) + x6 * y52 + x8 * y45 ) / 00931 12.0; 00932 00933 double y31 = y3 - y1; 00934 double y63 = y6 - y3; 00935 double y16 = y1 - y6; 00936 gradop[2][3] = 00937 ( x3 * ( y7 - y4 - y16 ) + x4 * y31 + x1 * ( y4 - y5 - y63 ) + x6 * ( y5 - y7 - y31 ) + x7 * y63 + x5 * y16 ) / 00938 12.0; 00939 00940 double y42 = y4 - y2; 00941 double y74 = y7 - y4; 00942 double y27 = y2 - y7; 00943 gradop[3][3] = 00944 ( x4 * ( y8 - y1 - y27 ) + x1 * y42 + x2 * ( y1 - y6 - y74 ) + x7 * ( y6 - y8 - y42 ) + x8 * y74 + x6 * y27 ) / 00945 12.0; 00946 00947 double y13 = y1 - y3; 00948 double y81 = y8 - y1; 00949 double y38 = y3 - y8; 00950 gradop[4][3] = 00951 ( x1 * ( y5 - y2 - y38 ) + x2 * y13 + x3 * ( y2 - y7 - y81 ) + x8 * ( y7 - y5 - y13 ) + x5 * y81 + x7 * y38 ) / 00952 12.0; 00953 00954 double y86 = y8 - y6; 00955 double y18 = y1 - y8; 00956 double y61 = y6 - y1; 00957 gradop[5][3] = 00958 ( x8 * ( y4 - y7 - y61 ) + x7 * y86 + x6 * ( y7 - y2 - y18 ) + x1 * ( y2 - y4 - y86 ) + x4 * y18 + x2 * y61 ) / 00959 12.0; 00960 00961 double y57 = y5 - y7; 00962 double y25 = y2 - y5; 00963 double y72 = y7 - y2; 00964 gradop[6][3] = 00965 ( x5 * ( y1 - y8 - y72 ) + x8 * y57 + x7 * ( y8 - y3 - y25 ) + x2 * ( y3 - y1 - y57 ) + x1 * y25 + x3 * y72 ) / 00966 12.0; 00967 00968 double y68 = y6 - y8; 00969 double y36 = y3 - y6; 00970 double y83 = y8 - y3; 00971 gradop[7][3] = 00972 ( x6 * ( y2 - y5 - y83 ) + x5 * y68 + x8 * ( y5 - y4 - y36 ) + x3 * ( y4 - y2 - y68 ) + x2 * y36 + x4 * y83 ) / 00973 12.0; 00974 00975 double y75 = y7 - y5; 00976 double y47 = y4 - y7; 00977 double y54 = y5 - y4; 00978 gradop[8][3] = 00979 ( x7 * ( y3 - y6 - y54 ) + x6 * y75 + x5 * ( y6 - y1 - y47 ) + x4 * ( y1 - y3 - y75 ) + x3 * y47 + x1 * y54 ) / 00980 12.0; 00981 00982 // calculate element volume and characteristic element aspect ratio 00983 // (used in time step and hourglass control) - 00984 00985 double volume = coordinates[0][0] * gradop[1][1] + coordinates[1][0] * gradop[2][1] + 00986 coordinates[2][0] * gradop[3][1] + coordinates[3][0] * gradop[4][1] + 00987 coordinates[4][0] * gradop[5][1] + coordinates[5][0] * gradop[6][1] + 00988 coordinates[6][0] * gradop[7][1] + coordinates[7][0] * gradop[8][1]; 00989 double aspect = 00990 .5 * SQR( volume ) / 00991 ( SQR( gradop[1][1] ) + SQR( gradop[2][1] ) + SQR( gradop[3][1] ) + SQR( gradop[4][1] ) + SQR( gradop[5][1] ) + 00992 SQR( gradop[6][1] ) + SQR( gradop[7][1] ) + SQR( gradop[8][1] ) + SQR( gradop[1][2] ) + SQR( gradop[2][2] ) + 00993 SQR( gradop[3][2] ) + SQR( gradop[4][2] ) + SQR( gradop[5][2] ) + SQR( gradop[6][2] ) + SQR( gradop[7][2] ) + 00994 SQR( gradop[8][2] ) + SQR( gradop[1][3] ) + SQR( gradop[2][3] ) + SQR( gradop[3][3] ) + SQR( gradop[4][3] ) + 00995 SQR( gradop[5][3] ) + SQR( gradop[6][3] ) + SQR( gradop[7][3] ) + SQR( gradop[8][3] ) ); 00996 00997 return (double)sqrt( aspect ); 00998 } 00999 01000 /*! 01001 oddy of a hex 01002 01003 General distortion measure based on left Cauchy-Green Tensor 01004 */ 01005 C_FUNC_DEF double v_hex_oddy( int /*num_nodes*/, double coordinates[][3] ) 01006 { 01007 01008 double oddy = 0.0, current_oddy; 01009 VerdictVector xxi, xet, xze; 01010 01011 VerdictVector node_pos[8]; 01012 make_hex_nodes( coordinates, node_pos ); 01013 01014 xxi = calc_hex_efg( 1, node_pos ); 01015 xet = calc_hex_efg( 2, node_pos ); 01016 xze = calc_hex_efg( 3, node_pos ); 01017 01018 current_oddy = oddy_comp( xxi, xet, xze ); 01019 if( current_oddy > oddy ) { oddy = current_oddy; } 01020 01021 xxi.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1], 01022 coordinates[1][2] - coordinates[0][2] ); 01023 01024 xet.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1], 01025 coordinates[3][2] - coordinates[0][2] ); 01026 01027 xze.set( coordinates[4][0] - coordinates[0][0], coordinates[4][1] - coordinates[0][1], 01028 coordinates[4][2] - coordinates[0][2] ); 01029 01030 current_oddy = oddy_comp( xxi, xet, xze ); 01031 if( current_oddy > oddy ) { oddy = current_oddy; } 01032 01033 xxi.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1], 01034 coordinates[2][2] - coordinates[1][2] ); 01035 01036 xet.set( coordinates[0][0] - coordinates[1][0], coordinates[0][1] - coordinates[1][1], 01037 coordinates[0][2] - coordinates[1][2] ); 01038 01039 xze.set( coordinates[5][0] - coordinates[1][0], coordinates[5][1] - coordinates[1][1], 01040 coordinates[5][2] - coordinates[1][2] ); 01041 01042 current_oddy = oddy_comp( xxi, xet, xze ); 01043 if( current_oddy > oddy ) { oddy = current_oddy; } 01044 01045 xxi.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1], 01046 coordinates[3][2] - coordinates[2][2] ); 01047 01048 xet.set( coordinates[1][0] - coordinates[2][0], coordinates[1][1] - coordinates[2][1], 01049 coordinates[1][2] - coordinates[2][2] ); 01050 01051 xze.set( coordinates[6][0] - coordinates[2][0], coordinates[6][1] - coordinates[2][1], 01052 coordinates[6][2] - coordinates[2][2] ); 01053 01054 current_oddy = oddy_comp( xxi, xet, xze ); 01055 if( current_oddy > oddy ) { oddy = current_oddy; } 01056 01057 xxi.set( coordinates[0][0] - coordinates[3][0], coordinates[0][1] - coordinates[3][1], 01058 coordinates[0][2] - coordinates[3][2] ); 01059 01060 xet.set( coordinates[2][0] - coordinates[3][0], coordinates[2][1] - coordinates[3][1], 01061 coordinates[2][2] - coordinates[3][2] ); 01062 01063 xze.set( coordinates[7][0] - coordinates[3][0], coordinates[7][1] - coordinates[3][1], 01064 coordinates[7][2] - coordinates[3][2] ); 01065 01066 current_oddy = oddy_comp( xxi, xet, xze ); 01067 if( current_oddy > oddy ) { oddy = current_oddy; } 01068 01069 xxi.set( coordinates[7][0] - coordinates[4][0], coordinates[7][1] - coordinates[4][1], 01070 coordinates[7][2] - coordinates[4][2] ); 01071 01072 xet.set( coordinates[5][0] - coordinates[4][0], coordinates[5][1] - coordinates[4][1], 01073 coordinates[5][2] - coordinates[4][2] ); 01074 01075 xze.set( coordinates[0][0] - coordinates[4][0], coordinates[0][1] - coordinates[4][1], 01076 coordinates[0][2] - coordinates[4][2] ); 01077 01078 current_oddy = oddy_comp( xxi, xet, xze ); 01079 if( current_oddy > oddy ) { oddy = current_oddy; } 01080 01081 xxi.set( coordinates[4][0] - coordinates[5][0], coordinates[4][1] - coordinates[5][1], 01082 coordinates[4][2] - coordinates[5][2] ); 01083 01084 xet.set( coordinates[6][0] - coordinates[5][0], coordinates[6][1] - coordinates[5][1], 01085 coordinates[6][2] - coordinates[5][2] ); 01086 01087 xze.set( coordinates[1][0] - coordinates[5][0], coordinates[1][1] - coordinates[5][1], 01088 coordinates[1][2] - coordinates[5][2] ); 01089 01090 current_oddy = oddy_comp( xxi, xet, xze ); 01091 if( current_oddy > oddy ) { oddy = current_oddy; } 01092 01093 xxi.set( coordinates[5][0] - coordinates[6][0], coordinates[5][1] - coordinates[6][1], 01094 coordinates[5][2] - coordinates[6][2] ); 01095 01096 xet.set( coordinates[7][0] - coordinates[6][0], coordinates[7][1] - coordinates[6][1], 01097 coordinates[7][2] - coordinates[6][2] ); 01098 01099 xze.set( coordinates[2][0] - coordinates[6][0], coordinates[2][1] - coordinates[6][1], 01100 coordinates[2][2] - coordinates[6][2] ); 01101 01102 current_oddy = oddy_comp( xxi, xet, xze ); 01103 if( current_oddy > oddy ) { oddy = current_oddy; } 01104 01105 xxi.set( coordinates[6][0] - coordinates[7][0], coordinates[6][1] - coordinates[7][1], 01106 coordinates[6][2] - coordinates[7][2] ); 01107 01108 xet.set( coordinates[4][0] - coordinates[7][0], coordinates[4][1] - coordinates[7][1], 01109 coordinates[4][2] - coordinates[7][2] ); 01110 01111 xze.set( coordinates[3][0] - coordinates[7][0], coordinates[3][1] - coordinates[7][1], 01112 coordinates[3][2] - coordinates[7][2] ); 01113 01114 current_oddy = oddy_comp( xxi, xet, xze ); 01115 if( current_oddy > oddy ) { oddy = current_oddy; } 01116 01117 if( oddy > 0 ) return (double)VERDICT_MIN( oddy, VERDICT_DBL_MAX ); 01118 return (double)VERDICT_MAX( oddy, -VERDICT_DBL_MAX ); 01119 } 01120 01121 /*! 01122 the average Frobenius aspect of a hex 01123 01124 NB (P. Pebay 01/20/07): 01125 this metric is calculated by averaging the 8 Frobenius aspects at 01126 each corner of the hex, when the reference corner is right isosceles. 01127 */ 01128 C_FUNC_DEF double v_hex_med_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] ) 01129 { 01130 01131 VerdictVector node_pos[8]; 01132 make_hex_nodes( coordinates, node_pos ); 01133 01134 double med_aspect_frobenius = 0.; 01135 VerdictVector xxi, xet, xze; 01136 01137 // J(0,0,0): 01138 01139 xxi = node_pos[1] - node_pos[0]; 01140 xet = node_pos[3] - node_pos[0]; 01141 xze = node_pos[4] - node_pos[0]; 01142 01143 med_aspect_frobenius += condition_comp( xxi, xet, xze ); 01144 // J(1,0,0): 01145 01146 xxi = node_pos[2] - node_pos[1]; 01147 xet = node_pos[0] - node_pos[1]; 01148 xze = node_pos[5] - node_pos[1]; 01149 01150 med_aspect_frobenius += condition_comp( xxi, xet, xze ); 01151 // J(1,1,0): 01152 01153 xxi = node_pos[3] - node_pos[2]; 01154 xet = node_pos[1] - node_pos[2]; 01155 xze = node_pos[6] - node_pos[2]; 01156 01157 med_aspect_frobenius += condition_comp( xxi, xet, xze ); 01158 // J(0,1,0): 01159 01160 xxi = node_pos[0] - node_pos[3]; 01161 xet = node_pos[2] - node_pos[3]; 01162 xze = node_pos[7] - node_pos[3]; 01163 01164 med_aspect_frobenius += condition_comp( xxi, xet, xze ); 01165 // J(0,0,1): 01166 01167 xxi = node_pos[7] - node_pos[4]; 01168 xet = node_pos[5] - node_pos[4]; 01169 xze = node_pos[0] - node_pos[4]; 01170 01171 med_aspect_frobenius += condition_comp( xxi, xet, xze ); 01172 // J(1,0,1): 01173 01174 xxi = node_pos[4] - node_pos[5]; 01175 xet = node_pos[6] - node_pos[5]; 01176 xze = node_pos[1] - node_pos[5]; 01177 01178 med_aspect_frobenius += condition_comp( xxi, xet, xze ); 01179 // J(1,1,1): 01180 01181 xxi = node_pos[5] - node_pos[6]; 01182 xet = node_pos[7] - node_pos[6]; 01183 xze = node_pos[2] - node_pos[6]; 01184 01185 med_aspect_frobenius += condition_comp( xxi, xet, xze ); 01186 // J(1,1,1): 01187 01188 xxi = node_pos[6] - node_pos[7]; 01189 xet = node_pos[4] - node_pos[7]; 01190 xze = node_pos[3] - node_pos[7]; 01191 01192 med_aspect_frobenius += condition_comp( xxi, xet, xze ); 01193 med_aspect_frobenius /= 24.; 01194 01195 if( med_aspect_frobenius > 0 ) return (double)VERDICT_MIN( med_aspect_frobenius, VERDICT_DBL_MAX ); 01196 return (double)VERDICT_MAX( med_aspect_frobenius, -VERDICT_DBL_MAX ); 01197 } 01198 01199 /*! 01200 maximum Frobenius condition number of a hex 01201 01202 Maximum Frobenius condition number of the Jacobian matrix at 8 corners 01203 NB (P. Pebay 01/25/07): 01204 this metric is calculated by taking the maximum of the 8 Frobenius aspects at 01205 each corner of the hex, when the reference corner is right isosceles. 01206 */ 01207 C_FUNC_DEF double v_hex_max_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] ) 01208 { 01209 01210 VerdictVector node_pos[8]; 01211 make_hex_nodes( coordinates, node_pos ); 01212 01213 double condition = 0.0, current_condition; 01214 VerdictVector xxi, xet, xze; 01215 01216 xxi = calc_hex_efg( 1, node_pos ); 01217 xet = calc_hex_efg( 2, node_pos ); 01218 xze = calc_hex_efg( 3, node_pos ); 01219 01220 current_condition = condition_comp( xxi, xet, xze ); 01221 if( current_condition > condition ) { condition = current_condition; } 01222 01223 // J(0,0,0): 01224 01225 xxi = node_pos[1] - node_pos[0]; 01226 xet = node_pos[3] - node_pos[0]; 01227 xze = node_pos[4] - node_pos[0]; 01228 01229 current_condition = condition_comp( xxi, xet, xze ); 01230 if( current_condition > condition ) { condition = current_condition; } 01231 01232 // J(1,0,0): 01233 01234 xxi = node_pos[2] - node_pos[1]; 01235 xet = node_pos[0] - node_pos[1]; 01236 xze = node_pos[5] - node_pos[1]; 01237 01238 current_condition = condition_comp( xxi, xet, xze ); 01239 if( current_condition > condition ) { condition = current_condition; } 01240 01241 // J(1,1,0): 01242 01243 xxi = node_pos[3] - node_pos[2]; 01244 xet = node_pos[1] - node_pos[2]; 01245 xze = node_pos[6] - node_pos[2]; 01246 01247 current_condition = condition_comp( xxi, xet, xze ); 01248 if( current_condition > condition ) { condition = current_condition; } 01249 01250 // J(0,1,0): 01251 01252 xxi = node_pos[0] - node_pos[3]; 01253 xet = node_pos[2] - node_pos[3]; 01254 xze = node_pos[7] - node_pos[3]; 01255 01256 current_condition = condition_comp( xxi, xet, xze ); 01257 if( current_condition > condition ) { condition = current_condition; } 01258 01259 // J(0,0,1): 01260 01261 xxi = node_pos[7] - node_pos[4]; 01262 xet = node_pos[5] - node_pos[4]; 01263 xze = node_pos[0] - node_pos[4]; 01264 01265 current_condition = condition_comp( xxi, xet, xze ); 01266 if( current_condition > condition ) { condition = current_condition; } 01267 01268 // J(1,0,1): 01269 01270 xxi = node_pos[4] - node_pos[5]; 01271 xet = node_pos[6] - node_pos[5]; 01272 xze = node_pos[1] - node_pos[5]; 01273 01274 current_condition = condition_comp( xxi, xet, xze ); 01275 if( current_condition > condition ) { condition = current_condition; } 01276 01277 // J(1,1,1): 01278 01279 xxi = node_pos[5] - node_pos[6]; 01280 xet = node_pos[7] - node_pos[6]; 01281 xze = node_pos[2] - node_pos[6]; 01282 01283 current_condition = condition_comp( xxi, xet, xze ); 01284 if( current_condition > condition ) { condition = current_condition; } 01285 01286 // J(1,1,1): 01287 01288 xxi = node_pos[6] - node_pos[7]; 01289 xet = node_pos[4] - node_pos[7]; 01290 xze = node_pos[3] - node_pos[7]; 01291 01292 current_condition = condition_comp( xxi, xet, xze ); 01293 if( current_condition > condition ) { condition = current_condition; } 01294 01295 condition /= 3.0; 01296 01297 if( condition > 0 ) return (double)VERDICT_MIN( condition, VERDICT_DBL_MAX ); 01298 return (double)VERDICT_MAX( condition, -VERDICT_DBL_MAX ); 01299 } 01300 01301 /*! 01302 The maximum Frobenius condition of a hex, a.k.a. condition 01303 NB (P. Pebay 01/25/07): 01304 this method is maintained for backwards compatibility only. 01305 It will become deprecated at some point. 01306 01307 */ 01308 C_FUNC_DEF double v_hex_condition( int /*num_nodes*/, double coordinates[][3] ) 01309 { 01310 01311 return v_hex_max_aspect_frobenius( 8, coordinates ); 01312 } 01313 01314 /*! 01315 jacobian of a hex 01316 01317 Minimum pointwise volume of local map at 8 corners & center of hex 01318 */ 01319 C_FUNC_DEF double v_hex_jacobian( int /*num_nodes*/, double coordinates[][3] ) 01320 { 01321 01322 VerdictVector node_pos[8]; 01323 make_hex_nodes( coordinates, node_pos ); 01324 01325 double jacobian = VERDICT_DBL_MAX; 01326 double current_jacobian; 01327 VerdictVector xxi, xet, xze; 01328 01329 xxi = calc_hex_efg( 1, node_pos ); 01330 xet = calc_hex_efg( 2, node_pos ); 01331 xze = calc_hex_efg( 3, node_pos ); 01332 01333 current_jacobian = xxi % ( xet * xze ) / 64.0; 01334 if( current_jacobian < jacobian ) { jacobian = current_jacobian; } 01335 01336 // J(0,0,0): 01337 01338 xxi = node_pos[1] - node_pos[0]; 01339 xet = node_pos[3] - node_pos[0]; 01340 xze = node_pos[4] - node_pos[0]; 01341 01342 current_jacobian = xxi % ( xet * xze ); 01343 if( current_jacobian < jacobian ) { jacobian = current_jacobian; } 01344 01345 // J(1,0,0): 01346 01347 xxi = node_pos[2] - node_pos[1]; 01348 xet = node_pos[0] - node_pos[1]; 01349 xze = node_pos[5] - node_pos[1]; 01350 01351 current_jacobian = xxi % ( xet * xze ); 01352 if( current_jacobian < jacobian ) { jacobian = current_jacobian; } 01353 01354 // J(1,1,0): 01355 01356 xxi = node_pos[3] - node_pos[2]; 01357 xet = node_pos[1] - node_pos[2]; 01358 xze = node_pos[6] - node_pos[2]; 01359 01360 current_jacobian = xxi % ( xet * xze ); 01361 if( current_jacobian < jacobian ) { jacobian = current_jacobian; } 01362 01363 // J(0,1,0): 01364 01365 xxi = node_pos[0] - node_pos[3]; 01366 xet = node_pos[2] - node_pos[3]; 01367 xze = node_pos[7] - node_pos[3]; 01368 01369 current_jacobian = xxi % ( xet * xze ); 01370 if( current_jacobian < jacobian ) { jacobian = current_jacobian; } 01371 01372 // J(0,0,1): 01373 01374 xxi = node_pos[7] - node_pos[4]; 01375 xet = node_pos[5] - node_pos[4]; 01376 xze = node_pos[0] - node_pos[4]; 01377 01378 current_jacobian = xxi % ( xet * xze ); 01379 if( current_jacobian < jacobian ) { jacobian = current_jacobian; } 01380 01381 // J(1,0,1): 01382 01383 xxi = node_pos[4] - node_pos[5]; 01384 xet = node_pos[6] - node_pos[5]; 01385 xze = node_pos[1] - node_pos[5]; 01386 01387 current_jacobian = xxi % ( xet * xze ); 01388 if( current_jacobian < jacobian ) { jacobian = current_jacobian; } 01389 01390 // J(1,1,1): 01391 01392 xxi = node_pos[5] - node_pos[6]; 01393 xet = node_pos[7] - node_pos[6]; 01394 xze = node_pos[2] - node_pos[6]; 01395 01396 current_jacobian = xxi % ( xet * xze ); 01397 if( current_jacobian < jacobian ) { jacobian = current_jacobian; } 01398 01399 // J(0,1,1): 01400 01401 xxi = node_pos[6] - node_pos[7]; 01402 xet = node_pos[4] - node_pos[7]; 01403 xze = node_pos[3] - node_pos[7]; 01404 01405 current_jacobian = xxi % ( xet * xze ); 01406 if( current_jacobian < jacobian ) { jacobian = current_jacobian; } 01407 01408 if( jacobian > 0 ) return (double)VERDICT_MIN( jacobian, VERDICT_DBL_MAX ); 01409 return (double)VERDICT_MAX( jacobian, -VERDICT_DBL_MAX ); 01410 } 01411 01412 /*! 01413 scaled jacobian of a hex 01414 01415 Minimum Jacobian divided by the lengths of the 3 edge vectors 01416 */ 01417 C_FUNC_DEF double v_hex_scaled_jacobian( int /*num_nodes*/, double coordinates[][3] ) 01418 { 01419 01420 double jacobi, min_norm_jac = VERDICT_DBL_MAX; 01421 // double min_jacobi = VERDICT_DBL_MAX; 01422 double temp_norm_jac, lengths; 01423 double len1_sq, len2_sq, len3_sq; 01424 VerdictVector xxi, xet, xze; 01425 01426 VerdictVector node_pos[8]; 01427 make_hex_nodes( coordinates, node_pos ); 01428 01429 xxi = calc_hex_efg( 1, node_pos ); 01430 xet = calc_hex_efg( 2, node_pos ); 01431 xze = calc_hex_efg( 3, node_pos ); 01432 01433 jacobi = xxi % ( xet * xze ); 01434 // if( jacobi < min_jacobi) { min_jacobi = jacobi; } 01435 01436 len1_sq = xxi.length_squared(); 01437 len2_sq = xet.length_squared(); 01438 len3_sq = xze.length_squared(); 01439 01440 if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) 01441 return (double)VERDICT_DBL_MAX; 01442 01443 lengths = sqrt( len1_sq * len2_sq * len3_sq ); 01444 temp_norm_jac = jacobi / lengths; 01445 01446 if( temp_norm_jac < min_norm_jac ) 01447 min_norm_jac = temp_norm_jac; 01448 else 01449 temp_norm_jac = jacobi; 01450 01451 // J(0,0,0): 01452 01453 xxi = node_pos[1] - node_pos[0]; 01454 xet = node_pos[3] - node_pos[0]; 01455 xze = node_pos[4] - node_pos[0]; 01456 01457 jacobi = xxi % ( xet * xze ); 01458 // if( jacobi < min_jacobi ) { min_jacobi = jacobi; } 01459 01460 len1_sq = xxi.length_squared(); 01461 len2_sq = xet.length_squared(); 01462 len3_sq = xze.length_squared(); 01463 01464 if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) 01465 return (double)VERDICT_DBL_MAX; 01466 01467 lengths = sqrt( len1_sq * len2_sq * len3_sq ); 01468 temp_norm_jac = jacobi / lengths; 01469 if( temp_norm_jac < min_norm_jac ) 01470 min_norm_jac = temp_norm_jac; 01471 else 01472 temp_norm_jac = jacobi; 01473 01474 // J(1,0,0): 01475 01476 xxi = node_pos[2] - node_pos[1]; 01477 xet = node_pos[0] - node_pos[1]; 01478 xze = node_pos[5] - node_pos[1]; 01479 01480 jacobi = xxi % ( xet * xze ); 01481 // if( jacobi < min_jacobi ) { min_jacobi = jacobi; } 01482 01483 len1_sq = xxi.length_squared(); 01484 len2_sq = xet.length_squared(); 01485 len3_sq = xze.length_squared(); 01486 01487 if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) 01488 return (double)VERDICT_DBL_MAX; 01489 01490 lengths = sqrt( len1_sq * len2_sq * len3_sq ); 01491 temp_norm_jac = jacobi / lengths; 01492 if( temp_norm_jac < min_norm_jac ) 01493 min_norm_jac = temp_norm_jac; 01494 else 01495 temp_norm_jac = jacobi; 01496 01497 // J(1,1,0): 01498 01499 xxi = node_pos[3] - node_pos[2]; 01500 xet = node_pos[1] - node_pos[2]; 01501 xze = node_pos[6] - node_pos[2]; 01502 01503 jacobi = xxi % ( xet * xze ); 01504 // if( jacobi < min_jacobi ) { min_jacobi = jacobi; } 01505 01506 len1_sq = xxi.length_squared(); 01507 len2_sq = xet.length_squared(); 01508 len3_sq = xze.length_squared(); 01509 01510 if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) 01511 return (double)VERDICT_DBL_MAX; 01512 01513 lengths = sqrt( len1_sq * len2_sq * len3_sq ); 01514 temp_norm_jac = jacobi / lengths; 01515 if( temp_norm_jac < min_norm_jac ) 01516 min_norm_jac = temp_norm_jac; 01517 else 01518 temp_norm_jac = jacobi; 01519 01520 // J(0,1,0): 01521 01522 xxi = node_pos[0] - node_pos[3]; 01523 xet = node_pos[2] - node_pos[3]; 01524 xze = node_pos[7] - node_pos[3]; 01525 01526 jacobi = xxi % ( xet * xze ); 01527 // if( jacobi < min_jacobi ) { min_jacobi = jacobi; } 01528 01529 len1_sq = xxi.length_squared(); 01530 len2_sq = xet.length_squared(); 01531 len3_sq = xze.length_squared(); 01532 01533 if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) 01534 return (double)VERDICT_DBL_MAX; 01535 01536 lengths = sqrt( len1_sq * len2_sq * len3_sq ); 01537 temp_norm_jac = jacobi / lengths; 01538 if( temp_norm_jac < min_norm_jac ) 01539 min_norm_jac = temp_norm_jac; 01540 else 01541 temp_norm_jac = jacobi; 01542 01543 // J(0,0,1): 01544 01545 xxi = node_pos[7] - node_pos[4]; 01546 xet = node_pos[5] - node_pos[4]; 01547 xze = node_pos[0] - node_pos[4]; 01548 01549 jacobi = xxi % ( xet * xze ); 01550 // if( jacobi < min_jacobi ) { min_jacobi = jacobi; } 01551 01552 len1_sq = xxi.length_squared(); 01553 len2_sq = xet.length_squared(); 01554 len3_sq = xze.length_squared(); 01555 01556 if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) 01557 return (double)VERDICT_DBL_MAX; 01558 01559 lengths = sqrt( len1_sq * len2_sq * len3_sq ); 01560 temp_norm_jac = jacobi / lengths; 01561 if( temp_norm_jac < min_norm_jac ) 01562 min_norm_jac = temp_norm_jac; 01563 else 01564 temp_norm_jac = jacobi; 01565 01566 // J(1,0,1): 01567 01568 xxi = node_pos[4] - node_pos[5]; 01569 xet = node_pos[6] - node_pos[5]; 01570 xze = node_pos[1] - node_pos[5]; 01571 01572 jacobi = xxi % ( xet * xze ); 01573 // if( jacobi < min_jacobi ) { min_jacobi = jacobi; } 01574 01575 len1_sq = xxi.length_squared(); 01576 len2_sq = xet.length_squared(); 01577 len3_sq = xze.length_squared(); 01578 01579 if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) 01580 return (double)VERDICT_DBL_MAX; 01581 01582 lengths = sqrt( len1_sq * len2_sq * len3_sq ); 01583 temp_norm_jac = jacobi / lengths; 01584 if( temp_norm_jac < min_norm_jac ) 01585 min_norm_jac = temp_norm_jac; 01586 else 01587 temp_norm_jac = jacobi; 01588 01589 // J(1,1,1): 01590 01591 xxi = node_pos[5] - node_pos[6]; 01592 xet = node_pos[7] - node_pos[6]; 01593 xze = node_pos[2] - node_pos[6]; 01594 01595 jacobi = xxi % ( xet * xze ); 01596 // if( jacobi < min_jacobi ) { min_jacobi = jacobi; } 01597 01598 len1_sq = xxi.length_squared(); 01599 len2_sq = xet.length_squared(); 01600 len3_sq = xze.length_squared(); 01601 01602 if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) 01603 return (double)VERDICT_DBL_MAX; 01604 01605 lengths = sqrt( len1_sq * len2_sq * len3_sq ); 01606 temp_norm_jac = jacobi / lengths; 01607 if( temp_norm_jac < min_norm_jac ) 01608 min_norm_jac = temp_norm_jac; 01609 else 01610 temp_norm_jac = jacobi; 01611 01612 // J(0,1,1): 01613 01614 xxi = node_pos[6] - node_pos[7]; 01615 xet = node_pos[4] - node_pos[7]; 01616 xze = node_pos[3] - node_pos[7]; 01617 01618 jacobi = xxi % ( xet * xze ); 01619 // if( jacobi < min_jacobi ) { min_jacobi = jacobi; } 01620 01621 len1_sq = xxi.length_squared(); 01622 len2_sq = xet.length_squared(); 01623 len3_sq = xze.length_squared(); 01624 01625 if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) 01626 return (double)VERDICT_DBL_MAX; 01627 01628 lengths = sqrt( len1_sq * len2_sq * len3_sq ); 01629 temp_norm_jac = jacobi / lengths; 01630 if( temp_norm_jac < min_norm_jac ) min_norm_jac = temp_norm_jac; 01631 // else 01632 // temp_norm_jac = jacobi; 01633 01634 if( min_norm_jac > 0 ) return (double)VERDICT_MIN( min_norm_jac, VERDICT_DBL_MAX ); 01635 return (double)VERDICT_MAX( min_norm_jac, -VERDICT_DBL_MAX ); 01636 } 01637 01638 /*! 01639 shear of a hex 01640 01641 3/Condition number of Jacobian Skew matrix 01642 */ 01643 C_FUNC_DEF double v_hex_shear( int /*num_nodes*/, double coordinates[][3] ) 01644 { 01645 01646 double shear; 01647 double min_shear = 1.0; 01648 VerdictVector xxi, xet, xze; 01649 double det, len1_sq, len2_sq, len3_sq, lengths; 01650 01651 VerdictVector node_pos[8]; 01652 make_hex_nodes( coordinates, node_pos ); 01653 01654 // J(0,0,0): 01655 01656 xxi = node_pos[1] - node_pos[0]; 01657 xet = node_pos[3] - node_pos[0]; 01658 xze = node_pos[4] - node_pos[0]; 01659 01660 len1_sq = xxi.length_squared(); 01661 len2_sq = xet.length_squared(); 01662 len3_sq = xze.length_squared(); 01663 01664 if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0; 01665 01666 lengths = sqrt( len1_sq * len2_sq * len3_sq ); 01667 det = xxi % ( xet * xze ); 01668 if( det < VERDICT_DBL_MIN ) { return 0; } 01669 01670 shear = det / lengths; 01671 min_shear = VERDICT_MIN( shear, min_shear ); 01672 01673 // J(1,0,0): 01674 01675 xxi = node_pos[2] - node_pos[1]; 01676 xet = node_pos[0] - node_pos[1]; 01677 xze = node_pos[5] - node_pos[1]; 01678 01679 len1_sq = xxi.length_squared(); 01680 len2_sq = xet.length_squared(); 01681 len3_sq = xze.length_squared(); 01682 01683 if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0; 01684 01685 lengths = sqrt( len1_sq * len2_sq * len3_sq ); 01686 det = xxi % ( xet * xze ); 01687 if( det < VERDICT_DBL_MIN ) { return 0; } 01688 01689 shear = det / lengths; 01690 min_shear = VERDICT_MIN( shear, min_shear ); 01691 01692 // J(1,1,0): 01693 01694 xxi = node_pos[3] - node_pos[2]; 01695 xet = node_pos[1] - node_pos[2]; 01696 xze = node_pos[6] - node_pos[2]; 01697 01698 len1_sq = xxi.length_squared(); 01699 len2_sq = xet.length_squared(); 01700 len3_sq = xze.length_squared(); 01701 01702 if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0; 01703 01704 lengths = sqrt( len1_sq * len2_sq * len3_sq ); 01705 det = xxi % ( xet * xze ); 01706 if( det < VERDICT_DBL_MIN ) { return 0; } 01707 01708 shear = det / lengths; 01709 min_shear = VERDICT_MIN( shear, min_shear ); 01710 01711 // J(0,1,0): 01712 01713 xxi = node_pos[0] - node_pos[3]; 01714 xet = node_pos[2] - node_pos[3]; 01715 xze = node_pos[7] - node_pos[3]; 01716 01717 len1_sq = xxi.length_squared(); 01718 len2_sq = xet.length_squared(); 01719 len3_sq = xze.length_squared(); 01720 01721 if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0; 01722 01723 lengths = sqrt( len1_sq * len2_sq * len3_sq ); 01724 det = xxi % ( xet * xze ); 01725 if( det < VERDICT_DBL_MIN ) { return 0; } 01726 01727 shear = det / lengths; 01728 min_shear = VERDICT_MIN( shear, min_shear ); 01729 01730 // J(0,0,1): 01731 01732 xxi = node_pos[7] - node_pos[4]; 01733 xet = node_pos[5] - node_pos[4]; 01734 xze = node_pos[0] - node_pos[4]; 01735 01736 len1_sq = xxi.length_squared(); 01737 len2_sq = xet.length_squared(); 01738 len3_sq = xze.length_squared(); 01739 01740 if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0; 01741 01742 lengths = sqrt( len1_sq * len2_sq * len3_sq ); 01743 det = xxi % ( xet * xze ); 01744 if( det < VERDICT_DBL_MIN ) { return 0; } 01745 01746 shear = det / lengths; 01747 min_shear = VERDICT_MIN( shear, min_shear ); 01748 01749 // J(1,0,1): 01750 01751 xxi = node_pos[4] - node_pos[5]; 01752 xet = node_pos[6] - node_pos[5]; 01753 xze = node_pos[1] - node_pos[5]; 01754 01755 len1_sq = xxi.length_squared(); 01756 len2_sq = xet.length_squared(); 01757 len3_sq = xze.length_squared(); 01758 01759 if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0; 01760 01761 lengths = sqrt( len1_sq * len2_sq * len3_sq ); 01762 det = xxi % ( xet * xze ); 01763 if( det < VERDICT_DBL_MIN ) { return 0; } 01764 01765 shear = det / lengths; 01766 min_shear = VERDICT_MIN( shear, min_shear ); 01767 01768 // J(1,1,1): 01769 01770 xxi = node_pos[5] - node_pos[6]; 01771 xet = node_pos[7] - node_pos[6]; 01772 xze = node_pos[2] - node_pos[6]; 01773 01774 len1_sq = xxi.length_squared(); 01775 len2_sq = xet.length_squared(); 01776 len3_sq = xze.length_squared(); 01777 01778 if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0; 01779 01780 lengths = sqrt( len1_sq * len2_sq * len3_sq ); 01781 det = xxi % ( xet * xze ); 01782 if( det < VERDICT_DBL_MIN ) { return 0; } 01783 01784 shear = det / lengths; 01785 min_shear = VERDICT_MIN( shear, min_shear ); 01786 01787 // J(0,1,1): 01788 01789 xxi = node_pos[6] - node_pos[7]; 01790 xet = node_pos[4] - node_pos[7]; 01791 xze = node_pos[3] - node_pos[7]; 01792 01793 len1_sq = xxi.length_squared(); 01794 len2_sq = xet.length_squared(); 01795 len3_sq = xze.length_squared(); 01796 01797 if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0; 01798 01799 lengths = sqrt( len1_sq * len2_sq * len3_sq ); 01800 det = xxi % ( xet * xze ); 01801 if( det < VERDICT_DBL_MIN ) { return 0; } 01802 01803 shear = det / lengths; 01804 min_shear = VERDICT_MIN( shear, min_shear ); 01805 01806 if( min_shear <= VERDICT_DBL_MIN ) min_shear = 0; 01807 01808 if( min_shear > 0 ) return (double)VERDICT_MIN( min_shear, VERDICT_DBL_MAX ); 01809 return (double)VERDICT_MAX( min_shear, -VERDICT_DBL_MAX ); 01810 } 01811 01812 /*! 01813 shape of a hex 01814 01815 3/Condition number of weighted Jacobian matrix 01816 */ 01817 C_FUNC_DEF double v_hex_shape( int /*num_nodes*/, double coordinates[][3] ) 01818 { 01819 01820 double det, shape; 01821 double min_shape = 1.0; 01822 static const double two_thirds = 2.0 / 3.0; 01823 01824 VerdictVector xxi, xet, xze; 01825 01826 VerdictVector node_pos[8]; 01827 make_hex_nodes( coordinates, node_pos ); 01828 01829 // J(0,0,0): 01830 01831 xxi = node_pos[1] - node_pos[0]; 01832 xet = node_pos[3] - node_pos[0]; 01833 xze = node_pos[4] - node_pos[0]; 01834 01835 det = xxi % ( xet * xze ); 01836 if( det > VERDICT_DBL_MIN ) 01837 shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze ); 01838 else 01839 return 0; 01840 01841 if( shape < min_shape ) { min_shape = shape; } 01842 01843 // J(1,0,0): 01844 01845 xxi = node_pos[2] - node_pos[1]; 01846 xet = node_pos[0] - node_pos[1]; 01847 xze = node_pos[5] - node_pos[1]; 01848 01849 det = xxi % ( xet * xze ); 01850 if( det > VERDICT_DBL_MIN ) 01851 shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze ); 01852 else 01853 return 0; 01854 01855 if( shape < min_shape ) { min_shape = shape; } 01856 01857 // J(1,1,0): 01858 01859 xxi = node_pos[3] - node_pos[2]; 01860 xet = node_pos[1] - node_pos[2]; 01861 xze = node_pos[6] - node_pos[2]; 01862 01863 det = xxi % ( xet * xze ); 01864 if( det > VERDICT_DBL_MIN ) 01865 shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze ); 01866 else 01867 return 0; 01868 01869 if( shape < min_shape ) { min_shape = shape; } 01870 01871 // J(0,1,0): 01872 01873 xxi = node_pos[0] - node_pos[3]; 01874 xet = node_pos[2] - node_pos[3]; 01875 xze = node_pos[7] - node_pos[3]; 01876 01877 det = xxi % ( xet * xze ); 01878 if( det > VERDICT_DBL_MIN ) 01879 shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze ); 01880 else 01881 return 0; 01882 01883 if( shape < min_shape ) { min_shape = shape; } 01884 01885 // J(0,0,1): 01886 01887 xxi = node_pos[7] - node_pos[4]; 01888 xet = node_pos[5] - node_pos[4]; 01889 xze = node_pos[0] - node_pos[4]; 01890 01891 det = xxi % ( xet * xze ); 01892 if( det > VERDICT_DBL_MIN ) 01893 shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze ); 01894 else 01895 return 0; 01896 01897 if( shape < min_shape ) { min_shape = shape; } 01898 01899 // J(1,0,1): 01900 01901 xxi = node_pos[4] - node_pos[5]; 01902 xet = node_pos[6] - node_pos[5]; 01903 xze = node_pos[1] - node_pos[5]; 01904 01905 det = xxi % ( xet * xze ); 01906 if( det > VERDICT_DBL_MIN ) 01907 shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze ); 01908 else 01909 return 0; 01910 01911 if( shape < min_shape ) { min_shape = shape; } 01912 01913 // J(1,1,1): 01914 01915 xxi = node_pos[5] - node_pos[6]; 01916 xet = node_pos[7] - node_pos[6]; 01917 xze = node_pos[2] - node_pos[6]; 01918 01919 det = xxi % ( xet * xze ); 01920 if( det > VERDICT_DBL_MIN ) 01921 shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze ); 01922 else 01923 return 0; 01924 01925 if( shape < min_shape ) { min_shape = shape; } 01926 01927 // J(1,1,1): 01928 01929 xxi = node_pos[6] - node_pos[7]; 01930 xet = node_pos[4] - node_pos[7]; 01931 xze = node_pos[3] - node_pos[7]; 01932 01933 det = xxi % ( xet * xze ); 01934 if( det > VERDICT_DBL_MIN ) 01935 shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze ); 01936 else 01937 return 0; 01938 01939 if( shape < min_shape ) { min_shape = shape; } 01940 01941 if( min_shape <= VERDICT_DBL_MIN ) min_shape = 0; 01942 01943 if( min_shape > 0 ) return (double)VERDICT_MIN( min_shape, VERDICT_DBL_MAX ); 01944 return (double)VERDICT_MAX( min_shape, -VERDICT_DBL_MAX ); 01945 } 01946 01947 /*! 01948 relative size of a hex 01949 01950 Min( J, 1/J ), where J is determinant of weighted Jacobian matrix 01951 */ 01952 C_FUNC_DEF double v_hex_relative_size_squared( int /*num_nodes*/, double coordinates[][3] ) 01953 { 01954 double size = 0; 01955 double tau; 01956 01957 VerdictVector xxi, xet, xze; 01958 double det, det_sum = 0; 01959 01960 v_hex_get_weight( xxi, xet, xze ); 01961 01962 // This is the average relative size 01963 double detw = xxi % ( xet * xze ); 01964 01965 if( detw < VERDICT_DBL_MIN ) return 0; 01966 01967 VerdictVector node_pos[8]; 01968 make_hex_nodes( coordinates, node_pos ); 01969 01970 // J(0,0,0): 01971 01972 xxi = node_pos[1] - node_pos[0]; 01973 xet = node_pos[3] - node_pos[0]; 01974 xze = node_pos[4] - node_pos[0]; 01975 01976 det = xxi % ( xet * xze ); 01977 det_sum += det; 01978 01979 // J(1,0,0): 01980 01981 xxi = node_pos[2] - node_pos[1]; 01982 xet = node_pos[0] - node_pos[1]; 01983 xze = node_pos[5] - node_pos[1]; 01984 01985 det = xxi % ( xet * xze ); 01986 det_sum += det; 01987 01988 // J(0,1,0): 01989 01990 xxi = node_pos[3] - node_pos[2]; 01991 xet = node_pos[1] - node_pos[2]; 01992 xze = node_pos[6] - node_pos[2]; 01993 01994 det = xxi % ( xet * xze ); 01995 det_sum += det; 01996 01997 // J(1,1,0): 01998 01999 xxi = node_pos[0] - node_pos[3]; 02000 xet = node_pos[2] - node_pos[3]; 02001 xze = node_pos[7] - node_pos[3]; 02002 02003 det = xxi % ( xet * xze ); 02004 det_sum += det; 02005 02006 // J(0,1,0): 02007 02008 xxi = node_pos[7] - node_pos[4]; 02009 xet = node_pos[5] - node_pos[4]; 02010 xze = node_pos[0] - node_pos[4]; 02011 02012 det = xxi % ( xet * xze ); 02013 det_sum += det; 02014 02015 // J(1,0,1): 02016 02017 xxi = node_pos[4] - node_pos[5]; 02018 xet = node_pos[6] - node_pos[5]; 02019 xze = node_pos[1] - node_pos[5]; 02020 02021 det = xxi % ( xet * xze ); 02022 det_sum += det; 02023 02024 // J(1,1,1): 02025 02026 xxi = node_pos[5] - node_pos[6]; 02027 xet = node_pos[7] - node_pos[6]; 02028 xze = node_pos[2] - node_pos[6]; 02029 02030 det = xxi % ( xet * xze ); 02031 det_sum += det; 02032 02033 // J(1,1,1): 02034 02035 xxi = node_pos[6] - node_pos[7]; 02036 xet = node_pos[4] - node_pos[7]; 02037 xze = node_pos[3] - node_pos[7]; 02038 02039 det = xxi % ( xet * xze ); 02040 det_sum += det; 02041 02042 if( det_sum > VERDICT_DBL_MIN ) 02043 { 02044 tau = det_sum / ( 8 * detw ); 02045 02046 tau = VERDICT_MIN( tau, 1.0 / tau ); 02047 02048 size = tau * tau; 02049 } 02050 02051 if( size > 0 ) return (double)VERDICT_MIN( size, VERDICT_DBL_MAX ); 02052 return (double)VERDICT_MAX( size, -VERDICT_DBL_MAX ); 02053 } 02054 02055 /*! 02056 shape and size of a hex 02057 02058 Product of Shape and Relative Size 02059 */ 02060 C_FUNC_DEF double v_hex_shape_and_size( int num_nodes, double coordinates[][3] ) 02061 { 02062 double size = v_hex_relative_size_squared( num_nodes, coordinates ); 02063 double shape = v_hex_shape( num_nodes, coordinates ); 02064 02065 double shape_size = size * shape; 02066 02067 if( shape_size > 0 ) return (double)VERDICT_MIN( shape_size, VERDICT_DBL_MAX ); 02068 return (double)VERDICT_MAX( shape_size, -VERDICT_DBL_MAX ); 02069 } 02070 02071 /*! 02072 shear and size of a hex 02073 02074 Product of Shear and Relative Size 02075 */ 02076 C_FUNC_DEF double v_hex_shear_and_size( int num_nodes, double coordinates[][3] ) 02077 { 02078 double size = v_hex_relative_size_squared( num_nodes, coordinates ); 02079 double shear = v_hex_shear( num_nodes, coordinates ); 02080 02081 double shear_size = shear * size; 02082 02083 if( shear_size > 0 ) return (double)VERDICT_MIN( shear_size, VERDICT_DBL_MAX ); 02084 return (double)VERDICT_MAX( shear_size, -VERDICT_DBL_MAX ); 02085 } 02086 02087 /*! 02088 distortion of a hex 02089 */ 02090 C_FUNC_DEF double v_hex_distortion( int num_nodes, double coordinates[][3] ) 02091 { 02092 02093 // use 2x2 gauss points for linear hex and 3x3 for 2nd order hex 02094 int number_of_gauss_points = 0; 02095 if( num_nodes == 8 ) 02096 // 2x2 quadrature rule 02097 number_of_gauss_points = 2; 02098 else if( num_nodes == 20 ) 02099 // 3x3 quadrature rule 02100 number_of_gauss_points = 3; 02101 02102 int number_dimension = 3; 02103 int total_number_of_gauss_points = number_of_gauss_points * number_of_gauss_points * number_of_gauss_points; 02104 double distortion = VERDICT_DBL_MAX; 02105 02106 // maxTotalNumberGaussPoints =27, maxNumberNodes = 20 02107 // they are defined in GaussIntegration.hpp 02108 // This is used to make these arrays static. 02109 // I tried dynamically allocated arrays but the new and delete 02110 // was very expensive 02111 02112 double shape_function[maxTotalNumberGaussPoints][maxNumberNodes]; 02113 double dndy1[maxTotalNumberGaussPoints][maxNumberNodes]; 02114 double dndy2[maxTotalNumberGaussPoints][maxNumberNodes]; 02115 double dndy3[maxTotalNumberGaussPoints][maxNumberNodes]; 02116 double weight[maxTotalNumberGaussPoints]; 02117 02118 // create an object of GaussIntegration 02119 GaussIntegration::initialize( number_of_gauss_points, num_nodes, number_dimension ); 02120 GaussIntegration::calculate_shape_function_3d_hex(); 02121 GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], dndy3[0], weight ); 02122 02123 VerdictVector xxi, xet, xze, xin; 02124 02125 double jacobian, minimum_jacobian; 02126 double element_volume = 0.0; 02127 minimum_jacobian = VERDICT_DBL_MAX; 02128 // calculate element volume 02129 int ife, ja; 02130 for( ife = 0; ife < total_number_of_gauss_points; ife++ ) 02131 { 02132 02133 xxi.set( 0.0, 0.0, 0.0 ); 02134 xet.set( 0.0, 0.0, 0.0 ); 02135 xze.set( 0.0, 0.0, 0.0 ); 02136 02137 for( ja = 0; ja < num_nodes; ja++ ) 02138 { 02139 xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] ); 02140 xxi += dndy1[ife][ja] * xin; 02141 xet += dndy2[ife][ja] * xin; 02142 xze += dndy3[ife][ja] * xin; 02143 } 02144 02145 jacobian = xxi % ( xet * xze ); 02146 if( minimum_jacobian > jacobian ) minimum_jacobian = jacobian; 02147 02148 element_volume += weight[ife] * jacobian; 02149 } 02150 02151 // loop through all nodes 02152 double dndy1_at_node[maxNumberNodes][maxNumberNodes]; 02153 double dndy2_at_node[maxNumberNodes][maxNumberNodes]; 02154 double dndy3_at_node[maxNumberNodes][maxNumberNodes]; 02155 02156 GaussIntegration::calculate_derivative_at_nodes_3d( dndy1_at_node, dndy2_at_node, dndy3_at_node ); 02157 int node_id; 02158 for( node_id = 0; node_id < num_nodes; node_id++ ) 02159 { 02160 02161 xxi.set( 0.0, 0.0, 0.0 ); 02162 xet.set( 0.0, 0.0, 0.0 ); 02163 xze.set( 0.0, 0.0, 0.0 ); 02164 02165 for( ja = 0; ja < num_nodes; ja++ ) 02166 { 02167 xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] ); 02168 xxi += dndy1_at_node[node_id][ja] * xin; 02169 xet += dndy2_at_node[node_id][ja] * xin; 02170 xze += dndy3_at_node[node_id][ja] * xin; 02171 } 02172 02173 jacobian = xxi % ( xet * xze ); 02174 if( minimum_jacobian > jacobian ) minimum_jacobian = jacobian; 02175 } 02176 distortion = minimum_jacobian / element_volume * 8.; 02177 return (double)distortion; 02178 } 02179 02180 /* 02181 C_FUNC_DEF double hex_jac_normjac_oddy_cond( int choices[], 02182 double coordinates[][3], 02183 double answers[4] ) 02184 { 02185 02186 //Define variables 02187 int i; 02188 02189 double xxi[3], xet[3], xze[3]; 02190 double norm_jacobian = 0.0, current_norm_jac = 0.0; 02191 double jacobian = 0.0, current_jacobian = 0.0; 02192 double oddy = 0.0, current_oddy = 0.0; 02193 double condition = 0.0, current_condition = 0.0; 02194 02195 02196 for( i=0; i<3; i++) 02197 xxi[i] = calc_hex_efg( 2, i, coordinates ); 02198 for( i=0; i<3; i++) 02199 xet[i] = calc_hex_efg( 3, i, coordinates ); 02200 for( i=0; i<3; i++) 02201 xze[i] = calc_hex_efg( 6, i, coordinates ); 02202 02203 // norm jacobian and jacobian 02204 if( choices[0] || choices[1] ) 02205 { 02206 current_jacobian = determinant( xxi, xet, xze ); 02207 current_norm_jac = normalize_jacobian( current_jacobian, 02208 xxi, xet, xze ); 02209 02210 if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 02211 02212 current_jacobian /= 64.0; 02213 if (current_jacobian < jacobian) { jacobian = current_jacobian; } 02214 } 02215 02216 // oddy 02217 if( choices[2] ) 02218 { 02219 current_oddy = oddy_comp( xxi, xet, xze ); 02220 if ( current_oddy > oddy ) { oddy = current_oddy; } 02221 } 02222 02223 // condition 02224 if( choices[3] ) 02225 { 02226 current_condition = condition_comp( xxi, xet, xze ); 02227 if ( current_condition > condition ) { condition = current_condition; } 02228 } 02229 02230 02231 for( i=0; i<3; i++) 02232 { 02233 xxi[i] = coordinates[1][i] - coordinates[0][i]; 02234 xet[i] = coordinates[3][i] - coordinates[0][i]; 02235 xze[i] = coordinates[4][i] - coordinates[0][i]; 02236 } 02237 02238 // norm jacobian and jacobian 02239 if( choices[0] || choices[1] ) 02240 { 02241 current_jacobian = determinant( xxi, xet, xze ); 02242 current_norm_jac = normalize_jacobian( current_jacobian, 02243 xxi, xet, xze ); 02244 02245 if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 02246 if (current_jacobian < jacobian) { jacobian = current_jacobian; } 02247 } 02248 02249 // oddy 02250 if( choices[2] ) 02251 { 02252 current_oddy = oddy_comp( xxi, xet, xze ); 02253 if ( current_oddy > oddy ) { oddy = current_oddy; } 02254 } 02255 02256 // condition 02257 if( choices[3] ) 02258 { 02259 current_condition = condition_comp( xxi, xet, xze ); 02260 if ( current_condition > condition ) { condition = current_condition; } 02261 } 02262 02263 02264 for( i=0; i<3; i++) 02265 { 02266 xxi[i] = coordinates[1][i] - coordinates[0][i]; 02267 xet[i] = coordinates[2][i] - coordinates[1][i]; 02268 xze[i] = coordinates[5][i] - coordinates[1][i]; 02269 } 02270 02271 // norm jacobian and jacobian 02272 if( choices[0] || choices[1] ) 02273 { 02274 current_jacobian = determinant( xxi, xet, xze ); 02275 current_norm_jac = normalize_jacobian( current_jacobian, 02276 xxi, xet, xze ); 02277 02278 if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 02279 if (current_jacobian < jacobian) { jacobian = current_jacobian; } 02280 } 02281 02282 // oddy 02283 if( choices[2] ) 02284 { 02285 current_oddy = oddy_comp( xxi, xet, xze ); 02286 if ( current_oddy > oddy ) { oddy = current_oddy; } 02287 } 02288 02289 // condition 02290 if( choices[3] ) 02291 { 02292 current_condition = condition_comp( xxi, xet, xze ); 02293 if ( current_condition > condition ) { condition = current_condition; } 02294 } 02295 02296 02297 for( i=0; i<3; i++) 02298 { 02299 xxi[i] = coordinates[2][i] - coordinates[3][i]; 02300 xet[i] = coordinates[3][i] - coordinates[0][i]; 02301 xze[i] = coordinates[7][i] - coordinates[3][i]; 02302 } 02303 02304 // norm jacobian and jacobian 02305 if( choices[0] || choices[1] ) 02306 { 02307 current_jacobian = determinant( xxi, xet, xze ); 02308 current_norm_jac = normalize_jacobian( current_jacobian, 02309 xxi, xet, xze ); 02310 02311 if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 02312 if (current_jacobian < jacobian) { jacobian = current_jacobian; } 02313 } 02314 02315 // oddy 02316 if( choices[2] ) 02317 { 02318 current_oddy = oddy_comp( xxi, xet, xze ); 02319 if ( current_oddy > oddy ) { oddy = current_oddy; } 02320 } 02321 02322 // condition 02323 if( choices[3] ) 02324 { 02325 current_condition = condition_comp( xxi, xet, xze ); 02326 if ( current_condition > condition ) { condition = current_condition; } 02327 } 02328 02329 02330 for( i=0; i<3; i++) 02331 { 02332 xxi[i] = coordinates[5][i] - coordinates[4][i]; 02333 xet[i] = coordinates[7][i] - coordinates[4][i]; 02334 xze[i] = coordinates[4][i] - coordinates[0][i]; 02335 } 02336 02337 02338 // norm jacobian and jacobian 02339 if( choices[0] || choices[1] ) 02340 { 02341 current_jacobian = determinant( xxi, xet, xze ); 02342 current_norm_jac = normalize_jacobian( current_jacobian, 02343 xxi, xet, xze ); 02344 02345 if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 02346 if (current_jacobian < jacobian) { jacobian = current_jacobian; } 02347 } 02348 02349 // oddy 02350 if( choices[2] ) 02351 { 02352 current_oddy = oddy_comp( xxi, xet, xze ); 02353 if ( current_oddy > oddy ) { oddy = current_oddy; } 02354 } 02355 02356 // condition 02357 if( choices[3] ) 02358 { 02359 current_condition = condition_comp( xxi, xet, xze ); 02360 if ( current_condition > condition ) { condition = current_condition; } 02361 } 02362 02363 02364 for( i=0; i<3; i++) 02365 { 02366 xxi[i] = coordinates[2][i] - coordinates[3][i]; 02367 xet[i] = coordinates[2][i] - coordinates[1][i]; 02368 xze[i] = coordinates[6][i] - coordinates[2][i]; 02369 } 02370 02371 // norm jacobian and jacobian 02372 if( choices[0] || choices[1] ) 02373 { 02374 current_jacobian = determinant( xxi, xet, xze ); 02375 current_norm_jac = normalize_jacobian( current_jacobian, 02376 xxi, xet, xze ); 02377 02378 if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 02379 if (current_jacobian < jacobian) { jacobian = current_jacobian; } 02380 } 02381 02382 // oddy 02383 if( choices[2] ) 02384 { 02385 current_oddy = oddy_comp( xxi, xet, xze ); 02386 if ( current_oddy > oddy ) { oddy = current_oddy; } 02387 } 02388 02389 // condition 02390 if( choices[3] ) 02391 { 02392 current_condition = condition_comp( xxi, xet, xze ); 02393 if ( current_condition > condition ) { condition = current_condition; } 02394 } 02395 02396 02397 for( i=0; i<3; i++) 02398 { 02399 xxi[i] = coordinates[5][i] - coordinates[4][i]; 02400 xet[i] = coordinates[6][i] - coordinates[5][i]; 02401 xze[i] = coordinates[5][i] - coordinates[1][i]; 02402 } 02403 02404 02405 // norm jacobian and jacobian 02406 if( choices[0] || choices[1] ) 02407 { 02408 current_jacobian = determinant( xxi, xet, xze ); 02409 current_norm_jac = normalize_jacobian( current_jacobian, 02410 xxi, xet, xze ); 02411 02412 if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 02413 if (current_jacobian < jacobian) { jacobian = current_jacobian; } 02414 } 02415 02416 // oddy 02417 if( choices[2] ) 02418 { 02419 current_oddy = oddy_comp( xxi, xet, xze ); 02420 if ( current_oddy > oddy ) { oddy = current_oddy; } 02421 } 02422 02423 // condition 02424 if( choices[3] ) 02425 { 02426 current_condition = condition_comp( xxi, xet, xze ); 02427 if ( current_condition > condition ) { condition = current_condition; } 02428 } 02429 02430 02431 for( i=0; i<3; i++) 02432 { 02433 xxi[i] = coordinates[6][i] - coordinates[7][i]; 02434 xet[i] = coordinates[7][i] - coordinates[4][i]; 02435 xze[i] = coordinates[7][i] - coordinates[3][i]; 02436 } 02437 02438 02439 // norm jacobian and jacobian 02440 if( choices[0] || choices[1] ) 02441 { 02442 current_jacobian = determinant( xxi, xet, xze ); 02443 current_norm_jac = normalize_jacobian( current_jacobian, 02444 xxi, xet, xze ); 02445 02446 if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 02447 if (current_jacobian < jacobian) { jacobian = current_jacobian; } 02448 } 02449 02450 // oddy 02451 if( choices[2] ) 02452 { 02453 current_oddy = oddy_comp( xxi, xet, xze ); 02454 if ( current_oddy > oddy ) { oddy = current_oddy; } 02455 } 02456 02457 // condition 02458 if( choices[3] ) 02459 { 02460 current_condition = condition_comp( xxi, xet, xze ); 02461 if ( current_condition > condition ) { condition = current_condition; } 02462 } 02463 02464 02465 for( i=0; i<3; i++) 02466 { 02467 xxi[i] = coordinates[6][i] - coordinates[7][i]; 02468 xet[i] = coordinates[6][i] - coordinates[5][i]; 02469 xze[i] = coordinates[6][i] - coordinates[2][i]; 02470 } 02471 02472 // norm jacobian and jacobian 02473 if( choices[0] || choices[1] ) 02474 { 02475 current_jacobian = determinant( xxi, xet, xze ); 02476 current_norm_jac = normalize_jacobian( current_jacobian, 02477 xxi, xet, xze ); 02478 02479 if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 02480 if (current_jacobian < jacobian) { jacobian = current_jacobian; } 02481 } 02482 02483 // oddy 02484 if( choices[2] ) 02485 { 02486 current_oddy = oddy_comp( xxi, xet, xze ); 02487 if ( current_oddy > oddy ) { oddy = current_oddy; } 02488 } 02489 02490 // condition 02491 if( choices[3] ) 02492 { 02493 current_condition = condition_comp( xxi, xet, xze ); 02494 if ( current_condition > condition ) { condition = current_condition; } 02495 02496 condition /= 3.0; 02497 } 02498 02499 02500 answers[0] = jacobian; 02501 answers[1] = norm_jacobian; 02502 answers[2] = oddy; 02503 answers[3] = condition; 02504 02505 return 1.0; 02506 02507 } 02508 */ 02509 02510 /*! 02511 multiple quality metrics of a hex 02512 */ 02513 C_FUNC_DEF void v_hex_quality( int num_nodes, double coordinates[][3], unsigned int metrics_request_flag, 02514 HexMetricVals* metric_vals ) 02515 { 02516 memset( metric_vals, 0, sizeof( HexMetricVals ) ); 02517 02518 // max edge ratio, skew, taper, and volume 02519 if( metrics_request_flag & ( V_HEX_MAX_EDGE_RATIO | V_HEX_SKEW | V_HEX_TAPER ) ) 02520 { 02521 VerdictVector node_pos[8]; 02522 make_hex_nodes( coordinates, node_pos ); 02523 02524 VerdictVector efg1, efg2, efg3; 02525 efg1 = calc_hex_efg( 1, node_pos ); 02526 efg2 = calc_hex_efg( 2, node_pos ); 02527 efg3 = calc_hex_efg( 3, node_pos ); 02528 02529 if( metrics_request_flag & V_HEX_MAX_EDGE_RATIO ) 02530 { 02531 double max_edge_ratio_12, max_edge_ratio_13, max_edge_ratio_23; 02532 02533 double mag_efg1 = efg1.length(); 02534 double mag_efg2 = efg2.length(); 02535 double mag_efg3 = efg3.length(); 02536 02537 max_edge_ratio_12 = safe_ratio( VERDICT_MAX( mag_efg1, mag_efg2 ), VERDICT_MIN( mag_efg1, mag_efg2 ) ); 02538 max_edge_ratio_13 = safe_ratio( VERDICT_MAX( mag_efg1, mag_efg3 ), VERDICT_MIN( mag_efg1, mag_efg3 ) ); 02539 max_edge_ratio_23 = safe_ratio( VERDICT_MAX( mag_efg2, mag_efg3 ), VERDICT_MIN( mag_efg2, mag_efg3 ) ); 02540 02541 metric_vals->max_edge_ratio = 02542 (double)VERDICT_MAX( max_edge_ratio_12, VERDICT_MAX( max_edge_ratio_13, max_edge_ratio_23 ) ); 02543 } 02544 02545 if( metrics_request_flag & V_HEX_SKEW ) 02546 { 02547 02548 VerdictVector vec1 = efg1; 02549 VerdictVector vec2 = efg2; 02550 VerdictVector vec3 = efg3; 02551 02552 if( vec1.normalize() <= VERDICT_DBL_MIN || vec2.normalize() <= VERDICT_DBL_MIN || 02553 vec3.normalize() <= VERDICT_DBL_MIN ) 02554 { metric_vals->skew = (double)VERDICT_DBL_MAX; } 02555 else 02556 { 02557 double skewx = fabs( vec1 % vec2 ); 02558 double skewy = fabs( vec1 % vec3 ); 02559 double skewz = fabs( vec2 % vec3 ); 02560 02561 metric_vals->skew = (double)( VERDICT_MAX( skewx, VERDICT_MAX( skewy, skewz ) ) ); 02562 } 02563 } 02564 02565 if( metrics_request_flag & V_HEX_TAPER ) 02566 { 02567 VerdictVector efg12 = calc_hex_efg( 12, node_pos ); 02568 VerdictVector efg13 = calc_hex_efg( 13, node_pos ); 02569 VerdictVector efg23 = calc_hex_efg( 23, node_pos ); 02570 02571 double taperx = fabs( safe_ratio( efg12.length(), VERDICT_MIN( efg1.length(), efg2.length() ) ) ); 02572 double tapery = fabs( safe_ratio( efg13.length(), VERDICT_MIN( efg1.length(), efg3.length() ) ) ); 02573 double taperz = fabs( safe_ratio( efg23.length(), VERDICT_MIN( efg2.length(), efg3.length() ) ) ); 02574 02575 metric_vals->taper = (double)VERDICT_MAX( taperx, VERDICT_MAX( tapery, taperz ) ); 02576 } 02577 } 02578 02579 if( metrics_request_flag & V_HEX_VOLUME ) { metric_vals->volume = v_hex_volume( 8, coordinates ); } 02580 02581 if( metrics_request_flag & 02582 ( V_HEX_JACOBIAN | V_HEX_SCALED_JACOBIAN | V_HEX_CONDITION | V_HEX_SHEAR | V_HEX_SHAPE | 02583 V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE | V_HEX_STRETCH ) ) 02584 { 02585 02586 static const double two_thirds = 2.0 / 3.0; 02587 VerdictVector edges[12]; 02588 // the length squares 02589 double length_squared[12]; 02590 // make vectors from coordinates 02591 make_hex_edges( coordinates, edges ); 02592 02593 // calculate the length squares if we need to 02594 if( metrics_request_flag & 02595 ( V_HEX_JACOBIAN | V_HEX_SHEAR | V_HEX_SCALED_JACOBIAN | V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE | 02596 V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHEAR_AND_SIZE | V_HEX_STRETCH ) ) 02597 make_edge_length_squares( edges, length_squared ); 02598 02599 double jacobian = VERDICT_DBL_MAX, scaled_jacobian = VERDICT_DBL_MAX, condition = 0.0, shear = 1.0, shape = 1.0, 02600 oddy = 0.0; 02601 double current_jacobian, current_scaled_jacobian, current_condition, current_shape, detw = 0, det_sum = 0, 02602 current_oddy; 02603 VerdictBoolean rel_size_error = VERDICT_FALSE; 02604 02605 VerdictVector xxi, xet, xze; 02606 02607 // get weights if we need based on average size of a hex 02608 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02609 { 02610 v_hex_get_weight( xxi, xet, xze ); 02611 detw = xxi % ( xet * xze ); 02612 if( detw < VERDICT_DBL_MIN ) rel_size_error = VERDICT_TRUE; 02613 } 02614 02615 xxi = edges[0] - edges[2] + edges[4] - edges[6]; 02616 xet = edges[1] - edges[3] + edges[5] - edges[7]; 02617 xze = edges[8] + edges[9] + edges[10] + edges[11]; 02618 02619 current_jacobian = xxi % ( xet * xze ) / 64.0; 02620 if( current_jacobian < jacobian ) jacobian = current_jacobian; 02621 02622 if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) ) 02623 { 02624 current_jacobian *= 64.0; 02625 current_scaled_jacobian = 02626 current_jacobian / sqrt( xxi.length_squared() * xet.length_squared() * xze.length_squared() ); 02627 if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian; 02628 } 02629 02630 if( metrics_request_flag & V_HEX_CONDITION ) 02631 { 02632 current_condition = condition_comp( xxi, xet, xze ); 02633 if( current_condition > condition ) { condition = current_condition; } 02634 } 02635 02636 if( metrics_request_flag & V_HEX_ODDY ) 02637 { 02638 current_oddy = oddy_comp( xxi, xet, xze ); 02639 if( current_oddy > oddy ) { oddy = current_oddy; } 02640 } 02641 02642 // J(0,0,0) 02643 current_jacobian = edges[0] % ( -edges[3] * edges[8] ); 02644 if( current_jacobian < jacobian ) jacobian = current_jacobian; 02645 02646 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02647 { det_sum += current_jacobian; } 02648 02649 if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) ) 02650 { 02651 if( length_squared[0] <= VERDICT_DBL_MIN || length_squared[3] <= VERDICT_DBL_MIN || 02652 length_squared[8] <= VERDICT_DBL_MIN ) 02653 { current_scaled_jacobian = VERDICT_DBL_MAX; } 02654 else 02655 { 02656 current_scaled_jacobian = 02657 current_jacobian / sqrt( length_squared[0] * length_squared[3] * length_squared[8] ); 02658 } 02659 if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian; 02660 } 02661 02662 if( metrics_request_flag & V_HEX_CONDITION ) 02663 { 02664 current_condition = condition_comp( edges[0], -edges[3], edges[8] ); 02665 if( current_condition > condition ) { condition = current_condition; } 02666 } 02667 02668 if( metrics_request_flag & V_HEX_ODDY ) 02669 { 02670 current_oddy = oddy_comp( edges[0], -edges[3], edges[8] ); 02671 if( current_oddy > oddy ) { oddy = current_oddy; } 02672 } 02673 02674 if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) ) 02675 { 02676 if( current_jacobian > VERDICT_DBL_MIN ) 02677 current_shape = 3 * pow( current_jacobian, two_thirds ) / 02678 ( length_squared[0] + length_squared[3] + length_squared[8] ); 02679 else 02680 current_shape = 0; 02681 02682 if( current_shape < shape ) { shape = current_shape; } 02683 } 02684 02685 // J(1,0,0) 02686 current_jacobian = edges[1] % ( -edges[0] * edges[9] ); 02687 if( current_jacobian < jacobian ) jacobian = current_jacobian; 02688 02689 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02690 { det_sum += current_jacobian; } 02691 02692 if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) ) 02693 { 02694 if( length_squared[1] <= VERDICT_DBL_MIN || length_squared[0] <= VERDICT_DBL_MIN || 02695 length_squared[9] <= VERDICT_DBL_MIN ) 02696 { current_scaled_jacobian = VERDICT_DBL_MAX; } 02697 else 02698 { 02699 current_scaled_jacobian = 02700 current_jacobian / sqrt( length_squared[1] * length_squared[0] * length_squared[9] ); 02701 } 02702 if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian; 02703 } 02704 02705 if( metrics_request_flag & V_HEX_CONDITION ) 02706 { 02707 current_condition = condition_comp( edges[1], -edges[0], edges[9] ); 02708 if( current_condition > condition ) { condition = current_condition; } 02709 } 02710 02711 if( metrics_request_flag & V_HEX_ODDY ) 02712 { 02713 current_oddy = oddy_comp( edges[1], -edges[0], edges[9] ); 02714 if( current_oddy > oddy ) { oddy = current_oddy; } 02715 } 02716 02717 if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) ) 02718 { 02719 if( current_jacobian > VERDICT_DBL_MIN ) 02720 current_shape = 3 * pow( current_jacobian, two_thirds ) / 02721 ( length_squared[1] + length_squared[0] + length_squared[9] ); 02722 else 02723 current_shape = 0; 02724 02725 if( current_shape < shape ) { shape = current_shape; } 02726 } 02727 02728 // J(1,1,0) 02729 current_jacobian = edges[2] % ( -edges[1] * edges[10] ); 02730 if( current_jacobian < jacobian ) jacobian = current_jacobian; 02731 02732 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02733 { det_sum += current_jacobian; } 02734 02735 if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) ) 02736 { 02737 if( length_squared[2] <= VERDICT_DBL_MIN || length_squared[1] <= VERDICT_DBL_MIN || 02738 length_squared[10] <= VERDICT_DBL_MIN ) 02739 { current_scaled_jacobian = VERDICT_DBL_MAX; } 02740 else 02741 { 02742 current_scaled_jacobian = 02743 current_jacobian / sqrt( length_squared[2] * length_squared[1] * length_squared[10] ); 02744 } 02745 if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian; 02746 } 02747 02748 if( metrics_request_flag & V_HEX_CONDITION ) 02749 { 02750 current_condition = condition_comp( edges[2], -edges[1], edges[10] ); 02751 if( current_condition > condition ) { condition = current_condition; } 02752 } 02753 02754 if( metrics_request_flag & V_HEX_ODDY ) 02755 { 02756 current_oddy = oddy_comp( edges[2], -edges[1], edges[10] ); 02757 if( current_oddy > oddy ) { oddy = current_oddy; } 02758 } 02759 02760 if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) ) 02761 { 02762 if( current_jacobian > VERDICT_DBL_MIN ) 02763 current_shape = 3 * pow( current_jacobian, two_thirds ) / 02764 ( length_squared[2] + length_squared[1] + length_squared[10] ); 02765 else 02766 current_shape = 0; 02767 02768 if( current_shape < shape ) { shape = current_shape; } 02769 } 02770 02771 // J(0,1,0) 02772 current_jacobian = edges[3] % ( -edges[2] * edges[11] ); 02773 if( current_jacobian < jacobian ) jacobian = current_jacobian; 02774 02775 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02776 { det_sum += current_jacobian; } 02777 02778 if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) ) 02779 { 02780 if( length_squared[3] <= VERDICT_DBL_MIN || length_squared[2] <= VERDICT_DBL_MIN || 02781 length_squared[11] <= VERDICT_DBL_MIN ) 02782 { current_scaled_jacobian = VERDICT_DBL_MAX; } 02783 else 02784 { 02785 current_scaled_jacobian = 02786 current_jacobian / sqrt( length_squared[3] * length_squared[2] * length_squared[11] ); 02787 } 02788 if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian; 02789 } 02790 02791 if( metrics_request_flag & V_HEX_CONDITION ) 02792 { 02793 current_condition = condition_comp( edges[3], -edges[2], edges[11] ); 02794 if( current_condition > condition ) { condition = current_condition; } 02795 } 02796 02797 if( metrics_request_flag & V_HEX_ODDY ) 02798 { 02799 current_oddy = oddy_comp( edges[3], -edges[2], edges[11] ); 02800 if( current_oddy > oddy ) { oddy = current_oddy; } 02801 } 02802 02803 if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) ) 02804 { 02805 if( current_jacobian > VERDICT_DBL_MIN ) 02806 current_shape = 3 * pow( current_jacobian, two_thirds ) / 02807 ( length_squared[3] + length_squared[2] + length_squared[11] ); 02808 else 02809 current_shape = 0; 02810 02811 if( current_shape < shape ) { shape = current_shape; } 02812 } 02813 02814 // J(0,0,1) 02815 current_jacobian = edges[4] % ( -edges[8] * -edges[7] ); 02816 if( current_jacobian < jacobian ) jacobian = current_jacobian; 02817 02818 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02819 { det_sum += current_jacobian; } 02820 02821 if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) ) 02822 { 02823 if( length_squared[4] <= VERDICT_DBL_MIN || length_squared[8] <= VERDICT_DBL_MIN || 02824 length_squared[7] <= VERDICT_DBL_MIN ) 02825 { current_scaled_jacobian = VERDICT_DBL_MAX; } 02826 else 02827 { 02828 current_scaled_jacobian = 02829 current_jacobian / sqrt( length_squared[4] * length_squared[8] * length_squared[7] ); 02830 } 02831 if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian; 02832 } 02833 02834 if( metrics_request_flag & V_HEX_CONDITION ) 02835 { 02836 current_condition = condition_comp( edges[4], -edges[8], -edges[7] ); 02837 if( current_condition > condition ) { condition = current_condition; } 02838 } 02839 02840 if( metrics_request_flag & V_HEX_ODDY ) 02841 { 02842 current_oddy = oddy_comp( edges[4], -edges[8], -edges[7] ); 02843 if( current_oddy > oddy ) { oddy = current_oddy; } 02844 } 02845 02846 if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) ) 02847 { 02848 if( current_jacobian > VERDICT_DBL_MIN ) 02849 current_shape = 3 * pow( current_jacobian, two_thirds ) / 02850 ( length_squared[4] + length_squared[8] + length_squared[7] ); 02851 else 02852 current_shape = 0; 02853 02854 if( current_shape < shape ) { shape = current_shape; } 02855 } 02856 02857 // J(1,0,1) 02858 current_jacobian = -edges[4] % ( edges[5] * -edges[9] ); 02859 if( current_jacobian < jacobian ) jacobian = current_jacobian; 02860 02861 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02862 { det_sum += current_jacobian; } 02863 02864 if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) ) 02865 { 02866 if( length_squared[4] <= VERDICT_DBL_MIN || length_squared[5] <= VERDICT_DBL_MIN || 02867 length_squared[9] <= VERDICT_DBL_MIN ) 02868 { current_scaled_jacobian = VERDICT_DBL_MAX; } 02869 else 02870 { 02871 current_scaled_jacobian = 02872 current_jacobian / sqrt( length_squared[4] * length_squared[5] * length_squared[9] ); 02873 } 02874 if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian; 02875 } 02876 02877 if( metrics_request_flag & V_HEX_CONDITION ) 02878 { 02879 current_condition = condition_comp( -edges[4], edges[5], -edges[9] ); 02880 if( current_condition > condition ) { condition = current_condition; } 02881 } 02882 02883 if( metrics_request_flag & V_HEX_ODDY ) 02884 { 02885 current_oddy = oddy_comp( -edges[4], edges[5], -edges[9] ); 02886 if( current_oddy > oddy ) { oddy = current_oddy; } 02887 } 02888 02889 if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) ) 02890 { 02891 if( current_jacobian > VERDICT_DBL_MIN ) 02892 current_shape = 3 * pow( current_jacobian, two_thirds ) / 02893 ( length_squared[4] + length_squared[5] + length_squared[9] ); 02894 else 02895 current_shape = 0; 02896 02897 if( current_shape < shape ) { shape = current_shape; } 02898 } 02899 02900 // J(1,1,1) 02901 current_jacobian = -edges[5] % ( edges[6] * -edges[10] ); 02902 if( current_jacobian < jacobian ) jacobian = current_jacobian; 02903 02904 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02905 { det_sum += current_jacobian; } 02906 02907 if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) ) 02908 { 02909 if( length_squared[5] <= VERDICT_DBL_MIN || length_squared[6] <= VERDICT_DBL_MIN || 02910 length_squared[10] <= VERDICT_DBL_MIN ) 02911 { current_scaled_jacobian = VERDICT_DBL_MAX; } 02912 else 02913 { 02914 current_scaled_jacobian = 02915 current_jacobian / sqrt( length_squared[5] * length_squared[6] * length_squared[10] ); 02916 } 02917 if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian; 02918 } 02919 02920 if( metrics_request_flag & V_HEX_CONDITION ) 02921 { 02922 current_condition = condition_comp( -edges[5], edges[6], -edges[10] ); 02923 if( current_condition > condition ) { condition = current_condition; } 02924 } 02925 02926 if( metrics_request_flag & V_HEX_ODDY ) 02927 { 02928 current_oddy = oddy_comp( -edges[5], edges[6], -edges[10] ); 02929 if( current_oddy > oddy ) { oddy = current_oddy; } 02930 } 02931 02932 if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) ) 02933 { 02934 if( current_jacobian > VERDICT_DBL_MIN ) 02935 current_shape = 3 * pow( current_jacobian, two_thirds ) / 02936 ( length_squared[5] + length_squared[6] + length_squared[10] ); 02937 else 02938 current_shape = 0; 02939 02940 if( current_shape < shape ) { shape = current_shape; } 02941 } 02942 02943 // J(0,1,1) 02944 current_jacobian = -edges[6] % ( edges[7] * -edges[11] ); 02945 if( current_jacobian < jacobian ) jacobian = current_jacobian; 02946 02947 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02948 { det_sum += current_jacobian; } 02949 02950 if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) ) 02951 { 02952 if( length_squared[6] <= VERDICT_DBL_MIN || length_squared[7] <= VERDICT_DBL_MIN || 02953 length_squared[11] <= VERDICT_DBL_MIN ) 02954 { current_scaled_jacobian = VERDICT_DBL_MAX; } 02955 else 02956 { 02957 current_scaled_jacobian = 02958 current_jacobian / sqrt( length_squared[6] * length_squared[7] * length_squared[11] ); 02959 } 02960 if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian; 02961 } 02962 02963 if( metrics_request_flag & V_HEX_CONDITION ) 02964 { 02965 current_condition = condition_comp( -edges[6], edges[7], -edges[11] ); 02966 if( current_condition > condition ) { condition = current_condition; } 02967 } 02968 02969 if( metrics_request_flag & V_HEX_ODDY ) 02970 { 02971 current_oddy = oddy_comp( -edges[6], edges[7], -edges[11] ); 02972 if( current_oddy > oddy ) { oddy = current_oddy; } 02973 } 02974 02975 if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) ) 02976 { 02977 if( current_jacobian > VERDICT_DBL_MIN ) 02978 current_shape = 3 * pow( current_jacobian, two_thirds ) / 02979 ( length_squared[6] + length_squared[7] + length_squared[11] ); 02980 else 02981 current_shape = 0; 02982 02983 if( current_shape < shape ) { shape = current_shape; } 02984 } 02985 02986 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02987 { 02988 if( det_sum > VERDICT_DBL_MIN && rel_size_error != VERDICT_TRUE ) 02989 { 02990 double tau = det_sum / ( 8 * detw ); 02991 metric_vals->relative_size_squared = (double)VERDICT_MIN( tau * tau, 1.0 / tau / tau ); 02992 } 02993 else 02994 metric_vals->relative_size_squared = 0.0; 02995 } 02996 02997 // set values from above calculations 02998 if( metrics_request_flag & V_HEX_JACOBIAN ) metric_vals->jacobian = (double)jacobian; 02999 03000 if( metrics_request_flag & V_HEX_SCALED_JACOBIAN ) metric_vals->scaled_jacobian = (double)scaled_jacobian; 03001 03002 if( metrics_request_flag & V_HEX_CONDITION ) metric_vals->condition = (double)( condition / 3.0 ); 03003 03004 if( metrics_request_flag & V_HEX_SHEAR ) 03005 { 03006 if( shear < VERDICT_DBL_MIN ) // shear has range 0 to +1 03007 shear = 0; 03008 metric_vals->shear = (double)shear; 03009 } 03010 03011 if( metrics_request_flag & V_HEX_SHAPE ) metric_vals->shape = (double)shape; 03012 03013 if( metrics_request_flag & V_HEX_SHAPE_AND_SIZE ) 03014 metric_vals->shape_and_size = (double)( shape * metric_vals->relative_size_squared ); 03015 03016 if( metrics_request_flag & V_HEX_SHEAR_AND_SIZE ) 03017 metric_vals->shear_and_size = (double)( shear * metric_vals->relative_size_squared ); 03018 03019 if( metrics_request_flag & V_HEX_ODDY ) metric_vals->oddy = (double)oddy; 03020 03021 if( metrics_request_flag & V_HEX_STRETCH ) 03022 { 03023 static const double HEX_STRETCH_SCALE_FACTOR = sqrt( 3.0 ); 03024 double min_edge = length_squared[0]; 03025 for( int j = 1; j < 12; j++ ) 03026 min_edge = VERDICT_MIN( min_edge, length_squared[j] ); 03027 03028 double max_diag = diag_length( 1, coordinates ); 03029 03030 metric_vals->stretch = (double)( HEX_STRETCH_SCALE_FACTOR * ( safe_ratio( sqrt( min_edge ), max_diag ) ) ); 03031 } 03032 } 03033 03034 if( metrics_request_flag & V_HEX_DIAGONAL ) metric_vals->diagonal = v_hex_diagonal( num_nodes, coordinates ); 03035 if( metrics_request_flag & V_HEX_DIMENSION ) metric_vals->dimension = v_hex_dimension( num_nodes, coordinates ); 03036 if( metrics_request_flag & V_HEX_DISTORTION ) metric_vals->distortion = v_hex_distortion( num_nodes, coordinates ); 03037 03038 // take care of any overflow problems 03039 // max_edge_ratio 03040 if( metric_vals->max_edge_ratio > 0 ) 03041 metric_vals->max_edge_ratio = (double)VERDICT_MIN( metric_vals->max_edge_ratio, VERDICT_DBL_MAX ); 03042 else 03043 metric_vals->max_edge_ratio = (double)VERDICT_MAX( metric_vals->max_edge_ratio, -VERDICT_DBL_MAX ); 03044 03045 // skew 03046 if( metric_vals->skew > 0 ) 03047 metric_vals->skew = (double)VERDICT_MIN( metric_vals->skew, VERDICT_DBL_MAX ); 03048 else 03049 metric_vals->skew = (double)VERDICT_MAX( metric_vals->skew, -VERDICT_DBL_MAX ); 03050 03051 // taper 03052 if( metric_vals->taper > 0 ) 03053 metric_vals->taper = (double)VERDICT_MIN( metric_vals->taper, VERDICT_DBL_MAX ); 03054 else 03055 metric_vals->taper = (double)VERDICT_MAX( metric_vals->taper, -VERDICT_DBL_MAX ); 03056 03057 // volume 03058 if( metric_vals->volume > 0 ) 03059 metric_vals->volume = (double)VERDICT_MIN( metric_vals->volume, VERDICT_DBL_MAX ); 03060 else 03061 metric_vals->volume = (double)VERDICT_MAX( metric_vals->volume, -VERDICT_DBL_MAX ); 03062 03063 // stretch 03064 if( metric_vals->stretch > 0 ) 03065 metric_vals->stretch = (double)VERDICT_MIN( metric_vals->stretch, VERDICT_DBL_MAX ); 03066 else 03067 metric_vals->stretch = (double)VERDICT_MAX( metric_vals->stretch, -VERDICT_DBL_MAX ); 03068 03069 // diagonal 03070 if( metric_vals->diagonal > 0 ) 03071 metric_vals->diagonal = (double)VERDICT_MIN( metric_vals->diagonal, VERDICT_DBL_MAX ); 03072 else 03073 metric_vals->diagonal = (double)VERDICT_MAX( metric_vals->diagonal, -VERDICT_DBL_MAX ); 03074 03075 // dimension 03076 if( metric_vals->dimension > 0 ) 03077 metric_vals->dimension = (double)VERDICT_MIN( metric_vals->dimension, VERDICT_DBL_MAX ); 03078 else 03079 metric_vals->dimension = (double)VERDICT_MAX( metric_vals->dimension, -VERDICT_DBL_MAX ); 03080 03081 // oddy 03082 if( metric_vals->oddy > 0 ) 03083 metric_vals->oddy = (double)VERDICT_MIN( metric_vals->oddy, VERDICT_DBL_MAX ); 03084 else 03085 metric_vals->oddy = (double)VERDICT_MAX( metric_vals->oddy, -VERDICT_DBL_MAX ); 03086 03087 // condition 03088 if( metric_vals->condition > 0 ) 03089 metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX ); 03090 else 03091 metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX ); 03092 03093 // jacobian 03094 if( metric_vals->jacobian > 0 ) 03095 metric_vals->jacobian = (double)VERDICT_MIN( metric_vals->jacobian, VERDICT_DBL_MAX ); 03096 else 03097 metric_vals->jacobian = (double)VERDICT_MAX( metric_vals->jacobian, -VERDICT_DBL_MAX ); 03098 03099 // scaled_jacobian 03100 if( metric_vals->scaled_jacobian > 0 ) 03101 metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX ); 03102 else 03103 metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX ); 03104 03105 // shear 03106 if( metric_vals->shear > 0 ) 03107 metric_vals->shear = (double)VERDICT_MIN( metric_vals->shear, VERDICT_DBL_MAX ); 03108 else 03109 metric_vals->shear = (double)VERDICT_MAX( metric_vals->shear, -VERDICT_DBL_MAX ); 03110 03111 // shape 03112 if( metric_vals->shape > 0 ) 03113 metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX ); 03114 else 03115 metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX ); 03116 03117 // relative_size_squared 03118 if( metric_vals->relative_size_squared > 0 ) 03119 metric_vals->relative_size_squared = (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX ); 03120 else 03121 metric_vals->relative_size_squared = 03122 (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX ); 03123 03124 // shape_and_size 03125 if( metric_vals->shape_and_size > 0 ) 03126 metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX ); 03127 else 03128 metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX ); 03129 03130 // shear_and_size 03131 if( metric_vals->shear_and_size > 0 ) 03132 metric_vals->shear_and_size = (double)VERDICT_MIN( metric_vals->shear_and_size, VERDICT_DBL_MAX ); 03133 else 03134 metric_vals->shear_and_size = (double)VERDICT_MAX( metric_vals->shear_and_size, -VERDICT_DBL_MAX ); 03135 03136 // distortion 03137 if( metric_vals->distortion > 0 ) 03138 metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX ); 03139 else 03140 metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX ); 03141 03142 if( metrics_request_flag & V_HEX_MED_ASPECT_FROBENIUS ) 03143 metric_vals->med_aspect_frobenius = v_hex_med_aspect_frobenius( 8, coordinates ); 03144 03145 if( metrics_request_flag & V_HEX_EDGE_RATIO ) metric_vals->edge_ratio = v_hex_edge_ratio( 8, coordinates ); 03146 }