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