MOAB: Mesh Oriented datABase  (version 5.2.1)
verdict_defines.hpp
Go to the documentation of this file.
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 <cmath>
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( const VerdictVector& v1, const VerdictVector& v2, const 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( const VerdictVector& x1, const VerdictVector& x2, const 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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines