 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 /*=========================================================================
2
3 Module: $RCSfile: V_TriMetric.cpp,v $
4
5 Copyright (c) 2007 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 * TriMetric.cpp contains quality calculations for Tris
18 *
19 * This file is part of VERDICT
20 *
21 */
22
23 #define VERDICT_EXPORTS
24
25 #include "moab/verdict.h"
26 #include "verdict_defines.hpp"
27 #include "V_GaussIntegration.hpp"
28 #include "VerdictVector.hpp"
29 #include <memory.h>
30 #include <cstddef>
31
32 // the average area of a tri
33 static double verdict_tri_size = 0;
34 static ComputeNormal compute_normal = NULL;
35
36 /*!
37 get weights based on the average area of a set of
38 tris
39 */
40 static int v_tri_get_weight( double& m11, double& m21, double& m12, double& m22 )
41 {
42 static const double rootOf3 = sqrt( 3.0 );
43 m11 = 1;
44 m21 = 0;
45 m12 = 0.5;
46 m22 = 0.5 * rootOf3;
47 double scale = sqrt( 2.0 * verdict_tri_size / ( m11 * m22 - m21 * m12 ) );
48
49 m11 *= scale;
50 m21 *= scale;
51 m12 *= scale;
52 m22 *= scale;
53
54 return 1;
55 }
56
57 /*! sets the average area of a tri */
58 C_FUNC_DEF void v_set_tri_size( double size )
59 {
60 verdict_tri_size = size;
61 }
62
63 C_FUNC_DEF void v_set_tri_normal_func( ComputeNormal func )
64 {
65 compute_normal = func;
66 }
67
68 /*!
69 the edge ratio of a triangle
70
71 NB (P. Pebay 01/14/07):
72 Hmax / Hmin where Hmax and Hmin are respectively the maximum and the
73 minimum edge lengths
74
75 */
76 C_FUNC_DEF double v_tri_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
77 {
78
79 // three vectors for each side
80 VerdictVector a( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
81 coordinates[1][2] - coordinates[0][2] );
82
83 VerdictVector b( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
84 coordinates[2][2] - coordinates[1][2] );
85
86 VerdictVector c( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
87 coordinates[0][2] - coordinates[2][2] );
88
89 double a2 = a.length_squared();
90 double b2 = b.length_squared();
91 double c2 = c.length_squared();
92
93 double m2, M2;
94 if( a2 < b2 )
95 {
96 if( b2 < c2 )
97 {
98 m2 = a2;
99 M2 = c2;
100 }
101 else // b2 <= a2
102 {
103 if( a2 < c2 )
104 {
105 m2 = a2;
106 M2 = b2;
107 }
108 else // c2 <= a2
109 {
110 m2 = c2;
111 M2 = b2;
112 }
113 }
114 }
115 else // b2 <= a2
116 {
117 if( a2 < c2 )
118 {
119 m2 = b2;
120 M2 = c2;
121 }
122 else // c2 <= a2
123 {
124 if( b2 < c2 )
125 {
126 m2 = b2;
127 M2 = a2;
128 }
129 else // c2 <= b2
130 {
131 m2 = c2;
132 M2 = a2;
133 }
134 }
135 }
136
137 if( m2 < VERDICT_DBL_MIN )
138 return (double)VERDICT_DBL_MAX;
139 else
140 {
141 double edge_ratio;
142 edge_ratio = sqrt( M2 / m2 );
143
144 if( edge_ratio > 0 ) return (double)VERDICT_MIN( edge_ratio, VERDICT_DBL_MAX );
145 return (double)VERDICT_MAX( edge_ratio, -VERDICT_DBL_MAX );
146 }
147 }
148
149 /*!
150 the aspect ratio of a triangle
151
152 NB (P. Pebay 01/14/07):
153 Hmax / ( 2.0 * sqrt(3.0) * IR) where Hmax is the maximum edge length
154 and IR is the inradius
155
156 note that previous incarnations of verdict used "v_tri_aspect_ratio" to denote
157 what is now called "v_tri_aspect_frobenius"
158
159 */
160 C_FUNC_DEF double v_tri_aspect_ratio( int /*num_nodes*/, double coordinates[][3] )
161 {
162 static const double normal_coeff = sqrt( 3. ) / 6.;
163
164 // three vectors for each side
165 VerdictVector a( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
166 coordinates[1][2] - coordinates[0][2] );
167
168 VerdictVector b( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
169 coordinates[2][2] - coordinates[1][2] );
170
171 VerdictVector c( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
172 coordinates[0][2] - coordinates[2][2] );
173
174 double a1 = a.length();
175 double b1 = b.length();
176 double c1 = c.length();
177
178 double hm = a1 > b1 ? a1 : b1;
179 hm = hm > c1 ? hm : c1;
180
181 VerdictVector ab = a * b;
182 double denominator = ab.length();
183
184 if( denominator < VERDICT_DBL_MIN )
185 return (double)VERDICT_DBL_MAX;
186 else
187 {
188 double aspect_ratio;
189 aspect_ratio = normal_coeff * hm * ( a1 + b1 + c1 ) / denominator;
190
191 if( aspect_ratio > 0 ) return (double)VERDICT_MIN( aspect_ratio, VERDICT_DBL_MAX );
192 return (double)VERDICT_MAX( aspect_ratio, -VERDICT_DBL_MAX );
193 }
194 }
195
196 /*!
197 the radius ratio of a triangle
198
199 NB (P. Pebay 01/13/07):
200 CR / (3.0*IR) where CR is the circumradius and IR is the inradius
201
202 this quality metric is also known to VERDICT, for tetrahedral elements only,
203 a the "aspect beta"
204
205 */
206 C_FUNC_DEF double v_tri_radius_ratio( int /*num_nodes*/, double coordinates[][3] )
207 {
208
209 // three vectors for each side
210 VerdictVector a( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
211 coordinates[1][2] - coordinates[0][2] );
212
213 VerdictVector b( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
214 coordinates[2][2] - coordinates[1][2] );
215
216 VerdictVector c( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
217 coordinates[0][2] - coordinates[2][2] );
218
219 double a2 = a.length_squared();
220 double b2 = b.length_squared();
221 double c2 = c.length_squared();
222
223 VerdictVector ab = a * b;
224 double denominator = ab.length_squared();
225
226 if( denominator < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
227
228 double radius_ratio;
229 radius_ratio = .25 * a2 * b2 * c2 * ( a2 + b2 + c2 ) / denominator;
230
231 if( radius_ratio > 0 ) return (double)VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
232 return (double)VERDICT_MAX( radius_ratio, -VERDICT_DBL_MAX );
233 }
234
235 /*!
236 the Frobenius aspect of a tri
237
238 srms^2/(2 * sqrt(3.0) * area)
239 where srms^2 is sum of the lengths squared
240
241 NB (P. Pebay 01/14/07):
242 this method was called "aspect ratio" in earlier incarnations of VERDICT
243
244 */
245
246 C_FUNC_DEF double v_tri_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
247 {
248 static const double two_times_root_of_3 = 2 * sqrt( 3.0 );
249
250 // three vectors for each side
251 VerdictVector side1( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
252 coordinates[1][2] - coordinates[0][2] );
253
254 VerdictVector side2( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
255 coordinates[2][2] - coordinates[1][2] );
256
257 VerdictVector side3( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
258 coordinates[0][2] - coordinates[2][2] );
259
260 // sum the lengths squared of each side
261 double srms = ( side1.length_squared() + side2.length_squared() + side3.length_squared() );
262
263 // find two times the area of the triangle by cross product
264 double areaX2 = ( ( side1 * ( -side3 ) ).length() );
265
266 if( areaX2 == 0.0 ) return (double)VERDICT_DBL_MAX;
267
268 double aspect = (double)( srms / ( two_times_root_of_3 * ( areaX2 ) ) );
269 if( aspect > 0 ) return (double)VERDICT_MIN( aspect, VERDICT_DBL_MAX );
270 return (double)VERDICT_MAX( aspect, -VERDICT_DBL_MAX );
271 }
272
273 /*!
274 The area of a tri
275
276 0.5 * jacobian at a node
277 */
278 C_FUNC_DEF double v_tri_area( int /*num_nodes*/, double coordinates[][3] )
279 {
280 // two vectors for two sides
281 VerdictVector side1( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
282 coordinates[1][2] - coordinates[0][2] );
283
284 VerdictVector side3( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
285 coordinates[2][2] - coordinates[0][2] );
286
287 // the cross product of the two vectors representing two sides of the
288 // triangle
289 VerdictVector tmp = side1 * side3;
290
291 // return the magnitude of the vector divided by two
292 double area = 0.5 * tmp.length();
293 if( area > 0 ) return (double)VERDICT_MIN( area, VERDICT_DBL_MAX );
294 return (double)VERDICT_MAX( area, -VERDICT_DBL_MAX );
295 }
296
297 /*!
298 The minimum angle of a tri
299
300 The minimum angle of a tri is the minimum angle between
301 two adjacents sides out of all three corners of the triangle.
302 */
303 C_FUNC_DEF double v_tri_minimum_angle( int /*num_nodes*/, double coordinates[][3] )
304 {
305
306 // vectors for all the sides
307 VerdictVector sides[4];
308 sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
309 coordinates[1][2] - coordinates[0][2] );
310 sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
311 coordinates[2][2] - coordinates[1][2] );
312 sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
313 coordinates[2][2] - coordinates[0][2] );
314
315 // in case we need to find the interior angle
316 // between sides 0 and 1
317 sides[3] = -sides[1];
318
319 // calculate the lengths squared of the sides
320 double sides_lengths[3];
321 sides_lengths[0] = sides[0].length_squared();
322 sides_lengths[1] = sides[1].length_squared();
323 sides_lengths[2] = sides[2].length_squared();
324
325 if( sides_lengths[0] == 0.0 || sides_lengths[1] == 0.0 || sides_lengths[2] == 0.0 ) return 0.0;
326
327 // using the law of sines, we know that the minimum
328 // angle is opposite of the shortest side
329
330 // find the shortest side
331 int short_side = 0;
332 if( sides_lengths[1] < sides_lengths[0] ) short_side = 1;
333 if( sides_lengths[2] < sides_lengths[short_side] ) short_side = 2;
334
335 // from the shortest side, calculate the angle of the
336 // opposite angle
337 double min_angle;
338 if( short_side == 0 )
339 {
340 min_angle = sides[2].interior_angle( sides[1] );
341 }
342 else if( short_side == 1 )
343 {
344 min_angle = sides[0].interior_angle( sides[2] );
345 }
346 else
347 {
348 min_angle = sides[0].interior_angle( sides[3] );
349 }
350
351 if( min_angle > 0 ) return (double)VERDICT_MIN( min_angle, VERDICT_DBL_MAX );
352 return (double)VERDICT_MAX( min_angle, -VERDICT_DBL_MAX );
353 }
354
355 /*!
356 The maximum angle of a tri
357
358 The maximum angle of a tri is the maximum angle between
359 two adjacents sides out of all three corners of the triangle.
360 */
361 C_FUNC_DEF double v_tri_maximum_angle( int /*num_nodes*/, double coordinates[][3] )
362 {
363
364 // vectors for all the sides
365 VerdictVector sides[4];
366 sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
367 coordinates[1][2] - coordinates[0][2] );
368 sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
369 coordinates[2][2] - coordinates[1][2] );
370 sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
371 coordinates[2][2] - coordinates[0][2] );
372
373 // in case we need to find the interior angle
374 // between sides 0 and 1
375 sides[3] = -sides[1];
376
377 // calculate the lengths squared of the sides
378 double sides_lengths[3];
379 sides_lengths[0] = sides[0].length_squared();
380 sides_lengths[1] = sides[1].length_squared();
381 sides_lengths[2] = sides[2].length_squared();
382
383 if( sides_lengths[0] == 0.0 || sides_lengths[1] == 0.0 || sides_lengths[2] == 0.0 )
384 {
385 return 0.0;
386 }
387
388 // using the law of sines, we know that the maximum
389 // angle is opposite of the longest side
390
391 // find the longest side
392 int short_side = 0;
393 if( sides_lengths[1] > sides_lengths[0] ) short_side = 1;
394 if( sides_lengths[2] > sides_lengths[short_side] ) short_side = 2;
395
396 // from the longest side, calculate the angle of the
397 // opposite angle
398 double max_angle;
399 if( short_side == 0 )
400 {
401 max_angle = sides[2].interior_angle( sides[1] );
402 }
403 else if( short_side == 1 )
404 {
405 max_angle = sides[0].interior_angle( sides[2] );
406 }
407 else
408 {
409 max_angle = sides[0].interior_angle( sides[3] );
410 }
411
412 if( max_angle > 0 ) return (double)VERDICT_MIN( max_angle, VERDICT_DBL_MAX );
413 return (double)VERDICT_MAX( max_angle, -VERDICT_DBL_MAX );
414 }
415
416 /*!
417 The condition of a tri
418
419 Condition number of the jacobian matrix at any corner
420 */
421 C_FUNC_DEF double v_tri_condition( int /*num_nodes*/, double coordinates[][3] )
422 {
423 static const double rt3 = sqrt( 3.0 );
424
425 VerdictVector v1( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
426 coordinates[1][2] - coordinates[0][2] );
427
428 VerdictVector v2( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
429 coordinates[2][2] - coordinates[0][2] );
430
431 VerdictVector tri_normal = v1 * v2;
432 double areax2 = tri_normal.length();
433
434 if( areax2 == 0.0 ) return (double)VERDICT_DBL_MAX;
435
436 double condition = (double)( ( ( v1 % v1 ) + ( v2 % v2 ) - ( v1 % v2 ) ) / ( areax2 * rt3 ) );
437
438 // check for inverted if we have access to the normal
439 if( compute_normal )
440 {
441 // center of tri
442 double point[3], surf_normal[3];
443 point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
444 point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
445 point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
446
447 // dot product
448 compute_normal( point, surf_normal );
449 if( ( tri_normal.x() * surf_normal[0] + tri_normal.y() * surf_normal[1] + tri_normal.z() * surf_normal[2] ) <
450 0 )
451 return (double)VERDICT_DBL_MAX;
452 }
453 return (double)VERDICT_MIN( condition, VERDICT_DBL_MAX );
454 }
455
456 /*!
457 The scaled jacobian of a tri
458
459 minimum of the jacobian divided by the lengths of 2 edge vectors
460 */
461 C_FUNC_DEF double v_tri_scaled_jacobian( int /*num_nodes*/, double coordinates[][3] )
462 {
463 static const double detw = 2. / sqrt( 3.0 );
464 VerdictVector first, second;
465 double jacobian;
466
467 VerdictVector edge[3];
468 edge[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
469 coordinates[1][2] - coordinates[0][2] );
470
471 edge[1].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
472 coordinates[2][2] - coordinates[0][2] );
473
474 edge[2].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
475 coordinates[2][2] - coordinates[1][2] );
476 first = edge[1] - edge[0];
477 second = edge[2] - edge[0];
478
479 VerdictVector cross = first * second;
480 jacobian = cross.length();
481
482 double max_edge_length_product;
483 max_edge_length_product =
484 VERDICT_MAX( edge[0].length() * edge[1].length(),
485 VERDICT_MAX( edge[1].length() * edge[2].length(), edge[0].length() * edge[2].length() ) );
486
487 if( max_edge_length_product < VERDICT_DBL_MIN ) return (double)0.0;
488
489 jacobian *= detw;
490 jacobian /= max_edge_length_product;
491
492 if( compute_normal )
493 {
494 // center of tri
495 double point[3], surf_normal[3];
496 point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
497 point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
498 point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
499
500 // dot product
501 compute_normal( point, surf_normal );
502 if( ( cross.x() * surf_normal[0] + cross.y() * surf_normal[1] + cross.z() * surf_normal[2] ) < 0 )
503 jacobian *= -1;
504 }
505
506 if( jacobian > 0 ) return (double)VERDICT_MIN( jacobian, VERDICT_DBL_MAX );
507 return (double)VERDICT_MAX( jacobian, -VERDICT_DBL_MAX );
508 }
509
510 /*!
511 The shape of a tri
512
513 2 / condition number of weighted jacobian matrix
514 */
515 C_FUNC_DEF double v_tri_shape( int num_nodes, double coordinates[][3] )
516 {
517 double condition = v_tri_condition( num_nodes, coordinates );
518
519 double shape;
520 if( condition <= VERDICT_DBL_MIN )
521 shape = VERDICT_DBL_MAX;
522 else
523 shape = ( 1 / condition );
524
525 if( shape > 0 ) return (double)VERDICT_MIN( shape, VERDICT_DBL_MAX );
526 return (double)VERDICT_MAX( shape, -VERDICT_DBL_MAX );
527 }
528
529 /*!
530 The relative size of a tri
531
532 Min(J,1/J) where J is the determinant of the weighted jacobian matrix.
533 */
534 C_FUNC_DEF double v_tri_relative_size_squared( int /*num_nodes*/, double coordinates[][3] )
535 {
536 double w11, w21, w12, w22;
537
538 VerdictVector xxi, xet, tri_normal;
539
540 v_tri_get_weight( w11, w21, w12, w22 );
541
542 double detw = determinant( w11, w21, w12, w22 );
543
544 if( detw == 0.0 ) return 0.0;
545
546 xxi.set( coordinates[0][0] - coordinates[1][0], coordinates[0][1] - coordinates[1][1],
547 coordinates[0][2] - coordinates[1][2] );
548
549 xet.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
550 coordinates[0][2] - coordinates[2][2] );
551
552 tri_normal = xxi * xet;
553
554 double deta = tri_normal.length();
555 if( deta == 0.0 || detw == 0.0 ) return 0.0;
556
557 double size = pow( deta / detw, 2 );
558
559 double rel_size = VERDICT_MIN( size, 1.0 / size );
560
561 if( rel_size > 0 ) return (double)VERDICT_MIN( rel_size, VERDICT_DBL_MAX );
562 return (double)VERDICT_MAX( rel_size, -VERDICT_DBL_MAX );
563 }
564
565 /*!
566 The shape and size of a tri
567
568 Product of the Shape and Relative Size
569 */
570 C_FUNC_DEF double v_tri_shape_and_size( int num_nodes, double coordinates[][3] )
571 {
572 double size, shape;
573
574 size = v_tri_relative_size_squared( num_nodes, coordinates );
575 shape = v_tri_shape( num_nodes, coordinates );
576
577 double shape_and_size = size * shape;
578
579 if( shape_and_size > 0 ) return (double)VERDICT_MIN( shape_and_size, VERDICT_DBL_MAX );
580 return (double)VERDICT_MAX( shape_and_size, -VERDICT_DBL_MAX );
581 }
582
583 /*!
584 The distortion of a tri
585
586 TODO: make a short definition of the distortion and comment below
587 */
588 C_FUNC_DEF double v_tri_distortion( int num_nodes, double coordinates[][3] )
589 {
590
591 double distortion;
592 int total_number_of_gauss_points = 0;
593 VerdictVector aa, bb, cc, normal_at_point, xin;
594 double element_area = 0.;
595
596 aa.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
597 coordinates[1][2] - coordinates[0][2] );
598
599 bb.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
600 coordinates[2][2] - coordinates[0][2] );
601
602 VerdictVector tri_normal = aa * bb;
603
604 int number_of_gauss_points = 0;
605 if( num_nodes == 3 )
606 {
607 distortion = 1.0;
608 return (double)distortion;
609 }
610
611 else if( num_nodes == 6 )
612 {
613 total_number_of_gauss_points = 6;
614 number_of_gauss_points = 6;
615 }
616
617 distortion = VERDICT_DBL_MAX;
618 double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
619 double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
620 double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
621 double weight[maxTotalNumberGaussPoints];
622
623 // create an object of GaussIntegration
624 int number_dims = 2;
625 int is_tri = 1;
626 GaussIntegration::initialize( number_of_gauss_points, num_nodes, number_dims, is_tri );
627 GaussIntegration::calculate_shape_function_2d_tri();
628 GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], weight );
629
630 // calculate element area
631 int ife, ja;
632 for( ife = 0; ife < total_number_of_gauss_points; ife++ )
633 {
634 aa.set( 0.0, 0.0, 0.0 );
635 bb.set( 0.0, 0.0, 0.0 );
636
637 for( ja = 0; ja < num_nodes; ja++ )
638 {
639 xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
640 aa += dndy1[ife][ja] * xin;
641 bb += dndy2[ife][ja] * xin;
642 }
643 normal_at_point = aa * bb;
644 double jacobian = normal_at_point.length();
645 element_area += weight[ife] * jacobian;
646 }
647
648 element_area *= 0.8660254;
649 double dndy1_at_node[maxNumberNodes][maxNumberNodes];
650 double dndy2_at_node[maxNumberNodes][maxNumberNodes];
651
652 GaussIntegration::calculate_derivative_at_nodes_2d_tri( dndy1_at_node, dndy2_at_node );
653
654 VerdictVector normal_at_nodes[7];
655
656 // evaluate normal at nodes and distortion values at nodes
657 int jai = 0;
658 for( ja = 0; ja < num_nodes; ja++ )
659 {
660 aa.set( 0.0, 0.0, 0.0 );
661 bb.set( 0.0, 0.0, 0.0 );
662 for( jai = 0; jai < num_nodes; jai++ )
663 {
664 xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
665 aa += dndy1_at_node[ja][jai] * xin;
666 bb += dndy2_at_node[ja][jai] * xin;
667 }
668 normal_at_nodes[ja] = aa * bb;
669 normal_at_nodes[ja].normalize();
670 }
671
672 // determine if element is flat
673 bool flat_element = true;
674 double dot_product;
675
676 for( ja = 0; ja < num_nodes; ja++ )
677 {
678 dot_product = normal_at_nodes[0] % normal_at_nodes[ja];
679 if( fabs( dot_product ) < 0.99 )
680 {
681 flat_element = false;
682 break;
683 }
684 }
685
686 // take into consideration of the thickness of the element
687 double thickness, thickness_gauss;
688 double distrt;
689 // get_tri_thickness(tri, element_area, thickness );
690 thickness = 0.001 * sqrt( element_area );
691
692 // set thickness gauss point location
693 double zl = 0.5773502691896;
694 if( flat_element ) zl = 0.0;
695
696 int no_gauss_pts_z = ( flat_element ) ? 1 : 2;
697 double thickness_z;
698
699 // loop on integration points
700 int igz;
701 for( ife = 0; ife < total_number_of_gauss_points; ife++ )
702 {
703 // loop on the thickness direction gauss points
704 for( igz = 0; igz < no_gauss_pts_z; igz++ )
705 {
706 zl = -zl;
707 thickness_z = zl * thickness / 2.0;
708
709 aa.set( 0.0, 0.0, 0.0 );
710 bb.set( 0.0, 0.0, 0.0 );
711 cc.set( 0.0, 0.0, 0.0 );
712
713 for( ja = 0; ja < num_nodes; ja++ )
714 {
715 xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
716 xin += thickness_z * normal_at_nodes[ja];
717 aa += dndy1[ife][ja] * xin;
718 bb += dndy2[ife][ja] * xin;
719 thickness_gauss = shape_function[ife][ja] * thickness / 2.0;
720 cc += thickness_gauss * normal_at_nodes[ja];
721 }
722
723 normal_at_point = aa * bb;
724 distrt = cc % normal_at_point;
725 if( distrt < distortion ) distortion = distrt;
726 }
727 }
728
729 // loop through nodal points
730 for( ja = 0; ja < num_nodes; ja++ )
731 {
732 for( igz = 0; igz < no_gauss_pts_z; igz++ )
733 {
734 zl = -zl;
735 thickness_z = zl * thickness / 2.0;
736
737 aa.set( 0.0, 0.0, 0.0 );
738 bb.set( 0.0, 0.0, 0.0 );
739 cc.set( 0.0, 0.0, 0.0 );
740
741 for( jai = 0; jai < num_nodes; jai++ )
742 {
743 xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
744 xin += thickness_z * normal_at_nodes[ja];
745 aa += dndy1_at_node[ja][jai] * xin;
746 bb += dndy2_at_node[ja][jai] * xin;
747 if( jai == ja )
748 thickness_gauss = thickness / 2.0;
749 else
750 thickness_gauss = 0.;
751 cc += thickness_gauss * normal_at_nodes[jai];
752 }
753 }
754
755 normal_at_point = aa * bb;
756 double sign_jacobian = ( tri_normal % normal_at_point ) > 0 ? 1. : -1.;
757 distrt = sign_jacobian * ( cc % normal_at_point );
758
759 if( distrt < distortion ) distortion = distrt;
760 }
761 if( element_area * thickness != 0 )
762 distortion *= 1. / ( element_area * thickness );
763 else
764 distortion *= 1.;
765
766 if( distortion > 0 ) return (double)VERDICT_MIN( distortion, VERDICT_DBL_MAX );
767 return (double)VERDICT_MAX( distortion, -VERDICT_DBL_MAX );
768 }
769
770 /*!
771 tri_quality for calculating multiple tri metrics at once
772
773 using this method is generally faster than using the individual
774 method multiple times.
775
776 */
777 C_FUNC_DEF void v_tri_quality( int num_nodes,
778 double coordinates[][3],
779 unsigned int metrics_request_flag,
780 TriMetricVals* metric_vals )
781 {
782
783 memset( metric_vals, 0, sizeof( TriMetricVals ) );
784
785 // for starts, lets set up some basic and common information
786
787 /* node numbers and side numbers used below
788
789 2
790 ++
791 / \
792 2 / \ 1
793 / \
794 / \
795 0 ---------+ 1
796 0
797 */
798
799 // vectors for each side
800 VerdictVector sides[3];
801 sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
802 coordinates[1][2] - coordinates[0][2] );
803 sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
804 coordinates[2][2] - coordinates[1][2] );
805 sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
806 coordinates[2][2] - coordinates[0][2] );
807 VerdictVector tri_normal = sides[0] * sides[2];
808 // if we have access to normal information, check to see if the
809 // element is inverted. If we don't have the normal information
810 // that we need for this, assume the element is not inverted.
811 // This flag will be used for condition number, jacobian, shape,
812 // and size and shape.
813 bool is_inverted = false;
814 if( compute_normal )
815 {
816 // center of tri
817 double point[3], surf_normal[3];
818 point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
819 point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
820 point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
821 // dot product
822 compute_normal( point, surf_normal );
823 if( ( tri_normal.x() * surf_normal[0] + tri_normal.y() * surf_normal[1] + tri_normal.z() * surf_normal[2] ) <
824 0 )
825 is_inverted = true;
826 }
827
828 // lengths squared of each side
829 double sides_lengths_squared[3];
830 sides_lengths_squared[0] = sides[0].length_squared();
831 sides_lengths_squared[1] = sides[1].length_squared();
832 sides_lengths_squared[2] = sides[2].length_squared();
833
834 // if we are doing angle calcuations
835 if( metrics_request_flag & ( V_TRI_MINIMUM_ANGLE | V_TRI_MAXIMUM_ANGLE ) )
836 {
837 // which is short and long side
838 int short_side = 0, long_side = 0;
839
840 if( sides_lengths_squared[1] < sides_lengths_squared[0] ) short_side = 1;
841 if( sides_lengths_squared[2] < sides_lengths_squared[short_side] ) short_side = 2;
842
843 if( sides_lengths_squared[1] > sides_lengths_squared[0] ) long_side = 1;
844 if( sides_lengths_squared[2] > sides_lengths_squared[long_side] ) long_side = 2;
845
846 // calculate the minimum angle of the tri
847 if( metrics_request_flag & V_TRI_MINIMUM_ANGLE )
848 {
849 if( sides_lengths_squared[0] == 0.0 || sides_lengths_squared[1] == 0.0 || sides_lengths_squared[2] == 0.0 )
850 {
851 metric_vals->minimum_angle = 0.0;
852 }
853 else if( short_side == 0 )
854 metric_vals->minimum_angle = (double)sides[2].interior_angle( sides[1] );
855 else if( short_side == 1 )
856 metric_vals->minimum_angle = (double)sides[0].interior_angle( sides[2] );
857 else
858 metric_vals->minimum_angle = (double)sides[0].interior_angle( -sides[1] );
859 }
860
861 // calculate the maximum angle of the tri
862 if( metrics_request_flag & V_TRI_MAXIMUM_ANGLE )
863 {
864 if( sides_lengths_squared[0] == 0.0 || sides_lengths_squared[1] == 0.0 || sides_lengths_squared[2] == 0.0 )
865 {
866 metric_vals->maximum_angle = 0.0;
867 }
868 else if( long_side == 0 )
869 metric_vals->maximum_angle = (double)sides[2].interior_angle( sides[1] );
870 else if( long_side == 1 )
871 metric_vals->maximum_angle = (double)sides[0].interior_angle( sides[2] );
872 else
873 metric_vals->maximum_angle = (double)sides[0].interior_angle( -sides[1] );
874 }
875 }
876
877 // calculate the area of the tri
878 // the following metrics depend on area
879 if( metrics_request_flag &
880 ( V_TRI_AREA | V_TRI_SCALED_JACOBIAN | V_TRI_SHAPE | V_TRI_RELATIVE_SIZE_SQUARED | V_TRI_SHAPE_AND_SIZE ) )
881 {
882 metric_vals->area = (double)( ( sides[0] * sides[2] ).length() * 0.5 );
883 }
884
885 // calculate the aspect ratio
886 if( metrics_request_flag & V_TRI_ASPECT_FROBENIUS )
887 {
888 // sum the lengths squared
889 double srms = sides_lengths_squared[0] + sides_lengths_squared[1] + sides_lengths_squared[2];
890
891 // calculate once and reuse
892 static const double twoTimesRootOf3 = 2 * sqrt( 3.0 );
893
894 double div = ( twoTimesRootOf3 * ( ( sides[0] * sides[2] ).length() ) );
895
896 if( div == 0.0 )
897 metric_vals->aspect_frobenius = (double)VERDICT_DBL_MAX;
898 else
899 metric_vals->aspect_frobenius = (double)( srms / div );
900 }
901
902 // calculate the scaled jacobian
903 if( metrics_request_flag & V_TRI_SCALED_JACOBIAN )
904 {
905 // calculate once and reuse
906 static const double twoOverRootOf3 = 2 / sqrt( 3.0 );
907 // use the area from above
908
909 double tmp = tri_normal.length() * twoOverRootOf3;
910
911 // now scale it by the lengths of the sides
912 double min_scaled_jac = VERDICT_DBL_MAX;
913 double temp_scaled_jac;
914 for( int i = 0; i < 3; i++ )
915 {
916 if( sides_lengths_squared[i % 3] == 0.0 || sides_lengths_squared[( i + 2 ) % 3] == 0.0 )
917 temp_scaled_jac = 0.0;
918 else
919 temp_scaled_jac =
920 tmp / sqrt( sides_lengths_squared[i % 3] ) / sqrt( sides_lengths_squared[( i + 2 ) % 3] );
921 if( temp_scaled_jac < min_scaled_jac ) min_scaled_jac = temp_scaled_jac;
922 }
923 // multiply by -1 if the normals are in opposite directions
924 if( is_inverted )
925 {
926 min_scaled_jac *= -1;
927 }
928 metric_vals->scaled_jacobian = (double)min_scaled_jac;
929 }
930
931 // calculate the condition number
932 if( metrics_request_flag & V_TRI_CONDITION )
933 {
934 // calculate once and reuse
935 static const double rootOf3 = sqrt( 3.0 );
936 // if it is inverted, the condition number is considered to be infinity.
937 if( is_inverted )
938 {
939 metric_vals->condition = VERDICT_DBL_MAX;
940 }
941 else
942 {
943 double area2x = ( sides[0] * sides[2] ).length();
944 if( area2x == 0.0 )
945 metric_vals->condition = (double)( VERDICT_DBL_MAX );
946 else
947 metric_vals->condition = (double)( ( sides[0] % sides[0] + sides[2] % sides[2] - sides[0] % sides[2] ) /
948 ( area2x * rootOf3 ) );
949 }
950 }
951
952 // calculate the shape
953 if( metrics_request_flag & V_TRI_SHAPE || metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
954 {
955 // if element is inverted, shape is zero. We don't need to
956 // calculate anything.
957 if( is_inverted )
958 {
959 metric_vals->shape = 0.0;
960 }
961 else
962 { // otherwise, we calculate the shape
963 // calculate once and reuse
964 static const double rootOf3 = sqrt( 3.0 );
965 // reuse area from before
966 double area2x = metric_vals->area * 2;
967 // dot products
968 double dots[3] = { sides[0] % sides[0], sides[2] % sides[2], sides[0] % sides[2] };
969
970 // add the dots
971 double sum_dots = dots[0] + dots[1] - dots[2];
972
973 // then the finale
974 if( sum_dots == 0.0 )
975 metric_vals->shape = 0.0;
976 else
977 metric_vals->shape = (double)( rootOf3 * area2x / sum_dots );
978 }
979 }
980
981 // calculate relative size squared
982 if( metrics_request_flag & V_TRI_RELATIVE_SIZE_SQUARED || metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
983 {
984 // get weights
985 double w11, w21, w12, w22;
986 v_tri_get_weight( w11, w21, w12, w22 );
987 // get the determinant
988 double detw = determinant( w11, w21, w12, w22 );
989 // use the area from above and divide with the determinant
990 if( metric_vals->area == 0.0 || detw == 0.0 )
991 metric_vals->relative_size_squared = 0.0;
992 else
993 {
994 double size = metric_vals->area * 2.0 / detw;
995 // square the size
996 size *= size;
997 // value ranges between 0 to 1
998 metric_vals->relative_size_squared = (double)VERDICT_MIN( size, 1.0 / size );
999 }
1000 }
1001
1002 // calculate shape and size
1003 if( metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
1004 {
1005 metric_vals->shape_and_size = metric_vals->relative_size_squared * metric_vals->shape;
1006 }
1007
1008 // calculate distortion
1009 if( metrics_request_flag & V_TRI_DISTORTION ) metric_vals->distortion = v_tri_distortion( num_nodes, coordinates );
1010
1011 // take care of any over-flow problems
1012 if( metric_vals->aspect_frobenius > 0 )
1013 metric_vals->aspect_frobenius = (double)VERDICT_MIN( metric_vals->aspect_frobenius, VERDICT_DBL_MAX );
1014 else
1015 metric_vals->aspect_frobenius = (double)VERDICT_MAX( metric_vals->aspect_frobenius, -VERDICT_DBL_MAX );
1016
1017 if( metric_vals->area > 0 )
1018 metric_vals->area = (double)VERDICT_MIN( metric_vals->area, VERDICT_DBL_MAX );
1019 else
1020 metric_vals->area = (double)VERDICT_MAX( metric_vals->area, -VERDICT_DBL_MAX );
1021
1022 if( metric_vals->minimum_angle > 0 )
1023 metric_vals->minimum_angle = (double)VERDICT_MIN( metric_vals->minimum_angle, VERDICT_DBL_MAX );
1024 else
1025 metric_vals->minimum_angle = (double)VERDICT_MAX( metric_vals->minimum_angle, -VERDICT_DBL_MAX );
1026
1027 if( metric_vals->maximum_angle > 0 )
1028 metric_vals->maximum_angle = (double)VERDICT_MIN( metric_vals->maximum_angle, VERDICT_DBL_MAX );
1029 else
1030 metric_vals->maximum_angle = (double)VERDICT_MAX( metric_vals->maximum_angle, -VERDICT_DBL_MAX );
1031
1032 if( metric_vals->condition > 0 )
1033 metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
1034 else
1035 metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
1036
1037 if( metric_vals->shape > 0 )
1038 metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
1039 else
1040 metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
1041
1042 if( metric_vals->scaled_jacobian > 0 )
1043 metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
1044 else
1045 metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
1046
1047 if( metric_vals->relative_size_squared > 0 )
1048 metric_vals->relative_size_squared = (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
1049 else
1050 metric_vals->relative_size_squared =
1051 (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
1052
1053 if( metric_vals->shape_and_size > 0 )
1054 metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
1055 else
1056 metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
1057
1058 if( metric_vals->distortion > 0 )
1059 metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
1060 else
1061 metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
1062
1063 if( metrics_request_flag & V_TRI_EDGE_RATIO )
1064 {
1065 metric_vals->edge_ratio = v_tri_edge_ratio( 3, coordinates );
1066 }
1067 if( metrics_request_flag & V_TRI_RADIUS_RATIO )
1068 {
1069 metric_vals->radius_ratio = v_tri_radius_ratio( 3, coordinates );
1070 }
1071 if( metrics_request_flag & V_TRI_ASPECT_FROBENIUS ) // there is no V_TRI_ASPECT_RATIO !
1072 {
1073 metric_vals->aspect_ratio = v_tri_radius_ratio( 3, coordinates );
1074 }
1075 }
1.9.1.