MOAB: Mesh Oriented datABase  (version 5.2.1)
LinearHex.cpp
Go to the documentation of this file.
00001 #include "moab/LocalDiscretization/LinearHex.hpp"
00002 #include "moab/Matrix3.hpp"
00003 #include "moab/Forward.hpp"
00004 #include <math.h>
00005 #include <limits>
00006 
00007 namespace moab
00008 {
00009 
00010 const double LinearHex::corner[8][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, { -1, 1, -1 },
00011                                          { -1, -1, 1 },  { 1, -1, 1 },  { 1, 1, 1 },  { -1, 1, 1 } };
00012 
00013 /* For each point, its weight and location are stored as an array.
00014    Hence, the inner dimension is 2, the outer dimension is gauss_count.
00015    We use a one-point Gaussian quadrature, since it integrates linear functions exactly.
00016 */
00017 const double LinearHex::gauss[1][2] = { { 2.0, 0.0 } };
00018 
00019 ErrorCode LinearHex::jacobianFcn( const double* params, const double* verts, const int /*nverts*/, const int ndim,
00020                                   double*, double* result )
00021 {
00022     assert( params && verts );
00023     Matrix3* J = reinterpret_cast< Matrix3* >( result );
00024     *J         = Matrix3( 0.0 );
00025     for( unsigned i = 0; i < 8; ++i )
00026     {
00027         const double params_p    = 1 + params[0] * corner[i][0];
00028         const double eta_p       = 1 + params[1] * corner[i][1];
00029         const double zeta_p      = 1 + params[2] * corner[i][2];
00030         const double dNi_dparams = corner[i][0] * eta_p * zeta_p;
00031         const double dNi_deta    = corner[i][1] * params_p * zeta_p;
00032         const double dNi_dzeta   = corner[i][2] * params_p * eta_p;
00033         ( *J )( 0, 0 ) += dNi_dparams * verts[i * ndim + 0];
00034         ( *J )( 1, 0 ) += dNi_dparams * verts[i * ndim + 1];
00035         ( *J )( 2, 0 ) += dNi_dparams * verts[i * ndim + 2];
00036         ( *J )( 0, 1 ) += dNi_deta * verts[i * ndim + 0];
00037         ( *J )( 1, 1 ) += dNi_deta * verts[i * ndim + 1];
00038         ( *J )( 2, 1 ) += dNi_deta * verts[i * ndim + 2];
00039         ( *J )( 0, 2 ) += dNi_dzeta * verts[i * ndim + 0];
00040         ( *J )( 1, 2 ) += dNi_dzeta * verts[i * ndim + 1];
00041         ( *J )( 2, 2 ) += dNi_dzeta * verts[i * ndim + 2];
00042     }
00043     ( *J ) *= 0.125;
00044     return MB_SUCCESS;
00045 }  // LinearHex::jacobian()
00046 
00047 ErrorCode LinearHex::evalFcn( const double* params, const double* field, const int /*ndim*/, const int num_tuples,
00048                               double*, double* result )
00049 {
00050     assert( params && field && num_tuples != -1 );
00051     for( int i = 0; i < num_tuples; i++ )
00052         result[i] = 0.0;
00053     for( unsigned i = 0; i < 8; ++i )
00054     {
00055         const double N_i =
00056             ( 1 + params[0] * corner[i][0] ) * ( 1 + params[1] * corner[i][1] ) * ( 1 + params[2] * corner[i][2] );
00057         for( int j = 0; j < num_tuples; j++ )
00058             result[j] += N_i * field[i * num_tuples + j];
00059     }
00060     for( int i = 0; i < num_tuples; i++ )
00061         result[i] *= 0.125;
00062 
00063     return MB_SUCCESS;
00064 }
00065 
00066 ErrorCode LinearHex::integrateFcn( const double* field, const double* verts, const int nverts, const int ndim,
00067                                    const int num_tuples, double* work, double* result )
00068 {
00069     assert( field && verts && num_tuples != -1 );
00070     double tmp_result[8];
00071     ErrorCode rval = MB_SUCCESS;
00072     for( int i = 0; i < num_tuples; i++ )
00073         result[i] = 0.0;
00074     CartVect x;
00075     Matrix3 J;
00076     for( unsigned int j1 = 0; j1 < LinearHex::gauss_count; ++j1 )
00077     {
00078         x[0]      = LinearHex::gauss[j1][1];
00079         double w1 = LinearHex::gauss[j1][0];
00080         for( unsigned int j2 = 0; j2 < LinearHex::gauss_count; ++j2 )
00081         {
00082             x[1]      = LinearHex::gauss[j2][1];
00083             double w2 = LinearHex::gauss[j2][0];
00084             for( unsigned int j3 = 0; j3 < LinearHex::gauss_count; ++j3 )
00085             {
00086                 x[2]      = LinearHex::gauss[j3][1];
00087                 double w3 = LinearHex::gauss[j3][0];
00088                 rval      = evalFcn( x.array(), field, ndim, num_tuples, NULL, tmp_result );
00089                 if( MB_SUCCESS != rval ) return rval;
00090                 rval = jacobianFcn( x.array(), verts, nverts, ndim, work, J[0] );
00091                 if( MB_SUCCESS != rval ) return rval;
00092                 double tmp_det = w1 * w2 * w3 * J.determinant();
00093                 for( int i = 0; i < num_tuples; i++ )
00094                     result[i] += tmp_result[i] * tmp_det;
00095             }
00096         }
00097     }
00098 
00099     return MB_SUCCESS;
00100 }  // LinearHex::integrate_vector()
00101 
00102 ErrorCode LinearHex::reverseEvalFcn( EvalFcn eval, JacobianFcn jacob, InsideFcn ins, const double* posn,
00103                                      const double* verts, const int nverts, const int ndim, const double iter_tol,
00104                                      const double inside_tol, double* work, double* params, int* is_inside )
00105 {
00106     assert( posn && verts );
00107     return EvalSet::evaluate_reverse( eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, params,
00108                                       is_inside );
00109 }
00110 
00111 int LinearHex::insideFcn( const double* params, const int ndim, const double tol )
00112 {
00113     return EvalSet::inside_function( params, ndim, tol );
00114 }
00115 
00116 ErrorCode LinearHex::normalFcn( const int ientDim, const int facet, const int nverts, const double* verts,
00117                                 double normal[3] )
00118 {
00119     // assert(facet < 6 && ientDim == 2 && nverts == 8);
00120     if( nverts != 8 ) MB_SET_ERR( MB_FAILURE, "Incorrect vertex count for passed hex :: expected value = 8 " );
00121     if( ientDim != 2 ) MB_SET_ERR( MB_FAILURE, "Requesting normal for unsupported dimension :: expected value = 2 " );
00122     if( facet > 6 || facet < 0 ) MB_SET_ERR( MB_FAILURE, "Incorrect local face id :: expected value = one of 0-5" );
00123 
00124     int id0 = CN::mConnectivityMap[MBHEX][ientDim - 1].conn[facet][0];
00125     int id1 = CN::mConnectivityMap[MBHEX][ientDim - 1].conn[facet][1];
00126     int id2 = CN::mConnectivityMap[MBHEX][ientDim - 1].conn[facet][3];
00127 
00128     double x0[3], x1[3];
00129 
00130     for( int i = 0; i < 3; i++ )
00131     {
00132         x0[i] = verts[3 * id1 + i] - verts[3 * id0 + i];
00133         x1[i] = verts[3 * id2 + i] - verts[3 * id0 + i];
00134     }
00135 
00136     double a   = x0[1] * x1[2] - x1[1] * x0[2];
00137     double b   = x1[0] * x0[2] - x0[0] * x1[2];
00138     double c   = x0[0] * x1[1] - x1[0] * x0[1];
00139     double nrm = sqrt( a * a + b * b + c * c );
00140 
00141     if( nrm > std::numeric_limits< double >::epsilon() )
00142     {
00143         normal[0] = a / nrm;
00144         normal[1] = b / nrm;
00145         normal[2] = c / nrm;
00146     }
00147     return MB_SUCCESS;
00148 }
00149 
00150 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines