verdict_defines.hpp

Go to the documentation of this file.

/*=========================================================================

  Module:    $RCSfile: verdict_defines.hpp,v $

  Copyright (c) 2006 Sandia Corporation.
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/

/*
 * verdict_defines.cpp contains common definitions
 *
 * This file is part of VERDICT
 *
 */

#ifndef VERDICT_DEFINES
#define VERDICT_DEFINES

#include <cmath>
#include "v_vector.h"
#include "VerdictVector.hpp"

enum VerdictBoolean
{
    VERDICT_FALSE = 0,
    VERDICT_TRUE  = 1
};

#define VERDICT_MIN( a, b ) ( ( a ) < ( b ) ? ( a ) : ( b ) )
#define VERDICT_MAX( a, b ) ( ( a ) > ( b ) ? ( a ) : ( b ) )

inline double determinant( double a, double b, double c, double d )
{
    return ( ( a ) * ( d ) - ( b ) * ( c ) );
}

inline double determinant( const VerdictVector& v1, const VerdictVector& v2, const VerdictVector& v3 )
{
    return v1 % ( v2 * v3 );
}

#define jacobian_matrix( a, b, c, d, e, f, g ) \
    double jac_mat_tmp;                        \
    jac_mat_tmp = sqrt( a );                   \
    if( jac_mat_tmp == 0 )                     \
    {                                          \
        ( d ) = 0;                             \
        ( e ) = 0;                             \
        ( f ) = 0;                             \
        ( g ) = 0;                             \
    }                                          \
    else                                       \
    {                                          \
        ( d ) = jac_mat_tmp;                   \
        ( e ) = 0;                             \
        ( f ) = ( b ) / jac_mat_tmp;           \
        ( g ) = ( c ) / jac_mat_tmp;           \
    }

// this assumes that detmw != 0
#define form_t( m11, m21, m12, m22, mw11, mw21, mw12, mw22, detmw, xm11, xm21, xm12, xm22 ) \
    xm11     = ( ( m11 ) * ( mw22 ) - ( m12 ) * ( mw21 ) ) / ( detmw );                     \
    ( xm21 ) = ( ( m21 ) * ( mw22 ) - ( m22 ) * ( mw21 ) ) / ( detmw );                     \
    ( xm12 ) = ( ( m12 ) * ( mw11 ) - ( m11 ) * ( mw12 ) ) / ( detmw );                     \
    ( xm22 ) = ( ( m22 ) * ( mw11 ) - ( m21 ) * ( mw12 ) ) / ( detmw );

extern double verdictSqrt2;

inline double normalize_jacobian( double jacobi,
                                  VerdictVector& v1,
                                  VerdictVector& v2,
                                  VerdictVector& v3,
                                  int tet_flag = 0 )
{
    double return_value = 0.0;

    if( jacobi != 0.0 )
    {

        double l1, l2, l3, length_product;
        // Note: there may be numerical problems if one is a lot shorter
        // than the others this way. But scaling each vector before the
        // triple product would involve 3 square roots instead of just
        // one.
        l1             = v1.length_squared();
        l2             = v2.length_squared();
        l3             = v3.length_squared();
        length_product = sqrt( l1 * l2 * l3 );

        // if some numerical scaling problem, or just plain roundoff,
        // then push back into range [-1,1].
        if( length_product < fabs( jacobi ) )
        {
            length_product = fabs( jacobi );
        }

        if( tet_flag == 1 )
            return_value = verdictSqrt2 * jacobi / length_product;
        else
            return_value = jacobi / length_product;
    }
    return return_value;
}

inline double norm_squared( double m11, double m21, double m12, double m22 )
{
    return m11 * m11 + m21 * m21 + m12 * m12 + m22 * m22;
}

#define metric_matrix( m11, m21, m12, m22, gm11, gm12, gm22 ) \
    gm11     = ( m11 ) * ( m11 ) + ( m21 ) * ( m21 );         \
    ( gm12 ) = ( m11 ) * ( m12 ) + ( m21 ) * ( m22 );         \
    ( gm22 ) = ( m12 ) * ( m12 ) + ( m22 ) * ( m22 );

inline int skew_matrix( double gm11,
                        double gm12,
                        double gm22,
                        double det,
                        double& qm11,
                        double& qm21,
                        double& qm12,
                        double& qm22 )
{
    double tmp = sqrt( gm11 * gm22 );
    if( tmp == 0 )
    {
        return false;
    }

    qm11 = 1;
    qm21 = 0;
    qm12 = gm12 / tmp;
    qm22 = det / tmp;
    return true;
}

inline void inverse( const VerdictVector& x1,
                     const VerdictVector& x2,
                     const VerdictVector& x3,
                     VerdictVector& u1,
                     VerdictVector& u2,
                     VerdictVector& u3 )
{
    double detx = determinant( x1, x2, x3 );
    VerdictVector rx1, rx2, rx3;

    rx1.set( x1.x(), x2.x(), x3.x() );
    rx2.set( x1.y(), x2.y(), x3.y() );
    rx3.set( x1.z(), x2.z(), x3.z() );

    u1 = rx2 * rx3;
    u2 = rx3 * rx1;
    u3 = rx1 * rx2;

    u1 /= detx;
    u2 /= detx;
    u3 /= detx;
}

/*
inline void form_T(double a1[3],
  double a2[3],
  double a3[3],
  double w1[3],
  double w2[3],
  double w3[3],
  double t1[3],
  double t2[3],
  double t3[3] )
{

  double x1[3], x2[3], x3[3];
  double ra1[3], ra2[3], ra3[3];

  x1[0] = a1[0];
  x1[1] = a2[0];
  x1[2] = a3[0];

  x2[0] = a1[1];
  x2[1] = a2[1];
  x2[2] = a3[1];

  x3[0] = a1[2];
  x3[1] = a2[2];
  x3[2] = a3[2];

  inverse(w1,w2,w3,x1,x2,x3);

  t1[0] = dot_product(ra1, x1);
  t1[1] = dot_product(ra1, x2);
  t1[2] = dot_product(ra1, x3);

  t2[0] = dot_product(ra2, x1);
  t2[0] = dot_product(ra2, x2);
  t2[0] = dot_product(ra2, x3);

  t3[0] = dot_product(ra3, x1);
  t3[0] = dot_product(ra3, x2);
  t3[0] = dot_product(ra3, x3);

}
*/

inline void form_Q( const VerdictVector& v1,
                    const VerdictVector& v2,
                    const VerdictVector& v3,
                    VerdictVector& q1,
                    VerdictVector& q2,
                    VerdictVector& q3 )
{

    double g11, g12, g13, g22, g23, g33;

    g11 = v1 % v1;
    g12 = v1 % v2;
    g13 = v1 % v3;
    g22 = v2 % v2;
    g23 = v2 % v3;
    g33 = v3 % v3;

    double rtg11 = sqrt( g11 );
    double rtg22 = sqrt( g22 );
    double rtg33 = sqrt( g33 );
    VerdictVector temp1;

    temp1 = v1 * v2;

    double cross = sqrt( temp1 % temp1 );

    double q11, q21, q31;
    double q12, q22, q32;
    double q13, q23, q33;

    q11 = 1;
    q21 = 0;
    q31 = 0;

    q12 = g12 / rtg11 / rtg22;
    q22 = cross / rtg11 / rtg22;
    q32 = 0;

    q13   = g13 / rtg11 / rtg33;
    q23   = ( g11 * g23 - g12 * g13 ) / rtg11 / rtg33 / cross;
    temp1 = v2 * v3;
    q33   = ( v1 % temp1 ) / rtg33 / cross;

    q1.set( q11, q21, q31 );
    q2.set( q12, q22, q32 );
    q3.set( q13, q23, q33 );
}

inline void product( VerdictVector& a1,
                     VerdictVector& a2,
                     VerdictVector& a3,
                     VerdictVector& b1,
                     VerdictVector& b2,
                     VerdictVector& b3,
                     VerdictVector& c1,
                     VerdictVector& c2,
                     VerdictVector& c3 )
{

    VerdictVector x1, x2, x3;

    x1.set( a1.x(), a2.x(), a3.x() );
    x2.set( a1.y(), a2.y(), a3.y() );
    x3.set( a1.z(), a2.z(), a3.z() );

    c1.set( x1 % b1, x2 % b1, x3 % b1 );
    c2.set( x1 % b2, x2 % b2, x3 % b2 );
    c3.set( x1 % b3, x2 % b3, x3 % b3 );
}

inline double norm_squared( VerdictVector& x1, VerdictVector& x2, VerdictVector& x3 )

{
    return ( x1 % x1 ) + ( x2 % x2 ) + ( x3 % x3 );
}

inline double skew_x( VerdictVector& q1,
                      VerdictVector& q2,
                      VerdictVector& q3,
                      VerdictVector& qw1,
                      VerdictVector& qw2,
                      VerdictVector& qw3 )
{
    double normsq1, normsq2, kappa;
    VerdictVector u1, u2, u3;
    VerdictVector x1, x2, x3;

    inverse( qw1, qw2, qw3, u1, u2, u3 );
    product( q1, q2, q3, u1, u2, u3, x1, x2, x3 );
    inverse( x1, x2, x3, u1, u2, u3 );
    normsq1 = norm_squared( x1, x2, x3 );
    normsq2 = norm_squared( u1, u2, u3 );
    kappa   = sqrt( normsq1 * normsq2 );

    double skew = 0;
    if( kappa > VERDICT_DBL_MIN ) skew = 3 / kappa;

    return skew;
}

#endif