MOAB: Mesh Oriented datABase  (version 5.2.1)
perftool.cpp
Go to the documentation of this file.
00001 /**
00002  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
00003  * storing and accessing finite element mesh data.
00004  *
00005  * Copyright 2004 Sandia Corporation.  Under the terms of Contract
00006  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
00007  * retains certain rights in this software.
00008  *
00009  * This library is free software; you can redistribute it and/or
00010  * modify it under the terms of the GNU Lesser General Public
00011  * License as published by the Free Software Foundation; either
00012  * version 2.1 of the License, or (at your option) any later version.
00013  *
00014  */
00015 
00016 // MOAB performance tests building mapped mesh with nodes and
00017 // hexes created one at a time.  This also creates the node to hex adjacencies.
00018 
00019 #include <math.h>
00020 #include <iostream>
00021 #include <iomanip>
00022 #include <time.h>
00023 #include <assert.h>
00024 #include <list>
00025 #include "moab/Core.hpp"
00026 #include "moab/Skinner.hpp"
00027 #include "moab/ReadUtilIface.hpp"
00028 
00029 using namespace moab;
00030 
00031 double LENGTH               = 1.0;
00032 const int DEFAULT_INTERVALS = 50;
00033 
00034 void create_regular_mesh( int interval, int dimension );
00035 void skin_common( int interval, int dim, int blocks, bool use_adj );
00036 void skin( int intervals, int dim, int num )
00037 {
00038     std::cout << "Skinning w/out adjacencies:" << std::endl;
00039     skin_common( intervals, dim, num, false );
00040 }
00041 void skin_adj( int intervals, int dim, int num )
00042 {
00043     std::cout << "Skinning with adjacencies:" << std::endl;
00044     skin_common( intervals, dim, num, true );
00045 }
00046 
00047 void tag_time( TagType storage, bool direct, int intervals, int dim, int blocks );
00048 
00049 void dense_tag( int intervals, int dim, int blocks )
00050 {
00051     std::cout << "Dense Tag Time:";
00052     tag_time( MB_TAG_DENSE, false, intervals, dim, blocks );
00053 }
00054 void sparse_tag( int intervals, int dim, int blocks )
00055 {
00056     std::cout << "Sparse Tag Time:";
00057     tag_time( MB_TAG_SPARSE, false, intervals, dim, blocks );
00058 }
00059 void direct_tag( int intervals, int dim, int blocks )
00060 {
00061     std::cout << "Direct Tag Time:";
00062     tag_time( MB_TAG_DENSE, true, intervals, dim, blocks );
00063 }
00064 
00065 typedef void ( *test_func_t )( int, int, int );
00066 const struct
00067 {
00068     std::string testName;
00069     test_func_t testFunc;
00070     std::string testDesc;
00071 } TestList[] = {
00072     { "skin", &skin, "Test time to get skin mesh w/out adjacencies" },
00073     { "skin_adj", &skin_adj, "Test time to get skin mesh with adjacencies" },
00074     { "sparse", &sparse_tag, "Sparse tag data manipulation" },
00075     { "dense", &dense_tag, "Dense tag data manipulation" },
00076     { "direct", &direct_tag, "Dense tag data manipulation using direct data access" },
00077 };
00078 const int TestListSize = sizeof( TestList ) / sizeof( TestList[0] );
00079 
00080 void usage( const char* argv0, bool error = true )
00081 {
00082     std::ostream& str = error ? std::cerr : std::cout;
00083     str << "Usage: " << argv0
00084         << " [-i <ints_per_side>] [-d <dimension>] [-n <test_specifc_int>] <test_name> "
00085            "[<test_name2> ...]"
00086         << std::endl;
00087     str << "       " << argv0 << " [-h|-l]" << std::endl;
00088     if( error ) return;
00089     str << "  -i  : specify interverals per side (num hex = ints^3, default: " << DEFAULT_INTERVALS << std::endl;
00090     str << "  -d  : specify element dimension, default: 3" << std::endl;
00091     str << "  -n  : specify an integer value that for which the meaning is test-specific" << std::endl;
00092     str << "  -h  : print this help text." << std::endl;
00093     str << "  -l  : list available tests" << std::endl;
00094 }
00095 
00096 void list_tests()
00097 {
00098     unsigned max_test_name = 0, max_test_desc = 0;
00099     for( int i = 0; i < TestListSize; ++i )
00100     {
00101         if( TestList[i].testName.size() > max_test_name ) max_test_name = TestList[i].testName.size();
00102         if( TestList[i].testDesc.size() > max_test_desc ) max_test_desc = TestList[i].testDesc.size();
00103     }
00104     std::cout << std::setw( max_test_name ) << "NAME"
00105               << "   " << std::setw( max_test_desc ) << std::left << "DESCRIPTION" << std::endl
00106               << std::setfill( '-' ) << std::setw( max_test_name ) << ""
00107               << "   " << std::setfill( '-' ) << std::setw( max_test_desc ) << "" << std::setfill( ' ' ) << std::endl;
00108     for( int i = 0; i < TestListSize; ++i )
00109         std::cout << std::setw( max_test_name ) << TestList[i].testName << " : " << std::setw( max_test_desc )
00110                   << std::left << TestList[i].testDesc << std::endl;
00111 }
00112 
00113 int main( int argc, char* argv[] )
00114 {
00115     int intervals = DEFAULT_INTERVALS;
00116     int dimension = 3;
00117     int number    = 0;
00118     std::vector< test_func_t > test_list;
00119     std::list< int* > expected_list;
00120     bool did_help = false;
00121     for( int i = 1; i < argc; ++i )
00122     {
00123         if( !expected_list.empty() )
00124         {
00125             int* ptr = expected_list.front();
00126             expected_list.pop_front();
00127             char* endptr;
00128             *ptr = strtol( argv[i], &endptr, 0 );
00129             if( !endptr )
00130             {
00131                 usage( argv[0] );
00132                 std::cerr << "Expected integer value, got \"" << argv[i] << '"' << std::endl;
00133                 return 1;
00134             }
00135         }
00136         else if( *argv[i] == '-' )
00137         {  // flag
00138             for( int j = 1; argv[i][j]; ++j )
00139             {
00140                 switch( argv[i][j] )
00141                 {
00142                     case 'i':
00143                         expected_list.push_back( &intervals );
00144                         break;
00145                     case 'd':
00146                         expected_list.push_back( &dimension );
00147                         break;
00148                     case 'n':
00149                         expected_list.push_back( &number );
00150                         break;
00151                     case 'h':
00152                         did_help = true;
00153                         usage( argv[0], false );
00154                         break;
00155                     case 'l':
00156                         did_help = true;
00157                         list_tests();
00158                         break;
00159                     default:
00160                         usage( argv[0] );
00161                         std::cerr << "Invalid option: -" << argv[i][j] << std::endl;
00162                         return 1;
00163                 }
00164             }
00165         }
00166         else
00167         {
00168             int j = -1;
00169             do
00170             {
00171                 ++j;
00172                 if( j >= TestListSize )
00173                 {
00174                     usage( argv[0] );
00175                     std::cerr << "Invalid test name: " << argv[i] << std::endl;
00176                     return 1;
00177                 }
00178             } while( TestList[j].testName != argv[i] );
00179             test_list.push_back( TestList[j].testFunc );
00180         }
00181     }
00182 
00183     if( !expected_list.empty() )
00184     {
00185         usage( argv[0] );
00186         std::cerr << "Missing final argument" << std::endl;
00187         return 1;
00188     }
00189 
00190     // error if no input
00191     if( test_list.empty() && !did_help )
00192     {
00193         usage( argv[0] );
00194         std::cerr << "No tests specified" << std::endl;
00195         return 1;
00196     }
00197 
00198     if( intervals < 1 )
00199     {
00200         std::cerr << "Invalid interval count: " << intervals << std::endl;
00201         return 1;
00202     }
00203 
00204     if( dimension < 1 || dimension > 3 )
00205     {
00206         std::cerr << "Invalid dimension: " << dimension << std::endl;
00207         return 1;
00208     }
00209 
00210     // now run the tests
00211     for( std::vector< test_func_t >::iterator i = test_list.begin(); i != test_list.end(); ++i )
00212     {
00213         test_func_t fptr = *i;
00214         fptr( intervals, dimension, number );
00215     }
00216 
00217     return 0;
00218 }
00219 
00220 void create_regular_mesh( Interface* gMB, int interval, int dim )
00221 {
00222     if( dim < 1 || dim > 3 || interval < 1 )
00223     {
00224         std::cerr << "Invalid arguments" << std::endl;
00225         exit( 1 );
00226     }
00227 
00228     const int nvi     = interval + 1;
00229     const int dims[3] = { nvi, dim > 1 ? nvi : 1, dim > 2 ? nvi : 1 };
00230     int num_vert      = dims[0] * dims[1] * dims[2];
00231 
00232     ReadUtilIface* readMeshIface;
00233     gMB->query_interface( readMeshIface );
00234 
00235     EntityHandle vstart;
00236     std::vector< double* > arrays;
00237     ErrorCode rval = readMeshIface->get_node_coords( 3, num_vert, 1, vstart, arrays );
00238     if( MB_SUCCESS != rval || arrays.size() < 3 )
00239     {
00240         std::cerr << "Vertex creation failed" << std::endl;
00241         exit( 2 );
00242     }
00243     double *x = arrays[0], *y = arrays[1], *z = arrays[2];
00244 
00245     // Calculate vertex coordinates
00246     for( int k = 0; k < dims[2]; ++k )
00247         for( int j = 0; j < dims[1]; ++j )
00248             for( int i = 0; i < dims[0]; ++i )
00249             {
00250                 *x = i;
00251                 ++x;
00252                 *y = j;
00253                 ++y;
00254                 *z = k;
00255                 ++z;
00256             }
00257 
00258     const long vert_per_elem = 1 << dim;  // 2^dim
00259     const long intervals[3]  = { interval, dim > 1 ? interval : 1, dim > 2 ? interval : 1 };
00260     const long num_elem      = intervals[0] * intervals[1] * intervals[2];
00261     const EntityType type    = ( dim == 1 ) ? MBEDGE : ( dim == 2 ) ? MBQUAD : MBHEX;
00262 
00263     EntityHandle estart, *conn = 0;
00264     rval = readMeshIface->get_element_connect( num_elem, vert_per_elem, type, 0, estart, conn );
00265     if( MB_SUCCESS != rval || !conn )
00266     {
00267         std::cerr << "Element creation failed" << std::endl;
00268         exit( 2 );
00269     }
00270 
00271     // Offsets of element vertices in grid relative to corner closest to origin
00272     long c                = dims[0] * dims[1];
00273     const long corners[8] = { 0, 1, 1 + dims[0], dims[0], c, c + 1, c + 1 + dims[0], c + dims[0] };
00274 
00275     // Populate element list
00276     EntityHandle* iter = conn;
00277     for( long z1 = 0; z1 < intervals[2]; ++z1 )
00278         for( long y1 = 0; y1 < intervals[1]; ++y1 )
00279             for( long x1 = 0; x1 < intervals[0]; ++x1 )
00280             {
00281                 const long index = x1 + y1 * dims[0] + z1 * ( dims[0] * dims[1] );
00282                 for( long j = 0; j < vert_per_elem; ++j, ++iter )
00283                     *iter = index + corners[j] + vstart;
00284             }
00285 
00286     // notify MOAB of the new elements
00287     rval = readMeshIface->update_adjacencies( estart, num_elem, vert_per_elem, conn );
00288     if( MB_SUCCESS != rval )
00289     {
00290         std::cerr << "Element update failed" << std::endl;
00291         exit( 2 );
00292     }
00293 }
00294 
00295 void skin_common( int interval, int dim, int num, bool use_adj )
00296 {
00297     Core moab;
00298     Interface* gMB = &moab;
00299     ErrorCode rval;
00300     double d;
00301     clock_t t, tt;
00302 
00303     create_regular_mesh( gMB, interval, dim );
00304 
00305     Range skin, verts, elems;
00306     rval = gMB->get_entities_by_dimension( 0, dim, elems );
00307     assert( MB_SUCCESS == rval );
00308     assert( !elems.empty() );
00309 
00310     Skinner tool( gMB );
00311 
00312     t    = clock();
00313     rval = tool.find_skin( 0, elems, true, verts, 0, use_adj, false );
00314     t    = clock() - t;
00315     if( MB_SUCCESS != rval )
00316     {
00317         std::cerr << "Search for skin vertices failed" << std::endl;
00318         exit( 2 );
00319     }
00320     d = ( (double)t ) / CLOCKS_PER_SEC;
00321     std::cout << "Got " << verts.size() << " skin vertices in " << d << " seconds." << std::endl;
00322 
00323     t = 0;
00324     if( num < 1 ) num = 1000;
00325     long blocksize = elems.size() / num;
00326     if( !blocksize ) blocksize = 1;
00327     long numblocks     = elems.size() / blocksize;
00328     Range::iterator it = elems.begin();
00329     for( long i = 0; i < numblocks; ++i )
00330     {
00331         verts.clear();
00332         Range::iterator end = it + blocksize;
00333         Range blockelems;
00334         blockelems.merge( it, end );
00335         it   = end;
00336         tt   = clock();
00337         rval = tool.find_skin( 0, blockelems, true, verts, 0, use_adj, false );
00338         t += clock() - tt;
00339         if( MB_SUCCESS != rval )
00340         {
00341             std::cerr << "Search for skin vertices failed" << std::endl;
00342             exit( 2 );
00343         }
00344     }
00345     d = ( (double)t ) / CLOCKS_PER_SEC;
00346     std::cout << "Got skin vertices for " << numblocks << " blocks of " << blocksize << " elements in " << d
00347               << " seconds." << std::endl;
00348 
00349     for( int e = 0; e < 2; ++e )
00350     {  // do this twice
00351         if( e == 1 )
00352         {
00353             // create all interior faces
00354             skin.clear();
00355             t = clock();
00356             gMB->get_adjacencies( elems, dim - 1, true, skin, Interface::UNION );
00357             t = clock() - t;
00358             d = ( (double)t ) / CLOCKS_PER_SEC;
00359             std::cout << "Created " << skin.size() << " entities of dimension-1 in " << d << " seconds" << std::endl;
00360         }
00361 
00362         skin.clear();
00363         t    = clock();
00364         rval = tool.find_skin( 0, elems, false, skin, 0, use_adj, true );
00365         t    = clock() - t;
00366         if( MB_SUCCESS != rval )
00367         {
00368             std::cerr << "Search for skin vertices failed" << std::endl;
00369             exit( 2 );
00370         }
00371         d = ( (double)t ) / CLOCKS_PER_SEC;
00372         std::cout << "Got " << skin.size() << " skin elements in " << d << " seconds." << std::endl;
00373 
00374         t  = 0;
00375         it = elems.begin();
00376         for( long i = 0; i < numblocks; ++i )
00377         {
00378             skin.clear();
00379             Range::iterator end = it + blocksize;
00380             Range blockelems;
00381             blockelems.merge( it, end );
00382             it   = end;
00383             tt   = clock();
00384             rval = tool.find_skin( 0, blockelems, false, skin, 0, use_adj, true );
00385             t += clock() - tt;
00386             if( MB_SUCCESS != rval )
00387             {
00388                 std::cerr << "Search for skin elements failed" << std::endl;
00389                 exit( 2 );
00390             }
00391         }
00392         d = ( (double)t ) / CLOCKS_PER_SEC;
00393         std::cout << "Got skin elements for " << numblocks << " blocks of " << blocksize << " elements in " << d
00394                   << " seconds." << std::endl;
00395     }
00396 }
00397 
00398 void tag_time( TagType storage, bool direct, int intervals, int dim, int blocks )
00399 {
00400     Core moab;
00401     Interface& mb = moab;
00402     create_regular_mesh( &mb, intervals, dim );
00403 
00404     // Create tag in which to store data
00405     Tag tag;
00406     mb.tag_get_handle( "data", 1, MB_TYPE_DOUBLE, tag, storage | MB_TAG_CREAT );
00407 
00408     // Make up some arbitrary iterative calculation for timing purposes:
00409     // set each value v_n = (V + v_n)/2 until all values are within
00410     // epsilon of V.
00411     std::vector< double > data;
00412     Range verts;
00413     mb.get_entities_by_type( 0, MBVERTEX, verts );
00414 
00415     clock_t t = clock();
00416 
00417     // initialize
00418     if( direct )
00419     {
00420         Range::iterator i, j = verts.begin();
00421         void* ptr;
00422         int count;
00423         while( j != verts.end() )
00424         {
00425             mb.tag_iterate( tag, i, verts.end(), count, ptr );
00426             double* arr = reinterpret_cast< double* >( ptr );
00427             for( j = i + count; i != j; ++i, ++arr )
00428                 *arr = ( 11.0 * *i + 7.0 ) / ( *i );
00429         }
00430     }
00431     else
00432     {
00433         data.resize( verts.size() );
00434         double* arr = &data[0];
00435         for( Range::iterator i = verts.begin(); i != verts.end(); ++i, ++arr )
00436             *arr = ( 11.0 * *i + 7.0 ) / ( *i );
00437         mb.tag_set_data( tag, verts, &data[0] );
00438     }
00439 
00440     // iterate
00441     const double v0   = acos( -1.0 );  // pi
00442     size_t iter_count = 0;
00443     double max_diff;
00444     const size_t num_verts = verts.size();
00445     do
00446     {
00447         if( direct )
00448         {
00449             max_diff = 0.0;
00450             Range::iterator i, j = verts.begin();
00451             void* ptr;
00452             while( j != verts.end() )
00453             {
00454                 int count;
00455                 mb.tag_iterate( tag, i, verts.end(), count, ptr );
00456                 double* arr = reinterpret_cast< double* >( ptr );
00457 
00458                 for( j = i + count; i != j; ++i, ++arr )
00459                 {
00460                     *arr        = 0.5 * ( *arr + v0 );
00461                     double diff = fabs( *arr - v0 );
00462                     if( diff > max_diff ) max_diff = diff;
00463                 }
00464             }
00465         }
00466         else if( blocks < 1 )
00467         {
00468             max_diff = 0.0;
00469             mb.tag_get_data( tag, verts, &data[0] );
00470             for( size_t i = 0; i < data.size(); ++i )
00471             {
00472                 data[i]     = 0.5 * ( v0 + data[i] );
00473                 double diff = fabs( data[i] - v0 );
00474                 if( diff > max_diff ) max_diff = diff;
00475             }
00476             mb.tag_set_data( tag, verts, &data[0] );
00477         }
00478         else
00479         {
00480             max_diff = 0.0;
00481             Range r;
00482             Range::iterator it = verts.begin();
00483             size_t step        = num_verts / blocks;
00484             for( int j = 0; j < blocks - 1; ++j )
00485             {
00486                 Range::iterator nx = it;
00487                 nx += step;
00488                 r.clear();
00489                 r.merge( it, nx );
00490                 mb.tag_get_data( tag, r, &data[0] );
00491                 it = nx;
00492 
00493                 for( size_t i = 0; i < step; ++i )
00494                 {
00495                     data[i]     = 0.5 * ( v0 + data[i] );
00496                     double diff = fabs( data[i] - v0 );
00497                     if( diff > max_diff ) max_diff = diff;
00498                 }
00499                 mb.tag_set_data( tag, r, &data[0] );
00500             }
00501 
00502             r.clear();
00503             r.merge( it, verts.end() );
00504             mb.tag_get_data( tag, r, &data[0] );
00505             for( size_t i = 0; i < ( num_verts - ( blocks - 1 ) * step ); ++i )
00506             {
00507                 data[i]     = 0.5 * ( v0 + data[i] );
00508                 double diff = fabs( data[i] - v0 );
00509                 if( diff > max_diff ) max_diff = diff;
00510             }
00511             mb.tag_set_data( tag, r, &data[0] );
00512         }
00513         ++iter_count;
00514         //    std::cout << iter_count << " " << max_diff << std::endl;
00515     } while( max_diff > 1e-6 );
00516 
00517     double secs = ( clock() - t ) / (double)CLOCKS_PER_SEC;
00518     std::cout << " " << iter_count << " iterations in " << secs << " seconds" << std::endl;
00519 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines