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