MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /*========================================================================= 00002 00003 Module: $RCSfile: verdict_defines.hpp,v $ 00004 00005 Copyright (c) 2006 Sandia Corporation. 00006 All rights reserved. 00007 See Copyright.txt or http://www.kitware.com/Copyright.htm for details. 00008 00009 This software is distributed WITHOUT ANY WARRANTY; without even 00010 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 00011 PURPOSE. See the above copyright notice for more information. 00012 00013 =========================================================================*/ 00014 00015 /* 00016 * verdict_defines.cpp contains common definitions 00017 * 00018 * This file is part of VERDICT 00019 * 00020 */ 00021 00022 #ifndef VERDICT_DEFINES 00023 #define VERDICT_DEFINES 00024 00025 #include <math.h> 00026 #include "v_vector.h" 00027 #include "VerdictVector.hpp" 00028 00029 enum VerdictBoolean 00030 { 00031 VERDICT_FALSE = 0, 00032 VERDICT_TRUE = 1 00033 }; 00034 00035 #define VERDICT_MIN( a, b ) ( ( a ) < ( b ) ? ( a ) : ( b ) ) 00036 #define VERDICT_MAX( a, b ) ( ( a ) > ( b ) ? ( a ) : ( b ) ) 00037 00038 inline double determinant( double a, double b, double c, double d ) 00039 { 00040 return ( ( a ) * ( d ) - ( b ) * ( c ) ); 00041 } 00042 00043 inline double determinant( VerdictVector v1, VerdictVector v2, VerdictVector v3 ) 00044 { 00045 return v1 % ( v2 * v3 ); 00046 } 00047 00048 #define jacobian_matrix( a, b, c, d, e, f, g ) \ 00049 double jac_mat_tmp; \ 00050 jac_mat_tmp = sqrt( a ); \ 00051 if( jac_mat_tmp == 0 ) \ 00052 { \ 00053 d = 0; \ 00054 e = 0; \ 00055 f = 0; \ 00056 g = 0; \ 00057 } \ 00058 else \ 00059 { \ 00060 d = jac_mat_tmp; \ 00061 e = 0; \ 00062 f = b / jac_mat_tmp; \ 00063 g = c / jac_mat_tmp; \ 00064 } 00065 00066 // this assumes that detmw != 0 00067 #define form_t( m11, m21, m12, m22, mw11, mw21, mw12, mw22, detmw, xm11, xm21, xm12, xm22 ) \ 00068 xm11 = ( m11 * mw22 - m12 * mw21 ) / detmw; \ 00069 xm21 = ( m21 * mw22 - m22 * mw21 ) / detmw; \ 00070 xm12 = ( m12 * mw11 - m11 * mw12 ) / detmw; \ 00071 xm22 = ( m22 * mw11 - m21 * mw12 ) / detmw; 00072 00073 extern double verdictSqrt2; 00074 00075 inline double normalize_jacobian( double jacobi, VerdictVector& v1, VerdictVector& v2, VerdictVector& v3, 00076 int tet_flag = 0 ) 00077 { 00078 double return_value = 0.0; 00079 00080 if( jacobi != 0.0 ) 00081 { 00082 00083 double l1, l2, l3, length_product; 00084 // Note: there may be numerical problems if one is a lot shorter 00085 // than the others this way. But scaling each vector before the 00086 // triple product would involve 3 square roots instead of just 00087 // one. 00088 l1 = v1.length_squared(); 00089 l2 = v2.length_squared(); 00090 l3 = v3.length_squared(); 00091 length_product = sqrt( l1 * l2 * l3 ); 00092 00093 // if some numerical scaling problem, or just plain roundoff, 00094 // then push back into range [-1,1]. 00095 if( length_product < fabs( jacobi ) ) { length_product = fabs( jacobi ); } 00096 00097 if( tet_flag == 1 ) 00098 return_value = verdictSqrt2 * jacobi / length_product; 00099 else 00100 return_value = jacobi / length_product; 00101 } 00102 return return_value; 00103 } 00104 00105 inline double norm_squared( double m11, double m21, double m12, double m22 ) 00106 { 00107 return m11 * m11 + m21 * m21 + m12 * m12 + m22 * m22; 00108 } 00109 00110 #define metric_matrix( m11, m21, m12, m22, gm11, gm12, gm22 ) \ 00111 gm11 = m11 * m11 + m21 * m21; \ 00112 gm12 = m11 * m12 + m21 * m22; \ 00113 gm22 = m12 * m12 + m22 * m22; 00114 00115 inline int skew_matrix( double gm11, double gm12, double gm22, double det, double& qm11, double& qm21, double& qm12, 00116 double& qm22 ) 00117 { 00118 double tmp = sqrt( gm11 * gm22 ); 00119 if( tmp == 0 ) { return false; } 00120 00121 qm11 = 1; 00122 qm21 = 0; 00123 qm12 = gm12 / tmp; 00124 qm22 = det / tmp; 00125 return true; 00126 } 00127 00128 inline void inverse( VerdictVector x1, VerdictVector x2, VerdictVector x3, VerdictVector& u1, VerdictVector& u2, 00129 VerdictVector& u3 ) 00130 { 00131 double detx = determinant( x1, x2, x3 ); 00132 VerdictVector rx1, rx2, rx3; 00133 00134 rx1.set( x1.x(), x2.x(), x3.x() ); 00135 rx2.set( x1.y(), x2.y(), x3.y() ); 00136 rx3.set( x1.z(), x2.z(), x3.z() ); 00137 00138 u1 = rx2 * rx3; 00139 u2 = rx3 * rx1; 00140 u3 = rx1 * rx2; 00141 00142 u1 /= detx; 00143 u2 /= detx; 00144 u3 /= detx; 00145 } 00146 00147 /* 00148 inline void form_T(double a1[3], 00149 double a2[3], 00150 double a3[3], 00151 double w1[3], 00152 double w2[3], 00153 double w3[3], 00154 double t1[3], 00155 double t2[3], 00156 double t3[3] ) 00157 { 00158 00159 double x1[3], x2[3], x3[3]; 00160 double ra1[3], ra2[3], ra3[3]; 00161 00162 x1[0] = a1[0]; 00163 x1[1] = a2[0]; 00164 x1[2] = a3[0]; 00165 00166 x2[0] = a1[1]; 00167 x2[1] = a2[1]; 00168 x2[2] = a3[1]; 00169 00170 x3[0] = a1[2]; 00171 x3[1] = a2[2]; 00172 x3[2] = a3[2]; 00173 00174 inverse(w1,w2,w3,x1,x2,x3); 00175 00176 t1[0] = dot_product(ra1, x1); 00177 t1[1] = dot_product(ra1, x2); 00178 t1[2] = dot_product(ra1, x3); 00179 00180 t2[0] = dot_product(ra2, x1); 00181 t2[0] = dot_product(ra2, x2); 00182 t2[0] = dot_product(ra2, x3); 00183 00184 t3[0] = dot_product(ra3, x1); 00185 t3[0] = dot_product(ra3, x2); 00186 t3[0] = dot_product(ra3, x3); 00187 00188 } 00189 */ 00190 00191 inline void form_Q( const VerdictVector& v1, const VerdictVector& v2, const VerdictVector& v3, VerdictVector& q1, 00192 VerdictVector& q2, VerdictVector& q3 ) 00193 { 00194 00195 double g11, g12, g13, g22, g23, g33; 00196 00197 g11 = v1 % v1; 00198 g12 = v1 % v2; 00199 g13 = v1 % v3; 00200 g22 = v2 % v2; 00201 g23 = v2 % v3; 00202 g33 = v3 % v3; 00203 00204 double rtg11 = sqrt( g11 ); 00205 double rtg22 = sqrt( g22 ); 00206 double rtg33 = sqrt( g33 ); 00207 VerdictVector temp1; 00208 00209 temp1 = v1 * v2; 00210 00211 double cross = sqrt( temp1 % temp1 ); 00212 00213 double q11, q21, q31; 00214 double q12, q22, q32; 00215 double q13, q23, q33; 00216 00217 q11 = 1; 00218 q21 = 0; 00219 q31 = 0; 00220 00221 q12 = g12 / rtg11 / rtg22; 00222 q22 = cross / rtg11 / rtg22; 00223 q32 = 0; 00224 00225 q13 = g13 / rtg11 / rtg33; 00226 q23 = ( g11 * g23 - g12 * g13 ) / rtg11 / rtg33 / cross; 00227 temp1 = v2 * v3; 00228 q33 = ( v1 % temp1 ) / rtg33 / cross; 00229 00230 q1.set( q11, q21, q31 ); 00231 q2.set( q12, q22, q32 ); 00232 q3.set( q13, q23, q33 ); 00233 } 00234 00235 inline void product( VerdictVector& a1, VerdictVector& a2, VerdictVector& a3, VerdictVector& b1, VerdictVector& b2, 00236 VerdictVector& b3, VerdictVector& c1, VerdictVector& c2, VerdictVector& c3 ) 00237 { 00238 00239 VerdictVector x1, x2, x3; 00240 00241 x1.set( a1.x(), a2.x(), a3.x() ); 00242 x2.set( a1.y(), a2.y(), a3.y() ); 00243 x3.set( a1.z(), a2.z(), a3.z() ); 00244 00245 c1.set( x1 % b1, x2 % b1, x3 % b1 ); 00246 c2.set( x1 % b2, x2 % b2, x3 % b2 ); 00247 c3.set( x1 % b3, x2 % b3, x3 % b3 ); 00248 } 00249 00250 inline double norm_squared( VerdictVector& x1, VerdictVector& x2, VerdictVector& x3 ) 00251 00252 { 00253 return ( x1 % x1 ) + ( x2 % x2 ) + ( x3 % x3 ); 00254 } 00255 00256 inline double skew_x( VerdictVector& q1, VerdictVector& q2, VerdictVector& q3, VerdictVector& qw1, VerdictVector& qw2, 00257 VerdictVector& qw3 ) 00258 { 00259 double normsq1, normsq2, kappa; 00260 VerdictVector u1, u2, u3; 00261 VerdictVector x1, x2, x3; 00262 00263 inverse( qw1, qw2, qw3, u1, u2, u3 ); 00264 product( q1, q2, q3, u1, u2, u3, x1, x2, x3 ); 00265 inverse( x1, x2, x3, u1, u2, u3 ); 00266 normsq1 = norm_squared( x1, x2, x3 ); 00267 normsq2 = norm_squared( u1, u2, u3 ); 00268 kappa = sqrt( normsq1 * normsq2 ); 00269 00270 double skew = 0; 00271 if( kappa > VERDICT_DBL_MIN ) skew = 3 / kappa; 00272 00273 return skew; 00274 } 00275 00276 #endif