MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #include "moab/DiscreteGeometry/DGMSolver.hpp" 00002 #include "moab/ErrorHandler.hpp" 00003 #include <iostream> 00004 #include <assert.h> 00005 #include <vector> 00006 #include <limits> 00007 #include <cmath> 00008 #include <algorithm> 00009 00010 namespace moab 00011 { 00012 00013 /* This class implements the lowest level solvers required by polynomial fitting for high-order 00014 * reconstruction. An underlying assumption of the matrices passed are that they are given in a 00015 * column-major form. So, the first mrows values of V is the first column, and so on. This 00016 * assumption is made because most of the operations are column-based in the current scenario. 00017 * */ 00018 00019 unsigned int DGMSolver::nchoosek( unsigned int n, unsigned int k ) 00020 { 00021 if( k > n ) { return 0; } 00022 unsigned long long ans = 1; 00023 if( k > ( n >> 1 ) ) { k = n - k; } 00024 for( unsigned int i = 1; i <= k; ++i ) 00025 { 00026 ans *= n--; 00027 ans /= i; 00028 if( ans > std::numeric_limits< unsigned int >::max() ) { return 0; } 00029 } 00030 return ans; 00031 } 00032 00033 unsigned int DGMSolver::compute_numcols_vander_multivar( unsigned int kvars, unsigned int degree ) 00034 { 00035 unsigned int mcols = 0; 00036 for( unsigned int i = 0; i <= degree; ++i ) 00037 { 00038 unsigned int temp = nchoosek( kvars - 1 + i, kvars - 1 ); 00039 if( !temp ) 00040 { 00041 std::cout << "overflow to compute nchoosek n= " << kvars - 1 + i << " k= " << kvars - 1 << std::endl; 00042 return 0; 00043 } 00044 mcols += temp; 00045 } 00046 return mcols; 00047 } 00048 00049 void DGMSolver::gen_multivar_monomial_basis( const int kvars, const double* vars, const int degree, 00050 std::vector< double >& basis ) 00051 { 00052 unsigned int len = compute_numcols_vander_multivar( kvars, degree ); 00053 basis.reserve( len - basis.capacity() + basis.size() ); 00054 size_t iend = basis.size(); 00055 #ifndef NDEBUG 00056 size_t istr = basis.size(); 00057 #endif 00058 basis.push_back( 1 ); 00059 ++iend; 00060 if( !degree ) { return; } 00061 std::vector< size_t > varspos( kvars ); 00062 // degree 1 00063 for( int ivar = 0; ivar < kvars; ++ivar ) 00064 { 00065 basis.push_back( vars[ivar] ); 00066 varspos[ivar] = iend++; 00067 } 00068 // degree 2 to degree 00069 for( int ideg = 2; ideg <= degree; ++ideg ) 00070 { 00071 size_t preend = iend; 00072 for( int ivar = 0; ivar < kvars; ++ivar ) 00073 { 00074 size_t varpreend = iend; 00075 for( size_t ilast = varspos[ivar]; ilast < preend; ++ilast ) 00076 { 00077 basis.push_back( vars[ivar] * basis[ilast] ); 00078 ++iend; 00079 } 00080 varspos[ivar] = varpreend; 00081 } 00082 } 00083 assert( len == iend - istr ); 00084 } 00085 00086 void DGMSolver::gen_vander_multivar( const int mrows, const int kvars, const double* us, const int degree, 00087 std::vector< double >& V ) 00088 { 00089 unsigned int ncols = compute_numcols_vander_multivar( kvars, degree ); 00090 V.reserve( mrows * ncols - V.capacity() + V.size() ); 00091 size_t istr = V.size(), icol = 0; 00092 // add ones, V is stored in an single array, elements placed in columnwise order 00093 for( int irow = 0; irow < mrows; ++irow ) 00094 { 00095 V.push_back( 1 ); 00096 } 00097 ++icol; 00098 if( !degree ) { return; } 00099 std::vector< size_t > varspos( kvars ); 00100 // degree 1 00101 for( int ivar = 0; ivar < kvars; ++ivar ) 00102 { 00103 for( int irow = 0; irow < mrows; ++irow ) 00104 { 00105 V.push_back( us[irow * kvars + ivar] ); // us stored in row-wise 00106 } 00107 varspos[ivar] = icol++; 00108 } 00109 // from 2 to degree 00110 for( int ideg = 2; ideg <= degree; ++ideg ) 00111 { 00112 size_t preendcol = icol; 00113 for( int ivar = 0; ivar < kvars; ++ivar ) 00114 { 00115 size_t varpreend = icol; 00116 for( size_t ilast = varspos[ivar]; ilast < preendcol; ++ilast ) 00117 { 00118 for( int irow = 0; irow < mrows; ++irow ) 00119 { 00120 V.push_back( us[irow * kvars + ivar] * V[istr + irow + ilast * mrows] ); 00121 } 00122 ++icol; 00123 } 00124 varspos[ivar] = varpreend; 00125 } 00126 } 00127 assert( icol == ncols ); 00128 } 00129 00130 void DGMSolver::rescale_matrix( int mrows, int ncols, double* V, double* ts ) 00131 { 00132 // This function rescales the input matrix using the norm of each column. 00133 double* v = new double[mrows]; 00134 for( int i = 0; i < ncols; i++ ) 00135 { 00136 for( int j = 0; j < mrows; j++ ) 00137 v[j] = V[mrows * i + j]; 00138 00139 // Compute norm of the column vector 00140 double w = vec_2norm( mrows, v ); 00141 00142 if( fabs( w ) == 0 ) 00143 ts[i] = 1; 00144 else 00145 { 00146 ts[i] = w; 00147 for( int j = 0; j < mrows; j++ ) 00148 V[mrows * i + j] = V[mrows * i + j] / ts[i]; 00149 } 00150 } 00151 delete[] v; 00152 } 00153 00154 void DGMSolver::compute_qtransposeB( int mrows, int ncols, const double* Q, int bncols, double* bs ) 00155 { 00156 for( int k = 0; k < ncols; k++ ) 00157 { 00158 for( int j = 0; j < bncols; j++ ) 00159 { 00160 double t2 = 0; 00161 for( int i = k; i < mrows; i++ ) 00162 t2 += Q[mrows * k + i] * bs[mrows * j + i]; 00163 t2 = t2 + t2; 00164 00165 for( int i = k; i < mrows; i++ ) 00166 bs[mrows * j + i] -= t2 * Q[mrows * k + i]; 00167 } 00168 } 00169 } 00170 00171 void DGMSolver::qr_polyfit_safeguarded( const int mrows, const int ncols, double* V, double* D, int* rank ) 00172 { 00173 double tol = 1e-8; 00174 *rank = ncols; 00175 double* v = new double[mrows]; 00176 00177 for( int k = 0; k < ncols; k++ ) 00178 { 00179 int nv = mrows - k; 00180 00181 for( int j = 0; j < nv; j++ ) 00182 v[j] = V[mrows * k + ( j + k )]; 00183 00184 double t2 = 0; 00185 00186 for( int j = 0; j < nv; j++ ) 00187 t2 = t2 + v[j] * v[j]; 00188 00189 double t = sqrt( t2 ); 00190 double vnrm = 0; 00191 00192 if( v[0] >= 0 ) 00193 { 00194 vnrm = sqrt( 2 * ( t2 + v[0] * t ) ); 00195 v[0] = v[0] + t; 00196 } 00197 else 00198 { 00199 vnrm = sqrt( 2 * ( t2 - v[0] * t ) ); 00200 v[0] = v[0] - t; 00201 } 00202 00203 if( vnrm > 0 ) 00204 { 00205 for( int j = 0; j < nv; j++ ) 00206 v[j] = v[j] / vnrm; 00207 } 00208 00209 for( int j = k; j < ncols; j++ ) 00210 { 00211 t2 = 0; 00212 for( int i = 0; i < nv; i++ ) 00213 t2 = t2 + v[i] * V[mrows * j + ( i + k )]; 00214 00215 t2 = t2 + t2; 00216 00217 for( int i = 0; i < nv; i++ ) 00218 V[mrows * j + ( i + k )] = V[mrows * j + ( i + k )] - t2 * v[i]; 00219 } 00220 00221 D[k] = V[mrows * k + k]; 00222 00223 for( int i = 0; i < nv; i++ ) 00224 V[mrows * k + ( i + k )] = v[i]; 00225 00226 if( ( fabs( D[k] ) ) < tol && ( *rank == ncols ) ) 00227 { 00228 *rank = k; 00229 break; 00230 } 00231 } 00232 00233 delete[] v; 00234 } 00235 00236 void DGMSolver::backsolve( int mrows, int ncols, double* R, int bncols, double* bs, double* ws ) 00237 { 00238 #if 0 00239 std::cout.precision(16); 00240 std::cout<<"Before backsolve "<<std::endl; 00241 std::cout<<" V = "<<std::endl; 00242 for (int k=0; k< ncols; k++){ 00243 for (int j=0; j<mrows; ++j){ 00244 std::cout<<R[mrows*k+j]<<std::endl; 00245 } 00246 std::cout<<std::endl; 00247 } 00248 std::cout<<std::endl; 00249 00250 //std::cout<<"#pnts = "<<mrows<<std::endl; 00251 std::cout<<"bs = "<<std::endl; 00252 std::cout<<std::endl; 00253 for (int k=0; k< bncols; k++){ 00254 for (int j=0; j<mrows; ++j){ 00255 std::cout<<" "<<bs[mrows*k+j]<<std::endl; 00256 } 00257 } 00258 std::cout<<std::endl; 00259 #endif 00260 00261 for( int k = 0; k < bncols; k++ ) 00262 { 00263 for( int j = ncols - 1; j >= 0; j-- ) 00264 { 00265 for( int i = j + 1; i < ncols; ++i ) 00266 bs[mrows * k + j] = bs[mrows * k + j] - R[mrows * i + j] * bs[mrows * k + i]; 00267 00268 assert( R[mrows * j + j] != 0 ); 00269 bs[mrows * k + j] = bs[mrows * k + j] / R[mrows * j + j]; 00270 } 00271 } 00272 00273 for( int k = 0; k < bncols; k++ ) 00274 { 00275 for( int j = 0; j < ncols; ++j ) 00276 bs[mrows * k + j] = bs[mrows * k + j] / ws[j]; 00277 } 00278 } 00279 00280 void DGMSolver::backsolve_polyfit_safeguarded( int dim, int degree, const bool interp, int mrows, int ncols, double* R, 00281 int bncols, double* bs, const double* ws, int* degree_out ) 00282 { 00283 #if 0 00284 std::cout.precision(12); 00285 std::cout<<"Before backsolve "<<std::endl; 00286 std::cout<<" V = "<<std::endl; 00287 for (int k=0; k< ncols; k++){ 00288 for (int j=0; j<mrows; ++j){ 00289 std::cout<<R[mrows*k+j]<<std::endl; 00290 } 00291 std::cout<<std::endl; 00292 } 00293 std::cout<<std::endl; 00294 00295 //std::cout<<"#pnts = "<<mrows<<std::endl; 00296 std::cout<<"bs = "<<std::endl; 00297 std::cout<<std::endl; 00298 for (int k=0; k< bncols; k++){ 00299 for (int j=0; j<mrows; ++j){ 00300 std::cout<<" "<<bs[mrows*k+j]<<std::endl; 00301 } 00302 } 00303 std::cout<<std::endl; 00304 00305 //std::cout<<" ] "<<std::endl; 00306 00307 std::cout<<"Input ws = [ "; 00308 for (int k=0; k< ncols; k++){ 00309 std::cout<<ws[k]<<", "; 00310 } 00311 std::cout<<" ] "<<std::endl; 00312 00313 std::cout << "R: " << R << "size: [" << mrows << "," << ncols << "]" << std::endl; 00314 std::cout << "bs: " << bs << "size: [" << mrows << "," << bncols << "]" << std::endl; 00315 std::cout << "ws: " << ws << "size: [" << ncols << "," << 1 << "]" << std::endl; 00316 std::cout << "degree_out: " << degree_out << std::endl; 00317 #endif 00318 00319 int deg, numcols = 0; 00320 00321 for( int k = 0; k < bncols; k++ ) 00322 { 00323 deg = degree; 00324 /* I think we should consider interp = true/false -Xinglin*/ 00325 if( dim == 1 ) 00326 numcols = deg + 1 - interp; 00327 else if( dim == 2 ) 00328 numcols = ( deg + 2 ) * ( deg + 1 ) / 2 - interp; 00329 00330 assert( numcols <= ncols ); 00331 00332 std::vector< double > bs_bak( numcols ); 00333 00334 if( deg >= 2 ) 00335 { 00336 for( int i = 0; i < numcols; i++ ) 00337 { 00338 assert( mrows * k + i < mrows * bncols ); 00339 bs_bak.at( i ) = bs[mrows * k + i]; 00340 } 00341 } 00342 00343 while( deg >= 1 ) 00344 { 00345 int cend = numcols - 1; 00346 bool downgrade = false; 00347 00348 // The reconstruction can be applied only on edges (2-d) or faces (3-d) 00349 assert( cend >= 0 ); 00350 assert( dim > 0 && dim < 3 ); 00351 00352 for( int d = deg; d >= 0; d-- ) 00353 { 00354 int cstart = 0; 00355 if( dim == 1 ) { cstart = d; } 00356 else if( dim == 2 ) 00357 { 00358 cstart = ( ( d + 1 ) * d ) / 2; 00359 // cstart = ((d*(d+1))>>1)-interp; 00360 } 00361 00362 // Solve for bs 00363 for( int j = cend; j >= cstart; j-- ) 00364 { 00365 assert( mrows * k + j < mrows * bncols ); 00366 for( int i = j + 1; i < numcols; ++i ) 00367 { 00368 assert( mrows * k + i < mrows * bncols ); 00369 assert( mrows * i + j < mrows * ncols ); // check R 00370 bs[mrows * k + j] = bs[mrows * k + j] - R[mrows * i + j] * bs[mrows * k + i]; 00371 } 00372 assert( mrows * j + j < mrows * ncols ); // check R 00373 bs[mrows * k + j] = bs[mrows * k + j] / R[mrows * j + j]; 00374 } 00375 00376 // Checking for change in the coefficient 00377 if( d >= 2 && d < deg ) 00378 { 00379 double tol; 00380 00381 if( dim == 1 ) 00382 { 00383 tol = 1e-06; 00384 assert( mrows * cstart + cstart < mrows * ncols ); // check R 00385 double tb = bs_bak.at( cstart ) / R[mrows * cstart + cstart]; 00386 assert( mrows * k + cstart < mrows * bncols ); 00387 if( fabs( bs[mrows * k + cstart] - tb ) > ( 1 + tol ) * fabs( tb ) ) 00388 { 00389 downgrade = true; 00390 break; 00391 } 00392 } 00393 00394 else if( dim == 2 ) 00395 { 00396 tol = 0.05; 00397 00398 std::vector< double > tb( cend - cstart + 1 ); 00399 for( int j = 0; j <= ( cend - cstart ); j++ ) 00400 { 00401 tb.at( j ) = bs_bak.at( cstart + j ); 00402 } 00403 00404 for( int j = cend; j >= cstart; j-- ) 00405 { 00406 int jind = j - cstart; 00407 00408 for( int i = j + 1; i <= cend; ++i ) 00409 { 00410 assert( mrows * i + j < mrows * ncols ); // check R 00411 tb.at( jind ) = tb.at( jind ) - R[mrows * i + j] * tb.at( i - cstart ); 00412 } 00413 assert( mrows * j + j < mrows * ncols ); // check R 00414 tb.at( jind ) = tb.at( jind ) / R[mrows * j + j]; 00415 assert( mrows * k + j < mrows * bncols ); 00416 double err = fabs( bs[mrows * k + j] - tb.at( jind ) ); 00417 00418 if( ( err > tol ) && ( err >= ( 1 + tol ) * fabs( tb.at( jind ) ) ) ) 00419 { 00420 downgrade = true; 00421 break; 00422 } 00423 } 00424 00425 if( downgrade ) break; 00426 } 00427 } 00428 00429 cend = cstart - 1; 00430 } 00431 00432 if( !downgrade ) 00433 break; 00434 else 00435 { 00436 deg = deg - 1; 00437 if( dim == 1 ) 00438 numcols = deg + 1; 00439 else if( dim == 2 ) 00440 numcols = ( deg + 2 ) * ( deg + 1 ) / 2; 00441 00442 for( int i = 0; i < numcols; i++ ) 00443 { 00444 assert( mrows * k + i < mrows * bncols ); 00445 bs[mrows * k + i] = bs_bak.at( i ); 00446 } 00447 } 00448 } 00449 assert( k < bncols ); 00450 degree_out[k] = deg; 00451 00452 for( int i = 0; i < numcols; i++ ) 00453 { 00454 // assert(mrows*k+i < mrows*bncols); 00455 // assert(i < ncols); 00456 bs[mrows * k + i] = bs[mrows * k + i] / ws[i]; 00457 } 00458 00459 for( int i = numcols; i < mrows; i++ ) 00460 { 00461 // assert(mrows*k+i < mrows*bncols); 00462 bs[mrows * k + i] = 0; 00463 } 00464 } 00465 } 00466 00467 void DGMSolver::vec_dotprod( const int len, const double* a, const double* b, double* c ) 00468 { 00469 if( !a || !b || !c ) { MB_SET_ERR_RET( "NULL Pointer" ); } 00470 for( int i = 0; i < len; ++i ) 00471 { 00472 c[i] = a[i] * b[i]; 00473 } 00474 } 00475 00476 void DGMSolver::vec_scalarprod( const int len, const double* a, const double c, double* b ) 00477 { 00478 if( !a || !b ) { MB_SET_ERR_RET( "NULL Pointer" ); } 00479 for( int i = 0; i < len; ++i ) 00480 { 00481 b[i] = c * a[i]; 00482 } 00483 } 00484 00485 void DGMSolver::vec_crossprod( const double a[3], const double b[3], double ( &c )[3] ) 00486 { 00487 c[0] = a[1] * b[2] - a[2] * b[1]; 00488 c[1] = a[2] * b[0] - a[0] * b[2]; 00489 c[2] = a[0] * b[1] - a[1] * b[0]; 00490 } 00491 00492 double DGMSolver::vec_innerprod( const int len, const double* a, const double* b ) 00493 { 00494 if( !a || !b ) { MB_SET_ERR_RET_VAL( "NULL Pointer", 0.0 ); } 00495 double ans = 0; 00496 for( int i = 0; i < len; ++i ) 00497 { 00498 ans += a[i] * b[i]; 00499 } 00500 return ans; 00501 } 00502 00503 double DGMSolver::vec_2norm( const int len, const double* a ) 00504 { 00505 if( !a ) { MB_SET_ERR_RET_VAL( "NULL Pointer", 0.0 ); } 00506 double w = 0, s = 0; 00507 for( int k = 0; k < len; k++ ) 00508 w = std::max( w, fabs( a[k] ) ); 00509 00510 if( w == 0 ) { return 0; } 00511 else 00512 { 00513 for( int k = 0; k < len; k++ ) 00514 { 00515 s += ( a[k] / w ) * ( a[k] / w ); 00516 } 00517 s = w * sqrt( s ); 00518 } 00519 return s; 00520 } 00521 00522 double DGMSolver::vec_normalize( const int len, const double* a, double* b ) 00523 { 00524 if( !a || !b ) { MB_SET_ERR_RET_VAL( "NULL Pointer", 0.0 ); } 00525 double nrm = 0, mx = 0; 00526 for( int i = 0; i < len; ++i ) 00527 { 00528 mx = std::max( fabs( a[i] ), mx ); 00529 } 00530 if( mx == 0 ) 00531 { 00532 for( int i = 0; i < len; ++i ) 00533 { 00534 b[i] = 0; 00535 } 00536 return 0; 00537 } 00538 for( int i = 0; i < len; ++i ) 00539 { 00540 nrm += ( a[i] / mx ) * ( a[i] / mx ); 00541 } 00542 nrm = mx * sqrt( nrm ); 00543 if( nrm == 0 ) { return nrm; } 00544 for( int i = 0; i < len; ++i ) 00545 { 00546 b[i] = a[i] / nrm; 00547 } 00548 return nrm; 00549 } 00550 00551 double DGMSolver::vec_distance( const int len, const double* a, const double* b ) 00552 { 00553 double res = 0; 00554 for( int i = 0; i < len; ++i ) 00555 { 00556 res += ( a[i] - b[i] ) * ( a[i] - b[i] ); 00557 } 00558 return sqrt( res ); 00559 } 00560 00561 void DGMSolver::vec_projoff( const int len, const double* a, const double* b, double* c ) 00562 { 00563 if( !a || !b || !c ) { MB_SET_ERR_RET( "NULL Pointer" ); } 00564 // c = a-<a,b>b/<b,b>; 00565 double bnrm = vec_2norm( len, b ); 00566 if( bnrm == 0 ) 00567 { 00568 for( int i = 0; i < len; ++i ) 00569 { 00570 c[i] = a[i]; 00571 } 00572 return; 00573 } 00574 double innerp = vec_innerprod( len, a, b ) / bnrm; 00575 00576 if( innerp == 0 ) 00577 { 00578 if( c != a ) 00579 { 00580 for( int i = 0; i < len; ++i ) 00581 { 00582 c[i] = a[i]; 00583 } 00584 } 00585 return; 00586 } 00587 00588 for( int i = 0; i < len; ++i ) 00589 { 00590 c[i] = a[i] - innerp * b[i] / bnrm; 00591 } 00592 } 00593 00594 void DGMSolver::vec_linear_operation( const int len, const double mu, const double* a, const double psi, 00595 const double* b, double* c ) 00596 { 00597 if( !a || !b || !c ) { MB_SET_ERR_RET( "NULL Pointer" ); } 00598 for( int i = 0; i < len; ++i ) 00599 { 00600 c[i] = mu * a[i] + psi * b[i]; 00601 } 00602 } 00603 00604 void DGMSolver::get_tri_natural_coords( const int dim, const double* cornercoords, const int npts, 00605 const double* currcoords, double* naturalcoords ) 00606 { 00607 assert( dim == 2 || dim == 3 ); 00608 double a = 0, b = 0, d = 0, tol = 1e-12; 00609 for( int i = 0; i < dim; ++i ) 00610 { 00611 a += ( cornercoords[dim + i] - cornercoords[i] ) * ( cornercoords[dim + i] - cornercoords[i] ); 00612 b += ( cornercoords[dim + i] - cornercoords[i] ) * ( cornercoords[2 * dim + i] - cornercoords[i] ); 00613 d += ( cornercoords[2 * dim + i] - cornercoords[i] ) * ( cornercoords[2 * dim + i] - cornercoords[i] ); 00614 } 00615 double det = a * d - b * b; 00616 assert( det > 0 ); 00617 for( int ipt = 0; ipt < npts; ++ipt ) 00618 { 00619 double e = 0, f = 0; 00620 for( int i = 0; i < dim; ++i ) 00621 { 00622 e += ( cornercoords[dim + i] - cornercoords[i] ) * ( currcoords[ipt * dim + i] - cornercoords[i] ); 00623 f += ( cornercoords[2 * dim + i] - cornercoords[i] ) * ( currcoords[ipt * dim + i] - cornercoords[i] ); 00624 } 00625 naturalcoords[ipt * 3 + 1] = ( e * d - b * f ) / det; 00626 naturalcoords[ipt * 3 + 2] = ( a * f - b * e ) / det; 00627 naturalcoords[ipt * 3] = 1 - naturalcoords[ipt * 3 + 1] - naturalcoords[ipt * 3 + 2]; 00628 if( naturalcoords[ipt * 3] < -tol || naturalcoords[ipt * 3 + 1] < -tol || naturalcoords[ipt * 3 + 2] < -tol ) 00629 { 00630 std::cout << "Corners: \n"; 00631 std::cout << cornercoords[0] << "\t" << cornercoords[1] << "\t" << cornercoords[3] << std::endl; 00632 std::cout << cornercoords[3] << "\t" << cornercoords[4] << "\t" << cornercoords[5] << std::endl; 00633 std::cout << cornercoords[6] << "\t" << cornercoords[7] << "\t" << cornercoords[8] << std::endl; 00634 std::cout << "Candidate: \n"; 00635 std::cout << currcoords[ipt * dim] << "\t" << currcoords[ipt * dim + 1] << "\t" << currcoords[ipt * dim + 2] 00636 << std::endl; 00637 exit( 0 ); 00638 } 00639 assert( fabs( naturalcoords[ipt * 3] + naturalcoords[ipt * 3 + 1] + naturalcoords[ipt * 3 + 2] - 1 ) < tol ); 00640 for( int i = 0; i < dim; ++i ) 00641 { 00642 assert( fabs( naturalcoords[ipt * 3] * cornercoords[i] + 00643 naturalcoords[ipt * 3 + 1] * cornercoords[dim + i] + 00644 naturalcoords[ipt * 3 + 2] * cornercoords[2 * dim + i] - currcoords[ipt * dim + i] ) < tol ); 00645 } 00646 } 00647 } 00648 } // namespace moab