Branch data Line data Source code
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 <math.h>
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 : 91 : inline double determinant( double a, double b, double c, double d )
39 : : {
40 : 91 : return ( ( a ) * ( d ) - ( b ) * ( c ) );
41 : : }
42 : :
43 : 46 : inline double determinant( VerdictVector v1, VerdictVector v2, VerdictVector v3 )
44 : : {
45 [ + - ]: 46 : 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, VerdictVector& v1, VerdictVector& v2, VerdictVector& v3,
76 : : int tet_flag = 0 )
77 : : {
78 : : double return_value = 0.0;
79 : :
80 : : if( jacobi != 0.0 )
81 : : {
82 : :
83 : : double l1, l2, l3, length_product;
84 : : // Note: there may be numerical problems if one is a lot shorter
85 : : // than the others this way. But scaling each vector before the
86 : : // triple product would involve 3 square roots instead of just
87 : : // one.
88 : : l1 = v1.length_squared();
89 : : l2 = v2.length_squared();
90 : : l3 = v3.length_squared();
91 : : length_product = sqrt( l1 * l2 * l3 );
92 : :
93 : : // if some numerical scaling problem, or just plain roundoff,
94 : : // then push back into range [-1,1].
95 : : if( length_product < fabs( jacobi ) ) { length_product = fabs( jacobi ); }
96 : :
97 : : if( tet_flag == 1 )
98 : : return_value = verdictSqrt2 * jacobi / length_product;
99 : : else
100 : : return_value = jacobi / length_product;
101 : : }
102 : : return return_value;
103 : : }
104 : :
105 : : inline double norm_squared( double m11, double m21, double m12, double m22 )
106 : : {
107 : : return m11 * m11 + m21 * m21 + m12 * m12 + m22 * m22;
108 : : }
109 : :
110 : : #define metric_matrix( m11, m21, m12, m22, gm11, gm12, gm22 ) \
111 : : gm11 = m11 * m11 + m21 * m21; \
112 : : gm12 = m11 * m12 + m21 * m22; \
113 : : gm22 = m12 * m12 + m22 * m22;
114 : :
115 : : inline int skew_matrix( double gm11, double gm12, double gm22, double det, double& qm11, double& qm21, double& qm12,
116 : : double& qm22 )
117 : : {
118 : : double tmp = sqrt( gm11 * gm22 );
119 : : if( tmp == 0 ) { return false; }
120 : :
121 : : qm11 = 1;
122 : : qm21 = 0;
123 : : qm12 = gm12 / tmp;
124 : : qm22 = det / tmp;
125 : : return true;
126 : : }
127 : :
128 : : inline void inverse( VerdictVector x1, VerdictVector x2, VerdictVector x3, VerdictVector& u1, VerdictVector& u2,
129 : : VerdictVector& u3 )
130 : : {
131 : : double detx = determinant( x1, x2, x3 );
132 : : VerdictVector rx1, rx2, rx3;
133 : :
134 : : rx1.set( x1.x(), x2.x(), x3.x() );
135 : : rx2.set( x1.y(), x2.y(), x3.y() );
136 : : rx3.set( x1.z(), x2.z(), x3.z() );
137 : :
138 : : u1 = rx2 * rx3;
139 : : u2 = rx3 * rx1;
140 : : u3 = rx1 * rx2;
141 : :
142 : : u1 /= detx;
143 : : u2 /= detx;
144 : : u3 /= detx;
145 : : }
146 : :
147 : : /*
148 : : inline void form_T(double a1[3],
149 : : double a2[3],
150 : : double a3[3],
151 : : double w1[3],
152 : : double w2[3],
153 : : double w3[3],
154 : : double t1[3],
155 : : double t2[3],
156 : : double t3[3] )
157 : : {
158 : :
159 : : double x1[3], x2[3], x3[3];
160 : : double ra1[3], ra2[3], ra3[3];
161 : :
162 : : x1[0] = a1[0];
163 : : x1[1] = a2[0];
164 : : x1[2] = a3[0];
165 : :
166 : : x2[0] = a1[1];
167 : : x2[1] = a2[1];
168 : : x2[2] = a3[1];
169 : :
170 : : x3[0] = a1[2];
171 : : x3[1] = a2[2];
172 : : x3[2] = a3[2];
173 : :
174 : : inverse(w1,w2,w3,x1,x2,x3);
175 : :
176 : : t1[0] = dot_product(ra1, x1);
177 : : t1[1] = dot_product(ra1, x2);
178 : : t1[2] = dot_product(ra1, x3);
179 : :
180 : : t2[0] = dot_product(ra2, x1);
181 : : t2[0] = dot_product(ra2, x2);
182 : : t2[0] = dot_product(ra2, x3);
183 : :
184 : : t3[0] = dot_product(ra3, x1);
185 : : t3[0] = dot_product(ra3, x2);
186 : : t3[0] = dot_product(ra3, x3);
187 : :
188 : : }
189 : : */
190 : :
191 : : inline void form_Q( const VerdictVector& v1, const VerdictVector& v2, const VerdictVector& v3, VerdictVector& q1,
192 : : VerdictVector& q2, VerdictVector& q3 )
193 : : {
194 : :
195 : : double g11, g12, g13, g22, g23, g33;
196 : :
197 : : g11 = v1 % v1;
198 : : g12 = v1 % v2;
199 : : g13 = v1 % v3;
200 : : g22 = v2 % v2;
201 : : g23 = v2 % v3;
202 : : g33 = v3 % v3;
203 : :
204 : : double rtg11 = sqrt( g11 );
205 : : double rtg22 = sqrt( g22 );
206 : : double rtg33 = sqrt( g33 );
207 : : VerdictVector temp1;
208 : :
209 : : temp1 = v1 * v2;
210 : :
211 : : double cross = sqrt( temp1 % temp1 );
212 : :
213 : : double q11, q21, q31;
214 : : double q12, q22, q32;
215 : : double q13, q23, q33;
216 : :
217 : : q11 = 1;
218 : : q21 = 0;
219 : : q31 = 0;
220 : :
221 : : q12 = g12 / rtg11 / rtg22;
222 : : q22 = cross / rtg11 / rtg22;
223 : : q32 = 0;
224 : :
225 : : q13 = g13 / rtg11 / rtg33;
226 : : q23 = ( g11 * g23 - g12 * g13 ) / rtg11 / rtg33 / cross;
227 : : temp1 = v2 * v3;
228 : : q33 = ( v1 % temp1 ) / rtg33 / cross;
229 : :
230 : : q1.set( q11, q21, q31 );
231 : : q2.set( q12, q22, q32 );
232 : : q3.set( q13, q23, q33 );
233 : : }
234 : :
235 : : inline void product( VerdictVector& a1, VerdictVector& a2, VerdictVector& a3, VerdictVector& b1, VerdictVector& b2,
236 : : VerdictVector& b3, VerdictVector& c1, VerdictVector& c2, VerdictVector& c3 )
237 : : {
238 : :
239 : : VerdictVector x1, x2, x3;
240 : :
241 : : x1.set( a1.x(), a2.x(), a3.x() );
242 : : x2.set( a1.y(), a2.y(), a3.y() );
243 : : x3.set( a1.z(), a2.z(), a3.z() );
244 : :
245 : : c1.set( x1 % b1, x2 % b1, x3 % b1 );
246 : : c2.set( x1 % b2, x2 % b2, x3 % b2 );
247 : : c3.set( x1 % b3, x2 % b3, x3 % b3 );
248 : : }
249 : :
250 : : inline double norm_squared( VerdictVector& x1, VerdictVector& x2, VerdictVector& x3 )
251 : :
252 : : {
253 : : return ( x1 % x1 ) + ( x2 % x2 ) + ( x3 % x3 );
254 : : }
255 : :
256 : : inline double skew_x( VerdictVector& q1, VerdictVector& q2, VerdictVector& q3, VerdictVector& qw1, VerdictVector& qw2,
257 : : VerdictVector& qw3 )
258 : : {
259 : : double normsq1, normsq2, kappa;
260 : : VerdictVector u1, u2, u3;
261 : : VerdictVector x1, x2, x3;
262 : :
263 : : inverse( qw1, qw2, qw3, u1, u2, u3 );
264 : : product( q1, q2, q3, u1, u2, u3, x1, x2, x3 );
265 : : inverse( x1, x2, x3, u1, u2, u3 );
266 : : normsq1 = norm_squared( x1, x2, x3 );
267 : : normsq2 = norm_squared( u1, u2, u3 );
268 : : kappa = sqrt( normsq1 * normsq2 );
269 : :
270 : : double skew = 0;
271 : : if( kappa > VERDICT_DBL_MIN ) skew = 3 / kappa;
272 : :
273 : : return skew;
274 : : }
275 : :
276 : : #endif
|
