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