MOAB: Mesh Oriented datABase  (version 5.3.1)
elem_eval_time.cpp
Go to the documentation of this file.
00001 #include "moab/Core.hpp"
00002 #include "moab/Range.hpp"
00003 #include "moab/ProgOptions.hpp"
00004 #include "moab/CartVect.hpp"
00005 #include "moab/ElemEvaluator.hpp"
00006 #include "moab/CN.hpp"
00007 #include "ElemUtil.hpp"
00008 #include <iostream>
00009 #include <ctime>
00010 #include <cstdlib>
00011 #include <cassert>
00012 #include <sstream>
00013 #include <string>
00014 
00015 // Different platforms follow different conventions for usage
00016 #if !defined( _MSC_VER ) && !defined( __MINGW32__ )
00017 #include <sys/resource.h>
00018 #endif
00019 #ifdef SOLARIS
00020 extern "C" int getrusage( int, struct rusage* );
00021 #ifndef RUSAGE_SELF
00022 #include </usr/ucbinclude/sys/rusage.h>
00023 #endif
00024 #endif
00025 
00026 using namespace moab;
00027 #define CHK( r, s )                                   \
00028     do                                                \
00029     {                                                 \
00030         if( MB_SUCCESS != ( r ) )                     \
00031         {                                             \
00032             fail( ( r ), ( s ), __FILE__, __LINE__ ); \
00033             return ( r );                             \
00034         }                                             \
00035     } while( false )
00036 
00037 const int ELEMEVAL = 0, ELEMUTIL = 1;
00038 
00039 static void fail( ErrorCode error_code, const char* str, const char* file_name, int line_number )
00040 {
00041     std::cerr << str << ", line " << line_number << " of file " << file_name << ", error code " << error_code
00042               << std::endl;
00043 }
00044 
00045 double mytime();
00046 
00047 ErrorCode get_ents( Interface& mbi, std::string& filename, int& dim, Range& elems, EntityType& tp, int& nv )
00048 {
00049     ErrorCode rval = mbi.load_file( filename.c_str(), 0 );
00050     while( elems.empty() && dim >= 1 )
00051     {
00052         rval = mbi.get_entities_by_dimension( 0, dim, elems );
00053         CHK( rval, "get_entities_by_dimension" );
00054         if( elems.empty() ) dim--;
00055     }
00056     if( elems.empty() ) { CHK( MB_FAILURE, "No elements in file" ); }
00057 
00058     // check to see they're all the same type & #vertices
00059     tp             = mbi.type_from_handle( *elems.begin() );
00060     EntityType tp2 = mbi.type_from_handle( *elems.rbegin() );
00061     if( tp != tp2 ) CHK( MB_FAILURE, "Elements must have same type" );
00062 
00063     int nv2;
00064     const EntityHandle* c1;
00065     rval = mbi.get_connectivity( *elems.begin(), c1, nv );
00066     CHK( rval, "get_connectivity" );
00067     rval = mbi.get_connectivity( *elems.rbegin(), c1, nv2 );
00068     CHK( rval, "get_connectivity" );
00069     if( nv != nv2 ) { CHK( MB_FAILURE, "Elements must have same #vertices" ); }
00070 
00071     return MB_SUCCESS;
00072 }
00073 
00074 void parse_options( ProgOptions& opts, int& dim, std::string& filename )
00075 {
00076     opts.addOpt< int >( std::string( "dim" ),
00077                         std::string( "Evaluate 1d/2d/3d elements (default: maximal dimension in mesh)" ), &dim,
00078                         ProgOptions::int_flag );
00079     opts.addOpt< std::string >( std::string( "filename,f" ), std::string( "Filename containing mesh" ), &filename );
00080 }
00081 
00082 ErrorCode get_elem_map( EntityType tp, std::vector< CartVect >& vcoords, int nconn, Element::Map*& elemmap )
00083 {
00084     switch( tp )
00085     {
00086         case MBHEX:
00087             if( nconn == 8 )
00088             {
00089                 elemmap = new Element::LinearHex( vcoords );
00090                 break;
00091             }
00092             else if( nconn == 27 )
00093             {
00094                 elemmap = new Element::QuadraticHex( vcoords );
00095                 break;
00096             }
00097             else
00098                 return MB_FAILURE;
00099         case MBTET:
00100             if( nconn == 4 )
00101             {
00102                 elemmap = new Element::LinearTet( vcoords );
00103                 break;
00104             }
00105             else
00106                 return MB_FAILURE;
00107         case MBQUAD:
00108             if( nconn == 4 )
00109             {
00110                 elemmap = new Element::LinearQuad( vcoords );
00111                 break;
00112             }
00113             else
00114                 return MB_FAILURE;
00115         default:
00116             return MB_FAILURE;
00117     }
00118 
00119     return MB_SUCCESS;
00120 }
00121 
00122 ErrorCode time_forward_eval( Interface* mbi, int method, Range& elems, std::vector< CartVect >& params,
00123                              std::vector< CartVect >& coords, double& evtime )
00124 {
00125     evtime = mytime();
00126     ErrorCode rval;
00127     Range::iterator rit;
00128     unsigned int i;
00129     if( ELEMEVAL == method )
00130     {
00131         // construct ElemEvaluator
00132         EvalSet eset;
00133         ElemEvaluator eeval( mbi );
00134         eeval.set_eval_set( *elems.begin() );
00135         eeval.set_tag_handle( 0, 0 );  // indicates coordinates as the field to evaluate
00136 
00137         for( rit = elems.begin(), i = 0; rit != elems.end(); ++rit, i++ )
00138         {
00139             eeval.set_ent_handle( *rit );
00140             rval = eeval.eval( params[i].array(), coords[i].array(), 3 );
00141 #ifndef NDEBUG
00142             if( MB_SUCCESS != rval ) return rval;
00143 #endif
00144         }
00145     }
00146     else if( ELEMUTIL == method )
00147     {
00148         std::vector< CartVect > vcoords( CN::MAX_NODES_PER_ELEMENT );
00149         const EntityHandle* connect;
00150         int nconn;
00151         Element::Map* elemmap = NULL;
00152         for( rit = elems.begin(), i = 0; rit != elems.end(); ++rit, i++ )
00153         {
00154             rval = mbi->get_connectivity( *rit, connect, nconn );
00155             CHK( rval, "get_connectivity" );
00156             rval = mbi->get_coords( connect, nconn, vcoords[0].array() );
00157             CHK( rval, "get_coords" );
00158             rval = get_elem_map( mbi->type_from_handle( *rit ), vcoords, nconn, elemmap );
00159             CHK( rval, "get_elem_map" );
00160             coords[i] = elemmap->evaluate( params[i] );
00161         }
00162     }
00163 
00164     evtime = mytime() - evtime;
00165     return MB_SUCCESS;
00166 }
00167 
00168 ErrorCode time_reverse_eval( Interface* mbi, int method, Range& elems, std::vector< CartVect >& coords,
00169                              std::vector< CartVect >& params, double& retime )
00170 {
00171     retime = mytime();
00172     ErrorCode rval;
00173     Range::iterator rit;
00174     unsigned int i;
00175     if( ELEMEVAL == method )
00176     {
00177         EvalSet eset;
00178         ElemEvaluator eeval( mbi );
00179         eeval.set_eval_set( *elems.begin() );
00180         eeval.set_tag_handle( 0, 0 );  // indicates coordinates as the field to evaluate
00181         int ins;
00182         for( rit = elems.begin(), i = 0; rit != elems.end(); ++rit, i++ )
00183         {
00184             eeval.set_ent_handle( *rit );
00185             rval = eeval.reverse_eval( coords[i].array(), 1.0e-10, 1.0e-6, params[i].array(), &ins );
00186             assert( ins );
00187 #ifndef NDEBUG
00188             if( MB_SUCCESS != rval ) return rval;
00189 #endif
00190         }
00191     }
00192     else if( ELEMUTIL == method )
00193     {
00194         std::vector< CartVect > vcoords( CN::MAX_NODES_PER_ELEMENT );
00195         const EntityHandle* connect;
00196         int nconn;
00197         Element::Map* elemmap = NULL;
00198         for( rit = elems.begin(), i = 0; rit != elems.end(); ++rit, i++ )
00199         {
00200             rval = mbi->get_connectivity( *rit, connect, nconn );
00201             CHK( rval, "get_connectivity" );
00202             rval = mbi->get_coords( connect, nconn, vcoords[0].array() );
00203             CHK( rval, "get_coords" );
00204             rval = get_elem_map( mbi->type_from_handle( *rit ), vcoords, nconn, elemmap );
00205             CHK( rval, "get_elem_map" );
00206             coords[i] = elemmap->ievaluate( coords[i], 1.0e-6 );
00207         }
00208     }
00209     retime = mytime() - retime;
00210     return MB_SUCCESS;
00211 }
00212 
00213 ErrorCode time_jacobian( Interface* mbi, int method, Range& elems, std::vector< CartVect >& params, double& jactime )
00214 {
00215     jactime = mytime();
00216     ErrorCode rval;
00217     Range::iterator rit;
00218     unsigned int i;
00219     Matrix3 jac;
00220     if( ELEMEVAL == method )
00221     {
00222         EvalSet eset;
00223         ElemEvaluator eeval( mbi );
00224         eeval.set_eval_set( *elems.begin() );
00225         eeval.set_tag_handle( 0, 0 );  // indicates coordinates as the field to evaluate
00226         for( rit = elems.begin(), i = 0; rit != elems.end(); ++rit, i++ )
00227         {
00228             eeval.set_ent_handle( *rit );
00229             rval = eeval.jacobian( params[i].array(), jac.array() );
00230 #ifndef NDEBUG
00231             if( MB_SUCCESS != rval ) return rval;
00232 #endif
00233         }
00234     }
00235     else if( ELEMUTIL == method )
00236     {
00237         std::vector< CartVect > vcoords( CN::MAX_NODES_PER_ELEMENT );
00238         const EntityHandle* connect;
00239         int nconn;
00240         Element::Map* elemmap = NULL;
00241         for( rit = elems.begin(), i = 0; rit != elems.end(); ++rit, i++ )
00242         {
00243             rval = mbi->get_connectivity( *rit, connect, nconn );
00244             CHK( rval, "get_connectivity" );
00245             rval = mbi->get_coords( connect, nconn, vcoords[0].array() );
00246             CHK( rval, "get_coords" );
00247             rval = get_elem_map( mbi->type_from_handle( *rit ), vcoords, nconn, elemmap );
00248             CHK( rval, "get_elem_map" );
00249             jac = elemmap->jacobian( params[i] );
00250         }
00251     }
00252     jactime = mytime() - jactime;
00253     return MB_SUCCESS;
00254 }
00255 
00256 ErrorCode time_integrate( Interface* mbi, int method, Tag tag, Range& elems, double& inttime )
00257 {
00258     inttime = mytime();
00259     ErrorCode rval;
00260     double integral = 0.0;
00261     if( ELEMEVAL == method )
00262     {
00263         EvalSet eset;
00264         ElemEvaluator eeval( mbi );
00265         eeval.set_eval_set( *elems.begin() );
00266         eeval.set_tag_handle( 0, 0 );  // indicates coordinates as the field to evaluate
00267         rval = eeval.set_tag_handle( tag, 0 );
00268         CHK( rval, "set_tag_handle" );
00269         for( Range::iterator rit = elems.begin(); rit != elems.end(); ++rit )
00270         {
00271             eeval.set_ent_handle( *rit );
00272             rval = eeval.integrate( &integral );
00273 #ifndef NDEBUG
00274             if( MB_SUCCESS != rval ) return rval;
00275 #endif
00276         }
00277     }
00278     else if( ELEMUTIL == method )
00279     {
00280         std::vector< CartVect > vcoords( CN::MAX_NODES_PER_ELEMENT );
00281         std::vector< double > tagval( CN::MAX_NODES_PER_ELEMENT );
00282         const EntityHandle* connect;
00283         int nconn;
00284         Element::Map* elemmap = NULL;
00285         Range::iterator rit;
00286         unsigned int i;
00287         for( rit = elems.begin(), i = 0; rit != elems.end(); ++rit, i++ )
00288         {
00289             rval = mbi->get_connectivity( *rit, connect, nconn );
00290             CHK( rval, "get_connectivity" );
00291             rval = mbi->get_coords( connect, nconn, vcoords[0].array() );
00292             CHK( rval, "get_coords" );
00293             rval = get_elem_map( mbi->type_from_handle( *rit ), vcoords, nconn, elemmap );
00294             CHK( rval, "get_elem_map" );
00295             rval = mbi->tag_get_data( tag, connect, nconn, &tagval[0] );
00296             CHK( rval, "tag_get_data" );
00297             integral = elemmap->integrate_scalar_field( &tagval[0] );
00298         }
00299     }
00300     inttime = mytime() - inttime;
00301     std::cout << "integral = " << integral << std::endl;
00302     return MB_SUCCESS;
00303 }
00304 
00305 ErrorCode put_random_field( Interface& mbi, Tag& tag, Range& elems )
00306 {
00307     Range verts;
00308     ErrorCode rval = mbi.get_adjacencies( elems, 0, false, verts, Interface::UNION );
00309     CHK( rval, "get_adjacencies" );
00310     rval = mbi.tag_get_handle( "mytag", 1, MB_TYPE_DOUBLE, tag, MB_TAG_CREAT | MB_TAG_DENSE );
00311     CHK( rval, "tag_get_handle" );
00312     std::vector< double > tag_vals( verts.size() );
00313     for( unsigned int i = 0; i < verts.size(); i++ )
00314         tag_vals[i] = ( (double)rand() ) / RAND_MAX;
00315     rval = mbi.tag_set_data( tag, verts, &tag_vals[0] );
00316     return rval;
00317 }
00318 
00319 ErrorCode elem_evals( Interface* mbi, int method, Range& elems, Tag tag, std::vector< CartVect >& params,
00320                       std::vector< CartVect >& coords, double& evtime, double& retime, double& jactime,
00321                       double& inttime )
00322 {
00323     evtime = 0, retime = 0, jactime = 0, inttime = 0;  // initializations to avoid compiler warning
00324 
00325     // time/test forward evaluation, putting results into vector
00326     ErrorCode rval = time_forward_eval( mbi, method, elems, params, coords, evtime );
00327     CHK( rval, "time_forward_eval" );
00328 
00329     // time/test reverse evaluation, putting results into vector
00330     rval = time_reverse_eval( mbi, method, elems, coords, params, retime );
00331     CHK( rval, "time_reverse_eval" );
00332 
00333     // time/test Jacobian evaluation
00334     rval = time_jacobian( mbi, method, elems, params, jactime );
00335     CHK( rval, "time_jacobian" );
00336 
00337     // time/test integration
00338     rval = time_integrate( mbi, method, tag, elems, inttime );
00339     CHK( rval, "time_integrate" );
00340 
00341     return rval;
00342 }
00343 
00344 int main( int argc, char* argv[] )
00345 {
00346     // parse options
00347     ProgOptions opts;
00348     std::string filename;
00349     int dim = 3;
00350     parse_options( opts, dim, filename );
00351     opts.parseCommandLine( argc, argv );
00352     if( filename.empty() )
00353         CHK( MB_FAILURE, "No file specified" );
00354     else if( dim < 1 || dim > 3 )
00355         CHK( MB_FAILURE, "Dimension must be > 0 and <= 3" );
00356 
00357     // read mesh file & gather element handles
00358     Core mbi;
00359     Range elems;
00360     int nv;
00361     EntityType tp;
00362     ErrorCode rval = get_ents( mbi, filename, dim, elems, tp, nv );
00363     if( MB_SUCCESS != rval ) return 1;
00364 
00365     // construct (random) parameter values for queries
00366     unsigned int ne = elems.size();
00367     std::vector< CartVect > params( ne ), coords( ne );
00368     srand( time( NULL ) );
00369     for( unsigned int i = 0; i < ne; i++ )
00370     {
00371         params[i][0] = -1 + 2 * ( (double)rand() ) / RAND_MAX;
00372         if( dim > 1 ) params[i][1] = -1 + 2 * ( (double)rand() ) / RAND_MAX;
00373         if( dim > 2 ) params[i][2] = -1 + 2 * ( (double)rand() ) / RAND_MAX;
00374     }
00375 
00376     // put random field on vertices
00377     Tag tag;
00378     rval = put_random_field( mbi, tag, elems );
00379     CHK( rval, "put_random_field" );
00380     double evtime[2], retime[2], jactime[2],
00381         inttime[2];  // initializations to avoid compiler warning
00382 
00383     rval = elem_evals( &mbi, ELEMEVAL, elems, tag, params, coords, evtime[0], retime[0], jactime[0], inttime[0] );
00384     CHK( rval, "new elem_evals" );
00385 
00386     rval = elem_evals( &mbi, ELEMUTIL, elems, tag, params, coords, evtime[1], retime[1], jactime[1], inttime[1] );
00387     CHK( rval, "old elem_evals" );
00388 
00389     std::cout << filename << ": " << elems.size() << " " << CN::EntityTypeName( tp ) << " elements, " << nv
00390               << " vertices per element." << std::endl
00391               << std::endl;
00392     std::cout << "New, old element evaluation code:" << std::endl;
00393     std::cout << "Evaluation type, time, time per element:" << std::endl;
00394     std::cout << "                             New                   Old            (New/Old)*100" << std::endl;
00395     std::cout << "Forward evaluation " << evtime[0] << ", " << evtime[0] / elems.size() << "    " << evtime[1] << ", "
00396               << evtime[1] / elems.size() << "    " << ( evtime[0] / ( evtime[1] ? evtime[1] : 1 ) ) * 100.0
00397               << std::endl;
00398     std::cout << "Reverse evaluation " << retime[0] << ", " << retime[0] / elems.size() << "    " << retime[1] << ", "
00399               << retime[1] / elems.size() << "    " << ( retime[0] / ( retime[1] ? retime[1] : 1 ) ) * 100.0
00400               << std::endl;
00401     std::cout << "Jacobian           " << jactime[0] << ", " << jactime[0] / elems.size() << "    " << jactime[1]
00402               << ", " << jactime[1] / elems.size() << "    " << ( jactime[0] / ( jactime[1] ? jactime[1] : 1 ) ) * 100.0
00403               << std::endl;
00404     std::cout << "Integration        " << inttime[0] << ", " << inttime[0] / elems.size() << "    " << inttime[1]
00405               << ", " << inttime[1] / elems.size() << "    " << ( inttime[0] / ( inttime[1] ? inttime[1] : 1 ) ) * 100.0
00406               << std::endl;
00407 }
00408 
00409 #if defined( _MSC_VER ) || defined( __MINGW32__ )
00410 double mytime2( double& tot_time, double& utime, double& stime, long& imem, long& rmem )
00411 {
00412     utime    = (double)clock() / CLOCKS_PER_SEC;
00413     tot_time = stime = 0;
00414     imem = rmem = 0;
00415     return tot_time;
00416 }
00417 #else
00418 double mytime2( double& tot_time, double& utime, double& stime, long& imem, long& rmem )
00419 {
00420     struct rusage r_usage;
00421     getrusage( RUSAGE_SELF, &r_usage );
00422     utime    = (double)r_usage.ru_utime.tv_sec + ( (double)r_usage.ru_utime.tv_usec / 1.e6 );
00423     stime    = (double)r_usage.ru_stime.tv_sec + ( (double)r_usage.ru_stime.tv_usec / 1.e6 );
00424     tot_time = utime + stime;
00425 #ifndef LINUX
00426     imem     = r_usage.ru_idrss;
00427     rmem     = r_usage.ru_maxrss;
00428 #else
00429     system( "ps o args,drs,rss | grep perf | grep -v grep" );  // RedHat 9.0 doesnt fill in actual
00430                                                                // memory data
00431     imem = rmem = 0;
00432 #endif
00433     return tot_time;
00434 }
00435 #endif
00436 double mytime()
00437 {
00438     double ttime, utime, stime;
00439     long imem, rmem;
00440     return mytime2( ttime, utime, stime, imem, rmem );
00441 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines