MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 <time.h> 00010 #include <stdlib.h> 00011 #include <assert.h> 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 }