MOAB: Mesh Oriented datABase  (version 5.4.0)
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 <cstdlib>
00037 #include <cstdio>
00038 #include <cassert>
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 )
00310     {
00311     }  // empty line to remove compiler warning
00312     rval = gMB->tag_get_handle( "position2_tag", 3, MB_TYPE_DOUBLE, pos2_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00313     assert( MB_SUCCESS == rval );
00314 }
00315 
00316 int main( int argc, char* argv[] )
00317 {
00318     int nelem = 20;
00319     if( argc < 3 )
00320     {
00321         std::cout << "Usage: " << argv[0] << " <ints_per_side> [A|B|C|D [1|2|3|4]|E]" << std::endl;
00322         return 1;
00323     }
00324 
00325     char which_test = '\0';
00326     int ver         = 0;
00327 
00328     sscanf( argv[1], "%d", &nelem );
00329     if( argc >= 3 ) sscanf( argv[2], "%c", &which_test );
00330     if( argc >= 4 ) sscanf( argv[3], "%d", &ver );
00331 
00332     if( 3 <= argc && which_test != 'A' && which_test != 'B' && which_test != 'C' && which_test != 'D' &&
00333         which_test != 'E' )
00334     {
00335         std::cout << "Must indicate A or B, C, D or E for test." << std::endl;
00336         return 1;
00337     }
00338 
00339     if( 4 <= argc && which_test == 'D' && ( ver < 1 || ver > 4 ) )
00340     {
00341         std::cout << "Must indicate version 1, 2, 3, or 4 for test D." << std::endl;
00342         return 1;
00343     }
00344 
00345     //  std::cout << "number of elements: " << nelem << "; test "
00346     //            << which_test << std::endl;
00347 
00348     // pre-build the coords array
00349     double* coords = NULL;
00350     build_coords( nelem, coords );
00351     assert( NULL != coords );
00352 
00353     int* connect = NULL;
00354 
00355     // test A: create structured mesh
00356     if( '\0' == which_test || 'A' == which_test ) testA( nelem, coords );
00357 
00358     build_connect( nelem, 1, connect );
00359 
00360     // test B: create mesh using bulk interface
00361     if( '\0' == which_test || 'B' == which_test ) testB( nelem, coords, connect );
00362 
00363     // test C: create mesh using individual interface
00364     if( '\0' == which_test || 'C' == which_test ) testC( nelem, coords );
00365 
00366     // test D: query mesh using iterators
00367     if( '\0' == which_test || 'D' == which_test ) testD( nelem, coords, ver );
00368 
00369     // test E: query mesh using direct access
00370     if( '\0' == which_test || 'E' == which_test ) testE( nelem, coords, connect );
00371 
00372     return 0;
00373 }
00374 
00375 void query_elem_to_vert()
00376 {
00377     Range all_hexes;
00378     ErrorCode result = gMB->get_entities_by_type( 0, MBHEX, all_hexes );
00379     RC( "query_elem_to_vert" );
00380     const EntityHandle* connect;
00381     int num_connect;
00382     double dum_coords[24];
00383     for( Range::iterator eit = all_hexes.begin(); eit != all_hexes.end(); ++eit )
00384     {
00385         result = gMB->get_connectivity( *eit, connect, num_connect );
00386         RC( "query_elem_to_vert" );
00387         result = gMB->get_coords( connect, num_connect, dum_coords );
00388         RC( "query_elem_to_vert" );
00389 
00390         // compute the centroid
00391         double centroid[3] = { 0.0, 0.0, 0.0 };
00392         for( int j = 0; j < 8; j++ )
00393         {
00394             centroid[0] += dum_coords[3 * j + 0];
00395             centroid[1] += dum_coords[3 * j + 1];
00396             centroid[2] += dum_coords[3 * j + 2];
00397         }
00398         centroid[0] *= 0.125;
00399         centroid[1] *= 0.125;
00400         centroid[2] *= 0.125;
00401         result = gMB->tag_set_data( pos_tag, &( *eit ), 1, centroid );
00402         RC( "query_elem_to_vert" );
00403     }
00404 }
00405 
00406 void query_vert_to_elem()
00407 {
00408     Range all_verts;
00409     std::vector< EntityHandle > neighbor_hexes;
00410     std::vector< double > neighbor_pos;
00411     double coords[3];
00412     neighbor_pos.resize( 3 * 8 );  // average vertex will have 8 adjacent hexes
00413     ErrorCode result = gMB->get_entities_by_type( 0, MBVERTEX, all_verts );
00414     RC( "query_vert_to_elem" );
00415     for( Range::iterator vit = all_verts.begin(); vit != all_verts.end(); ++vit )
00416     {
00417         neighbor_hexes.clear();
00418         result = gMB->get_coords( &( *vit ), 1, coords );
00419         RC( "query_vert_to_elem" );
00420         result = gMB->get_adjacencies( &( *vit ), 1, 3, false, neighbor_hexes );
00421         RC( "query_vert_to_elem" );
00422         assert( neighbor_pos.size() >= 3 * neighbor_hexes.size() );
00423         result = gMB->tag_get_data( pos2_tag, &neighbor_hexes[0], neighbor_hexes.size(), &neighbor_pos[0] );
00424         RC( "query_vert_to_elem" );
00425         for( unsigned int i = 0; i < neighbor_hexes.size(); i++ )
00426         {
00427             neighbor_pos[3 * i] += coords[0];
00428             neighbor_pos[3 * i + 1] += coords[1];
00429             neighbor_pos[3 * i + 2] += coords[2];
00430         }
00431 
00432         result = gMB->tag_set_data( pos2_tag, &neighbor_hexes[0], neighbor_hexes.size(), &neighbor_pos[0] );
00433         RC( "query_vert_to_elem" );
00434     }
00435 
00436     // get all hexes and divide pos_tag by 8; reuse all_verts
00437     result = normalize_elems( coords );
00438     RC( "query_vert_to_elem" );
00439 }
00440 
00441 void query_elem_to_vert_iters( int chunk_size,
00442                                bool check_valid,
00443                                std::vector< EntityHandle >& connect,
00444                                double* dum_coords,
00445                                double* dum_pos )
00446 {
00447     std::vector< EntityHandle > hexes;
00448     SetIterator* iter;
00449     ErrorCode result = gMB->create_set_iterator( 0, MBHEX, -1, chunk_size, check_valid, iter );
00450     RC( "query_elem_to_vert_iters" );
00451     bool atend = false;
00452     while( !atend )
00453     {
00454         hexes.clear();
00455         result = iter->get_next_arr( hexes, atend );
00456         RC( "query_elem_to_vert_iters" );
00457         result = gMB->get_connectivity( &hexes[0], hexes.size(), connect );
00458         RC( "query_elem_to_vert_iters" );
00459         result = gMB->get_coords( &connect[0], connect.size(), dum_coords );
00460         RC( "query_elem_to_vert_iters" );
00461         result = gMB->tag_get_data( pos_tag, &hexes[0], hexes.size(), dum_pos );
00462         RC( "query_elem_to_vert_iters" );
00463         for( unsigned int i = 0; i < hexes.size(); i++ )
00464         {
00465             // compute the centroid
00466             for( int j = 0; j < 8; j++ )
00467             {
00468                 dum_pos[3 * i + 0] += dum_coords[24 * i + 3 * j];
00469                 dum_pos[3 * i + 1] += dum_coords[24 * i + 3 * j + 1];
00470                 dum_pos[3 * i + 2] += dum_coords[24 * i + 3 * j + 2];
00471             }
00472             dum_pos[3 * i + 0] *= 0.125;
00473             dum_pos[3 * i + 1] *= 0.125;
00474             dum_pos[3 * i + 2] *= 0.125;
00475         }
00476         result = gMB->tag_set_data( pos_tag, &hexes[0], hexes.size(), dum_pos );
00477         RC( "query_elem_to_vert_iters" );
00478     }
00479 
00480     delete iter;
00481 }
00482 
00483 void query_vert_to_elem_iters( int chunk_size,
00484                                bool check_valid,
00485                                std::vector< EntityHandle >& /*connect*/,
00486                                double* dum_coords,
00487                                double* dum_pos )
00488 {
00489     std::vector< EntityHandle > verts, neighbor_hexes;
00490     SetIterator* iter;
00491     ErrorCode result = gMB->create_set_iterator( 0, MBVERTEX, -1, chunk_size, check_valid, iter );
00492     RC( "query_vert_to_elem_iters" );
00493     assert( MB_SUCCESS == result );
00494     bool atend = false;
00495     while( !atend )
00496     {
00497         verts.clear();
00498         result = iter->get_next_arr( verts, atend );
00499         RC( "query_vert_to_elem_iters" );
00500         result = gMB->get_coords( &verts[0], verts.size(), dum_coords );
00501         RC( "query_vert_to_elem_iters" );
00502         chunk_size = std::min( (int)verts.size(), chunk_size );
00503         for( int i = 0; i < chunk_size; i++ )
00504         {
00505             neighbor_hexes.clear();
00506             result = gMB->get_adjacencies( &verts[i], 1, 3, false, neighbor_hexes );
00507             RC( "query_vert_to_elem_iters" );
00508             result = gMB->tag_get_data( pos2_tag, &neighbor_hexes[0], neighbor_hexes.size(), dum_pos );
00509             RC( "query_vert_to_elem_iters" );
00510             for( unsigned int j = 0; j < neighbor_hexes.size(); j++ )
00511             {
00512                 dum_pos[3 * j + 0] += dum_coords[3 * i + 0];
00513                 dum_pos[3 * j + 1] += dum_coords[3 * i + 1];
00514                 dum_pos[3 * j + 2] += dum_coords[3 * i + 2];
00515             }
00516             result = gMB->tag_set_data( pos2_tag, &neighbor_hexes[0], neighbor_hexes.size(), dum_pos );
00517             RC( "query_vert_to_elem_iters" );
00518         }
00519     }
00520 
00521     result = normalize_elems( dum_coords );
00522     RC( "query_vert_to_elem_iters" );
00523 }
00524 
00525 // normalize pos_tag on all elems by 1/8; coords assumed large enough to hold 3 doubles
00526 ErrorCode normalize_elems( double* coords )
00527 {
00528     // get all hexes and divide pos_tag by 8
00529     Range elems;
00530     ErrorCode result = gMB->get_entities_by_type( 0, MBHEX, elems );
00531     RR( "normalize" );
00532 
00533     for( Range::iterator vit = elems.begin(); vit != elems.end(); ++vit )
00534     {
00535         result = gMB->tag_get_data( pos2_tag, &( *vit ), 1, coords );
00536         RR( "normalize" );
00537         coords[0] *= 0.125;
00538         coords[1] *= 0.125;
00539         coords[2] *= 0.125;
00540         result = gMB->tag_set_data( pos2_tag, &( *vit ), 1, coords );
00541         RR( "normalize" );
00542     }
00543 
00544     return result;
00545 }
00546 
00547 void query_struct_elem_to_vert()
00548 {
00549     // assumes brick mapped mesh with handles starting at zero
00550     Range all_hexes;
00551     ErrorCode result = gMB->get_entities_by_type( 0, MBHEX, all_hexes );
00552     RC( "query_struct_elem_to_vert" );
00553     double dum_coords[24];
00554     std::vector< EntityHandle > connect;
00555     for( Range::iterator eit = all_hexes.begin(); eit != all_hexes.end(); ++eit )
00556     {
00557         result = gMB->get_connectivity( &( *eit ), 1, connect );
00558         RC( "query_struct_elem_to_vert" );
00559         result = gMB->get_coords( &connect[0], connect.size(), dum_coords );
00560         RC( "query_struct_elem_to_vert" );
00561 
00562         double centroid[3] = { 0.0, 0.0, 0.0 };
00563         for( int j = 0; j < 8; j++ )
00564         {
00565             centroid[0] += dum_coords[3 * j + 0];
00566             centroid[1] += dum_coords[3 * j + 1];
00567             centroid[2] += dum_coords[3 * j + 2];
00568         }
00569         centroid[0] *= 0.125;
00570         centroid[1] *= 0.125;
00571         centroid[2] *= 0.125;
00572         result = gMB->tag_set_data( pos_tag, &( *eit ), 1, centroid );
00573         RC( "query_struct_elem_to_vert" );
00574     }
00575 }
00576 
00577 #if defined( _MSC_VER ) || defined( __MINGW32__ )
00578 void print_time( const bool print_em, double& tot_time, double& utime, double& stime, long& imem, long& rmem )
00579 {
00580     utime = (double)clock() / CLOCKS_PER_SEC;
00581     if( print_em ) std::cout << "Total wall time = " << utime << std::endl;
00582     tot_time = stime = 0;
00583     imem = rmem = 0;
00584 }
00585 #else
00586 void print_time( const bool print_em, double& tot_time, double& utime, double& stime, long& imem, long& rmem )
00587 {
00588     struct rusage r_usage;
00589     getrusage( RUSAGE_SELF, &r_usage );
00590     utime = (double)r_usage.ru_utime.tv_sec + ( (double)r_usage.ru_utime.tv_usec / 1.e6 );
00591     stime = (double)r_usage.ru_stime.tv_sec + ( (double)r_usage.ru_stime.tv_usec / 1.e6 );
00592     tot_time = utime + stime;
00593     if( print_em )
00594         std::cout << "User, system, total time = " << utime << ", " << stime << ", " << tot_time << std::endl;
00595 #ifndef LINUX
00596     if( print_em )
00597     {
00598         std::cout << "Max resident set size = " << r_usage.ru_maxrss << " kbytes" << std::endl;
00599         std::cout << "Int resident set size = " << r_usage.ru_idrss << std::endl;
00600     }
00601     imem = r_usage.ru_idrss;
00602     rmem = r_usage.ru_maxrss;
00603 #else
00604     system( "ps o args,drs,rss | grep perf | grep -v grep" );  // RedHat 9.0 doesnt fill in actual
00605                                                                // memory data
00606     imem = rmem = 0;
00607 #endif
00608 }
00609 #endif
00610 
00611 void testA( const int nelem, const double* coords )
00612 {
00613     double ttime0 = 0.0, ttime1 = 0.0, ttime2 = 0.0, ttime3 = 0.0, ttime4 = 0.0, utime = 0.0, stime = 0.0;
00614     long imem0 = 0, rmem0 = 0, imem1 = 0, rmem1 = 0, imem2 = 0, rmem2 = 0, imem3 = 0, rmem3 = 0, imem4 = 0, rmem4 = 0;
00615 
00616     print_time( false, ttime0, utime, stime, imem0, rmem0 );
00617 
00618     // make a 3d block of vertices
00619     EntitySequence* dum_seq    = NULL;
00620     ScdVertexData* vseq        = NULL;
00621     StructuredElementSeq* eseq = NULL;
00622     init();
00623     Core* mbcore = dynamic_cast< Core* >( gMB );
00624     assert( mbcore != NULL );
00625     SequenceManager* seq_mgr = mbcore->sequence_manager();
00626     HomCoord vseq_minmax[2]  = { HomCoord( 0, 0, 0 ), HomCoord( nelem, nelem, nelem ) };
00627     EntityHandle vstart, estart;
00628 
00629     ErrorCode result = seq_mgr->create_scd_sequence( vseq_minmax[0], vseq_minmax[1], MBVERTEX, 1, vstart, dum_seq );
00630     RC( "testA" );
00631     if( NULL != dum_seq ) vseq = dynamic_cast< ScdVertexData* >( dum_seq->data() );
00632     assert( MB_FAILURE != result && vstart != 0 && dum_seq != NULL && vseq != NULL );
00633     // now the element sequence
00634     result = seq_mgr->create_scd_sequence( vseq_minmax[0], vseq_minmax[1], MBHEX, 1, estart, dum_seq );
00635     if( NULL != dum_seq ) eseq = dynamic_cast< StructuredElementSeq* >( dum_seq );
00636     assert( MB_FAILURE != result && estart != 0 && dum_seq != NULL && eseq != NULL );
00637 
00638     // only need to add one vseq to this, unity transform
00639     // trick: if I know it's going to be unity, just input 3 sets of equivalent points
00640     result = eseq->sdata()->add_vsequence( vseq, vseq_minmax[0], vseq_minmax[0], vseq_minmax[0], vseq_minmax[0],
00641                                            vseq_minmax[0], vseq_minmax[0] );
00642     assert( MB_SUCCESS == result );
00643 
00644     // set the coordinates of the vertices
00645     EntityHandle handle;
00646     int i;
00647     double dumv[3];
00648     int num_verts = ( nelem + 1 ) * ( nelem + 1 ) * ( nelem + 1 );
00649     for( i = 0, handle = vstart; i < num_verts; i++, handle++ )
00650     {
00651         dumv[0] = coords[i];
00652         dumv[1] = coords[num_verts + i];
00653         dumv[2] = coords[2 * num_verts + i];
00654         result  = gMB->set_coords( &handle, 1, dumv );
00655         assert( MB_SUCCESS == result );
00656     }
00657 
00658     print_time( false, ttime1, utime, stime, imem1, rmem1 );
00659 
00660     // query the mesh 2 ways
00661     query_struct_elem_to_vert();
00662 
00663     print_time( false, ttime2, utime, stime, imem2, rmem2 );
00664 
00665     query_vert_to_elem();
00666 
00667     print_time( false, ttime3, utime, stime, imem3, rmem3 );
00668 
00669 #ifndef NDEBUG
00670     check_answers( "A" );
00671 #endif
00672 
00673     delete gMB;
00674 
00675     print_time( false, ttime4, utime, stime, imem4, rmem4 );
00676 
00677     std::cout << "MOAB_scd:nelem,construct,e_to_v,v_to_e,after_dtor,total= " << nelem << " " << ttime1 - ttime0 << " "
00678               << ttime2 - ttime1 << " " << ttime3 - ttime2 << " " << ttime4 - ttime3 << " " << ttime4 - ttime0
00679               << " seconds" << std::endl;
00680     std::cout << "MOAB_scd_memory(rss):initial,after_construction,e-v,v-e,after_dtor= " << rmem0 << " " << rmem1 << " "
00681               << rmem2 << " " << rmem3 << " " << rmem4 << " kb" << std::endl;
00682 }
00683 
00684 void testB( const int nelem, const double* coords, int* connect )
00685 {
00686     double ttime0 = 0.0, ttime1 = 0.0, ttime2 = 0.0, ttime3 = 0.0, ttime4 = 0.0, utime = 0.0, stime = 0.0;
00687     long imem0 = 0, rmem0 = 0, imem1 = 0, rmem1 = 0, imem2 = 0, rmem2 = 0, imem3 = 0, rmem3 = 0, imem4 = 0, rmem4 = 0;
00688 
00689     print_time( false, ttime0, utime, stime, imem0, rmem0 );
00690 
00691     int num_verts = ( nelem + 1 ) * ( nelem + 1 ) * ( nelem + 1 );
00692     int num_elems = nelem * nelem * nelem;
00693     EntityHandle vstart, estart;
00694 
00695     // get the read interface
00696     ReadUtilIface* readMeshIface;
00697     init();
00698     gMB->query_interface( readMeshIface );
00699 
00700     // create a sequence to hold the node coordinates
00701     // get the current number of entities and start at the next slot
00702     std::vector< double* > coord_arrays;
00703     ErrorCode result = readMeshIface->get_node_coords( 3, num_verts, 1, vstart, coord_arrays );
00704     RC( "testB" );
00705     assert( MB_SUCCESS == result && 1 == vstart && coord_arrays[0] && coord_arrays[1] && coord_arrays[2] );
00706     // memcpy the coordinate data into place
00707     memcpy( coord_arrays[0], coords, sizeof( double ) * num_verts );
00708     memcpy( coord_arrays[1], &coords[num_verts], sizeof( double ) * num_verts );
00709     memcpy( coord_arrays[2], &coords[2 * num_verts], sizeof( double ) * num_verts );
00710 
00711     EntityHandle* conn = 0;
00712     result             = readMeshIface->get_element_connect( num_elems, 8, MBHEX, 1, estart, conn );
00713     assert( MB_SUCCESS == result );
00714     for( int i = 0; i < num_elems * 8; i++ )
00715         conn[i] = vstart + connect[i];
00716 
00717     delete[] connect;
00718 
00719     result = readMeshIface->update_adjacencies( estart, num_elems, 8, conn );
00720     assert( MB_SUCCESS == result );
00721 
00722     print_time( false, ttime1, utime, stime, imem1, rmem1 );
00723 
00724     // query the mesh 2 ways
00725     query_elem_to_vert();
00726 
00727     print_time( false, ttime2, utime, stime, imem2, rmem2 );
00728 
00729     query_vert_to_elem();
00730 
00731     print_time( false, ttime3, utime, stime, imem3, rmem3 );
00732 
00733 #ifndef NDEBUG
00734     check_answers( "B" );
00735 #endif
00736 
00737     delete gMB;
00738 
00739     print_time( false, ttime4, utime, stime, imem4, rmem4 );
00740 
00741     std::cout << "MOAB_ucd_blocked:nelem,construct,e_to_v,v_to_e,after_dtor,total= " << nelem << " " << ttime1 - ttime0
00742               << " " << ttime2 - ttime1 << " " << ttime3 - ttime2 << " " << ttime4 - ttime3 << " " << ttime4 - ttime0
00743               << " seconds" << std::endl;
00744     std::cout << "MOAB_ucdblocked_memory_(rss):initial,after_construction,e-v,v-e,after_dtor= " << rmem0 << " " << rmem1
00745               << " " << rmem2 << " " << rmem3 << " " << rmem4 << " kb" << std::endl;
00746 }
00747 
00748 void testC( const int nelem, const double* coords )
00749 {
00750     double ttime0 = 0.0, ttime1 = 0.0, ttime2 = 0.0, ttime3 = 0.0, ttime4 = 0.0, utime = 0.0, stime = 0.0;
00751     long imem0 = 0, rmem0 = 0, imem1 = 0, rmem1 = 0, imem2 = 0, rmem2 = 0, imem3 = 0, rmem3 = 0, imem4 = 0, rmem4 = 0;
00752 
00753     print_time( false, ttime0, utime, stime, imem0, rmem0 );
00754 
00755     // create the vertices; assume we don't need to keep a list of vertex handles, since they'll
00756     // be created in sequence
00757     int numv             = nelem + 1;
00758     int numv_sq          = numv * numv;
00759     int num_verts        = numv * numv * numv;
00760     double dum_coords[3] = { coords[0], coords[num_verts], coords[2 * num_verts] };
00761     EntityHandle vstart;
00762 
00763     init();
00764     ErrorCode result = gMB->create_vertex( dum_coords, vstart );
00765     RC( "testC" );
00766     assert( MB_SUCCESS == result && 1 == vstart );
00767 
00768     EntityHandle dum_vert, vijk;
00769     int i;
00770     for( i = 1; i < num_verts; i++ )
00771     {
00772         dum_coords[0] = coords[i];
00773         dum_coords[1] = coords[num_verts + i];
00774         dum_coords[2] = coords[2 * num_verts + i];
00775         result        = gMB->create_vertex( dum_coords, dum_vert );
00776         assert( MB_SUCCESS == result );
00777     }
00778 
00779     EntityHandle dum_conn[8];
00780     for( i = 0; i < nelem; i++ )
00781     {
00782         for( int j = 0; j < nelem; j++ )
00783         {
00784             for( int k = 0; k < nelem; k++ )
00785             {
00786                 vijk        = vstart + VINDEX( i, j, k );
00787                 dum_conn[0] = vijk;
00788                 dum_conn[1] = vijk + 1;
00789                 dum_conn[2] = vijk + 1 + numv;
00790                 dum_conn[3] = vijk + numv;
00791                 dum_conn[4] = vijk + numv * numv;
00792                 dum_conn[5] = vijk + 1 + numv * numv;
00793                 dum_conn[6] = vijk + 1 + numv + numv * numv;
00794                 dum_conn[7] = vijk + numv + numv * numv;
00795                 result      = gMB->create_element( MBHEX, dum_conn, 8, dum_vert );
00796                 assert( MB_SUCCESS == result );
00797             }
00798         }
00799     }
00800 
00801     print_time( false, ttime1, utime, stime, imem1, rmem1 );
00802 
00803     // query the mesh 2 ways
00804     query_elem_to_vert();
00805 
00806     print_time( false, ttime2, utime, stime, imem2, rmem2 );
00807 
00808     query_vert_to_elem();
00809 
00810     print_time( false, ttime3, utime, stime, imem3, rmem3 );
00811 
00812 #ifndef NDEBUG
00813     check_answers( "C" );
00814 #endif
00815 
00816     delete gMB;
00817 
00818     print_time( false, ttime4, utime, stime, imem4, rmem4 );
00819 
00820     std::cout << "MOAB_ucd_indiv:nelem,construct,e_to_v,v_to_e,after_dtor,total= " << nelem << " " << ttime1 - ttime0
00821               << " " << ttime2 - ttime1 << " " << ttime3 - ttime2 << " " << ttime4 - ttime3 << " " << ttime4 - ttime0
00822               << " seconds" << std::endl;
00823     std::cout << "MOAB_ucd_indiv_memory_(rss):initial,after_construction,e-v,v-e,after_dtor= " << rmem0 << " " << rmem1
00824               << " " << rmem2 << " " << rmem3 << " " << rmem4 << " kb" << std::endl;
00825 }
00826 
00827 void testD( const int nelem, const double* coords, int ver )
00828 {
00829     double ttime0 = 0.0, ttime1 = 0.0, ttime2 = 0.0, ttime3 = 0.0, ttime4 = 0.0, ttime5 = 0.0, ttime6 = 0.0,
00830            ttime7 = 0.0, ttime8 = 0.0, ttime9 = 0.0, ttime10 = 0.0, utime = 0.0, stime = 0.0;
00831     long imem0 = 0, rmem0 = 0, imem1 = 0, rmem1 = 0, imem2 = 0, rmem2 = 0, imem3 = 0, rmem3 = 0, imem4 = 0, rmem4 = 0,
00832          imem5 = 0, rmem5 = 0, imem6 = 0, rmem6 = 0, imem7 = 0, rmem7 = 0, imem8 = 0, rmem8 = 0, imem9 = 0, rmem9 = 0,
00833          imem10 = 0, rmem10 = 0;
00834 
00835     print_time( false, ttime0, utime, stime, imem0, rmem0 );
00836 
00837     // create the vertices; assume we don't need to keep a list of vertex handles, since they'll
00838     // be created in sequence
00839     int numv      = nelem + 1;
00840     int numv_sq   = numv * numv;
00841     int num_verts = numv * numv * numv;
00842     std::vector< double > dum_coords( 24 ), dum_pos( 24 );
00843     dum_coords[0] = coords[0];
00844     dum_coords[1] = coords[num_verts];
00845     dum_coords[2] = coords[2 * num_verts];
00846     EntityHandle vstart;
00847 
00848     init();
00849     ErrorCode result = gMB->create_vertex( &dum_coords[0], vstart );
00850     RC( "testD" );
00851     assert( MB_SUCCESS == result && 1 == vstart );
00852 
00853     EntityHandle dum_vert, vijk;
00854     int i;
00855     for( i = 1; i < num_verts; i++ )
00856     {
00857         dum_coords[0] = coords[i];
00858         dum_coords[1] = coords[num_verts + i];
00859         dum_coords[2] = coords[2 * num_verts + i];
00860         result        = gMB->create_vertex( &dum_coords[0], dum_vert );
00861         assert( MB_SUCCESS == result );
00862     }
00863 
00864     EntityHandle dum_conn[8];
00865     for( i = 0; i < nelem; i++ )
00866     {
00867         for( int j = 0; j < nelem; j++ )
00868         {
00869             for( int k = 0; k < nelem; k++ )
00870             {
00871                 vijk        = vstart + VINDEX( i, j, k );
00872                 dum_conn[0] = vijk;
00873                 dum_conn[1] = vijk + 1;
00874                 dum_conn[2] = vijk + 1 + numv;
00875                 dum_conn[3] = vijk + numv;
00876                 dum_conn[4] = vijk + numv * numv;
00877                 dum_conn[5] = vijk + 1 + numv * numv;
00878                 dum_conn[6] = vijk + 1 + numv + numv * numv;
00879                 dum_conn[7] = vijk + numv + numv * numv;
00880                 result      = gMB->create_element( MBHEX, dum_conn, 8, dum_vert );
00881                 assert( MB_SUCCESS == result );
00882             }
00883         }
00884     }
00885 
00886     print_time( false, ttime1, utime, stime, imem1, rmem1 );
00887 
00888     // query the mesh 2 ways with !check_valid
00889     std::vector< EntityHandle > connect( 8 );
00890 #ifndef NDEBUG
00891     // used only in debug mode
00892     double def_val[3] = { 0.0, 0.0, 0.0 };
00893 #endif
00894     if( ver == 0 || ver == 1 )
00895     {
00896         query_elem_to_vert_iters( 1, false, connect, &dum_coords[0], &dum_pos[0] );
00897         print_time( false, ttime2, utime, stime, imem2, rmem2 );
00898         query_vert_to_elem_iters( 1, false, connect, &dum_coords[0], &dum_pos[0] );
00899         print_time( false, ttime3, utime, stime, imem3, rmem3 );
00900 #ifndef NDEBUG
00901         check_answers( "D" );
00902         result = gMB->tag_delete( pos_tag );
00903         assert( MB_SUCCESS == result );
00904         result =
00905             gMB->tag_get_handle( "position_tag", 3, MB_TYPE_DOUBLE, pos_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00906         assert( MB_SUCCESS == result );
00907         result = gMB->tag_delete( pos2_tag );
00908         assert( MB_SUCCESS == result );
00909         result =
00910             gMB->tag_get_handle( "position2_tag", 3, MB_TYPE_DOUBLE, pos2_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00911         assert( MB_SUCCESS == result );
00912 #endif
00913     }
00914 
00915     if( ver == 0 || ver == 2 )
00916     {
00917         if( ver != 0 ) print_time( false, ttime3, utime, stime, imem3, rmem3 );
00918         query_elem_to_vert_iters( 1, true, connect, &dum_coords[0], &dum_pos[0] );
00919         print_time( false, ttime4, utime, stime, imem4, rmem4 );
00920         query_vert_to_elem_iters( 1, true, connect, &dum_coords[0], &dum_pos[0] );
00921         print_time( false, ttime5, utime, stime, imem5, rmem5 );
00922 #ifndef NDEBUG
00923         check_answers( "D" );
00924         result = gMB->tag_delete( pos_tag );
00925         assert( MB_SUCCESS == result );
00926         result =
00927             gMB->tag_get_handle( "position_tag", 3, MB_TYPE_DOUBLE, pos_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00928         assert( MB_SUCCESS == result );
00929         result = gMB->tag_delete( pos2_tag );
00930         assert( MB_SUCCESS == result );
00931         result =
00932             gMB->tag_get_handle( "position2_tag", 3, MB_TYPE_DOUBLE, pos2_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00933         assert( MB_SUCCESS == result );
00934 #endif
00935     }
00936 
00937     if( ver == 0 || ver >= 3 )
00938     {
00939         dum_coords.resize( 2400 );
00940         dum_pos.resize( 300 );
00941     }
00942     if( ver == 0 || ver == 3 )
00943     {
00944         if( ver != 0 ) print_time( false, ttime5, utime, stime, imem3, rmem3 );
00945         query_elem_to_vert_iters( 100, false, connect, &dum_coords[0], &dum_pos[0] );
00946         print_time( false, ttime6, utime, stime, imem6, rmem6 );
00947         query_vert_to_elem_iters( 100, false, connect, &dum_coords[0], &dum_pos[0] );
00948         print_time( false, ttime7, utime, stime, imem7, rmem7 );
00949 #ifndef NDEBUG
00950         check_answers( "D" );
00951         result = gMB->tag_delete( pos_tag );
00952         assert( MB_SUCCESS == result );
00953         result =
00954             gMB->tag_get_handle( "position_tag", 3, MB_TYPE_DOUBLE, pos_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00955         assert( MB_SUCCESS == result );
00956         result = gMB->tag_delete( pos2_tag );
00957         assert( MB_SUCCESS == result );
00958         result =
00959             gMB->tag_get_handle( "position2_tag", 3, MB_TYPE_DOUBLE, pos2_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00960         assert( MB_SUCCESS == result );
00961 #endif
00962     }
00963 
00964     if( ver == 0 || ver == 4 )
00965     {
00966         if( ver != 0 ) print_time( false, ttime7, utime, stime, imem3, rmem3 );
00967         query_elem_to_vert_iters( 100, true, connect, &dum_coords[0], &dum_pos[0] );
00968         print_time( false, ttime8, utime, stime, imem8, rmem8 );
00969         query_vert_to_elem_iters( 100, true, connect, &dum_coords[0], &dum_pos[0] );
00970         print_time( false, ttime9, utime, stime, imem9, rmem9 );
00971 #ifndef NDEBUG
00972         check_answers( "D" );
00973         result = gMB->tag_delete( pos_tag );
00974         assert( MB_SUCCESS == result );
00975         result =
00976             gMB->tag_get_handle( "position_tag", 3, MB_TYPE_DOUBLE, pos_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00977         assert( MB_SUCCESS == result );
00978         result = gMB->tag_delete( pos2_tag );
00979         assert( MB_SUCCESS == result );
00980         result =
00981             gMB->tag_get_handle( "position2_tag", 3, MB_TYPE_DOUBLE, pos2_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
00982         assert( MB_SUCCESS == result );
00983 #endif
00984     }
00985 
00986     if( ver > 0 && ver < 4 ) print_time( false, ttime9, utime, stime, imem9, rmem9 );
00987     delete gMB;
00988 
00989     print_time( false, ttime10, utime, stime, imem10, rmem10 );
00990 
00991     if( ver == 0 || ver == 1 )
00992     {
00993         std::cout << "MOAB_ucd_iters_!check_valid_1:nelem,construct,e-v,v-e,after_dtor,total= " << nelem << " "
00994                   << ttime1 - ttime0 << " " << ttime2 - ttime1 << " " << ttime3 - ttime2 << " " << ttime10 - ttime9
00995                   << " " << ttime3 - ttime0 + ttime10 - ttime9 << std::endl;
00996         std::cout << "MOAB_ucd_iters_memory_(rss)_!check_valid_1:initial,after_construction,e-v,v-"
00997                      "e,after_dtor= "
00998                   << rmem0 << " " << rmem1 << " " << rmem2 << " " << rmem3 << " " << rmem10 << " kb" << std::endl;
00999     }
01000     if( ver == 0 || ver == 2 )
01001     {
01002         std::cout << "MOAB_ucd_iters_check_valid_1:nelem,construct,e-v,v-e,after_dtor,total= " << nelem << " "
01003                   << ttime1 - ttime0 << " " << ttime4 - ttime3 << " " << ttime5 - ttime4 << " " << ttime10 - ttime9
01004                   << " " << ttime1 - ttime0 + ttime5 - ttime3 + ttime10 - ttime9 << std::endl;
01005         std::cout << "MOAB_ucd_iters_memory_(rss)_check_valid_1:initial,after_construction,e-v,v-e,"
01006                      "after_dtor= "
01007                   << rmem0 << " " << rmem1 << " " << rmem2 << " " << rmem3 << " " << rmem10 << " kb" << std::endl;
01008     }
01009     if( ver == 0 || ver == 3 )
01010     {
01011         std::cout << "MOAB_ucd_iters_!check_valid_100:nelem,construct,e-v,v-e,after_dtor,total= " << nelem << " "
01012                   << ttime1 - ttime0 << " " << ttime6 - ttime5 << " " << ttime7 - ttime6 << " " << ttime10 - ttime9
01013                   << " " << ttime1 - ttime0 + ttime7 - ttime5 + ttime10 - ttime9 << std::endl;
01014         std::cout << "MOAB_ucd_iters_memory_(rss)_!check_valid_100:initial,after_construction,e-v,"
01015                      "v-e,after_dtor= "
01016                   << rmem0 << " " << rmem1 << " " << rmem6 << " " << rmem7 << " " << rmem10 << " kb" << std::endl;
01017     }
01018     if( ver == 0 || ver == 4 )
01019     {
01020         std::cout << "MOAB_ucd_iters_check_valid_100:nelem,construct,e-v,v-e,after_dtor,total= " << nelem << " "
01021                   << ttime1 - ttime0 << " " << ttime8 - ttime7 << " " << ttime9 - ttime8 << " " << ttime10 - ttime9
01022                   << " " << ttime1 - ttime0 + ttime10 - ttime7 << std::endl;
01023         std::cout << "MOAB_ucd_iters_memory_(rss)_check_valid_100:initial,after_construction,e-v,v-"
01024                      "e,after_dtor= "
01025                   << rmem0 << " " << rmem1 << " " << rmem8 << " " << rmem9 << " " << rmem10 << " kb" << std::endl;
01026     }
01027 }
01028 
01029 void testE( const int nelem, const double* coords, int* connect )
01030 {
01031     double ttime0 = 0.0, ttime1 = 0.0, ttime2 = 0.0, ttime3 = 0.0, ttime4 = 0.0, ttime5 = 0.0, ttime6 = 0.0,
01032            utime = 0.0, stime = 0.0;
01033     long imem0 = 0, rmem0 = 0, imem1 = 0, rmem1 = 0, imem2 = 0, rmem2 = 0, imem3 = 0, rmem3 = 0, imem4 = 0, rmem4 = 0,
01034          imem5 = 0, rmem5 = 0, imem6 = 0, rmem6 = 0;
01035 
01036     print_time( false, ttime0, utime, stime, imem0, rmem0 );
01037 
01038     int num_verts = ( nelem + 1 ) * ( nelem + 1 ) * ( nelem + 1 );
01039     int num_elems = nelem * nelem * nelem;
01040     EntityHandle vstart, estart;
01041 
01042     // get the read interface
01043     ReadUtilIface* readMeshIface;
01044     init();
01045     gMB->query_interface( readMeshIface );
01046 
01047     // create a sequence to hold the node coordinates
01048     // get the current number of entities and start at the next slot
01049     std::vector< double* > coord_arrays;
01050     ErrorCode result = readMeshIface->get_node_coords( 3, num_verts, 1, vstart, coord_arrays );
01051     RC( "testE" );
01052     assert( MB_SUCCESS == result && 1 == vstart && coord_arrays[0] && coord_arrays[1] && coord_arrays[2] );
01053     // memcpy the coordinate data into place
01054     memcpy( coord_arrays[0], coords, sizeof( double ) * num_verts );
01055     memcpy( coord_arrays[1], &coords[num_verts], sizeof( double ) * num_verts );
01056     memcpy( coord_arrays[2], &coords[2 * num_verts], sizeof( double ) * num_verts );
01057 
01058     EntityHandle* conn = 0;
01059     result             = readMeshIface->get_element_connect( num_elems, 8, MBHEX, 1, estart, conn );
01060     assert( MB_SUCCESS == result );
01061     for( int i = 0; i < num_elems * 8; i++ )
01062         conn[i] = vstart + connect[i];
01063 
01064     delete[] connect;
01065 
01066     Range verts( vstart, vstart + num_verts - 1 ), elems( estart, estart + num_elems - 1 );
01067 
01068     result = readMeshIface->update_adjacencies( estart, num_elems, 8, conn );
01069     assert( MB_SUCCESS == result );
01070 
01071     print_time( false, ttime1, utime, stime, imem1, rmem1 );
01072 
01073     // query the mesh 2 ways
01074     query_elem_to_vert_direct();
01075 
01076     print_time( false, ttime2, utime, stime, imem2, rmem2 );
01077 
01078     query_vert_to_elem_direct();
01079 
01080     print_time( false, ttime3, utime, stime, imem3, rmem3 );
01081 
01082 #ifndef NDEBUG
01083     check_answers( "E" );
01084 #endif
01085 
01086     query_elem_to_vert_direct();
01087 
01088     print_time( false, ttime4, utime, stime, imem4, rmem4 );
01089 
01090     query_vert_to_elem_direct();
01091 
01092     print_time( false, ttime5, utime, stime, imem5, rmem5 );
01093 
01094 #ifndef NDEBUG
01095     check_answers( "E" );
01096 #endif
01097 
01098     delete gMB;
01099 
01100     print_time( false, ttime6, utime, stime, imem6, rmem6 );
01101 
01102     std::cout << "MOAB_ucd_direct:nelem,construct,e_to_v,v_to_e,after_dtor,total= " << nelem << " " << ttime1 - ttime0
01103               << " " << ttime2 - ttime1 << " " << ttime3 - ttime2 << " " << ttime6 - ttime5 << " "
01104               << ttime3 - ttime0 + ttime6 - ttime5 << " seconds" << std::endl;
01105     std::cout << "MOAB_ucd_direct_memory_(rss):initial,after_construction,e-v,v-e,after_dtor= " << rmem0 << " " << rmem1
01106               << " " << rmem2 << " " << rmem3 << " " << rmem6 << " kb" << std::endl;
01107 
01108     std::cout << "MOAB_ucd_direct2:nelem,construct,e_to_v,v_to_e,after_dtor,total= " << nelem << " " << ttime1 - ttime0
01109               << " " << ttime4 - ttime3 << " " << ttime5 - ttime4 << " " << ttime6 - ttime5 << " "
01110               << ttime1 - ttime0 + ttime6 - ttime3 << " seconds" << std::endl;
01111     std::cout << "MOAB_ucd_direct2_memory_(rss):initial,after_construction,e-v,v-e,after_dtor= " << rmem0 << " "
01112               << rmem1 << " " << rmem4 << " " << rmem5 << " " << rmem6 << " kb" << std::endl;
01113 }
01114 
01115 void query_elem_to_vert_direct()
01116 {
01117     Range all_hexes, all_verts;
01118     ErrorCode result = gMB->get_entities_by_type( 0, MBHEX, all_hexes );
01119     assert( MB_SUCCESS == result );
01120     result = gMB->get_entities_by_type( 0, MBVERTEX, all_verts );
01121     RC( "query_elem_to_vert_direct" );
01122     EntityHandle* connect;
01123     int ecount, vcount, vpere;
01124     double* coords[3];
01125 
01126     result = gMB->connect_iterate( all_hexes.begin(), all_hexes.end(), connect, vpere, ecount );
01127     if( MB_SUCCESS != result || ecount != (int)all_hexes.size() )
01128     {
01129         std::cout << "FAILED in connect_iterate!" << std::endl;
01130         return;
01131     }
01132     result = gMB->coords_iterate( all_verts.begin(), all_verts.end(), coords[0], coords[1], coords[2], vcount );
01133     if( MB_SUCCESS != result || vcount != (int)all_verts.size() )
01134     {
01135         std::cout << "FAILED in coords_iterate!" << std::endl;
01136         return;
01137     }
01138     double* centroid;
01139     result = gMB->tag_iterate( pos_tag, all_hexes.begin(), all_hexes.end(), ecount, (void*&)centroid );
01140     if( MB_SUCCESS != result || ecount != (int)all_hexes.size() )
01141     {
01142         std::cout << "FAILED in connect_iterate!" << std::endl;
01143         return;
01144     }
01145 
01146     EntityHandle vstart = *all_verts.begin();
01147     for( int i = 0; i < ecount; i++ )
01148     {
01149         // compute the centroid
01150         for( int j = 0; j < vpere; j++ )
01151         {
01152             int vind = *connect - vstart;
01153             connect++;
01154             centroid[3 * i + 0] += coords[0][vind];
01155             centroid[3 * i + 1] += coords[1][vind];
01156             centroid[3 * i + 2] += coords[2][vind];
01157         }
01158 
01159         // now normalize
01160         centroid[3 * i + 0] *= 0.125;
01161         centroid[3 * i + 1] *= 0.125;
01162         centroid[3 * i + 2] *= 0.125;
01163     }
01164 }
01165 
01166 void query_vert_to_elem_direct()
01167 {
01168     Range all_verts, tmp_ents;
01169     ErrorCode result = gMB->get_entities_by_type( 0, MBVERTEX, all_verts );
01170     assert( MB_SUCCESS == result );
01171 
01172     // make sure vertex-element adjacencies are created
01173     result = gMB->get_adjacencies( &( *all_verts.begin() ), 1, 3, false, tmp_ents );
01174     assert( MB_SUCCESS == result );
01175 
01176     const std::vector< EntityHandle >** adjs;
01177     int count;
01178     result = gMB->adjacencies_iterate( all_verts.begin(), all_verts.end(), adjs, count );
01179     if( MB_SUCCESS != result && count != (int)all_verts.size() )
01180     {
01181         std::cout << "FAILED:adjacencies_iterate." << std::endl;
01182         return;
01183     }
01184 
01185     double* coords[3];
01186     result = gMB->coords_iterate( all_verts.begin(), all_verts.end(), coords[0], coords[1], coords[2], count );
01187     if( MB_SUCCESS != result || count != (int)all_verts.size() )
01188     {
01189         std::cout << "FAILED in coords_iterate!" << std::endl;
01190         return;
01191     }
01192     // get all hexes, then iterator over pos2_tag
01193     double* centroid;
01194     int ecount;
01195     tmp_ents.clear();
01196     result = gMB->get_entities_by_type( 0, MBHEX, tmp_ents );
01197     assert( MB_SUCCESS == result );
01198     result = gMB->tag_iterate( pos2_tag, tmp_ents.begin(), tmp_ents.end(), ecount, (void*&)centroid );
01199     if( MB_SUCCESS != result || ecount != (int)tmp_ents.size() )
01200     {
01201         std::cout << "FAILED in tag_iterate!" << std::endl;
01202         return;
01203     }
01204 
01205     int i;
01206     Range::iterator vit;
01207     EntityHandle estart = *tmp_ents.begin();
01208     for( vit = all_verts.begin(), i = 0; vit != all_verts.end(); ++vit, i++ )
01209     {
01210         assert( adjs[i] );
01211         for( std::vector< EntityHandle >::const_iterator vit2 = adjs[i]->begin(); vit2 != adjs[i]->end(); ++vit2 )
01212             if( *vit >= estart )
01213             {
01214                 int eind = *vit2 - estart;
01215                 centroid[3 * eind + 0] += coords[0][i];
01216                 centroid[3 * eind + 1] += coords[1][i];
01217                 centroid[3 * eind + 2] += coords[2][i];
01218             }
01219     }
01220 
01221     // now normalize
01222     for( i = 0; i < (int)tmp_ents.size(); i++ )
01223     {
01224         centroid[3 * i + 0] *= 0.125;
01225         centroid[3 * i + 1] *= 0.125;
01226         centroid[3 * i + 2] *= 0.125;
01227     }
01228 }
01229 
01230 void check_answers( const char* /*test_name*/ )
01231 {
01232     Range elems;
01233     ErrorCode result = gMB->get_entities_by_type( 0, MBHEX, elems );
01234     RC( "check_answers" );
01235 
01236     double coords1[3], coords2[3], del[3];
01237     double diff = 0.0;
01238     for( Range::iterator vit = elems.begin(); vit != elems.end(); ++vit )
01239     {
01240         result = gMB->tag_get_data( pos_tag, &( *vit ), 1, coords1 );
01241         RC( "check_answers" );
01242         result = gMB->tag_get_data( pos2_tag, &( *vit ), 1, coords2 );
01243         RC( "check_answers" );
01244         for( int i = 0; i < 3; i++ )
01245             del[i] = fabs( coords1[i] - coords2[i] );
01246         if( del[0] || del[1] || del[2] )
01247         {
01248             double len2 = std::max( coords1[0] * coords1[0] + coords1[1] * coords1[1] + coords1[2] * coords1[2],
01249                                     coords2[0] * coords2[0] + coords2[1] * coords2[1] + coords2[2] * coords2[2] ),
01250                    num  = del[0] * del[0] + del[1] * del[1] + del[2] * del[2];
01251             if( len2 > 0.0 )
01252                 diff = std::max( diff, num / sqrt( len2 ) );
01253             else if( num > 0.0 )
01254                 diff = sqrt( num );
01255         }
01256     }
01257     if( diff > 0.0 ) std::cout << "Max relative difference = " << diff << std::endl;
01258 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines