MOAB: Mesh Oriented datABase
(version 5.3.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 { 02555 metric_vals->skew = (double)VERDICT_DBL_MAX; 02556 } 02557 else 02558 { 02559 double skewx = fabs( vec1 % vec2 ); 02560 double skewy = fabs( vec1 % vec3 ); 02561 double skewz = fabs( vec2 % vec3 ); 02562 02563 metric_vals->skew = (double)( VERDICT_MAX( skewx, VERDICT_MAX( skewy, skewz ) ) ); 02564 } 02565 } 02566 02567 if( metrics_request_flag & V_HEX_TAPER ) 02568 { 02569 VerdictVector efg12 = calc_hex_efg( 12, node_pos ); 02570 VerdictVector efg13 = calc_hex_efg( 13, node_pos ); 02571 VerdictVector efg23 = calc_hex_efg( 23, node_pos ); 02572 02573 double taperx = fabs( safe_ratio( efg12.length(), VERDICT_MIN( efg1.length(), efg2.length() ) ) ); 02574 double tapery = fabs( safe_ratio( efg13.length(), VERDICT_MIN( efg1.length(), efg3.length() ) ) ); 02575 double taperz = fabs( safe_ratio( efg23.length(), VERDICT_MIN( efg2.length(), efg3.length() ) ) ); 02576 02577 metric_vals->taper = (double)VERDICT_MAX( taperx, VERDICT_MAX( tapery, taperz ) ); 02578 } 02579 } 02580 02581 if( metrics_request_flag & V_HEX_VOLUME ) { metric_vals->volume = v_hex_volume( 8, coordinates ); } 02582 02583 if( metrics_request_flag & 02584 ( V_HEX_JACOBIAN | V_HEX_SCALED_JACOBIAN | V_HEX_CONDITION | V_HEX_SHEAR | V_HEX_SHAPE | 02585 V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE | V_HEX_STRETCH ) ) 02586 { 02587 02588 static const double two_thirds = 2.0 / 3.0; 02589 VerdictVector edges[12]; 02590 // the length squares 02591 double length_squared[12]; 02592 // make vectors from coordinates 02593 make_hex_edges( coordinates, edges ); 02594 02595 // calculate the length squares if we need to 02596 if( metrics_request_flag & 02597 ( V_HEX_JACOBIAN | V_HEX_SHEAR | V_HEX_SCALED_JACOBIAN | V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE | 02598 V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHEAR_AND_SIZE | V_HEX_STRETCH ) ) 02599 make_edge_length_squares( edges, length_squared ); 02600 02601 double jacobian = VERDICT_DBL_MAX, scaled_jacobian = VERDICT_DBL_MAX, condition = 0.0, shear = 1.0, shape = 1.0, 02602 oddy = 0.0; 02603 double current_jacobian, current_scaled_jacobian, current_condition, current_shape, detw = 0, det_sum = 0, 02604 current_oddy; 02605 VerdictBoolean rel_size_error = VERDICT_FALSE; 02606 02607 VerdictVector xxi, xet, xze; 02608 02609 // get weights if we need based on average size of a hex 02610 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02611 { 02612 v_hex_get_weight( xxi, xet, xze ); 02613 detw = xxi % ( xet * xze ); 02614 if( detw < VERDICT_DBL_MIN ) rel_size_error = VERDICT_TRUE; 02615 } 02616 02617 xxi = edges[0] - edges[2] + edges[4] - edges[6]; 02618 xet = edges[1] - edges[3] + edges[5] - edges[7]; 02619 xze = edges[8] + edges[9] + edges[10] + edges[11]; 02620 02621 current_jacobian = xxi % ( xet * xze ) / 64.0; 02622 if( current_jacobian < jacobian ) jacobian = current_jacobian; 02623 02624 if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) ) 02625 { 02626 current_jacobian *= 64.0; 02627 current_scaled_jacobian = 02628 current_jacobian / sqrt( xxi.length_squared() * xet.length_squared() * xze.length_squared() ); 02629 if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian; 02630 } 02631 02632 if( metrics_request_flag & V_HEX_CONDITION ) 02633 { 02634 current_condition = condition_comp( xxi, xet, xze ); 02635 if( current_condition > condition ) { condition = current_condition; } 02636 } 02637 02638 if( metrics_request_flag & V_HEX_ODDY ) 02639 { 02640 current_oddy = oddy_comp( xxi, xet, xze ); 02641 if( current_oddy > oddy ) { oddy = current_oddy; } 02642 } 02643 02644 // J(0,0,0) 02645 current_jacobian = edges[0] % ( -edges[3] * edges[8] ); 02646 if( current_jacobian < jacobian ) jacobian = current_jacobian; 02647 02648 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02649 { 02650 det_sum += current_jacobian; 02651 } 02652 02653 if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) ) 02654 { 02655 if( length_squared[0] <= VERDICT_DBL_MIN || length_squared[3] <= VERDICT_DBL_MIN || 02656 length_squared[8] <= VERDICT_DBL_MIN ) 02657 { 02658 current_scaled_jacobian = VERDICT_DBL_MAX; 02659 } 02660 else 02661 { 02662 current_scaled_jacobian = 02663 current_jacobian / sqrt( length_squared[0] * length_squared[3] * length_squared[8] ); 02664 } 02665 if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian; 02666 } 02667 02668 if( metrics_request_flag & V_HEX_CONDITION ) 02669 { 02670 current_condition = condition_comp( edges[0], -edges[3], edges[8] ); 02671 if( current_condition > condition ) { condition = current_condition; } 02672 } 02673 02674 if( metrics_request_flag & V_HEX_ODDY ) 02675 { 02676 current_oddy = oddy_comp( edges[0], -edges[3], edges[8] ); 02677 if( current_oddy > oddy ) { oddy = current_oddy; } 02678 } 02679 02680 if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) ) 02681 { 02682 if( current_jacobian > VERDICT_DBL_MIN ) 02683 current_shape = 3 * pow( current_jacobian, two_thirds ) / 02684 ( length_squared[0] + length_squared[3] + length_squared[8] ); 02685 else 02686 current_shape = 0; 02687 02688 if( current_shape < shape ) { shape = current_shape; } 02689 } 02690 02691 // J(1,0,0) 02692 current_jacobian = edges[1] % ( -edges[0] * edges[9] ); 02693 if( current_jacobian < jacobian ) jacobian = current_jacobian; 02694 02695 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02696 { 02697 det_sum += current_jacobian; 02698 } 02699 02700 if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) ) 02701 { 02702 if( length_squared[1] <= VERDICT_DBL_MIN || length_squared[0] <= VERDICT_DBL_MIN || 02703 length_squared[9] <= VERDICT_DBL_MIN ) 02704 { 02705 current_scaled_jacobian = VERDICT_DBL_MAX; 02706 } 02707 else 02708 { 02709 current_scaled_jacobian = 02710 current_jacobian / sqrt( length_squared[1] * length_squared[0] * length_squared[9] ); 02711 } 02712 if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian; 02713 } 02714 02715 if( metrics_request_flag & V_HEX_CONDITION ) 02716 { 02717 current_condition = condition_comp( edges[1], -edges[0], edges[9] ); 02718 if( current_condition > condition ) { condition = current_condition; } 02719 } 02720 02721 if( metrics_request_flag & V_HEX_ODDY ) 02722 { 02723 current_oddy = oddy_comp( edges[1], -edges[0], edges[9] ); 02724 if( current_oddy > oddy ) { oddy = current_oddy; } 02725 } 02726 02727 if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) ) 02728 { 02729 if( current_jacobian > VERDICT_DBL_MIN ) 02730 current_shape = 3 * pow( current_jacobian, two_thirds ) / 02731 ( length_squared[1] + length_squared[0] + length_squared[9] ); 02732 else 02733 current_shape = 0; 02734 02735 if( current_shape < shape ) { shape = current_shape; } 02736 } 02737 02738 // J(1,1,0) 02739 current_jacobian = edges[2] % ( -edges[1] * edges[10] ); 02740 if( current_jacobian < jacobian ) jacobian = current_jacobian; 02741 02742 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02743 { 02744 det_sum += current_jacobian; 02745 } 02746 02747 if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) ) 02748 { 02749 if( length_squared[2] <= VERDICT_DBL_MIN || length_squared[1] <= VERDICT_DBL_MIN || 02750 length_squared[10] <= VERDICT_DBL_MIN ) 02751 { 02752 current_scaled_jacobian = VERDICT_DBL_MAX; 02753 } 02754 else 02755 { 02756 current_scaled_jacobian = 02757 current_jacobian / sqrt( length_squared[2] * length_squared[1] * length_squared[10] ); 02758 } 02759 if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian; 02760 } 02761 02762 if( metrics_request_flag & V_HEX_CONDITION ) 02763 { 02764 current_condition = condition_comp( edges[2], -edges[1], edges[10] ); 02765 if( current_condition > condition ) { condition = current_condition; } 02766 } 02767 02768 if( metrics_request_flag & V_HEX_ODDY ) 02769 { 02770 current_oddy = oddy_comp( edges[2], -edges[1], edges[10] ); 02771 if( current_oddy > oddy ) { oddy = current_oddy; } 02772 } 02773 02774 if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) ) 02775 { 02776 if( current_jacobian > VERDICT_DBL_MIN ) 02777 current_shape = 3 * pow( current_jacobian, two_thirds ) / 02778 ( length_squared[2] + length_squared[1] + length_squared[10] ); 02779 else 02780 current_shape = 0; 02781 02782 if( current_shape < shape ) { shape = current_shape; } 02783 } 02784 02785 // J(0,1,0) 02786 current_jacobian = edges[3] % ( -edges[2] * edges[11] ); 02787 if( current_jacobian < jacobian ) jacobian = current_jacobian; 02788 02789 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02790 { 02791 det_sum += current_jacobian; 02792 } 02793 02794 if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) ) 02795 { 02796 if( length_squared[3] <= VERDICT_DBL_MIN || length_squared[2] <= VERDICT_DBL_MIN || 02797 length_squared[11] <= VERDICT_DBL_MIN ) 02798 { 02799 current_scaled_jacobian = VERDICT_DBL_MAX; 02800 } 02801 else 02802 { 02803 current_scaled_jacobian = 02804 current_jacobian / sqrt( length_squared[3] * length_squared[2] * length_squared[11] ); 02805 } 02806 if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian; 02807 } 02808 02809 if( metrics_request_flag & V_HEX_CONDITION ) 02810 { 02811 current_condition = condition_comp( edges[3], -edges[2], edges[11] ); 02812 if( current_condition > condition ) { condition = current_condition; } 02813 } 02814 02815 if( metrics_request_flag & V_HEX_ODDY ) 02816 { 02817 current_oddy = oddy_comp( edges[3], -edges[2], edges[11] ); 02818 if( current_oddy > oddy ) { oddy = current_oddy; } 02819 } 02820 02821 if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) ) 02822 { 02823 if( current_jacobian > VERDICT_DBL_MIN ) 02824 current_shape = 3 * pow( current_jacobian, two_thirds ) / 02825 ( length_squared[3] + length_squared[2] + length_squared[11] ); 02826 else 02827 current_shape = 0; 02828 02829 if( current_shape < shape ) { shape = current_shape; } 02830 } 02831 02832 // J(0,0,1) 02833 current_jacobian = edges[4] % ( -edges[8] * -edges[7] ); 02834 if( current_jacobian < jacobian ) jacobian = current_jacobian; 02835 02836 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02837 { 02838 det_sum += current_jacobian; 02839 } 02840 02841 if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) ) 02842 { 02843 if( length_squared[4] <= VERDICT_DBL_MIN || length_squared[8] <= VERDICT_DBL_MIN || 02844 length_squared[7] <= VERDICT_DBL_MIN ) 02845 { 02846 current_scaled_jacobian = VERDICT_DBL_MAX; 02847 } 02848 else 02849 { 02850 current_scaled_jacobian = 02851 current_jacobian / sqrt( length_squared[4] * length_squared[8] * length_squared[7] ); 02852 } 02853 if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian; 02854 } 02855 02856 if( metrics_request_flag & V_HEX_CONDITION ) 02857 { 02858 current_condition = condition_comp( edges[4], -edges[8], -edges[7] ); 02859 if( current_condition > condition ) { condition = current_condition; } 02860 } 02861 02862 if( metrics_request_flag & V_HEX_ODDY ) 02863 { 02864 current_oddy = oddy_comp( edges[4], -edges[8], -edges[7] ); 02865 if( current_oddy > oddy ) { oddy = current_oddy; } 02866 } 02867 02868 if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) ) 02869 { 02870 if( current_jacobian > VERDICT_DBL_MIN ) 02871 current_shape = 3 * pow( current_jacobian, two_thirds ) / 02872 ( length_squared[4] + length_squared[8] + length_squared[7] ); 02873 else 02874 current_shape = 0; 02875 02876 if( current_shape < shape ) { shape = current_shape; } 02877 } 02878 02879 // J(1,0,1) 02880 current_jacobian = -edges[4] % ( edges[5] * -edges[9] ); 02881 if( current_jacobian < jacobian ) jacobian = current_jacobian; 02882 02883 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02884 { 02885 det_sum += current_jacobian; 02886 } 02887 02888 if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) ) 02889 { 02890 if( length_squared[4] <= VERDICT_DBL_MIN || length_squared[5] <= VERDICT_DBL_MIN || 02891 length_squared[9] <= VERDICT_DBL_MIN ) 02892 { 02893 current_scaled_jacobian = VERDICT_DBL_MAX; 02894 } 02895 else 02896 { 02897 current_scaled_jacobian = 02898 current_jacobian / sqrt( length_squared[4] * length_squared[5] * length_squared[9] ); 02899 } 02900 if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian; 02901 } 02902 02903 if( metrics_request_flag & V_HEX_CONDITION ) 02904 { 02905 current_condition = condition_comp( -edges[4], edges[5], -edges[9] ); 02906 if( current_condition > condition ) { condition = current_condition; } 02907 } 02908 02909 if( metrics_request_flag & V_HEX_ODDY ) 02910 { 02911 current_oddy = oddy_comp( -edges[4], edges[5], -edges[9] ); 02912 if( current_oddy > oddy ) { oddy = current_oddy; } 02913 } 02914 02915 if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) ) 02916 { 02917 if( current_jacobian > VERDICT_DBL_MIN ) 02918 current_shape = 3 * pow( current_jacobian, two_thirds ) / 02919 ( length_squared[4] + length_squared[5] + length_squared[9] ); 02920 else 02921 current_shape = 0; 02922 02923 if( current_shape < shape ) { shape = current_shape; } 02924 } 02925 02926 // J(1,1,1) 02927 current_jacobian = -edges[5] % ( edges[6] * -edges[10] ); 02928 if( current_jacobian < jacobian ) jacobian = current_jacobian; 02929 02930 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02931 { 02932 det_sum += current_jacobian; 02933 } 02934 02935 if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) ) 02936 { 02937 if( length_squared[5] <= VERDICT_DBL_MIN || length_squared[6] <= VERDICT_DBL_MIN || 02938 length_squared[10] <= VERDICT_DBL_MIN ) 02939 { 02940 current_scaled_jacobian = VERDICT_DBL_MAX; 02941 } 02942 else 02943 { 02944 current_scaled_jacobian = 02945 current_jacobian / sqrt( length_squared[5] * length_squared[6] * length_squared[10] ); 02946 } 02947 if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian; 02948 } 02949 02950 if( metrics_request_flag & V_HEX_CONDITION ) 02951 { 02952 current_condition = condition_comp( -edges[5], edges[6], -edges[10] ); 02953 if( current_condition > condition ) { condition = current_condition; } 02954 } 02955 02956 if( metrics_request_flag & V_HEX_ODDY ) 02957 { 02958 current_oddy = oddy_comp( -edges[5], edges[6], -edges[10] ); 02959 if( current_oddy > oddy ) { oddy = current_oddy; } 02960 } 02961 02962 if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) ) 02963 { 02964 if( current_jacobian > VERDICT_DBL_MIN ) 02965 current_shape = 3 * pow( current_jacobian, two_thirds ) / 02966 ( length_squared[5] + length_squared[6] + length_squared[10] ); 02967 else 02968 current_shape = 0; 02969 02970 if( current_shape < shape ) { shape = current_shape; } 02971 } 02972 02973 // J(0,1,1) 02974 current_jacobian = -edges[6] % ( edges[7] * -edges[11] ); 02975 if( current_jacobian < jacobian ) jacobian = current_jacobian; 02976 02977 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 02978 { 02979 det_sum += current_jacobian; 02980 } 02981 02982 if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) ) 02983 { 02984 if( length_squared[6] <= VERDICT_DBL_MIN || length_squared[7] <= VERDICT_DBL_MIN || 02985 length_squared[11] <= VERDICT_DBL_MIN ) 02986 { 02987 current_scaled_jacobian = VERDICT_DBL_MAX; 02988 } 02989 else 02990 { 02991 current_scaled_jacobian = 02992 current_jacobian / sqrt( length_squared[6] * length_squared[7] * length_squared[11] ); 02993 } 02994 if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian; 02995 } 02996 02997 if( metrics_request_flag & V_HEX_CONDITION ) 02998 { 02999 current_condition = condition_comp( -edges[6], edges[7], -edges[11] ); 03000 if( current_condition > condition ) { condition = current_condition; } 03001 } 03002 03003 if( metrics_request_flag & V_HEX_ODDY ) 03004 { 03005 current_oddy = oddy_comp( -edges[6], edges[7], -edges[11] ); 03006 if( current_oddy > oddy ) { oddy = current_oddy; } 03007 } 03008 03009 if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) ) 03010 { 03011 if( current_jacobian > VERDICT_DBL_MIN ) 03012 current_shape = 3 * pow( current_jacobian, two_thirds ) / 03013 ( length_squared[6] + length_squared[7] + length_squared[11] ); 03014 else 03015 current_shape = 0; 03016 03017 if( current_shape < shape ) { shape = current_shape; } 03018 } 03019 03020 if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) ) 03021 { 03022 if( det_sum > VERDICT_DBL_MIN && rel_size_error != VERDICT_TRUE ) 03023 { 03024 double tau = det_sum / ( 8 * detw ); 03025 metric_vals->relative_size_squared = (double)VERDICT_MIN( tau * tau, 1.0 / tau / tau ); 03026 } 03027 else 03028 metric_vals->relative_size_squared = 0.0; 03029 } 03030 03031 // set values from above calculations 03032 if( metrics_request_flag & V_HEX_JACOBIAN ) metric_vals->jacobian = (double)jacobian; 03033 03034 if( metrics_request_flag & V_HEX_SCALED_JACOBIAN ) metric_vals->scaled_jacobian = (double)scaled_jacobian; 03035 03036 if( metrics_request_flag & V_HEX_CONDITION ) metric_vals->condition = (double)( condition / 3.0 ); 03037 03038 if( metrics_request_flag & V_HEX_SHEAR ) 03039 { 03040 if( shear < VERDICT_DBL_MIN ) // shear has range 0 to +1 03041 shear = 0; 03042 metric_vals->shear = (double)shear; 03043 } 03044 03045 if( metrics_request_flag & V_HEX_SHAPE ) metric_vals->shape = (double)shape; 03046 03047 if( metrics_request_flag & V_HEX_SHAPE_AND_SIZE ) 03048 metric_vals->shape_and_size = (double)( shape * metric_vals->relative_size_squared ); 03049 03050 if( metrics_request_flag & V_HEX_SHEAR_AND_SIZE ) 03051 metric_vals->shear_and_size = (double)( shear * metric_vals->relative_size_squared ); 03052 03053 if( metrics_request_flag & V_HEX_ODDY ) metric_vals->oddy = (double)oddy; 03054 03055 if( metrics_request_flag & V_HEX_STRETCH ) 03056 { 03057 static const double HEX_STRETCH_SCALE_FACTOR = sqrt( 3.0 ); 03058 double min_edge = length_squared[0]; 03059 for( int j = 1; j < 12; j++ ) 03060 min_edge = VERDICT_MIN( min_edge, length_squared[j] ); 03061 03062 double max_diag = diag_length( 1, coordinates ); 03063 03064 metric_vals->stretch = (double)( HEX_STRETCH_SCALE_FACTOR * ( safe_ratio( sqrt( min_edge ), max_diag ) ) ); 03065 } 03066 } 03067 03068 if( metrics_request_flag & V_HEX_DIAGONAL ) metric_vals->diagonal = v_hex_diagonal( num_nodes, coordinates ); 03069 if( metrics_request_flag & V_HEX_DIMENSION ) metric_vals->dimension = v_hex_dimension( num_nodes, coordinates ); 03070 if( metrics_request_flag & V_HEX_DISTORTION ) metric_vals->distortion = v_hex_distortion( num_nodes, coordinates ); 03071 03072 // take care of any overflow problems 03073 // max_edge_ratio 03074 if( metric_vals->max_edge_ratio > 0 ) 03075 metric_vals->max_edge_ratio = (double)VERDICT_MIN( metric_vals->max_edge_ratio, VERDICT_DBL_MAX ); 03076 else 03077 metric_vals->max_edge_ratio = (double)VERDICT_MAX( metric_vals->max_edge_ratio, -VERDICT_DBL_MAX ); 03078 03079 // skew 03080 if( metric_vals->skew > 0 ) 03081 metric_vals->skew = (double)VERDICT_MIN( metric_vals->skew, VERDICT_DBL_MAX ); 03082 else 03083 metric_vals->skew = (double)VERDICT_MAX( metric_vals->skew, -VERDICT_DBL_MAX ); 03084 03085 // taper 03086 if( metric_vals->taper > 0 ) 03087 metric_vals->taper = (double)VERDICT_MIN( metric_vals->taper, VERDICT_DBL_MAX ); 03088 else 03089 metric_vals->taper = (double)VERDICT_MAX( metric_vals->taper, -VERDICT_DBL_MAX ); 03090 03091 // volume 03092 if( metric_vals->volume > 0 ) 03093 metric_vals->volume = (double)VERDICT_MIN( metric_vals->volume, VERDICT_DBL_MAX ); 03094 else 03095 metric_vals->volume = (double)VERDICT_MAX( metric_vals->volume, -VERDICT_DBL_MAX ); 03096 03097 // stretch 03098 if( metric_vals->stretch > 0 ) 03099 metric_vals->stretch = (double)VERDICT_MIN( metric_vals->stretch, VERDICT_DBL_MAX ); 03100 else 03101 metric_vals->stretch = (double)VERDICT_MAX( metric_vals->stretch, -VERDICT_DBL_MAX ); 03102 03103 // diagonal 03104 if( metric_vals->diagonal > 0 ) 03105 metric_vals->diagonal = (double)VERDICT_MIN( metric_vals->diagonal, VERDICT_DBL_MAX ); 03106 else 03107 metric_vals->diagonal = (double)VERDICT_MAX( metric_vals->diagonal, -VERDICT_DBL_MAX ); 03108 03109 // dimension 03110 if( metric_vals->dimension > 0 ) 03111 metric_vals->dimension = (double)VERDICT_MIN( metric_vals->dimension, VERDICT_DBL_MAX ); 03112 else 03113 metric_vals->dimension = (double)VERDICT_MAX( metric_vals->dimension, -VERDICT_DBL_MAX ); 03114 03115 // oddy 03116 if( metric_vals->oddy > 0 ) 03117 metric_vals->oddy = (double)VERDICT_MIN( metric_vals->oddy, VERDICT_DBL_MAX ); 03118 else 03119 metric_vals->oddy = (double)VERDICT_MAX( metric_vals->oddy, -VERDICT_DBL_MAX ); 03120 03121 // condition 03122 if( metric_vals->condition > 0 ) 03123 metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX ); 03124 else 03125 metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX ); 03126 03127 // jacobian 03128 if( metric_vals->jacobian > 0 ) 03129 metric_vals->jacobian = (double)VERDICT_MIN( metric_vals->jacobian, VERDICT_DBL_MAX ); 03130 else 03131 metric_vals->jacobian = (double)VERDICT_MAX( metric_vals->jacobian, -VERDICT_DBL_MAX ); 03132 03133 // scaled_jacobian 03134 if( metric_vals->scaled_jacobian > 0 ) 03135 metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX ); 03136 else 03137 metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX ); 03138 03139 // shear 03140 if( metric_vals->shear > 0 ) 03141 metric_vals->shear = (double)VERDICT_MIN( metric_vals->shear, VERDICT_DBL_MAX ); 03142 else 03143 metric_vals->shear = (double)VERDICT_MAX( metric_vals->shear, -VERDICT_DBL_MAX ); 03144 03145 // shape 03146 if( metric_vals->shape > 0 ) 03147 metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX ); 03148 else 03149 metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX ); 03150 03151 // relative_size_squared 03152 if( metric_vals->relative_size_squared > 0 ) 03153 metric_vals->relative_size_squared = (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX ); 03154 else 03155 metric_vals->relative_size_squared = 03156 (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX ); 03157 03158 // shape_and_size 03159 if( metric_vals->shape_and_size > 0 ) 03160 metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX ); 03161 else 03162 metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX ); 03163 03164 // shear_and_size 03165 if( metric_vals->shear_and_size > 0 ) 03166 metric_vals->shear_and_size = (double)VERDICT_MIN( metric_vals->shear_and_size, VERDICT_DBL_MAX ); 03167 else 03168 metric_vals->shear_and_size = (double)VERDICT_MAX( metric_vals->shear_and_size, -VERDICT_DBL_MAX ); 03169 03170 // distortion 03171 if( metric_vals->distortion > 0 ) 03172 metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX ); 03173 else 03174 metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX ); 03175 03176 if( metrics_request_flag & V_HEX_MED_ASPECT_FROBENIUS ) 03177 metric_vals->med_aspect_frobenius = v_hex_med_aspect_frobenius( 8, coordinates ); 03178 03179 if( metrics_request_flag & V_HEX_EDGE_RATIO ) metric_vals->edge_ratio = v_hex_edge_ratio( 8, coordinates ); 03180 }