Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
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,
00076                                   VerdictVector& v1,
00077                                   VerdictVector& v2,
00078                                   VerdictVector& v3,
00079                                   int tet_flag = 0 )
00080 {
00081     double return_value = 0.0;
00082 
00083     if( jacobi != 0.0 )
00084     {
00085 
00086         double l1, l2, l3, length_product;
00087         // Note: there may be numerical problems if one is a lot shorter
00088         // than the others this way. But scaling each vector before the
00089         // triple product would involve 3 square roots instead of just
00090         // one.
00091         l1             = v1.length_squared();
00092         l2             = v2.length_squared();
00093         l3             = v3.length_squared();
00094         length_product = sqrt( l1 * l2 * l3 );
00095 
00096         // if some numerical scaling problem, or just plain roundoff,
00097         // then push back into range [-1,1].
00098         if( length_product < fabs( jacobi ) )
00099         {
00100             length_product = fabs( jacobi );
00101         }
00102 
00103         if( tet_flag == 1 )
00104             return_value = verdictSqrt2 * jacobi / length_product;
00105         else
00106             return_value = jacobi / length_product;
00107     }
00108     return return_value;
00109 }
00110 
00111 inline double norm_squared( double m11, double m21, double m12, double m22 )
00112 {
00113     return m11 * m11 + m21 * m21 + m12 * m12 + m22 * m22;
00114 }
00115 
00116 #define metric_matrix( m11, m21, m12, m22, gm11, gm12, gm22 ) \
00117     gm11     = ( m11 ) * ( m11 ) + ( m21 ) * ( m21 );         \
00118     ( gm12 ) = ( m11 ) * ( m12 ) + ( m21 ) * ( m22 );         \
00119     ( gm22 ) = ( m12 ) * ( m12 ) + ( m22 ) * ( m22 );
00120 
00121 inline int skew_matrix( double gm11,
00122                         double gm12,
00123                         double gm22,
00124                         double det,
00125                         double& qm11,
00126                         double& qm21,
00127                         double& qm12,
00128                         double& qm22 )
00129 {
00130     double tmp = sqrt( gm11 * gm22 );
00131     if( tmp == 0 )
00132     {
00133         return false;
00134     }
00135 
00136     qm11 = 1;
00137     qm21 = 0;
00138     qm12 = gm12 / tmp;
00139     qm22 = det / tmp;
00140     return true;
00141 }
00142 
00143 inline void inverse( const VerdictVector& x1,
00144                      const VerdictVector& x2,
00145                      const VerdictVector& x3,
00146                      VerdictVector& u1,
00147                      VerdictVector& u2,
00148                      VerdictVector& u3 )
00149 {
00150     double detx = determinant( x1, x2, x3 );
00151     VerdictVector rx1, rx2, rx3;
00152 
00153     rx1.set( x1.x(), x2.x(), x3.x() );
00154     rx2.set( x1.y(), x2.y(), x3.y() );
00155     rx3.set( x1.z(), x2.z(), x3.z() );
00156 
00157     u1 = rx2 * rx3;
00158     u2 = rx3 * rx1;
00159     u3 = rx1 * rx2;
00160 
00161     u1 /= detx;
00162     u2 /= detx;
00163     u3 /= detx;
00164 }
00165 
00166 /*
00167 inline void form_T(double a1[3],
00168   double a2[3],
00169   double a3[3],
00170   double w1[3],
00171   double w2[3],
00172   double w3[3],
00173   double t1[3],
00174   double t2[3],
00175   double t3[3] )
00176 {
00177 
00178   double x1[3], x2[3], x3[3];
00179   double ra1[3], ra2[3], ra3[3];
00180 
00181   x1[0] = a1[0];
00182   x1[1] = a2[0];
00183   x1[2] = a3[0];
00184 
00185   x2[0] = a1[1];
00186   x2[1] = a2[1];
00187   x2[2] = a3[1];
00188 
00189   x3[0] = a1[2];
00190   x3[1] = a2[2];
00191   x3[2] = a3[2];
00192 
00193   inverse(w1,w2,w3,x1,x2,x3);
00194 
00195   t1[0] = dot_product(ra1, x1);
00196   t1[1] = dot_product(ra1, x2);
00197   t1[2] = dot_product(ra1, x3);
00198 
00199   t2[0] = dot_product(ra2, x1);
00200   t2[0] = dot_product(ra2, x2);
00201   t2[0] = dot_product(ra2, x3);
00202 
00203   t3[0] = dot_product(ra3, x1);
00204   t3[0] = dot_product(ra3, x2);
00205   t3[0] = dot_product(ra3, x3);
00206 
00207 }
00208 */
00209 
00210 inline void form_Q( const VerdictVector& v1,
00211                     const VerdictVector& v2,
00212                     const VerdictVector& v3,
00213                     VerdictVector& q1,
00214                     VerdictVector& q2,
00215                     VerdictVector& q3 )
00216 {
00217 
00218     double g11, g12, g13, g22, g23, g33;
00219 
00220     g11 = v1 % v1;
00221     g12 = v1 % v2;
00222     g13 = v1 % v3;
00223     g22 = v2 % v2;
00224     g23 = v2 % v3;
00225     g33 = v3 % v3;
00226 
00227     double rtg11 = sqrt( g11 );
00228     double rtg22 = sqrt( g22 );
00229     double rtg33 = sqrt( g33 );
00230     VerdictVector temp1;
00231 
00232     temp1 = v1 * v2;
00233 
00234     double cross = sqrt( temp1 % temp1 );
00235 
00236     double q11, q21, q31;
00237     double q12, q22, q32;
00238     double q13, q23, q33;
00239 
00240     q11 = 1;
00241     q21 = 0;
00242     q31 = 0;
00243 
00244     q12 = g12 / rtg11 / rtg22;
00245     q22 = cross / rtg11 / rtg22;
00246     q32 = 0;
00247 
00248     q13   = g13 / rtg11 / rtg33;
00249     q23   = ( g11 * g23 - g12 * g13 ) / rtg11 / rtg33 / cross;
00250     temp1 = v2 * v3;
00251     q33   = ( v1 % temp1 ) / rtg33 / cross;
00252 
00253     q1.set( q11, q21, q31 );
00254     q2.set( q12, q22, q32 );
00255     q3.set( q13, q23, q33 );
00256 }
00257 
00258 inline void product( VerdictVector& a1,
00259                      VerdictVector& a2,
00260                      VerdictVector& a3,
00261                      VerdictVector& b1,
00262                      VerdictVector& b2,
00263                      VerdictVector& b3,
00264                      VerdictVector& c1,
00265                      VerdictVector& c2,
00266                      VerdictVector& c3 )
00267 {
00268 
00269     VerdictVector x1, x2, x3;
00270 
00271     x1.set( a1.x(), a2.x(), a3.x() );
00272     x2.set( a1.y(), a2.y(), a3.y() );
00273     x3.set( a1.z(), a2.z(), a3.z() );
00274 
00275     c1.set( x1 % b1, x2 % b1, x3 % b1 );
00276     c2.set( x1 % b2, x2 % b2, x3 % b2 );
00277     c3.set( x1 % b3, x2 % b3, x3 % b3 );
00278 }
00279 
00280 inline double norm_squared( VerdictVector& x1, VerdictVector& x2, VerdictVector& x3 )
00281 
00282 {
00283     return ( x1 % x1 ) + ( x2 % x2 ) + ( x3 % x3 );
00284 }
00285 
00286 inline double skew_x( VerdictVector& q1,
00287                       VerdictVector& q2,
00288                       VerdictVector& q3,
00289                       VerdictVector& qw1,
00290                       VerdictVector& qw2,
00291                       VerdictVector& qw3 )
00292 {
00293     double normsq1, normsq2, kappa;
00294     VerdictVector u1, u2, u3;
00295     VerdictVector x1, x2, x3;
00296 
00297     inverse( qw1, qw2, qw3, u1, u2, u3 );
00298     product( q1, q2, q3, u1, u2, u3, x1, x2, x3 );
00299     inverse( x1, x2, x3, u1, u2, u3 );
00300     normsq1 = norm_squared( x1, x2, x3 );
00301     normsq2 = norm_squared( u1, u2, u3 );
00302     kappa   = sqrt( normsq1 * normsq2 );
00303 
00304     double skew = 0;
00305     if( kappa > VERDICT_DBL_MIN ) skew = 3 / kappa;
00306 
00307     return skew;
00308 }
00309 
00310 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines