MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #include "moab/LocalDiscretization/LinearQuad.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 LinearQuad::corner[4][2] = { { -1, -1 }, { 1, -1 }, { 1, 1 }, { -1, 1 } }; 00011 00012 /* For each point, its weight and location are stored as an array. 00013 Hence, the inner dimension is 2, the outer dimension is gauss_count. 00014 We use a one-point Gaussian quadrature, since it integrates linear functions exactly. 00015 */ 00016 const double LinearQuad::gauss[1][2] = { { 2.0, 0.0 } }; 00017 00018 ErrorCode LinearQuad::jacobianFcn( const double* params, const double* verts, const int /*nverts*/, const int /*ndim*/, 00019 double*, double* result ) 00020 { 00021 Matrix3* J = reinterpret_cast< Matrix3* >( result ); 00022 *J = Matrix3( 0.0 ); 00023 for( unsigned i = 0; i < 4; ++i ) 00024 { 00025 const double xi_p = 1 + params[0] * corner[i][0]; 00026 const double eta_p = 1 + params[1] * corner[i][1]; 00027 const double dNi_dxi = corner[i][0] * eta_p; 00028 const double dNi_deta = corner[i][1] * xi_p; 00029 ( *J )( 0, 0 ) += dNi_dxi * verts[i * 3 + 0]; 00030 ( *J )( 1, 0 ) += dNi_dxi * verts[i * 3 + 1]; 00031 ( *J )( 0, 1 ) += dNi_deta * verts[i * 3 + 0]; 00032 ( *J )( 1, 1 ) += dNi_deta * verts[i * 3 + 1]; 00033 } 00034 ( *J ) *= 0.25; 00035 ( *J )( 2, 2 ) = 1.0; /* to make sure the Jacobian determinant is non-zero */ 00036 return MB_SUCCESS; 00037 } // LinearQuad::jacobian() 00038 00039 ErrorCode LinearQuad::evalFcn( const double* params, const double* field, const int /*ndim*/, const int num_tuples, 00040 double*, double* result ) 00041 { 00042 for( int i = 0; i < num_tuples; i++ ) 00043 result[i] = 0.0; 00044 for( unsigned i = 0; i < 4; ++i ) 00045 { 00046 const double N_i = ( 1 + params[0] * corner[i][0] ) * ( 1 + params[1] * corner[i][1] ); 00047 for( int j = 0; j < num_tuples; j++ ) 00048 result[j] += N_i * field[i * num_tuples + j]; 00049 } 00050 for( int i = 0; i < num_tuples; i++ ) 00051 result[i] *= 0.25; 00052 00053 return MB_SUCCESS; 00054 } 00055 00056 ErrorCode LinearQuad::integrateFcn( const double* field, const double* verts, const int nverts, const int ndim, 00057 const int num_tuples, double* work, double* result ) 00058 { 00059 double tmp_result[4]; 00060 ErrorCode rval = MB_SUCCESS; 00061 for( int i = 0; i < num_tuples; i++ ) 00062 result[i] = 0.0; 00063 CartVect x; 00064 Matrix3 J; 00065 for( unsigned int j1 = 0; j1 < LinearQuad::gauss_count; ++j1 ) 00066 { 00067 x[0] = LinearQuad::gauss[j1][1]; 00068 double w1 = LinearQuad::gauss[j1][0]; 00069 for( unsigned int j2 = 0; j2 < LinearQuad::gauss_count; ++j2 ) 00070 { 00071 x[1] = LinearQuad::gauss[j2][1]; 00072 double w2 = LinearQuad::gauss[j2][0]; 00073 rval = evalFcn( x.array(), field, ndim, num_tuples, NULL, tmp_result ); 00074 if( MB_SUCCESS != rval ) return rval; 00075 rval = jacobianFcn( x.array(), verts, nverts, ndim, work, J[0] ); 00076 if( MB_SUCCESS != rval ) return rval; 00077 double tmp_det = w1 * w2 * J.determinant(); 00078 for( int i = 0; i < num_tuples; i++ ) 00079 result[i] += tmp_result[i] * tmp_det; 00080 } 00081 } 00082 return MB_SUCCESS; 00083 } // LinearHex::integrate_vector() 00084 00085 ErrorCode LinearQuad::reverseEvalFcn( EvalFcn eval, JacobianFcn jacob, InsideFcn ins, const double* posn, 00086 const double* verts, const int nverts, const int ndim, const double iter_tol, 00087 const double inside_tol, double* work, double* params, int* is_inside ) 00088 { 00089 return EvalSet::evaluate_reverse( eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, params, 00090 is_inside ); 00091 } 00092 00093 int LinearQuad::insideFcn( const double* params, const int ndim, const double tol ) 00094 { 00095 return EvalSet::inside_function( params, ndim, tol ); 00096 } 00097 00098 ErrorCode LinearQuad::normalFcn( const int ientDim, const int facet, const int nverts, const double* verts, 00099 double normal[3] ) 00100 { 00101 // assert(facet <4 && ientDim == 1 && nverts==4); 00102 if( nverts != 4 ) MB_SET_ERR( MB_FAILURE, "Incorrect vertex count for passed quad :: expected value = 4" ); 00103 if( ientDim != 1 ) MB_SET_ERR( MB_FAILURE, "Requesting normal for unsupported dimension :: expected value = 1 " ); 00104 if( facet > 4 || facet < 0 ) MB_SET_ERR( MB_FAILURE, "Incorrect local edge id :: expected value = one of 0-3" ); 00105 00106 // Get the local vertex ids of local edge 00107 int id0 = CN::mConnectivityMap[MBQUAD][ientDim - 1].conn[facet][0]; 00108 int id1 = CN::mConnectivityMap[MBQUAD][ientDim - 1].conn[facet][1]; 00109 00110 // Find a vector along the edge 00111 double edge[3]; 00112 for( int i = 0; i < 3; i++ ) 00113 { 00114 edge[i] = verts[3 * id1 + i] - verts[3 * id0 + i]; 00115 } 00116 // Find the normal of the face 00117 double x0[3], x1[3], fnrm[3]; 00118 for( int i = 0; i < 3; i++ ) 00119 { 00120 x0[i] = verts[3 * 1 + i] - verts[3 * 0 + i]; 00121 x1[i] = verts[3 * 3 + i] - verts[3 * 0 + i]; 00122 } 00123 fnrm[0] = x0[1] * x1[2] - x1[1] * x0[2]; 00124 fnrm[1] = x1[0] * x0[2] - x0[0] * x1[2]; 00125 fnrm[2] = x0[0] * x1[1] - x1[0] * x0[1]; 00126 00127 // Find the normal of the edge as the cross product of edge and face normal 00128 00129 double a = edge[1] * fnrm[2] - fnrm[1] * edge[2]; 00130 double b = edge[2] * fnrm[0] - fnrm[2] * edge[0]; 00131 double c = edge[0] * fnrm[1] - fnrm[0] * edge[1]; 00132 double nrm = sqrt( a * a + b * b + c * c ); 00133 00134 if( nrm > std::numeric_limits< double >::epsilon() ) 00135 { 00136 normal[0] = a / nrm; 00137 normal[1] = b / nrm; 00138 normal[2] = c / nrm; 00139 } 00140 return MB_SUCCESS; 00141 } 00142 00143 } // namespace moab