MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #include "moab/LocalDiscretization/LinearTet.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 LinearTet::corner[4][3] = { { 0, 0, 0 }, { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } }; 00011 00012 ErrorCode LinearTet::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 == 4 && 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], 00023 verts[3 * 3 + 0] - verts[0 * 3 + 0], verts[1 * 3 + 1] - verts[0 * 3 + 1], 00024 verts[2 * 3 + 1] - verts[0 * 3 + 1], verts[3 * 3 + 1] - verts[0 * 3 + 1], 00025 verts[1 * 3 + 2] - verts[0 * 3 + 2], verts[2 * 3 + 2] - verts[0 * 3 + 2], 00026 verts[3 * 3 + 2] - verts[0 * 3 + 2] ); 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 LinearTet::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 std::vector< double > f0( num_tuples ); 00040 std::copy( field, field + num_tuples, f0.begin() ); 00041 std::copy( field, field + num_tuples, result ); 00042 00043 for( unsigned i = 1; i < 4; ++i ) 00044 { 00045 double p = 0.5 * ( params[i - 1] + 1 ); // transform from -1 <= p <= 1 to 0 <= p <= 1 00046 for( int j = 0; j < num_tuples; j++ ) 00047 result[j] += ( field[i * num_tuples + j] - f0[j] ) * p; 00048 } 00049 00050 return MB_SUCCESS; 00051 } 00052 00053 ErrorCode LinearTet::integrateFcn( const double* field, const double* /*verts*/, const int nverts, const int /*ndim*/, 00054 const int num_tuples, double* work, double* result ) 00055 { 00056 assert( field && num_tuples > 0 ); 00057 std::fill( result, result + num_tuples, 0.0 ); 00058 for( int i = 0; i < nverts; ++i ) 00059 { 00060 for( int j = 0; j < num_tuples; j++ ) 00061 result[j] += field[i * num_tuples + j]; 00062 } 00063 double tmp = work[18] / 24.0; 00064 for( int i = 0; i < num_tuples; i++ ) 00065 result[i] *= tmp; 00066 00067 return MB_SUCCESS; 00068 } 00069 00070 ErrorCode LinearTet::jacobianFcn( const double*, const double*, const int, const int, double* work, double* result ) 00071 { 00072 // jacobian is cached in work array 00073 assert( work ); 00074 std::copy( work, work + 9, result ); 00075 return MB_SUCCESS; 00076 } 00077 00078 ErrorCode LinearTet::reverseEvalFcn( EvalFcn eval, JacobianFcn jacob, InsideFcn ins, const double* posn, 00079 const double* verts, const int nverts, const int ndim, const double iter_tol, 00080 const double inside_tol, double* work, double* params, int* is_inside ) 00081 { 00082 assert( posn && verts ); 00083 return evaluate_reverse( eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, params, 00084 is_inside ); 00085 } 00086 00087 int LinearTet::insideFcn( const double* params, const int, const double tol ) 00088 { 00089 return ( params[0] >= -1.0 - tol && params[1] >= -1.0 - tol && params[2] >= -1.0 - tol && 00090 params[0] + params[1] + params[2] <= 1.0 + tol ); 00091 } 00092 00093 ErrorCode LinearTet::evaluate_reverse( EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f, const double* posn, 00094 const double* verts, const int nverts, const int ndim, const double iter_tol, 00095 const double inside_tol, double* work, double* params, int* inside ) 00096 { 00097 // TODO: should differentiate between epsilons used for 00098 // Newton Raphson iteration, and epsilons used for curved boundary geometry errors 00099 // right now, fix the tolerance used for NR 00100 const double error_tol_sqr = iter_tol * iter_tol; 00101 CartVect* cvparams = reinterpret_cast< CartVect* >( params ); 00102 const CartVect* cvposn = reinterpret_cast< const CartVect* >( posn ); 00103 00104 // find best initial guess to improve convergence 00105 CartVect tmp_params[] = { CartVect( -1, -1, -1 ), CartVect( 1, -1, -1 ), CartVect( -1, 1, -1 ), 00106 CartVect( -1, -1, 1 ) }; 00107 double resl = std::numeric_limits< double >::max(); 00108 CartVect new_pos, tmp_pos; 00109 ErrorCode rval; 00110 for( unsigned int i = 0; i < 4; i++ ) 00111 { 00112 rval = ( *eval )( tmp_params[i].array(), verts, ndim, ndim, work, tmp_pos.array() ); 00113 if( MB_SUCCESS != rval ) return rval; 00114 double tmp_resl = ( tmp_pos - *cvposn ).length_squared(); 00115 if( tmp_resl < resl ) 00116 { 00117 *cvparams = tmp_params[i]; 00118 new_pos = tmp_pos; 00119 resl = tmp_resl; 00120 } 00121 } 00122 00123 // residual is diff between old and new pos; need to minimize that 00124 CartVect res = new_pos - *cvposn; 00125 Matrix3 J; 00126 rval = ( *jacob )( cvparams->array(), verts, nverts, ndim, work, J.array() ); 00127 #ifndef NDEBUG 00128 double det = J.determinant(); 00129 assert( det > std::numeric_limits< double >::epsilon() ); 00130 #endif 00131 Matrix3 Ji = J.inverse(); 00132 00133 int iters = 0; 00134 // while |res| larger than tol 00135 int dum, *tmp_inside = ( inside ? inside : &dum ); 00136 while( res % res > error_tol_sqr ) 00137 { 00138 if( ++iters > 25 ) 00139 { 00140 // if we haven't converged but we're outside, that's defined as success 00141 *tmp_inside = ( *inside_f )( params, ndim, inside_tol ); 00142 if( !( *tmp_inside ) ) 00143 return MB_SUCCESS; 00144 else 00145 return MB_INDEX_OUT_OF_RANGE; 00146 } 00147 00148 // new params tries to eliminate residual 00149 *cvparams -= Ji * res; 00150 00151 // get the new forward-evaluated position, and its difference from the target pt 00152 rval = ( *eval )( params, verts, ndim, ndim, work, new_pos.array() ); 00153 if( MB_SUCCESS != rval ) return rval; 00154 res = new_pos - *cvposn; 00155 } 00156 00157 if( inside ) *inside = ( *inside_f )( params, ndim, inside_tol ); 00158 00159 return MB_SUCCESS; 00160 } // Map::evaluate_reverse() 00161 00162 ErrorCode LinearTet::normalFcn( const int ientDim, const int facet, const int nverts, const double* verts, 00163 double normal[3] ) 00164 { 00165 // assert(facet < 4 && ientDim == 2 && nverts == 4); 00166 if( nverts != 4 ) MB_SET_ERR( MB_FAILURE, "Incorrect vertex count for passed tet :: expected value = 4 " ); 00167 if( ientDim != 2 ) MB_SET_ERR( MB_FAILURE, "Requesting normal for unsupported dimension :: expected value = 2 " ); 00168 if( facet > 4 || facet < 0 ) MB_SET_ERR( MB_FAILURE, "Incorrect local face id :: expected value = one of 0-3" ); 00169 00170 int id0 = CN::mConnectivityMap[MBTET][ientDim - 1].conn[facet][0]; 00171 int id1 = CN::mConnectivityMap[MBTET][ientDim - 1].conn[facet][1]; 00172 int id2 = CN::mConnectivityMap[MBTET][ientDim - 1].conn[facet][2]; 00173 00174 double x0[3], x1[3]; 00175 00176 for( int i = 0; i < 3; i++ ) 00177 { 00178 x0[i] = verts[3 * id1 + i] - verts[3 * id0 + i]; 00179 x1[i] = verts[3 * id2 + i] - verts[3 * id0 + i]; 00180 } 00181 00182 double a = x0[1] * x1[2] - x1[1] * x0[2]; 00183 double b = x1[0] * x0[2] - x0[0] * x1[2]; 00184 double c = x0[0] * x1[1] - x1[0] * x0[1]; 00185 double nrm = sqrt( a * a + b * b + c * c ); 00186 00187 if( nrm > std::numeric_limits< double >::epsilon() ) 00188 { 00189 normal[0] = a / nrm; 00190 normal[1] = b / nrm; 00191 normal[2] = c / nrm; 00192 } 00193 return MB_SUCCESS; 00194 } 00195 00196 } // namespace moab