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