MOAB: Mesh Oriented datABase  (version 5.2.1)
perf.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 // Different platforms follow different conventions for usage
00020 #if !defined( _MSC_VER ) && !defined( __MINGW32__ )
00021 #include <sys/resource.h>
00022 #else
00023 #include <time.h>
00024 #endif
00025 #ifdef SOLARIS
00026 extern "C" int getrusage( int, struct rusage* );
00027 #ifndef RUSAGE_SELF
00028 #include </usr/ucbinclude/sys/rusage.h>
00029 #endif
00030 #endif
00031 
00032 #ifndef IS_BUILDING_MB
00033 #define IS_BUILDING_MB
00034 #endif
00035 
00036 #include <stdlib.h>
00037 #include <stdio.h>
00038 #include <assert.h>
00039 #include <iostream>
00040 #include "moab/Core.hpp"
00041 #include "moab/ReadUtilIface.hpp"
00042 #include "VertexSequence.hpp"
00043 #include "StructuredElementSeq.hpp"
00044 #include "EntitySequence.hpp"
00045 #include "SequenceManager.hpp"
00046 #include "moab/HomXform.hpp"
00047 #include "moab/SetIterator.hpp"
00048 
00049 using namespace moab;
00050 
00051 double LENGTH = 1.0;
00052 
00053 void testA( const int nelem, const double* coords );
00054 void testB( const int nelem, const double* coords, int* connect );
00055 void testC( const int nelem, const double* coords );
00056 void testD( const int nelem, const double* coords, int ver );
00057 void testE( const int nelem, const double* coords, int* connect );
00058 void print_time( const bool print_em, double& tot_time, double& utime, double& stime, long& imem, long& rmem );
00059 void query_vert_to_elem();
00060 void query_elem_to_vert();
00061 void query_struct_elem_to_vert();
00062 void query_elem_to_vert_direct();
00063 void query_vert_to_elem_direct();
00064 ErrorCode normalize_elems( double* coords );
00065 void check_answers( const char* );
00066 
00067 #define RC( msg )                                        \
00068     if( MB_SUCCESS != result ) do                        \
00069         {                                                \
00070             std::cout << "FAIL in " << msg << std::endl; \
00071             return;                                      \
00072     } while( true )
00073 #define RR( msg )                                        \
00074     if( MB_SUCCESS != result ) do                        \
00075         {                                                \
00076             std::cout << "FAIL in " << msg << std::endl; \
00077             return result;                               \
00078     } while( true )
00079 
00080 void compute_edge( double* start, const int nelem, const double xint, const int stride )
00081 {
00082     for( int i = 1; i < nelem; i++ )
00083     {
00084         start[i * stride]                     = start[0] + i * xint;
00085         start[nelem + 1 + i * stride]         = start[nelem + 1] + i * xint;
00086         start[2 * ( nelem + 1 ) + i * stride] = start[2 * ( nelem + 1 )] + i * xint;
00087     }
00088 }
00089 
00090 void compute_face( double* a, const int nelem, const double xint, const int stride1, const int stride2 )
00091 {
00092     // 2D TFI on a face starting at a, with strides stride1 in ada and stride2 in tse
00093     for( int j = 1; j < nelem; j++ )
00094     {
00095         double tse = j * xint;
00096         for( int i = 1; i < nelem; i++ )
00097         {
00098             double ada = i * xint;
00099 
00100             a[i * stride1 + j * stride2] =
00101                 ( 1.0 - ada ) * a[i * stride1] + ada * a[i * stride1 + nelem * stride2] +
00102                 ( 1.0 - tse ) * a[j * stride2] + tse * a[j * stride2 + nelem * stride1] -
00103                 ( 1.0 - tse ) * ( 1.0 - ada ) * a[0] - ( 1.0 - tse ) * ada * a[nelem * stride1] -
00104                 tse * ( 1.0 - ada ) * a[nelem * stride2] - tse * ada * a[nelem * ( stride1 + stride2 )];
00105             a[nelem + 1 + i * stride1 + j * stride2] =
00106                 ( 1.0 - ada ) * a[nelem + 1 + i * stride1] + ada * a[nelem + 1 + i * stride1 + nelem * stride2] +
00107                 ( 1.0 - tse ) * a[nelem + 1 + j * stride2] + tse * a[nelem + 1 + j * stride2 + nelem * stride1] -
00108                 ( 1.0 - tse ) * ( 1.0 - ada ) * a[nelem + 1 + 0] -
00109                 ( 1.0 - tse ) * ada * a[nelem + 1 + nelem * stride1] -
00110                 tse * ( 1.0 - ada ) * a[nelem + 1 + nelem * stride2] -
00111                 tse * ada * a[nelem + 1 + nelem * ( stride1 + stride2 )];
00112             a[2 * ( nelem + 1 ) + i * stride1 + j * stride2] =
00113                 ( 1.0 - ada ) * a[2 * ( nelem + 1 ) + i * stride1] +
00114                 ada * a[2 * ( nelem + 1 ) + i * stride1 + nelem * stride2] +
00115                 ( 1.0 - tse ) * a[2 * ( nelem + 1 ) + j * stride2] +
00116                 tse * a[2 * ( nelem + 1 ) + j * stride2 + nelem * stride1] -
00117                 ( 1.0 - tse ) * ( 1.0 - ada ) * a[2 * ( nelem + 1 ) + 0] -
00118                 ( 1.0 - tse ) * ada * a[2 * ( nelem + 1 ) + nelem * stride1] -
00119                 tse * ( 1.0 - ada ) * a[2 * ( nelem + 1 ) + nelem * stride2] -
00120                 tse * ada * a[2 * ( nelem + 1 ) + nelem * ( stride1 + stride2 )];
00121         }
00122     }
00123 }
00124 
00125 void build_coords( const int nelem, double*& coords )
00126 {
00127     double ttime0 = 0.0, ttime1 = 0.0, utime1 = 0.0, stime1 = 0.0;
00128     long imem = 0, rmem = 0;
00129     print_time( false, ttime0, utime1, stime1, imem, rmem );
00130     // allocate the memory
00131     int numv     = nelem + 1;
00132     int numv_sq  = numv * numv;
00133     int tot_numv = numv * numv * numv;
00134     coords       = new double[3 * tot_numv];
00135 
00136 // use FORTRAN-like indexing
00137 #define VINDEX( i, j, k ) ( i + ( j * numv ) + ( k * numv_sq ) )
00138     int idx;
00139     double scale1, scale2, scale3;
00140     // use these to prevent optimization on 1-scale, etc (real map wouldn't have
00141     // all these equal)
00142     scale1 = LENGTH / nelem;
00143     scale2 = LENGTH / nelem;
00144     scale3 = LENGTH / nelem;
00145 
00146 #ifdef REALTFI
00147     // use a real TFI xform to compute coordinates
00148     // compute edges
00149     // i (stride=1)
00150     compute_edge( &coords[VINDEX( 0, 0, 0 )], nelem, scale1, 1 );
00151     compute_edge( &coords[VINDEX( 0, nelem, 0 )], nelem, scale1, 1 );
00152     compute_edge( &coords[VINDEX( 0, 0, nelem )], nelem, scale1, 1 );
00153     compute_edge( &coords[VINDEX( 0, nelem, nelem )], nelem, scale1, 1 );
00154     // j (stride=numv)
00155     compute_edge( &coords[VINDEX( 0, 0, 0 )], nelem, scale1, numv );
00156     compute_edge( &coords[VINDEX( nelem, 0, 0 )], nelem, scale1, numv );
00157     compute_edge( &coords[VINDEX( 0, 0, nelem )], nelem, scale1, numv );
00158     compute_edge( &coords[VINDEX( nelem, 0, nelem )], nelem, scale1, numv );
00159     // k (stride=numv^2)
00160     compute_edge( &coords[VINDEX( 0, 0, 0 )], nelem, scale1, numv_sq );
00161     compute_edge( &coords[VINDEX( nelem, 0, 0 )], nelem, scale1, numv_sq );
00162     compute_edge( &coords[VINDEX( 0, nelem, 0 )], nelem, scale1, numv_sq );
00163     compute_edge( &coords[VINDEX( nelem, nelem, 0 )], nelem, scale1, numv_sq );
00164 
00165     // compute faces
00166     // i=0, nelem
00167     compute_face( &coords[VINDEX( 0, 0, 0 )], nelem, scale1, numv, numv_sq );
00168     compute_face( &coords[VINDEX( nelem, 0, 0 )], nelem, scale1, numv, numv_sq );
00169     // j=0, nelem
00170     compute_face( &coords[VINDEX( 0, 0, 0 )], nelem, scale1, 1, numv_sq );
00171     compute_face( &coords[VINDEX( 0, nelem, 0 )], nelem, scale1, 1, numv_sq );
00172     // k=0, nelem
00173     compute_face( &coords[VINDEX( 0, 0, 0 )], nelem, scale1, 1, numv );
00174     compute_face( &coords[VINDEX( 0, 0, nelem )], nelem, scale1, 1, numv );
00175 
00176     // initialize corner indices
00177     int i000 = VINDEX( 0, 0, 0 );
00178     int ia00 = VINDEX( nelem, 0, 0 );
00179     int i0t0 = VINDEX( 0, nelem, 0 );
00180     int iat0 = VINDEX( nelem, nelem, 0 );
00181     int i00g = VINDEX( 0, 0, nelem );
00182     int ia0g = VINDEX( nelem, 0, nelem );
00183     int i0tg = VINDEX( 0, nelem, nelem );
00184     int iatg = VINDEX( nelem, nelem, nelem );
00185     double cX, cY, cZ;
00186     int adaInts   = nelem;
00187     int tseInts   = nelem;
00188     int gammaInts = nelem;
00189 
00190     for( int i = 1; i < nelem; i++ )
00191     {
00192         for( int j = 1; j < nelem; j++ )
00193         {
00194             for( int k = 1; k < nelem; k++ )
00195             {
00196                 // idx = VINDEX(i,j,k);
00197                 double tse   = i * scale1;
00198                 double ada   = j * scale2;
00199                 double gamma = k * scale3;
00200                 double tm1   = 1.0 - tse;
00201                 double am1   = 1.0 - ada;
00202                 double gm1   = 1.0 - gamma;
00203 
00204                 cX = gm1 * ( am1 * ( tm1 * coords[i000] + tse * coords[i0t0] ) +
00205                              ada * ( tm1 * coords[ia00] + tse * coords[iat0] ) ) +
00206                      gamma * ( am1 * ( tm1 * coords[i00g] + tse * coords[i0tg] ) +
00207                                ada * ( tm1 * coords[ia0g] + tse * coords[iatg] ) );
00208 
00209                 cY = gm1 * ( am1 * ( tm1 * coords[i000] + tse * coords[i0t0] ) +
00210                              ada * ( tm1 * coords[ia00] + tse * coords[iat0] ) ) +
00211                      gamma * ( am1 * ( tm1 * coords[i00g] + tse * coords[i0tg] ) +
00212                                ada * ( tm1 * coords[ia0g] + tse * coords[iatg] ) );
00213 
00214                 cZ = gm1 * ( am1 * ( tm1 * coords[i000] + tse * coords[i0t0] ) +
00215                              ada * ( tm1 * coords[ia00] + tse * coords[iat0] ) ) +
00216                      gamma * ( am1 * ( tm1 * coords[i00g] + tse * coords[i0tg] ) +
00217                                ada * ( tm1 * coords[ia0g] + tse * coords[iatg] ) );
00218 
00219                 double* ai0k = &coords[VINDEX( k, 0, i )];
00220                 double* aiak = &coords[VINDEX( k, adaInts, i )];
00221                 double* a0jk = &coords[VINDEX( k, j, 0 )];
00222                 double* atjk = &coords[VINDEX( k, j, tseInts )];
00223                 double* aij0 = &coords[VINDEX( 0, j, i )];
00224                 double* aijg = &coords[VINDEX( gammaInts, j, i )];
00225 
00226                 coords[VINDEX( i, j, k )] = ( am1 * ai0k[0] + ada * aiak[0] + tm1 * a0jk[0] + tse * atjk[0] +
00227                                               gm1 * aij0[0] + gamma * aijg[0] ) /
00228                                                 2.0 -
00229                                             cX / 2.0;
00230 
00231                 coords[nelem + 1 + VINDEX( i, j, k )] =
00232                     ( am1 * ai0k[nelem + 1] + ada * aiak[nelem + 1] + tm1 * a0jk[nelem + 1] + tse * atjk[nelem + 1] +
00233                       gm1 * aij0[nelem + 1] + gamma * aijg[nelem + 1] ) /
00234                         2.0 -
00235                     cY / 2.0;
00236 
00237                 coords[2 * ( nelem + 1 ) + VINDEX( i, j, k )] =
00238                     ( am1 * ai0k[2 * ( nelem + 1 )] + ada * aiak[2 * ( nelem + 1 )] + tm1 * a0jk[2 * ( nelem + 1 )] +
00239                       tse * atjk[2 * ( nelem + 1 )] + gm1 * aij0[2 * ( nelem + 1 )] +
00240                       gamma * aijg[2 * ( nelem + 1 )] ) /
00241                         2.0 -
00242                     cZ / 2.0;
00243             }
00244         }
00245     }
00246 
00247 #else
00248     for( int i = 0; i < numv; i++ )
00249     {
00250         for( int j = 0; j < numv; j++ )
00251         {
00252             for( int k = 0; k < numv; k++ )
00253             {
00254                 idx = VINDEX( i, j, k );
00255                 // blocked coordinate ordering
00256                 coords[idx]                = i * scale1;
00257                 coords[tot_numv + idx]     = j * scale2;
00258                 coords[2 * tot_numv + idx] = k * scale3;
00259             }
00260         }
00261     }
00262 #endif
00263     print_time( false, ttime1, utime1, stime1, imem, rmem );
00264     //  std::cout << "MOAB: TFI time = " << ttime1-ttime0 << " sec"
00265     //            << std::endl;
00266 }
00267 
00268 void build_connect( const int nelem, const EntityHandle /*vstart*/, int*& connect )
00269 {
00270     // allocate the memory
00271     int nume_tot = nelem * nelem * nelem;
00272     connect      = new int[8 * nume_tot];
00273 
00274     EntityHandle vijk;
00275     int numv    = nelem + 1;
00276     int numv_sq = numv * numv;
00277     int idx     = 0;
00278     for( int i = 0; i < nelem; i++ )
00279     {
00280         for( int j = 0; j < nelem; j++ )
00281         {
00282             for( int k = 0; k < nelem; k++ )
00283             {
00284                 vijk           = VINDEX( i, j, k );
00285                 connect[idx++] = vijk;
00286                 connect[idx++] = vijk + 1;
00287                 connect[idx++] = vijk + 1 + numv;
00288                 connect[idx++] = vijk + numv;
00289                 connect[idx++] = vijk + numv * numv;
00290                 connect[idx++] = vijk + 1 + numv * numv;
00291                 connect[idx++] = vijk + 1 + numv + numv * numv;
00292                 connect[idx++] = vijk + numv + numv * numv;
00293                 assert( i <= numv * numv * numv );
00294             }
00295         }
00296     }
00297 }
00298 
00299 Interface* gMB;
00300 Tag pos_tag, pos2_tag;
00301 
00302 void init()
00303 {
00304     gMB               = new Core();
00305     double def_val[3] = { 0.0, 0.0, 0.0 };
00306     ErrorCode rval =
00307         gMB->tag_get_handle( "position_tag", 3, MB_TYPE_DOUBLE, pos_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00308     assert( MB_SUCCESS == rval );
00309     if( rval ) {}  // empty line to remove compiler warning
00310     rval = gMB->tag_get_handle( "position2_tag", 3, MB_TYPE_DOUBLE, pos2_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00311     assert( MB_SUCCESS == rval );
00312 }
00313 
00314 int main( int argc, char* argv[] )
00315 {
00316     int nelem = 20;
00317     if( argc < 3 )
00318     {
00319         std::cout << "Usage: " << argv[0] << " <ints_per_side> [A|B|C|D [1|2|3|4]|E]" << std::endl;
00320         return 1;
00321     }
00322 
00323     char which_test = '\0';
00324     int ver         = 0;
00325 
00326     sscanf( argv[1], "%d", &nelem );
00327     if( argc >= 3 ) sscanf( argv[2], "%c", &which_test );
00328     if( argc >= 4 ) sscanf( argv[3], "%d", &ver );
00329 
00330     if( 3 <= argc && which_test != 'A' && which_test != 'B' && which_test != 'C' && which_test != 'D' &&
00331         which_test != 'E' )
00332     {
00333         std::cout << "Must indicate A or B, C, D or E for test." << std::endl;
00334         return 1;
00335     }
00336 
00337     if( 4 <= argc && which_test == 'D' && ( ver < 1 || ver > 4 ) )
00338     {
00339         std::cout << "Must indicate version 1, 2, 3, or 4 for test D." << std::endl;
00340         return 1;
00341     }
00342 
00343     //  std::cout << "number of elements: " << nelem << "; test "
00344     //            << which_test << std::endl;
00345 
00346     // pre-build the coords array
00347     double* coords = NULL;
00348     build_coords( nelem, coords );
00349     assert( NULL != coords );
00350 
00351     int* connect = NULL;
00352 
00353     // test A: create structured mesh
00354     if( '\0' == which_test || 'A' == which_test ) testA( nelem, coords );
00355 
00356     build_connect( nelem, 1, connect );
00357 
00358     // test B: create mesh using bulk interface
00359     if( '\0' == which_test || 'B' == which_test ) testB( nelem, coords, connect );
00360 
00361     // test C: create mesh using individual interface
00362     if( '\0' == which_test || 'C' == which_test ) testC( nelem, coords );
00363 
00364     // test D: query mesh using iterators
00365     if( '\0' == which_test || 'D' == which_test ) testD( nelem, coords, ver );
00366 
00367     // test E: query mesh using direct access
00368     if( '\0' == which_test || 'E' == which_test ) testE( nelem, coords, connect );
00369 
00370     return 0;
00371 }
00372 
00373 void query_elem_to_vert()
00374 {
00375     Range all_hexes;
00376     ErrorCode result = gMB->get_entities_by_type( 0, MBHEX, all_hexes );
00377     RC( "query_elem_to_vert" );
00378     const EntityHandle* connect;
00379     int num_connect;
00380     double dum_coords[24];
00381     for( Range::iterator eit = all_hexes.begin(); eit != all_hexes.end(); ++eit )
00382     {
00383         result = gMB->get_connectivity( *eit, connect, num_connect );
00384         RC( "query_elem_to_vert" );
00385         result = gMB->get_coords( connect, num_connect, dum_coords );
00386         RC( "query_elem_to_vert" );
00387 
00388         // compute the centroid
00389         double centroid[3] = { 0.0, 0.0, 0.0 };
00390         for( int j = 0; j < 8; j++ )
00391         {
00392             centroid[0] += dum_coords[3 * j + 0];
00393             centroid[1] += dum_coords[3 * j + 1];
00394             centroid[2] += dum_coords[3 * j + 2];
00395         }
00396         centroid[0] *= 0.125;
00397         centroid[1] *= 0.125;
00398         centroid[2] *= 0.125;
00399         result = gMB->tag_set_data( pos_tag, &( *eit ), 1, centroid );
00400         RC( "query_elem_to_vert" );
00401     }
00402 }
00403 
00404 void query_vert_to_elem()
00405 {
00406     Range all_verts;
00407     std::vector< EntityHandle > neighbor_hexes;
00408     std::vector< double > neighbor_pos;
00409     double coords[3];
00410     neighbor_pos.resize( 3 * 8 );  // average vertex will have 8 adjacent hexes
00411     ErrorCode result = gMB->get_entities_by_type( 0, MBVERTEX, all_verts );
00412     RC( "query_vert_to_elem" );
00413     for( Range::iterator vit = all_verts.begin(); vit != all_verts.end(); ++vit )
00414     {
00415         neighbor_hexes.clear();
00416         result = gMB->get_coords( &( *vit ), 1, coords );
00417         RC( "query_vert_to_elem" );
00418         result = gMB->get_adjacencies( &( *vit ), 1, 3, false, neighbor_hexes );
00419         RC( "query_vert_to_elem" );
00420         assert( neighbor_pos.size() >= 3 * neighbor_hexes.size() );
00421         result = gMB->tag_get_data( pos2_tag, &neighbor_hexes[0], neighbor_hexes.size(), &neighbor_pos[0] );
00422         RC( "query_vert_to_elem" );
00423         for( unsigned int i = 0; i < neighbor_hexes.size(); i++ )
00424         {
00425             neighbor_pos[3 * i] += coords[0];
00426             neighbor_pos[3 * i + 1] += coords[1];
00427             neighbor_pos[3 * i + 2] += coords[2];
00428         }
00429 
00430         result = gMB->tag_set_data( pos2_tag, &neighbor_hexes[0], neighbor_hexes.size(), &neighbor_pos[0] );
00431         RC( "query_vert_to_elem" );
00432     }
00433 
00434     // get all hexes and divide pos_tag by 8; reuse all_verts
00435     result = normalize_elems( coords );
00436     RC( "query_vert_to_elem" );
00437 }
00438 
00439 void query_elem_to_vert_iters( int chunk_size, bool check_valid, std::vector< EntityHandle >& connect,
00440                                double* dum_coords, double* dum_pos )
00441 {
00442     std::vector< EntityHandle > hexes;
00443     SetIterator* iter;
00444     ErrorCode result = gMB->create_set_iterator( 0, MBHEX, -1, chunk_size, check_valid, iter );
00445     RC( "query_elem_to_vert_iters" );
00446     bool atend = false;
00447     while( !atend )
00448     {
00449         hexes.clear();
00450         result = iter->get_next_arr( hexes, atend );
00451         RC( "query_elem_to_vert_iters" );
00452         result = gMB->get_connectivity( &hexes[0], hexes.size(), connect );
00453         RC( "query_elem_to_vert_iters" );
00454         result = gMB->get_coords( &connect[0], connect.size(), dum_coords );
00455         RC( "query_elem_to_vert_iters" );
00456         result = gMB->tag_get_data( pos_tag, &hexes[0], hexes.size(), dum_pos );
00457         RC( "query_elem_to_vert_iters" );
00458         for( unsigned int i = 0; i < hexes.size(); i++ )
00459         {
00460             // compute the centroid
00461             for( int j = 0; j < 8; j++ )
00462             {
00463                 dum_pos[3 * i + 0] += dum_coords[24 * i + 3 * j];
00464                 dum_pos[3 * i + 1] += dum_coords[24 * i + 3 * j + 1];
00465                 dum_pos[3 * i + 2] += dum_coords[24 * i + 3 * j + 2];
00466             }
00467             dum_pos[3 * i + 0] *= 0.125;
00468             dum_pos[3 * i + 1] *= 0.125;
00469             dum_pos[3 * i + 2] *= 0.125;
00470         }
00471         result = gMB->tag_set_data( pos_tag, &hexes[0], hexes.size(), dum_pos );
00472         RC( "query_elem_to_vert_iters" );
00473     }
00474 
00475     delete iter;
00476 }
00477 
00478 void query_vert_to_elem_iters( int chunk_size, bool check_valid, std::vector< EntityHandle >& /*connect*/,
00479                                double* dum_coords, double* dum_pos )
00480 {
00481     std::vector< EntityHandle > verts, neighbor_hexes;
00482     SetIterator* iter;
00483     ErrorCode result = gMB->create_set_iterator( 0, MBVERTEX, -1, chunk_size, check_valid, iter );
00484     RC( "query_vert_to_elem_iters" );
00485     assert( MB_SUCCESS == result );
00486     bool atend = false;
00487     while( !atend )
00488     {
00489         verts.clear();
00490         result = iter->get_next_arr( verts, atend );
00491         RC( "query_vert_to_elem_iters" );
00492         result = gMB->get_coords( &verts[0], verts.size(), dum_coords );
00493         RC( "query_vert_to_elem_iters" );
00494         chunk_size = std::min( (int)verts.size(), chunk_size );
00495         for( int i = 0; i < chunk_size; i++ )
00496         {
00497             neighbor_hexes.clear();
00498             result = gMB->get_adjacencies( &verts[i], 1, 3, false, neighbor_hexes );
00499             RC( "query_vert_to_elem_iters" );
00500             result = gMB->tag_get_data( pos2_tag, &neighbor_hexes[0], neighbor_hexes.size(), dum_pos );
00501             RC( "query_vert_to_elem_iters" );
00502             for( unsigned int j = 0; j < neighbor_hexes.size(); j++ )
00503             {
00504                 dum_pos[3 * j + 0] += dum_coords[3 * i + 0];
00505                 dum_pos[3 * j + 1] += dum_coords[3 * i + 1];
00506                 dum_pos[3 * j + 2] += dum_coords[3 * i + 2];
00507             }
00508             result = gMB->tag_set_data( pos2_tag, &neighbor_hexes[0], neighbor_hexes.size(), dum_pos );
00509             RC( "query_vert_to_elem_iters" );
00510         }
00511     }
00512 
00513     result = normalize_elems( dum_coords );
00514     RC( "query_vert_to_elem_iters" );
00515 }
00516 
00517 // normalize pos_tag on all elems by 1/8; coords assumed large enough to hold 3 doubles
00518 ErrorCode normalize_elems( double* coords )
00519 {
00520     // get all hexes and divide pos_tag by 8
00521     Range elems;
00522     ErrorCode result = gMB->get_entities_by_type( 0, MBHEX, elems );
00523     RR( "normalize" );
00524 
00525     for( Range::iterator vit = elems.begin(); vit != elems.end(); ++vit )
00526     {
00527         result = gMB->tag_get_data( pos2_tag, &( *vit ), 1, coords );
00528         RR( "normalize" );
00529         coords[0] *= 0.125;
00530         coords[1] *= 0.125;
00531         coords[2] *= 0.125;
00532         result = gMB->tag_set_data( pos2_tag, &( *vit ), 1, coords );
00533         RR( "normalize" );
00534     }
00535 
00536     return result;
00537 }
00538 
00539 void query_struct_elem_to_vert()
00540 {
00541     // assumes brick mapped mesh with handles starting at zero
00542     Range all_hexes;
00543     ErrorCode result = gMB->get_entities_by_type( 0, MBHEX, all_hexes );
00544     RC( "query_struct_elem_to_vert" );
00545     double dum_coords[24];
00546     std::vector< EntityHandle > connect;
00547     for( Range::iterator eit = all_hexes.begin(); eit != all_hexes.end(); ++eit )
00548     {
00549         result = gMB->get_connectivity( &( *eit ), 1, connect );
00550         RC( "query_struct_elem_to_vert" );
00551         result = gMB->get_coords( &connect[0], connect.size(), dum_coords );
00552         RC( "query_struct_elem_to_vert" );
00553 
00554         double centroid[3] = { 0.0, 0.0, 0.0 };
00555         for( int j = 0; j < 8; j++ )
00556         {
00557             centroid[0] += dum_coords[3 * j + 0];
00558             centroid[1] += dum_coords[3 * j + 1];
00559             centroid[2] += dum_coords[3 * j + 2];
00560         }
00561         centroid[0] *= 0.125;
00562         centroid[1] *= 0.125;
00563         centroid[2] *= 0.125;
00564         result = gMB->tag_set_data( pos_tag, &( *eit ), 1, centroid );
00565         RC( "query_struct_elem_to_vert" );
00566     }
00567 }
00568 
00569 #if defined( _MSC_VER ) || defined( __MINGW32__ )
00570 void print_time( const bool print_em, double& tot_time, double& utime, double& stime, long& imem, long& rmem )
00571 {
00572     utime = (double)clock() / CLOCKS_PER_SEC;
00573     if( print_em ) std::cout << "Total wall time = " << utime << std::endl;
00574     tot_time = stime = 0;
00575     imem = rmem = 0;
00576 }
00577 #else
00578 void print_time( const bool print_em, double& tot_time, double& utime, double& stime, long& imem, long& rmem )
00579 {
00580     struct rusage r_usage;
00581     getrusage( RUSAGE_SELF, &r_usage );
00582     utime = (double)r_usage.ru_utime.tv_sec + ( (double)r_usage.ru_utime.tv_usec / 1.e6 );
00583     stime = (double)r_usage.ru_stime.tv_sec + ( (double)r_usage.ru_stime.tv_usec / 1.e6 );
00584     tot_time = utime + stime;
00585     if( print_em )
00586         std::cout << "User, system, total time = " << utime << ", " << stime << ", " << tot_time << std::endl;
00587 #ifndef LINUX
00588     if( print_em )
00589     {
00590         std::cout << "Max resident set size = " << r_usage.ru_maxrss << " kbytes" << std::endl;
00591         std::cout << "Int resident set size = " << r_usage.ru_idrss << std::endl;
00592     }
00593     imem = r_usage.ru_idrss;
00594     rmem = r_usage.ru_maxrss;
00595 #else
00596     system( "ps o args,drs,rss | grep perf | grep -v grep" );  // RedHat 9.0 doesnt fill in actual
00597                                                                // memory data
00598     imem = rmem = 0;
00599 #endif
00600 }
00601 #endif
00602 
00603 void testA( const int nelem, const double* coords )
00604 {
00605     double ttime0 = 0.0, ttime1 = 0.0, ttime2 = 0.0, ttime3 = 0.0, ttime4 = 0.0, utime = 0.0, stime = 0.0;
00606     long imem0 = 0, rmem0 = 0, imem1 = 0, rmem1 = 0, imem2 = 0, rmem2 = 0, imem3 = 0, rmem3 = 0, imem4 = 0, rmem4 = 0;
00607 
00608     print_time( false, ttime0, utime, stime, imem0, rmem0 );
00609 
00610     // make a 3d block of vertices
00611     EntitySequence* dum_seq    = NULL;
00612     ScdVertexData* vseq        = NULL;
00613     StructuredElementSeq* eseq = NULL;
00614     init();
00615     Core* mbcore = dynamic_cast< Core* >( gMB );
00616     assert( mbcore != NULL );
00617     SequenceManager* seq_mgr = mbcore->sequence_manager();
00618     HomCoord vseq_minmax[2]  = { HomCoord( 0, 0, 0 ), HomCoord( nelem, nelem, nelem ) };
00619     EntityHandle vstart, estart;
00620 
00621     ErrorCode result = seq_mgr->create_scd_sequence( vseq_minmax[0], vseq_minmax[1], MBVERTEX, 1, vstart, dum_seq );
00622     RC( "testA" );
00623     if( NULL != dum_seq ) vseq = dynamic_cast< ScdVertexData* >( dum_seq->data() );
00624     assert( MB_FAILURE != result && vstart != 0 && dum_seq != NULL && vseq != NULL );
00625     // now the element sequence
00626     result = seq_mgr->create_scd_sequence( vseq_minmax[0], vseq_minmax[1], MBHEX, 1, estart, dum_seq );
00627     if( NULL != dum_seq ) eseq = dynamic_cast< StructuredElementSeq* >( dum_seq );
00628     assert( MB_FAILURE != result && estart != 0 && dum_seq != NULL && eseq != NULL );
00629 
00630     // only need to add one vseq to this, unity transform
00631     // trick: if I know it's going to be unity, just input 3 sets of equivalent points
00632     result = eseq->sdata()->add_vsequence( vseq, vseq_minmax[0], vseq_minmax[0], vseq_minmax[0], vseq_minmax[0],
00633                                            vseq_minmax[0], vseq_minmax[0] );
00634     assert( MB_SUCCESS == result );
00635 
00636     // set the coordinates of the vertices
00637     EntityHandle handle;
00638     int i;
00639     double dumv[3];
00640     int num_verts = ( nelem + 1 ) * ( nelem + 1 ) * ( nelem + 1 );
00641     for( i = 0, handle = vstart; i < num_verts; i++, handle++ )
00642     {
00643         dumv[0] = coords[i];
00644         dumv[1] = coords[num_verts + i];
00645         dumv[2] = coords[2 * num_verts + i];
00646         result  = gMB->set_coords( &handle, 1, dumv );
00647         assert( MB_SUCCESS == result );
00648     }
00649 
00650     print_time( false, ttime1, utime, stime, imem1, rmem1 );
00651 
00652     // query the mesh 2 ways
00653     query_struct_elem_to_vert();
00654 
00655     print_time( false, ttime2, utime, stime, imem2, rmem2 );
00656 
00657     query_vert_to_elem();
00658 
00659     print_time( false, ttime3, utime, stime, imem3, rmem3 );
00660 
00661 #ifndef NDEBUG
00662     check_answers( "A" );
00663 #endif
00664 
00665     delete gMB;
00666 
00667     print_time( false, ttime4, utime, stime, imem4, rmem4 );
00668 
00669     std::cout << "MOAB_scd:nelem,construct,e_to_v,v_to_e,after_dtor,total= " << nelem << " " << ttime1 - ttime0 << " "
00670               << ttime2 - ttime1 << " " << ttime3 - ttime2 << " " << ttime4 - ttime3 << " " << ttime4 - ttime0
00671               << " seconds" << std::endl;
00672     std::cout << "MOAB_scd_memory(rss):initial,after_construction,e-v,v-e,after_dtor= " << rmem0 << " " << rmem1 << " "
00673               << rmem2 << " " << rmem3 << " " << rmem4 << " kb" << std::endl;
00674 }
00675 
00676 void testB( const int nelem, const double* coords, int* connect )
00677 {
00678     double ttime0 = 0.0, ttime1 = 0.0, ttime2 = 0.0, ttime3 = 0.0, ttime4 = 0.0, utime = 0.0, stime = 0.0;
00679     long imem0 = 0, rmem0 = 0, imem1 = 0, rmem1 = 0, imem2 = 0, rmem2 = 0, imem3 = 0, rmem3 = 0, imem4 = 0, rmem4 = 0;
00680 
00681     print_time( false, ttime0, utime, stime, imem0, rmem0 );
00682 
00683     int num_verts = ( nelem + 1 ) * ( nelem + 1 ) * ( nelem + 1 );
00684     int num_elems = nelem * nelem * nelem;
00685     EntityHandle vstart, estart;
00686 
00687     // get the read interface
00688     ReadUtilIface* readMeshIface;
00689     init();
00690     gMB->query_interface( readMeshIface );
00691 
00692     // create a sequence to hold the node coordinates
00693     // get the current number of entities and start at the next slot
00694     std::vector< double* > coord_arrays;
00695     ErrorCode result = readMeshIface->get_node_coords( 3, num_verts, 1, vstart, coord_arrays );
00696     RC( "testB" );
00697     assert( MB_SUCCESS == result && 1 == vstart && coord_arrays[0] && coord_arrays[1] && coord_arrays[2] );
00698     // memcpy the coordinate data into place
00699     memcpy( coord_arrays[0], coords, sizeof( double ) * num_verts );
00700     memcpy( coord_arrays[1], &coords[num_verts], sizeof( double ) * num_verts );
00701     memcpy( coord_arrays[2], &coords[2 * num_verts], sizeof( double ) * num_verts );
00702 
00703     EntityHandle* conn = 0;
00704     result             = readMeshIface->get_element_connect( num_elems, 8, MBHEX, 1, estart, conn );
00705     assert( MB_SUCCESS == result );
00706     for( int i = 0; i < num_elems * 8; i++ )
00707         conn[i] = vstart + connect[i];
00708 
00709     delete[] connect;
00710 
00711     result = readMeshIface->update_adjacencies( estart, num_elems, 8, conn );
00712     assert( MB_SUCCESS == result );
00713 
00714     print_time( false, ttime1, utime, stime, imem1, rmem1 );
00715 
00716     // query the mesh 2 ways
00717     query_elem_to_vert();
00718 
00719     print_time( false, ttime2, utime, stime, imem2, rmem2 );
00720 
00721     query_vert_to_elem();
00722 
00723     print_time( false, ttime3, utime, stime, imem3, rmem3 );
00724 
00725 #ifndef NDEBUG
00726     check_answers( "B" );
00727 #endif
00728 
00729     delete gMB;
00730 
00731     print_time( false, ttime4, utime, stime, imem4, rmem4 );
00732 
00733     std::cout << "MOAB_ucd_blocked:nelem,construct,e_to_v,v_to_e,after_dtor,total= " << nelem << " " << ttime1 - ttime0
00734               << " " << ttime2 - ttime1 << " " << ttime3 - ttime2 << " " << ttime4 - ttime3 << " " << ttime4 - ttime0
00735               << " seconds" << std::endl;
00736     std::cout << "MOAB_ucdblocked_memory_(rss):initial,after_construction,e-v,v-e,after_dtor= " << rmem0 << " " << rmem1
00737               << " " << rmem2 << " " << rmem3 << " " << rmem4 << " kb" << std::endl;
00738 }
00739 
00740 void testC( const int nelem, const double* coords )
00741 {
00742     double ttime0 = 0.0, ttime1 = 0.0, ttime2 = 0.0, ttime3 = 0.0, ttime4 = 0.0, utime = 0.0, stime = 0.0;
00743     long imem0 = 0, rmem0 = 0, imem1 = 0, rmem1 = 0, imem2 = 0, rmem2 = 0, imem3 = 0, rmem3 = 0, imem4 = 0, rmem4 = 0;
00744 
00745     print_time( false, ttime0, utime, stime, imem0, rmem0 );
00746 
00747     // create the vertices; assume we don't need to keep a list of vertex handles, since they'll
00748     // be created in sequence
00749     int numv             = nelem + 1;
00750     int numv_sq          = numv * numv;
00751     int num_verts        = numv * numv * numv;
00752     double dum_coords[3] = { coords[0], coords[num_verts], coords[2 * num_verts] };
00753     EntityHandle vstart;
00754 
00755     init();
00756     ErrorCode result = gMB->create_vertex( dum_coords, vstart );
00757     RC( "testC" );
00758     assert( MB_SUCCESS == result && 1 == vstart );
00759 
00760     EntityHandle dum_vert, vijk;
00761     int i;
00762     for( i = 1; i < num_verts; i++ )
00763     {
00764         dum_coords[0] = coords[i];
00765         dum_coords[1] = coords[num_verts + i];
00766         dum_coords[2] = coords[2 * num_verts + i];
00767         result        = gMB->create_vertex( dum_coords, dum_vert );
00768         assert( MB_SUCCESS == result );
00769     }
00770 
00771     EntityHandle dum_conn[8];
00772     for( i = 0; i < nelem; i++ )
00773     {
00774         for( int j = 0; j < nelem; j++ )
00775         {
00776             for( int k = 0; k < nelem; k++ )
00777             {
00778                 vijk        = vstart + VINDEX( i, j, k );
00779                 dum_conn[0] = vijk;
00780                 dum_conn[1] = vijk + 1;
00781                 dum_conn[2] = vijk + 1 + numv;
00782                 dum_conn[3] = vijk + numv;
00783                 dum_conn[4] = vijk + numv * numv;
00784                 dum_conn[5] = vijk + 1 + numv * numv;
00785                 dum_conn[6] = vijk + 1 + numv + numv * numv;
00786                 dum_conn[7] = vijk + numv + numv * numv;
00787                 result      = gMB->create_element( MBHEX, dum_conn, 8, dum_vert );
00788                 assert( MB_SUCCESS == result );
00789             }
00790         }
00791     }
00792 
00793     print_time( false, ttime1, utime, stime, imem1, rmem1 );
00794 
00795     // query the mesh 2 ways
00796     query_elem_to_vert();
00797 
00798     print_time( false, ttime2, utime, stime, imem2, rmem2 );
00799 
00800     query_vert_to_elem();
00801 
00802     print_time( false, ttime3, utime, stime, imem3, rmem3 );
00803 
00804 #ifndef NDEBUG
00805     check_answers( "C" );
00806 #endif
00807 
00808     delete gMB;
00809 
00810     print_time( false, ttime4, utime, stime, imem4, rmem4 );
00811 
00812     std::cout << "MOAB_ucd_indiv:nelem,construct,e_to_v,v_to_e,after_dtor,total= " << nelem << " " << ttime1 - ttime0
00813               << " " << ttime2 - ttime1 << " " << ttime3 - ttime2 << " " << ttime4 - ttime3 << " " << ttime4 - ttime0
00814               << " seconds" << std::endl;
00815     std::cout << "MOAB_ucd_indiv_memory_(rss):initial,after_construction,e-v,v-e,after_dtor= " << rmem0 << " " << rmem1
00816               << " " << rmem2 << " " << rmem3 << " " << rmem4 << " kb" << std::endl;
00817 }
00818 
00819 void testD( const int nelem, const double* coords, int ver )
00820 {
00821     double ttime0 = 0.0, ttime1 = 0.0, ttime2 = 0.0, ttime3 = 0.0, ttime4 = 0.0, ttime5 = 0.0, ttime6 = 0.0,
00822            ttime7 = 0.0, ttime8 = 0.0, ttime9 = 0.0, ttime10 = 0.0, utime = 0.0, stime = 0.0;
00823     long imem0 = 0, rmem0 = 0, imem1 = 0, rmem1 = 0, imem2 = 0, rmem2 = 0, imem3 = 0, rmem3 = 0, imem4 = 0, rmem4 = 0,
00824          imem5 = 0, rmem5 = 0, imem6 = 0, rmem6 = 0, imem7 = 0, rmem7 = 0, imem8 = 0, rmem8 = 0, imem9 = 0, rmem9 = 0,
00825          imem10 = 0, rmem10 = 0;
00826 
00827     print_time( false, ttime0, utime, stime, imem0, rmem0 );
00828 
00829     // create the vertices; assume we don't need to keep a list of vertex handles, since they'll
00830     // be created in sequence
00831     int numv      = nelem + 1;
00832     int numv_sq   = numv * numv;
00833     int num_verts = numv * numv * numv;
00834     std::vector< double > dum_coords( 24 ), dum_pos( 24 );
00835     dum_coords[0] = coords[0];
00836     dum_coords[1] = coords[num_verts];
00837     dum_coords[2] = coords[2 * num_verts];
00838     EntityHandle vstart;
00839 
00840     init();
00841     ErrorCode result = gMB->create_vertex( &dum_coords[0], vstart );
00842     RC( "testD" );
00843     assert( MB_SUCCESS == result && 1 == vstart );
00844 
00845     EntityHandle dum_vert, vijk;
00846     int i;
00847     for( i = 1; i < num_verts; i++ )
00848     {
00849         dum_coords[0] = coords[i];
00850         dum_coords[1] = coords[num_verts + i];
00851         dum_coords[2] = coords[2 * num_verts + i];
00852         result        = gMB->create_vertex( &dum_coords[0], dum_vert );
00853         assert( MB_SUCCESS == result );
00854     }
00855 
00856     EntityHandle dum_conn[8];
00857     for( i = 0; i < nelem; i++ )
00858     {
00859         for( int j = 0; j < nelem; j++ )
00860         {
00861             for( int k = 0; k < nelem; k++ )
00862             {
00863                 vijk        = vstart + VINDEX( i, j, k );
00864                 dum_conn[0] = vijk;
00865                 dum_conn[1] = vijk + 1;
00866                 dum_conn[2] = vijk + 1 + numv;
00867                 dum_conn[3] = vijk + numv;
00868                 dum_conn[4] = vijk + numv * numv;
00869                 dum_conn[5] = vijk + 1 + numv * numv;
00870                 dum_conn[6] = vijk + 1 + numv + numv * numv;
00871                 dum_conn[7] = vijk + numv + numv * numv;
00872                 result      = gMB->create_element( MBHEX, dum_conn, 8, dum_vert );
00873                 assert( MB_SUCCESS == result );
00874             }
00875         }
00876     }
00877 
00878     print_time( false, ttime1, utime, stime, imem1, rmem1 );
00879 
00880     // query the mesh 2 ways with !check_valid
00881     std::vector< EntityHandle > connect( 8 );
00882 #ifndef NDEBUG
00883     // used only in debug mode
00884     double def_val[3] = { 0.0, 0.0, 0.0 };
00885 #endif
00886     if( ver == 0 || ver == 1 )
00887     {
00888         query_elem_to_vert_iters( 1, false, connect, &dum_coords[0], &dum_pos[0] );
00889         print_time( false, ttime2, utime, stime, imem2, rmem2 );
00890         query_vert_to_elem_iters( 1, false, connect, &dum_coords[0], &dum_pos[0] );
00891         print_time( false, ttime3, utime, stime, imem3, rmem3 );
00892 #ifndef NDEBUG
00893         check_answers( "D" );
00894         result = gMB->tag_delete( pos_tag );
00895         assert( MB_SUCCESS == result );
00896         result =
00897             gMB->tag_get_handle( "position_tag", 3, MB_TYPE_DOUBLE, pos_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00898         assert( MB_SUCCESS == result );
00899         result = gMB->tag_delete( pos2_tag );
00900         assert( MB_SUCCESS == result );
00901         result =
00902             gMB->tag_get_handle( "position2_tag", 3, MB_TYPE_DOUBLE, pos2_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00903         assert( MB_SUCCESS == result );
00904 #endif
00905     }
00906 
00907     if( ver == 0 || ver == 2 )
00908     {
00909         if( ver != 0 ) print_time( false, ttime3, utime, stime, imem3, rmem3 );
00910         query_elem_to_vert_iters( 1, true, connect, &dum_coords[0], &dum_pos[0] );
00911         print_time( false, ttime4, utime, stime, imem4, rmem4 );
00912         query_vert_to_elem_iters( 1, true, connect, &dum_coords[0], &dum_pos[0] );
00913         print_time( false, ttime5, utime, stime, imem5, rmem5 );
00914 #ifndef NDEBUG
00915         check_answers( "D" );
00916         result = gMB->tag_delete( pos_tag );
00917         assert( MB_SUCCESS == result );
00918         result =
00919             gMB->tag_get_handle( "position_tag", 3, MB_TYPE_DOUBLE, pos_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00920         assert( MB_SUCCESS == result );
00921         result = gMB->tag_delete( pos2_tag );
00922         assert( MB_SUCCESS == result );
00923         result =
00924             gMB->tag_get_handle( "position2_tag", 3, MB_TYPE_DOUBLE, pos2_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00925         assert( MB_SUCCESS == result );
00926 #endif
00927     }
00928 
00929     if( ver == 0 || ver >= 3 )
00930     {
00931         dum_coords.resize( 2400 );
00932         dum_pos.resize( 300 );
00933     }
00934     if( ver == 0 || ver == 3 )
00935     {
00936         if( ver != 0 ) print_time( false, ttime5, utime, stime, imem3, rmem3 );
00937         query_elem_to_vert_iters( 100, false, connect, &dum_coords[0], &dum_pos[0] );
00938         print_time( false, ttime6, utime, stime, imem6, rmem6 );
00939         query_vert_to_elem_iters( 100, false, connect, &dum_coords[0], &dum_pos[0] );
00940         print_time( false, ttime7, utime, stime, imem7, rmem7 );
00941 #ifndef NDEBUG
00942         check_answers( "D" );
00943         result = gMB->tag_delete( pos_tag );
00944         assert( MB_SUCCESS == result );
00945         result =
00946             gMB->tag_get_handle( "position_tag", 3, MB_TYPE_DOUBLE, pos_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00947         assert( MB_SUCCESS == result );
00948         result = gMB->tag_delete( pos2_tag );
00949         assert( MB_SUCCESS == result );
00950         result =
00951             gMB->tag_get_handle( "position2_tag", 3, MB_TYPE_DOUBLE, pos2_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00952         assert( MB_SUCCESS == result );
00953 #endif
00954     }
00955 
00956     if( ver == 0 || ver == 4 )
00957     {
00958         if( ver != 0 ) print_time( false, ttime7, utime, stime, imem3, rmem3 );
00959         query_elem_to_vert_iters( 100, true, connect, &dum_coords[0], &dum_pos[0] );
00960         print_time( false, ttime8, utime, stime, imem8, rmem8 );
00961         query_vert_to_elem_iters( 100, true, connect, &dum_coords[0], &dum_pos[0] );
00962         print_time( false, ttime9, utime, stime, imem9, rmem9 );
00963 #ifndef NDEBUG
00964         check_answers( "D" );
00965         result = gMB->tag_delete( pos_tag );
00966         assert( MB_SUCCESS == result );
00967         result =
00968             gMB->tag_get_handle( "position_tag", 3, MB_TYPE_DOUBLE, pos_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00969         assert( MB_SUCCESS == result );
00970         result = gMB->tag_delete( pos2_tag );
00971         assert( MB_SUCCESS == result );
00972         result =
00973             gMB->tag_get_handle( "position2_tag", 3, MB_TYPE_DOUBLE, pos2_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00974         assert( MB_SUCCESS == result );
00975 #endif
00976     }
00977 
00978     if( ver > 0 && ver < 4 ) print_time( false, ttime9, utime, stime, imem9, rmem9 );
00979     delete gMB;
00980 
00981     print_time( false, ttime10, utime, stime, imem10, rmem10 );
00982 
00983     if( ver == 0 || ver == 1 )
00984     {
00985         std::cout << "MOAB_ucd_iters_!check_valid_1:nelem,construct,e-v,v-e,after_dtor,total= " << nelem << " "
00986                   << ttime1 - ttime0 << " " << ttime2 - ttime1 << " " << ttime3 - ttime2 << " " << ttime10 - ttime9
00987                   << " " << ttime3 - ttime0 + ttime10 - ttime9 << std::endl;
00988         std::cout << "MOAB_ucd_iters_memory_(rss)_!check_valid_1:initial,after_construction,e-v,v-"
00989                      "e,after_dtor= "
00990                   << rmem0 << " " << rmem1 << " " << rmem2 << " " << rmem3 << " " << rmem10 << " kb" << std::endl;
00991     }
00992     if( ver == 0 || ver == 2 )
00993     {
00994         std::cout << "MOAB_ucd_iters_check_valid_1:nelem,construct,e-v,v-e,after_dtor,total= " << nelem << " "
00995                   << ttime1 - ttime0 << " " << ttime4 - ttime3 << " " << ttime5 - ttime4 << " " << ttime10 - ttime9
00996                   << " " << ttime1 - ttime0 + ttime5 - ttime3 + ttime10 - ttime9 << std::endl;
00997         std::cout << "MOAB_ucd_iters_memory_(rss)_check_valid_1:initial,after_construction,e-v,v-e,"
00998                      "after_dtor= "
00999                   << rmem0 << " " << rmem1 << " " << rmem2 << " " << rmem3 << " " << rmem10 << " kb" << std::endl;
01000     }
01001     if( ver == 0 || ver == 3 )
01002     {
01003         std::cout << "MOAB_ucd_iters_!check_valid_100:nelem,construct,e-v,v-e,after_dtor,total= " << nelem << " "
01004                   << ttime1 - ttime0 << " " << ttime6 - ttime5 << " " << ttime7 - ttime6 << " " << ttime10 - ttime9
01005                   << " " << ttime1 - ttime0 + ttime7 - ttime5 + ttime10 - ttime9 << std::endl;
01006         std::cout << "MOAB_ucd_iters_memory_(rss)_!check_valid_100:initial,after_construction,e-v,"
01007                      "v-e,after_dtor= "
01008                   << rmem0 << " " << rmem1 << " " << rmem6 << " " << rmem7 << " " << rmem10 << " kb" << std::endl;
01009     }
01010     if( ver == 0 || ver == 4 )
01011     {
01012         std::cout << "MOAB_ucd_iters_check_valid_100:nelem,construct,e-v,v-e,after_dtor,total= " << nelem << " "
01013                   << ttime1 - ttime0 << " " << ttime8 - ttime7 << " " << ttime9 - ttime8 << " " << ttime10 - ttime9
01014                   << " " << ttime1 - ttime0 + ttime10 - ttime7 << std::endl;
01015         std::cout << "MOAB_ucd_iters_memory_(rss)_check_valid_100:initial,after_construction,e-v,v-"
01016                      "e,after_dtor= "
01017                   << rmem0 << " " << rmem1 << " " << rmem8 << " " << rmem9 << " " << rmem10 << " kb" << std::endl;
01018     }
01019 }
01020 
01021 void testE( const int nelem, const double* coords, int* connect )
01022 {
01023     double ttime0 = 0.0, ttime1 = 0.0, ttime2 = 0.0, ttime3 = 0.0, ttime4 = 0.0, ttime5 = 0.0, ttime6 = 0.0,
01024            utime = 0.0, stime = 0.0;
01025     long imem0 = 0, rmem0 = 0, imem1 = 0, rmem1 = 0, imem2 = 0, rmem2 = 0, imem3 = 0, rmem3 = 0, imem4 = 0, rmem4 = 0,
01026          imem5 = 0, rmem5 = 0, imem6 = 0, rmem6 = 0;
01027 
01028     print_time( false, ttime0, utime, stime, imem0, rmem0 );
01029 
01030     int num_verts = ( nelem + 1 ) * ( nelem + 1 ) * ( nelem + 1 );
01031     int num_elems = nelem * nelem * nelem;
01032     EntityHandle vstart, estart;
01033 
01034     // get the read interface
01035     ReadUtilIface* readMeshIface;
01036     init();
01037     gMB->query_interface( readMeshIface );
01038 
01039     // create a sequence to hold the node coordinates
01040     // get the current number of entities and start at the next slot
01041     std::vector< double* > coord_arrays;
01042     ErrorCode result = readMeshIface->get_node_coords( 3, num_verts, 1, vstart, coord_arrays );
01043     RC( "testE" );
01044     assert( MB_SUCCESS == result && 1 == vstart && coord_arrays[0] && coord_arrays[1] && coord_arrays[2] );
01045     // memcpy the coordinate data into place
01046     memcpy( coord_arrays[0], coords, sizeof( double ) * num_verts );
01047     memcpy( coord_arrays[1], &coords[num_verts], sizeof( double ) * num_verts );
01048     memcpy( coord_arrays[2], &coords[2 * num_verts], sizeof( double ) * num_verts );
01049 
01050     EntityHandle* conn = 0;
01051     result             = readMeshIface->get_element_connect( num_elems, 8, MBHEX, 1, estart, conn );
01052     assert( MB_SUCCESS == result );
01053     for( int i = 0; i < num_elems * 8; i++ )
01054         conn[i] = vstart + connect[i];
01055 
01056     delete[] connect;
01057 
01058     Range verts( vstart, vstart + num_verts - 1 ), elems( estart, estart + num_elems - 1 );
01059 
01060     result = readMeshIface->update_adjacencies( estart, num_elems, 8, conn );
01061     assert( MB_SUCCESS == result );
01062 
01063     print_time( false, ttime1, utime, stime, imem1, rmem1 );
01064 
01065     // query the mesh 2 ways
01066     query_elem_to_vert_direct();
01067 
01068     print_time( false, ttime2, utime, stime, imem2, rmem2 );
01069 
01070     query_vert_to_elem_direct();
01071 
01072     print_time( false, ttime3, utime, stime, imem3, rmem3 );
01073 
01074 #ifndef NDEBUG
01075     check_answers( "E" );
01076 #endif
01077 
01078     query_elem_to_vert_direct();
01079 
01080     print_time( false, ttime4, utime, stime, imem4, rmem4 );
01081 
01082     query_vert_to_elem_direct();
01083 
01084     print_time( false, ttime5, utime, stime, imem5, rmem5 );
01085 
01086 #ifndef NDEBUG
01087     check_answers( "E" );
01088 #endif
01089 
01090     delete gMB;
01091 
01092     print_time( false, ttime6, utime, stime, imem6, rmem6 );
01093 
01094     std::cout << "MOAB_ucd_direct:nelem,construct,e_to_v,v_to_e,after_dtor,total= " << nelem << " " << ttime1 - ttime0
01095               << " " << ttime2 - ttime1 << " " << ttime3 - ttime2 << " " << ttime6 - ttime5 << " "
01096               << ttime3 - ttime0 + ttime6 - ttime5 << " seconds" << std::endl;
01097     std::cout << "MOAB_ucd_direct_memory_(rss):initial,after_construction,e-v,v-e,after_dtor= " << rmem0 << " " << rmem1
01098               << " " << rmem2 << " " << rmem3 << " " << rmem6 << " kb" << std::endl;
01099 
01100     std::cout << "MOAB_ucd_direct2:nelem,construct,e_to_v,v_to_e,after_dtor,total= " << nelem << " " << ttime1 - ttime0
01101               << " " << ttime4 - ttime3 << " " << ttime5 - ttime4 << " " << ttime6 - ttime5 << " "
01102               << ttime1 - ttime0 + ttime6 - ttime3 << " seconds" << std::endl;
01103     std::cout << "MOAB_ucd_direct2_memory_(rss):initial,after_construction,e-v,v-e,after_dtor= " << rmem0 << " "
01104               << rmem1 << " " << rmem4 << " " << rmem5 << " " << rmem6 << " kb" << std::endl;
01105 }
01106 
01107 void query_elem_to_vert_direct()
01108 {
01109     Range all_hexes, all_verts;
01110     ErrorCode result = gMB->get_entities_by_type( 0, MBHEX, all_hexes );
01111     assert( MB_SUCCESS == result );
01112     result = gMB->get_entities_by_type( 0, MBVERTEX, all_verts );
01113     RC( "query_elem_to_vert_direct" );
01114     EntityHandle* connect;
01115     int ecount, vcount, vpere;
01116     double* coords[3];
01117 
01118     result = gMB->connect_iterate( all_hexes.begin(), all_hexes.end(), connect, vpere, ecount );
01119     if( MB_SUCCESS != result || ecount != (int)all_hexes.size() )
01120     {
01121         std::cout << "FAILED in connect_iterate!" << std::endl;
01122         return;
01123     }
01124     result = gMB->coords_iterate( all_verts.begin(), all_verts.end(), coords[0], coords[1], coords[2], vcount );
01125     if( MB_SUCCESS != result || vcount != (int)all_verts.size() )
01126     {
01127         std::cout << "FAILED in coords_iterate!" << std::endl;
01128         return;
01129     }
01130     double* centroid;
01131     result = gMB->tag_iterate( pos_tag, all_hexes.begin(), all_hexes.end(), ecount, (void*&)centroid );
01132     if( MB_SUCCESS != result || ecount != (int)all_hexes.size() )
01133     {
01134         std::cout << "FAILED in connect_iterate!" << std::endl;
01135         return;
01136     }
01137 
01138     EntityHandle vstart = *all_verts.begin();
01139     for( int i = 0; i < ecount; i++ )
01140     {
01141         // compute the centroid
01142         for( int j = 0; j < vpere; j++ )
01143         {
01144             int vind = *connect - vstart;
01145             connect++;
01146             centroid[3 * i + 0] += coords[0][vind];
01147             centroid[3 * i + 1] += coords[1][vind];
01148             centroid[3 * i + 2] += coords[2][vind];
01149         }
01150 
01151         // now normalize
01152         centroid[3 * i + 0] *= 0.125;
01153         centroid[3 * i + 1] *= 0.125;
01154         centroid[3 * i + 2] *= 0.125;
01155     }
01156 }
01157 
01158 void query_vert_to_elem_direct()
01159 {
01160     Range all_verts, tmp_ents;
01161     ErrorCode result = gMB->get_entities_by_type( 0, MBVERTEX, all_verts );
01162     assert( MB_SUCCESS == result );
01163 
01164     // make sure vertex-element adjacencies are created
01165     result = gMB->get_adjacencies( &( *all_verts.begin() ), 1, 3, false, tmp_ents );
01166     assert( MB_SUCCESS == result );
01167 
01168     const std::vector< EntityHandle >** adjs;
01169     int count;
01170     result = gMB->adjacencies_iterate( all_verts.begin(), all_verts.end(), adjs, count );
01171     if( MB_SUCCESS != result && count != (int)all_verts.size() )
01172     {
01173         std::cout << "FAILED:adjacencies_iterate." << std::endl;
01174         return;
01175     }
01176 
01177     double* coords[3];
01178     result = gMB->coords_iterate( all_verts.begin(), all_verts.end(), coords[0], coords[1], coords[2], count );
01179     if( MB_SUCCESS != result || count != (int)all_verts.size() )
01180     {
01181         std::cout << "FAILED in coords_iterate!" << std::endl;
01182         return;
01183     }
01184     // get all hexes, then iterator over pos2_tag
01185     double* centroid;
01186     int ecount;
01187     tmp_ents.clear();
01188     result = gMB->get_entities_by_type( 0, MBHEX, tmp_ents );
01189     assert( MB_SUCCESS == result );
01190     result = gMB->tag_iterate( pos2_tag, tmp_ents.begin(), tmp_ents.end(), ecount, (void*&)centroid );
01191     if( MB_SUCCESS != result || ecount != (int)tmp_ents.size() )
01192     {
01193         std::cout << "FAILED in tag_iterate!" << std::endl;
01194         return;
01195     }
01196 
01197     int i;
01198     Range::iterator vit;
01199     EntityHandle estart = *tmp_ents.begin();
01200     for( vit = all_verts.begin(), i = 0; vit != all_verts.end(); ++vit, i++ )
01201     {
01202         assert( adjs[i] );
01203         for( std::vector< EntityHandle >::const_iterator vit2 = adjs[i]->begin(); vit2 != adjs[i]->end(); ++vit2 )
01204             if( *vit >= estart )
01205             {
01206                 int eind = *vit2 - estart;
01207                 centroid[3 * eind + 0] += coords[0][i];
01208                 centroid[3 * eind + 1] += coords[1][i];
01209                 centroid[3 * eind + 2] += coords[2][i];
01210             }
01211     }
01212 
01213     // now normalize
01214     for( i = 0; i < (int)tmp_ents.size(); i++ )
01215     {
01216         centroid[3 * i + 0] *= 0.125;
01217         centroid[3 * i + 1] *= 0.125;
01218         centroid[3 * i + 2] *= 0.125;
01219     }
01220 }
01221 
01222 void check_answers( const char* /*test_name*/ )
01223 {
01224     Range elems;
01225     ErrorCode result = gMB->get_entities_by_type( 0, MBHEX, elems );
01226     RC( "check_answers" );
01227 
01228     double coords1[3], coords2[3], del[3];
01229     double diff = 0.0;
01230     for( Range::iterator vit = elems.begin(); vit != elems.end(); ++vit )
01231     {
01232         result = gMB->tag_get_data( pos_tag, &( *vit ), 1, coords1 );
01233         RC( "check_answers" );
01234         result = gMB->tag_get_data( pos2_tag, &( *vit ), 1, coords2 );
01235         RC( "check_answers" );
01236         for( int i = 0; i < 3; i++ )
01237             del[i] = fabs( coords1[i] - coords2[i] );
01238         if( del[0] || del[1] || del[2] )
01239         {
01240             double len2 = std::max( coords1[0] * coords1[0] + coords1[1] * coords1[1] + coords1[2] * coords1[2],
01241                                     coords2[0] * coords2[0] + coords2[1] * coords2[1] + coords2[2] * coords2[2] ),
01242                    num  = del[0] * del[0] + del[1] * del[1] + del[2] * del[2];
01243             if( len2 > 0.0 )
01244                 diff = std::max( diff, num / sqrt( len2 ) );
01245             else if( num > 0.0 )
01246                 diff = sqrt( num );
01247         }
01248     }
01249     if( diff > 0.0 ) std::cout << "Max relative difference = " << diff << std::endl;
01250 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines