 |
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: verdict_defines.hpp,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 * verdict_defines.cpp contains common definitions
17 *
18 * This file is part of VERDICT
19 *
20 */
21
22 #ifndef VERDICT_DEFINES
23 #define VERDICT_DEFINES
24
25 #include <cmath>
26 #include "v_vector.h"
27 #include "VerdictVector.hpp"
28
29 enum VerdictBoolean
30 {
31 VERDICT_FALSE = 0,
32 VERDICT_TRUE = 1
33 };
34
35 #define VERDICT_MIN( a, b ) ( ( a ) < ( b ) ? ( a ) : ( b ) )
36 #define VERDICT_MAX( a, b ) ( ( a ) > ( b ) ? ( a ) : ( b ) )
37
38 inline double determinant( double a, double b, double c, double d )
39 {
40 return ( ( a ) * ( d ) - ( b ) * ( c ) );
41 }
42
43 inline double determinant( const VerdictVector& v1, const VerdictVector& v2, const VerdictVector& v3 )
44 {
45 return v1 % ( v2 * v3 );
46 }
47
48 #define jacobian_matrix( a, b, c, d, e, f, g ) \
49 double jac_mat_tmp; \
50 jac_mat_tmp = sqrt( a ); \
51 if( jac_mat_tmp == 0 ) \
52 { \
53 ( d ) = 0; \
54 ( e ) = 0; \
55 ( f ) = 0; \
56 ( g ) = 0; \
57 } \
58 else \
59 { \
60 ( d ) = jac_mat_tmp; \
61 ( e ) = 0; \
62 ( f ) = ( b ) / jac_mat_tmp; \
63 ( g ) = ( c ) / jac_mat_tmp; \
64 }
65
66 // this assumes that detmw != 0
67 #define form_t( m11, m21, m12, m22, mw11, mw21, mw12, mw22, detmw, xm11, xm21, xm12, xm22 ) \
68 xm11 = ( ( m11 ) * ( mw22 ) - ( m12 ) * ( mw21 ) ) / ( detmw ); \
69 ( xm21 ) = ( ( m21 ) * ( mw22 ) - ( m22 ) * ( mw21 ) ) / ( detmw ); \
70 ( xm12 ) = ( ( m12 ) * ( mw11 ) - ( m11 ) * ( mw12 ) ) / ( detmw ); \
71 ( xm22 ) = ( ( m22 ) * ( mw11 ) - ( m21 ) * ( mw12 ) ) / ( detmw );
72
73 extern double verdictSqrt2;
74
75 inline double normalize_jacobian( double jacobi,
76 VerdictVector& v1,
77 VerdictVector& v2,
78 VerdictVector& v3,
79 int tet_flag = 0 )
80 {
81 double return_value = 0.0;
82
83 if( jacobi != 0.0 )
84 {
85
86 double l1, l2, l3, length_product;
87 // Note: there may be numerical problems if one is a lot shorter
88 // than the others this way. But scaling each vector before the
89 // triple product would involve 3 square roots instead of just
90 // one.
91 l1 = v1.length_squared();
92 l2 = v2.length_squared();
93 l3 = v3.length_squared();
94 length_product = sqrt( l1 * l2 * l3 );
95
96 // if some numerical scaling problem, or just plain roundoff,
97 // then push back into range [-1,1].
98 if( length_product < fabs( jacobi ) )
99 {
100 length_product = fabs( jacobi );
101 }
102
103 if( tet_flag == 1 )
104 return_value = verdictSqrt2 * jacobi / length_product;
105 else
106 return_value = jacobi / length_product;
107 }
108 return return_value;
109 }
110
111 inline double norm_squared( double m11, double m21, double m12, double m22 )
112 {
113 return m11 * m11 + m21 * m21 + m12 * m12 + m22 * m22;
114 }
115
116 #define metric_matrix( m11, m21, m12, m22, gm11, gm12, gm22 ) \
117 gm11 = ( m11 ) * ( m11 ) + ( m21 ) * ( m21 ); \
118 ( gm12 ) = ( m11 ) * ( m12 ) + ( m21 ) * ( m22 ); \
119 ( gm22 ) = ( m12 ) * ( m12 ) + ( m22 ) * ( m22 );
120
121 inline int skew_matrix( double gm11,
122 double gm12,
123 double gm22,
124 double det,
125 double& qm11,
126 double& qm21,
127 double& qm12,
128 double& qm22 )
129 {
130 double tmp = sqrt( gm11 * gm22 );
131 if( tmp == 0 )
132 {
133 return false;
134 }
135
136 qm11 = 1;
137 qm21 = 0;
138 qm12 = gm12 / tmp;
139 qm22 = det / tmp;
140 return true;
141 }
142
143 inline void inverse( const VerdictVector& x1,
144 const VerdictVector& x2,
145 const VerdictVector& x3,
146 VerdictVector& u1,
147 VerdictVector& u2,
148 VerdictVector& u3 )
149 {
150 double detx = determinant( x1, x2, x3 );
151 VerdictVector rx1, rx2, rx3;
152
153 rx1.set( x1.x(), x2.x(), x3.x() );
154 rx2.set( x1.y(), x2.y(), x3.y() );
155 rx3.set( x1.z(), x2.z(), x3.z() );
156
157 u1 = rx2 * rx3;
158 u2 = rx3 * rx1;
159 u3 = rx1 * rx2;
160
161 u1 /= detx;
162 u2 /= detx;
163 u3 /= detx;
164 }
165
166 /*
167 inline void form_T(double a1[3],
168 double a2[3],
169 double a3[3],
170 double w1[3],
171 double w2[3],
172 double w3[3],
173 double t1[3],
174 double t2[3],
175 double t3[3] )
176 {
177
178 double x1[3], x2[3], x3[3];
179 double ra1[3], ra2[3], ra3[3];
180
181 x1[0] = a1[0];
182 x1[1] = a2[0];
183 x1[2] = a3[0];
184
185 x2[0] = a1[1];
186 x2[1] = a2[1];
187 x2[2] = a3[1];
188
189 x3[0] = a1[2];
190 x3[1] = a2[2];
191 x3[2] = a3[2];
192
193 inverse(w1,w2,w3,x1,x2,x3);
194
195 t1[0] = dot_product(ra1, x1);
196 t1[1] = dot_product(ra1, x2);
197 t1[2] = dot_product(ra1, x3);
198
199 t2[0] = dot_product(ra2, x1);
200 t2[0] = dot_product(ra2, x2);
201 t2[0] = dot_product(ra2, x3);
202
203 t3[0] = dot_product(ra3, x1);
204 t3[0] = dot_product(ra3, x2);
205 t3[0] = dot_product(ra3, x3);
206
207 }
208 */
209
210 inline void form_Q( const VerdictVector& v1,
211 const VerdictVector& v2,
212 const VerdictVector& v3,
213 VerdictVector& q1,
214 VerdictVector& q2,
215 VerdictVector& q3 )
216 {
217
218 double g11, g12, g13, g22, g23, g33;
219
220 g11 = v1 % v1;
221 g12 = v1 % v2;
222 g13 = v1 % v3;
223 g22 = v2 % v2;
224 g23 = v2 % v3;
225 g33 = v3 % v3;
226
227 double rtg11 = sqrt( g11 );
228 double rtg22 = sqrt( g22 );
229 double rtg33 = sqrt( g33 );
230 VerdictVector temp1;
231
232 temp1 = v1 * v2;
233
234 double cross = sqrt( temp1 % temp1 );
235
236 double q11, q21, q31;
237 double q12, q22, q32;
238 double q13, q23, q33;
239
240 q11 = 1;
241 q21 = 0;
242 q31 = 0;
243
244 q12 = g12 / rtg11 / rtg22;
245 q22 = cross / rtg11 / rtg22;
246 q32 = 0;
247
248 q13 = g13 / rtg11 / rtg33;
249 q23 = ( g11 * g23 - g12 * g13 ) / rtg11 / rtg33 / cross;
250 temp1 = v2 * v3;
251 q33 = ( v1 % temp1 ) / rtg33 / cross;
252
253 q1.set( q11, q21, q31 );
254 q2.set( q12, q22, q32 );
255 q3.set( q13, q23, q33 );
256 }
257
258 inline void product( VerdictVector& a1,
259 VerdictVector& a2,
260 VerdictVector& a3,
261 VerdictVector& b1,
262 VerdictVector& b2,
263 VerdictVector& b3,
264 VerdictVector& c1,
265 VerdictVector& c2,
266 VerdictVector& c3 )
267 {
268
269 VerdictVector x1, x2, x3;
270
271 x1.set( a1.x(), a2.x(), a3.x() );
272 x2.set( a1.y(), a2.y(), a3.y() );
273 x3.set( a1.z(), a2.z(), a3.z() );
274
275 c1.set( x1 % b1, x2 % b1, x3 % b1 );
276 c2.set( x1 % b2, x2 % b2, x3 % b2 );
277 c3.set( x1 % b3, x2 % b3, x3 % b3 );
278 }
279
280 inline double norm_squared( VerdictVector& x1, VerdictVector& x2, VerdictVector& x3 )
281
282 {
283 return ( x1 % x1 ) + ( x2 % x2 ) + ( x3 % x3 );
284 }
285
286 inline double skew_x( VerdictVector& q1,
287 VerdictVector& q2,
288 VerdictVector& q3,
289 VerdictVector& qw1,
290 VerdictVector& qw2,
291 VerdictVector& qw3 )
292 {
293 double normsq1, normsq2, kappa;
294 VerdictVector u1, u2, u3;
295 VerdictVector x1, x2, x3;
296
297 inverse( qw1, qw2, qw3, u1, u2, u3 );
298 product( q1, q2, q3, u1, u2, u3, x1, x2, x3 );
299 inverse( x1, x2, x3, u1, u2, u3 );
300 normsq1 = norm_squared( x1, x2, x3 );
301 normsq2 = norm_squared( u1, u2, u3 );
302 kappa = sqrt( normsq1 * normsq2 );
303
304 double skew = 0;
305 if( kappa > VERDICT_DBL_MIN ) skew = 3 / kappa;
306
307 return skew;
308 }
309
310 #endif
1.9.1.