MOAB: Mesh Oriented datABase  (version 5.3.1)
DGMSolver.cpp
Go to the documentation of this file.
00001 #include "moab/DiscreteGeometry/DGMSolver.hpp"
00002 #include "moab/ErrorHandler.hpp"
00003 #include <iostream>
00004 #include <cassert>
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines