MOAB: Mesh Oriented datABase  (version 5.4.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 <cmath>
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,
00020                                   const double* verts,
00021                                   const int /*nverts*/,
00022                                   const int ndim,
00023                                   double*,
00024                                   double* result )
00025 {
00026     assert( params && verts );
00027     Matrix3* J = reinterpret_cast< Matrix3* >( result );
00028     *J         = Matrix3( 0.0 );
00029     for( unsigned i = 0; i < 8; ++i )
00030     {
00031         const double params_p    = 1 + params[0] * corner[i][0];
00032         const double eta_p       = 1 + params[1] * corner[i][1];
00033         const double zeta_p      = 1 + params[2] * corner[i][2];
00034         const double dNi_dparams = corner[i][0] * eta_p * zeta_p;
00035         const double dNi_deta    = corner[i][1] * params_p * zeta_p;
00036         const double dNi_dzeta   = corner[i][2] * params_p * eta_p;
00037         ( *J )( 0, 0 ) += dNi_dparams * verts[i * ndim + 0];
00038         ( *J )( 1, 0 ) += dNi_dparams * verts[i * ndim + 1];
00039         ( *J )( 2, 0 ) += dNi_dparams * verts[i * ndim + 2];
00040         ( *J )( 0, 1 ) += dNi_deta * verts[i * ndim + 0];
00041         ( *J )( 1, 1 ) += dNi_deta * verts[i * ndim + 1];
00042         ( *J )( 2, 1 ) += dNi_deta * verts[i * ndim + 2];
00043         ( *J )( 0, 2 ) += dNi_dzeta * verts[i * ndim + 0];
00044         ( *J )( 1, 2 ) += dNi_dzeta * verts[i * ndim + 1];
00045         ( *J )( 2, 2 ) += dNi_dzeta * verts[i * ndim + 2];
00046     }
00047     ( *J ) *= 0.125;
00048     return MB_SUCCESS;
00049 }  // LinearHex::jacobian()
00050
00051 ErrorCode LinearHex::evalFcn( const double* params,
00052                               const double* field,
00053                               const int /*ndim*/,
00054                               const int num_tuples,
00055                               double*,
00056                               double* result )
00057 {
00058     assert( params && field && num_tuples != -1 );
00059     for( int i = 0; i < num_tuples; i++ )
00060         result[i] = 0.0;
00061     for( unsigned i = 0; i < 8; ++i )
00062     {
00063         const double N_i =
00064             ( 1 + params[0] * corner[i][0] ) * ( 1 + params[1] * corner[i][1] ) * ( 1 + params[2] * corner[i][2] );
00065         for( int j = 0; j < num_tuples; j++ )
00066             result[j] += N_i * field[i * num_tuples + j];
00067     }
00068     for( int i = 0; i < num_tuples; i++ )
00069         result[i] *= 0.125;
00070
00071     return MB_SUCCESS;
00072 }
00073
00074 ErrorCode LinearHex::integrateFcn( const double* field,
00075                                    const double* verts,
00076                                    const int nverts,
00077                                    const int ndim,
00078                                    const int num_tuples,
00079                                    double* work,
00080                                    double* result )
00081 {
00082     assert( field && verts && num_tuples != -1 );
00083     double tmp_result[8];
00084     ErrorCode rval = MB_SUCCESS;
00085     for( int i = 0; i < num_tuples; i++ )
00086         result[i] = 0.0;
00087     CartVect x;
00088     Matrix3 J;
00089     for( unsigned int j1 = 0; j1 < LinearHex::gauss_count; ++j1 )
00090     {
00091         x[0]      = LinearHex::gauss[j1][1];
00092         double w1 = LinearHex::gauss[j1][0];
00093         for( unsigned int j2 = 0; j2 < LinearHex::gauss_count; ++j2 )
00094         {
00095             x[1]      = LinearHex::gauss[j2][1];
00096             double w2 = LinearHex::gauss[j2][0];
00097             for( unsigned int j3 = 0; j3 < LinearHex::gauss_count; ++j3 )
00098             {
00099                 x[2]      = LinearHex::gauss[j3][1];
00100                 double w3 = LinearHex::gauss[j3][0];
00101                 rval      = evalFcn( x.array(), field, ndim, num_tuples, NULL, tmp_result );
00102                 if( MB_SUCCESS != rval ) return rval;
00103                 rval = jacobianFcn( x.array(), verts, nverts, ndim, work, J[0] );
00104                 if( MB_SUCCESS != rval ) return rval;
00105                 double tmp_det = w1 * w2 * w3 * J.determinant();
00106                 for( int i = 0; i < num_tuples; i++ )
00107                     result[i] += tmp_result[i] * tmp_det;
00108             }
00109         }
00110     }
00111
00112     return MB_SUCCESS;
00113 }  // LinearHex::integrate_vector()
00114
00115 ErrorCode LinearHex::reverseEvalFcn( EvalFcn eval,
00116                                      JacobianFcn jacob,
00117                                      InsideFcn ins,
00118                                      const double* posn,
00119                                      const double* verts,
00120                                      const int nverts,
00121                                      const int ndim,
00122                                      const double iter_tol,
00123                                      const double inside_tol,
00124                                      double* work,
00125                                      double* params,
00126                                      int* is_inside )
00127 {
00128     assert( posn && verts );
00129     return EvalSet::evaluate_reverse( eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, params,
00130                                       is_inside );
00131 }
00132
00133 int LinearHex::insideFcn( const double* params, const int ndim, const double tol )
00134 {
00135     return EvalSet::inside_function( params, ndim, tol );
00136 }
00137
00138 ErrorCode LinearHex::normalFcn( const int ientDim,
00139                                 const int facet,
00140                                 const int nverts,
00141                                 const double* verts,
00142                                 double normal[3] )
00143 {
00144     // assert(facet < 6 && ientDim == 2 && nverts == 8);
00145     if( nverts != 8 ) MB_SET_ERR( MB_FAILURE, "Incorrect vertex count for passed hex :: expected value = 8 " );
00146     if( ientDim != 2 ) MB_SET_ERR( MB_FAILURE, "Requesting normal for unsupported dimension :: expected value = 2 " );
00147     if( facet > 6 || facet < 0 ) MB_SET_ERR( MB_FAILURE, "Incorrect local face id :: expected value = one of 0-5" );
00148
00149     int id0 = CN::mConnectivityMap[MBHEX][ientDim - 1].conn[facet][0];
00150     int id1 = CN::mConnectivityMap[MBHEX][ientDim - 1].conn[facet][1];
00151     int id2 = CN::mConnectivityMap[MBHEX][ientDim - 1].conn[facet][3];
00152
00153     double x0[3], x1[3];
00154
00155     for( int i = 0; i < 3; i++ )
00156     {
00157         x0[i] = verts[3 * id1 + i] - verts[3 * id0 + i];
00158         x1[i] = verts[3 * id2 + i] - verts[3 * id0 + i];
00159     }
00160
00161     double a   = x0[1] * x1[2] - x1[1] * x0[2];
00162     double b   = x1[0] * x0[2] - x0[0] * x1[2];
00163     double c   = x0[0] * x1[1] - x1[0] * x0[1];
00164     double nrm = sqrt( a * a + b * b + c * c );
00165
00166     if( nrm > std::numeric_limits< double >::epsilon() )
00167     {
00168         normal[0] = a / nrm;
00169         normal[1] = b / nrm;
00170         normal[2] = c / nrm;
00171     }
00172     return MB_SUCCESS;
00173 }
00174
00175 }  // namespace moab