MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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