MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #include "moab/LocalDiscretization/LinearTri.hpp" 00002 #include "moab/Forward.hpp" 00003 #include <algorithm> 00004 #include <math.h> 00005 #include <limits> 00006 00007 namespace moab 00008 { 00009 00010 const double LinearTri::corner[3][2] = { { 0, 0 }, { 1, 0 }, { 0, 1 } }; 00011 00012 ErrorCode LinearTri::initFcn( const double* verts, const int nverts, double*& work ) 00013 { 00014 // allocate work array as: 00015 // work[0..8] = T 00016 // work[9..17] = Tinv 00017 // work[18] = detT 00018 // work[19] = detTinv 00019 assert( nverts == 3 && verts ); 00020 if( !work ) work = new double[20]; 00021 00022 Matrix3 J( verts[1 * 3 + 0] - verts[0 * 3 + 0], verts[2 * 3 + 0] - verts[0 * 3 + 0], 0.0, 00023 verts[1 * 3 + 1] - verts[0 * 3 + 1], verts[2 * 3 + 1] - verts[0 * 3 + 1], 0.0, 00024 verts[1 * 3 + 2] - verts[0 * 3 + 2], verts[2 * 3 + 2] - verts[0 * 3 + 2], 1.0 ); 00025 J *= 0.5; 00026 00027 J.copyto( work ); 00028 J.inverse().copyto( work + Matrix3::size ); 00029 work[18] = J.determinant(); 00030 work[19] = ( work[18] < 1e-12 ? std::numeric_limits< double >::max() : 1.0 / work[18] ); 00031 00032 return MB_SUCCESS; 00033 } 00034 00035 ErrorCode LinearTri::evalFcn( const double* params, const double* field, const int /*ndim*/, const int num_tuples, 00036 double* /*work*/, double* result ) 00037 { 00038 assert( params && field && num_tuples > 0 ); 00039 // convert to [0,1] 00040 double p1 = 0.5 * ( 1.0 + params[0] ), p2 = 0.5 * ( 1.0 + params[1] ), p0 = 1.0 - p1 - p2; 00041 00042 for( int j = 0; j < num_tuples; j++ ) 00043 result[j] = p0 * field[0 * num_tuples + j] + p1 * field[1 * num_tuples + j] + p2 * field[2 * num_tuples + j]; 00044 00045 return MB_SUCCESS; 00046 } 00047 00048 ErrorCode LinearTri::integrateFcn( const double* field, const double* /*verts*/, const int nverts, const int /*ndim*/, 00049 const int num_tuples, double* work, double* result ) 00050 { 00051 assert( field && num_tuples > 0 ); 00052 std::fill( result, result + num_tuples, 0.0 ); 00053 for( int i = 0; i < nverts; ++i ) 00054 { 00055 for( int j = 0; j < num_tuples; j++ ) 00056 result[j] += field[i * num_tuples + j]; 00057 } 00058 double tmp = work[18] / 6.0; 00059 for( int i = 0; i < num_tuples; i++ ) 00060 result[i] *= tmp; 00061 00062 return MB_SUCCESS; 00063 } 00064 00065 ErrorCode LinearTri::jacobianFcn( const double*, const double*, const int, const int, double* work, double* result ) 00066 { 00067 // jacobian is cached in work array 00068 assert( work ); 00069 std::copy( work, work + 9, result ); 00070 return MB_SUCCESS; 00071 } 00072 00073 ErrorCode LinearTri::reverseEvalFcn( EvalFcn eval, JacobianFcn jacob, InsideFcn ins, const double* posn, 00074 const double* verts, const int nverts, const int ndim, const double iter_tol, 00075 const double inside_tol, double* work, double* params, int* is_inside ) 00076 { 00077 assert( posn && verts ); 00078 return evaluate_reverse( eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, params, 00079 is_inside ); 00080 } 00081 00082 int LinearTri::insideFcn( const double* params, const int, const double tol ) 00083 { 00084 return ( params[0] >= -1.0 - tol && params[1] >= -1.0 - tol && params[0] + params[1] <= 1.0 + tol ); 00085 } 00086 00087 ErrorCode LinearTri::evaluate_reverse( EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f, const double* posn, 00088 const double* verts, const int nverts, const int ndim, const double iter_tol, 00089 const double inside_tol, double* work, double* params, int* inside ) 00090 { 00091 // TODO: should differentiate between epsilons used for 00092 // Newton Raphson iteration, and epsilons used for curved boundary geometry errors 00093 // right now, fix the tolerance used for NR 00094 const double error_tol_sqr = iter_tol * iter_tol; 00095 CartVect* cvparams = reinterpret_cast< CartVect* >( params ); 00096 const CartVect* cvposn = reinterpret_cast< const CartVect* >( posn ); 00097 00098 // find best initial guess to improve convergence 00099 CartVect tmp_params[] = { CartVect( -1, -1, -1 ), CartVect( 1, -1, -1 ), CartVect( -1, 1, -1 ) }; 00100 double resl = std::numeric_limits< double >::max(); 00101 CartVect new_pos, tmp_pos; 00102 ErrorCode rval; 00103 for( unsigned int i = 0; i < 3; i++ ) 00104 { 00105 rval = ( *eval )( tmp_params[i].array(), verts, ndim, 3, work, tmp_pos.array() ); 00106 if( MB_SUCCESS != rval ) return rval; 00107 double tmp_resl = ( tmp_pos - *cvposn ).length_squared(); 00108 if( tmp_resl < resl ) 00109 { 00110 *cvparams = tmp_params[i]; 00111 new_pos = tmp_pos; 00112 resl = tmp_resl; 00113 } 00114 } 00115 00116 // residual is diff between old and new pos; need to minimize that 00117 CartVect res = new_pos - *cvposn; 00118 Matrix3 J; 00119 rval = ( *jacob )( cvparams->array(), verts, nverts, ndim, work, J[0] ); 00120 #ifndef NDEBUG 00121 double det = J.determinant(); 00122 assert( det > std::numeric_limits< double >::epsilon() ); 00123 #endif 00124 Matrix3 Ji = J.inverse(); 00125 00126 int iters = 0; 00127 // while |res| larger than tol 00128 while( res % res > error_tol_sqr ) 00129 { 00130 if( ++iters > 25 ) return MB_FAILURE; 00131 00132 // new params tries to eliminate residual 00133 *cvparams -= Ji * res; 00134 00135 // get the new forward-evaluated position, and its difference from the target pt 00136 rval = ( *eval )( params, verts, ndim, 3, work, new_pos.array() ); 00137 if( MB_SUCCESS != rval ) return rval; 00138 res = new_pos - *cvposn; 00139 } 00140 00141 if( inside ) *inside = ( *inside_f )( params, ndim, inside_tol ); 00142 00143 return MB_SUCCESS; 00144 } // Map::evaluate_reverse() 00145 00146 /* ErrorCode LinearTri::get_normal( int facet, double *work, double *normal) 00147 { 00148 ErrorCode error; 00149 //Get the local vertex ids of local edge 00150 int id1 = ledges[facet][0]; 00151 int id2 = ledges[facet][1]; 00152 00153 //Find the normal to the face 00154 double face_normal[3]; 00155 00156 00157 }*/ 00158 00159 ErrorCode LinearTri::normalFcn( const int ientDim, const int facet, const int nverts, const double* verts, 00160 double normal[3] ) 00161 { 00162 // assert(facet < 3 && ientDim == 1 && nverts==3); 00163 if( nverts != 3 ) MB_SET_ERR( MB_FAILURE, "Incorrect vertex count for passed triangle :: expected value = 3 " ); 00164 if( ientDim != 1 ) MB_SET_ERR( MB_FAILURE, "Requesting normal for unsupported dimension :: expected value = 1 " ); 00165 if( facet > 3 || facet < 0 ) MB_SET_ERR( MB_FAILURE, "Incorrect local edge id :: expected value = one of 0-2" ); 00166 00167 // Get the local vertex ids of local edge 00168 int id0 = CN::mConnectivityMap[MBTRI][ientDim - 1].conn[facet][0]; 00169 int id1 = CN::mConnectivityMap[MBTRI][ientDim - 1].conn[facet][1]; 00170 00171 // Find a vector along the edge 00172 double edge[3]; 00173 for( int i = 0; i < 3; i++ ) 00174 { 00175 edge[i] = verts[3 * id1 + i] - verts[3 * id0 + i]; 00176 } 00177 // Find the normal of the face 00178 double x0[3], x1[3], fnrm[3]; 00179 for( int i = 0; i < 3; i++ ) 00180 { 00181 x0[i] = verts[3 * 1 + i] - verts[3 * 0 + i]; 00182 x1[i] = verts[3 * 2 + i] - verts[3 * 0 + i]; 00183 } 00184 fnrm[0] = x0[1] * x1[2] - x1[1] * x0[2]; 00185 fnrm[1] = x1[0] * x0[2] - x0[0] * x1[2]; 00186 fnrm[2] = x0[0] * x1[1] - x1[0] * x0[1]; 00187 00188 // Find the normal of the edge as the cross product of edge and face normal 00189 00190 double a = edge[1] * fnrm[2] - fnrm[1] * edge[2]; 00191 double b = edge[2] * fnrm[0] - fnrm[2] * edge[0]; 00192 double c = edge[0] * fnrm[1] - fnrm[0] * edge[1]; 00193 double nrm = sqrt( a * a + b * b + c * c ); 00194 00195 if( nrm > std::numeric_limits< double >::epsilon() ) 00196 { 00197 normal[0] = a / nrm; 00198 normal[1] = b / nrm; 00199 normal[2] = c / nrm; 00200 } 00201 return MB_SUCCESS; 00202 } 00203 00204 } // namespace moab