MOAB: Mesh Oriented datABase  (version 5.4.1)
tstt_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 /** TSTT Mesh Interface brick mesh performance test
00017  *
00018  * This program tests TSTT mesh interface functions used to create a square structured
00019  * mesh.  Boilerplate taken from tstt mesh interface test in MOAB and performance test in MOAB
00020  *
00021  */
00022 
00023 // Different platforms follow different conventions for usage
00024 #if !defined( _NT ) && !defined( _MSC_VER ) && !defined( __MINGW32__ )
00025 #include <sys/resource.h>
00026 #endif
00027 #ifdef SOLARIS
00028 extern "C" int getrusage( int, struct rusage* );
00029 #ifndef RUSAGE_SELF
00030 #include </usr/ucbinclude/sys/rusage.h>
00031 #endif
00032 #endif
00033 
00034 #include <cstdlib>
00035 #include <cstdio>
00036 #include <iostream>
00037 
00038 //-------------------------------------------------------
00039 // CHANGE THESE DEFINITIONS TO SUIT YOUR IMPLEMENTATION
00040 
00041 // #include here any files with declarations needed to instantiate your
00042 // implementation instance
00043 #include "TSTT_MOAB.hh"
00044 
00045 // define IMPLEMENTATION_CLASS to be the namespace::class of your implementation
00046 #define IMPLEMENTATION_CLASS TSTT_MOAB::MoabMesh
00047 //-------------------------------------------------------
00048 
00049 #define CAST_INTERFACE( var_in, var_out, iface )                           \
00050     TSTT::iface var_out;                                                   \
00051     try                                                                    \
00052     {                                                                      \
00053         ( var_out ) = var_in;                                              \
00054     }                                                                      \
00055     catch( TSTT::Error )                                                   \
00056     {                                                                      \
00057         cerr << "Error: current interface doesn't support iface." << endl; \
00058         return;                                                            \
00059     }
00060 
00061 #define CAST_MINTERFACE( var_in, var_out, iface )                          \
00062     TSTTM::iface var_out;                                                  \
00063     try                                                                    \
00064     {                                                                      \
00065         ( var_out ) = var_in;                                              \
00066     }                                                                      \
00067     catch( TSTT::Error )                                                   \
00068     {                                                                      \
00069         cerr << "Error: current interface doesn't support iface." << endl; \
00070         return;                                                            \
00071     }
00072 
00073 #include <iostream>
00074 #include "TSTT.hh"
00075 #include "TSTTM.hh"
00076 
00077 // needed to get the proper size for handles
00078 typedef void* Entity_Handle;
00079 
00080 using namespace std;
00081 
00082 #define ARRAY_PTR( array, type )  ( reinterpret_cast< ( type )* >( ( array )._get_ior()->d_firstElement ) )
00083 #define HANDLE_ARRAY_PTR( array ) ( reinterpret_cast< Entity_Handle* >( ( array )._get_ior()->d_firstElement ) )
00084 #define ARRAY_SIZE( array )       ( ( array )._is_nil() ? 0 : ( array ).upper( 0 ) - ( array ).lower( 0 ) + 1 )
00085 #define CHECK_SIZE( array, size )                                     \
00086     if( ( array )._is_nil() || ARRAY_SIZE( array ) == 0 )             \
00087         ( array ) = ( array ).create1d( size );                       \
00088     else if( ARRAY_SIZE( array ) < ( size ) )                         \
00089     {                                                                 \
00090         cerr << "Array passed in is non-zero but too short." << endl; \
00091         assert( false );                                              \
00092     }
00093 
00094 double LENGTH = 1.0;
00095 
00096 // forward declare some functions
00097 void query_vert_to_elem( TSTTM::Mesh& mesh );
00098 void query_elem_to_vert( TSTTM::Mesh& mesh );
00099 void print_time( const bool print_em, double& tot_time, double& utime, double& stime );
00100 void build_connect( const int nelem, const int vstart, int*& connect );
00101 void testB( TSTTM::Mesh& mesh, const int nelem, sidl::array< double >& coords, const int* connect );
00102 void testC( TSTTM::Mesh& mesh, const int nelem, sidl::array< double >& coords );
00103 void compute_edge( double* start, const int nelem, const double xint, const int stride );
00104 void compute_face( double* a, const int nelem, const double xint, const int stride1, const int stride2 );
00105 void build_coords( const int nelem, sidl::array< double >& coords );
00106 void build_connect( const int nelem, const int vstart, int*& connect );
00107 
00108 int main( int argc, char* argv[] )
00109 {
00110     int nelem = 20;
00111     if( argc < 3 )
00112     {
00113         std::cout << "Usage: " << argv[0] << " <ints_per_side> <A|B|C>" << std::endl;
00114         return 1;
00115     }
00116 
00117     char which_test = '\0';
00118 
00119     sscanf( argv[1], "%d", &nelem );
00120     sscanf( argv[2], "%c", &which_test );
00121 
00122     if( which_test != 'B' && which_test != 'C' )
00123     {
00124         std::cout << "Must indicate B or C for test." << std::endl;
00125         return 1;
00126     }
00127 
00128     std::cout << "number of elements: " << nelem << "; test " << which_test << std::endl;
00129 
00130     // pre-build the coords array
00131     sidl::array< double > coords;
00132     build_coords( nelem, coords );
00133     assert( NULL != coords );
00134 
00135     int* connect = NULL;
00136     build_connect( nelem, 1, connect );
00137 
00138     // create an implementation
00139     TSTTM::Mesh mesh = IMPLEMENTATION_CLASS::_create();
00140 
00141     switch( which_test )
00142     {
00143         case 'B':
00144             // test B: create mesh using bulk interface
00145             testB( mesh, nelem, coords, connect );
00146             break;
00147 
00148         case 'C':
00149             // test C: create mesh using individual interface
00150             testC( mesh, nelem, coords );
00151             break;
00152     }
00153 
00154     return 0;
00155 }
00156 
00157 void testB( TSTTM::Mesh& mesh, const int nelem, sidl::array< double >& coords, const int* connect )
00158 {
00159     double utime, stime, ttime0, ttime1, ttime2, ttime3;
00160 
00161     print_time( false, ttime0, utime, stime );
00162     int num_verts = ( nelem + 1 ) * ( nelem + 1 ) * ( nelem + 1 );
00163     int num_elems = nelem * nelem * nelem;
00164 
00165     // create vertices as a block
00166     CAST_MINTERFACE( mesh, mesh_arrmod, ArrMod );
00167     sidl::array< Entity_Handle > sidl_vertices, dum_handles;
00168     CHECK_SIZE( dum_handles, num_verts );
00169     int sidl_vertices_size;
00170 
00171     try
00172     {
00173         mesh_arrmod.createVtxArr( num_verts, TSTTM::StorageOrder_BLOCKED, coords, 3 * num_verts, sidl_vertices,
00174                                   sidl_vertices_size );
00175     }
00176     catch( TSTT::Error& /* err */ )
00177     {
00178         cerr << "Couldn't create vertices in bulk call" << endl;
00179         return;
00180     }
00181 
00182     // need to explicitly fill connectivity array, since we don't know
00183     // the format of entity handles
00184     sidl::array< Entity_Handle > sidl_connect;
00185     int sidl_connect_size = 8 * num_elems;
00186     CHECK_SIZE( sidl_connect, 8 * num_elems );
00187     Entity_Handle* sidl_connect_ptr = HANDLE_ARRAY_PTR( sidl_connect );
00188 
00189     Entity_Handle* sidl_vertices_ptr = HANDLE_ARRAY_PTR( sidl_vertices );
00190     for( int i = 0; i < sidl_connect_size; i++ )
00191     {
00192         // use connect[i]-1 because we used starting vertex index (vstart) of 1
00193         assert( connect[i] - 1 < num_verts );
00194         sidl_connect_ptr[i] = sidl_vertices_ptr[connect[i] - 1];
00195     }
00196 
00197     // create the entities
00198     sidl::array< Entity_Handle > new_hexes;
00199     int new_hexes_size;
00200     sidl::array< TSTTM::CreationStatus > status;
00201     int status_size;
00202 
00203     try
00204     {
00205         mesh_arrmod.createEntArr( TSTTM::EntityTopology_HEXAHEDRON, sidl_connect, sidl_connect_size, new_hexes,
00206                                   new_hexes_size, status, status_size );
00207     }
00208     catch( TSTT::Error& /* err */ )
00209     {
00210         cerr << "Couldn't create hex elements in bulk call" << endl;
00211         return;
00212     }
00213 
00214     print_time( false, ttime1, utime, stime );
00215 
00216     // query the mesh 2 ways
00217     query_elem_to_vert( mesh );
00218 
00219     print_time( false, ttime2, utime, stime );
00220 
00221     query_vert_to_elem( mesh );
00222 
00223     print_time( false, ttime3, utime, stime );
00224 
00225     std::cout << "TSTT/MOAB ucd blocked: nelem, construct, e_to_v query, v_to_e query = " << nelem << ", "
00226               << ttime1 - ttime0 << ", " << ttime2 - ttime1 << ", " << ttime3 - ttime2 << " seconds" << std::endl;
00227 }
00228 
00229 void testC( TSTTM::Mesh& mesh, const int nelem, sidl::array< double >& coords )
00230 {
00231     double utime, stime, ttime0, ttime1, ttime2, ttime3;
00232     print_time( false, ttime0, utime, stime );
00233 
00234     CAST_MINTERFACE( mesh, mesh_arrmod, ArrMod );
00235 
00236     // need some dimensions
00237     int numv      = nelem + 1;
00238     int numv_sq   = numv * numv;
00239     int num_verts = numv * numv * numv;
00240 #define VINDEX( i, j, k ) ( ( i ) + ( (j)*numv ) + ( (k)*numv_sq ) )
00241 
00242     // array to hold vertices created individually
00243     sidl::array< Entity_Handle > sidl_vertices;
00244     // int sidl_vertices_size = num_verts;
00245     CHECK_SIZE( sidl_vertices, num_verts );
00246 
00247     // temporary array to hold vertex positions for single vertex
00248     sidl::array< double > tmp_coords;
00249     int tmp_coords_size = 3;
00250     CHECK_SIZE( tmp_coords, 3 );
00251     double* dum_coords = ARRAY_PTR( tmp_coords, double );
00252 
00253     // get direct pointer to coordinate array
00254     double* coords_ptr = ARRAY_PTR( coords, double );
00255 
00256     for( int i = 0; i < num_verts; i++ )
00257     {
00258 
00259         // temporary array to hold (single) vertices
00260         sidl::array< Entity_Handle > tmp_vertices;
00261         int tmp_vertices_size = 0;
00262 
00263         // create the vertex
00264         dum_coords[0] = coords_ptr[i];
00265         dum_coords[1] = coords_ptr[num_verts + i];
00266         dum_coords[2] = coords_ptr[2 * num_verts + i];
00267         try
00268         {
00269             mesh_arrmod.createVtxArr( 1, TSTTM::StorageOrder_BLOCKED, tmp_coords, tmp_coords_size, tmp_vertices,
00270                                       tmp_vertices_size );
00271         }
00272         catch( TSTT::Error& /* err */ )
00273         {
00274             cerr << "Couldn't create vertex in individual call" << endl;
00275             return;
00276         }
00277 
00278         // assign into permanent vertex array
00279         sidl_vertices.set( i, tmp_vertices.get( 0 ) );
00280     }
00281 
00282     // get vertex array pointer for reading into tmp_conn
00283     Entity_Handle* tmp_sidl_vertices = HANDLE_ARRAY_PTR( sidl_vertices );
00284 
00285     for( int i = 0; i < nelem; i++ )
00286     {
00287         for( int j = 0; j < nelem; j++ )
00288         {
00289             for( int k = 0; k < nelem; k++ )
00290             {
00291 
00292                 sidl::array< Entity_Handle > tmp_conn;
00293                 int tmp_conn_size = 8;
00294                 CHECK_SIZE( tmp_conn, 8 );
00295 
00296                 int vijk = VINDEX( i, j, k );
00297                 tmp_conn.set( 0, tmp_sidl_vertices[vijk] );
00298                 tmp_conn.set( 1, tmp_sidl_vertices[vijk + 1] );
00299                 tmp_conn.set( 2, tmp_sidl_vertices[vijk + 1 + numv] );
00300                 tmp_conn.set( 3, tmp_sidl_vertices[vijk + numv] );
00301                 tmp_conn.set( 4, tmp_sidl_vertices[vijk + numv * numv] );
00302                 tmp_conn.set( 5, tmp_sidl_vertices[vijk + 1 + numv * numv] );
00303                 tmp_conn.set( 6, tmp_sidl_vertices[vijk + 1 + numv + numv * numv] );
00304                 tmp_conn.set( 7, tmp_sidl_vertices[vijk + numv + numv * numv] );
00305 
00306                 // create the entity
00307 
00308                 sidl::array< Entity_Handle > new_hex;
00309                 int new_hex_size = 0;
00310                 sidl::array< TSTTM::CreationStatus > status;
00311                 int status_size = 0;
00312 
00313                 try
00314                 {
00315                     mesh_arrmod.createEntArr( TSTTM::EntityTopology_HEXAHEDRON, tmp_conn, tmp_conn_size, new_hex,
00316                                               new_hex_size, status, status_size );
00317                 }
00318                 catch( TSTT::Error& /* err */ )
00319                 {
00320                     cerr << "Couldn't create hex element in individual call" << endl;
00321                     return;
00322                 }
00323             }
00324         }
00325     }
00326 
00327     print_time( false, ttime1, utime, stime );
00328 
00329     // query the mesh 2 ways
00330     query_elem_to_vert( mesh );
00331 
00332     print_time( false, ttime2, utime, stime );
00333 
00334     query_vert_to_elem( mesh );
00335 
00336     print_time( false, ttime3, utime, stime );
00337 
00338     std::cout << "TSTT/MOAB ucd indiv: nelem, construct, e_to_v query, v_to_e query = " << nelem << ", "
00339               << ttime1 - ttime0 << ", " << ttime2 - ttime1 << ", " << ttime3 - ttime2 << " seconds" << std::endl;
00340 }
00341 
00342 void query_elem_to_vert( TSTTM::Mesh& mesh )
00343 {
00344     sidl::array< Entity_Handle > all_hexes;
00345     int all_hexes_size;
00346     CAST_MINTERFACE( mesh, mesh_ent, Entity );
00347 
00348     // get all the hex elements
00349     try
00350     {
00351         mesh.getEntities( 0, TSTTM::EntityType_REGION, TSTTM::EntityTopology_HEXAHEDRON, all_hexes, all_hexes_size );
00352     }
00353     catch( TSTT::Error& /* err */ )
00354     {
00355         cerr << "Couldn't get all hex elements in query_mesh" << endl;
00356         return;
00357     }
00358 
00359     try
00360     {
00361         // set up some tmp arrays and array ptrs
00362         Entity_Handle* all_hexes_ptr = HANDLE_ARRAY_PTR( all_hexes );
00363 
00364         // now loop over elements
00365         for( int i = 0; i < all_hexes_size; i++ )
00366         {
00367             sidl::array< int > dum_offsets;
00368             sidl::array< Entity_Handle > dum_connect;
00369             int dum_connect_size = 0;
00370             // get the connectivity of this element; will allocate space on 1st iteration,
00371             // but will have correct size on subsequent ones
00372             mesh_ent.getEntAdj( all_hexes_ptr[i], TSTTM::EntityType_VERTEX, dum_connect, dum_connect_size );
00373 
00374             // get vertex coordinates; ; will allocate space on 1st iteration,
00375             // but will have correct size on subsequent ones
00376             sidl::array< double > dum_coords;
00377             int dum_coords_size       = 0;
00378             TSTTM::StorageOrder order = TSTTM::StorageOrder_UNDETERMINED;
00379             mesh.getVtxArrCoords( dum_connect, dum_connect_size, order, dum_coords, dum_coords_size );
00380 
00381             assert( 24 == dum_coords_size && ARRAY_SIZE( dum_coords ) == 24 );
00382             double* dum_coords_ptr = ARRAY_PTR( dum_coords, double );
00383             double centroid[3]     = { 0.0, 0.0, 0.0 };
00384             if( order == TSTTM::StorageOrder_BLOCKED )
00385             {
00386                 for( int j = 0; j < 8; j++ )
00387                 {
00388                     centroid[0] += dum_coords_ptr[j];
00389                     centroid[1] += dum_coords_ptr[8 + j];
00390                     centroid[2] += dum_coords_ptr[16 + j];
00391                     centroid[0] += dum_coords.get( j );
00392                     centroid[1] += dum_coords.get( 8 + j );
00393                     centroid[2] += dum_coords.get( 16 + j );
00394                 }
00395             }
00396             else
00397             {
00398                 for( int j = 0; j < 8; j++ )
00399                 {
00400                     centroid[0] += dum_coords_ptr[3 * j];
00401                     centroid[1] += dum_coords_ptr[3 * j + 1];
00402                     centroid[2] += dum_coords_ptr[3 * j + 2];
00403                 }
00404             }
00405         }
00406     }
00407     catch( TSTT::Error& /* err */ )
00408     {
00409         cerr << "Problem getting connectivity or vertex coords." << endl;
00410         return;
00411     }
00412 }
00413 
00414 void query_vert_to_elem( TSTTM::Mesh& mesh )
00415 {
00416     sidl::array< Entity_Handle > all_verts;
00417     int all_verts_size;
00418     CAST_MINTERFACE( mesh, mesh_ent, Entity );
00419 
00420     // get all the vertices elements
00421     try
00422     {
00423         mesh.getEntities( 0, TSTTM::EntityType_VERTEX, TSTTM::EntityTopology_POINT, all_verts, all_verts_size );
00424     }
00425     catch( TSTT::Error& /* err */ )
00426     {
00427         cerr << "Couldn't get all vertices in query_vert_to_elem" << endl;
00428         return;
00429     }
00430 
00431     try
00432     {
00433         // set up some tmp arrays and array ptrs
00434         Entity_Handle* all_verts_ptr = HANDLE_ARRAY_PTR( all_verts );
00435 
00436         // now loop over vertices
00437         for( int i = 0; i < all_verts_size; i++ )
00438         {
00439             sidl::array< Entity_Handle > dum_hexes;
00440             int dum_hexes_size;
00441 
00442             // get the connectivity of this element; will have to allocate space on every
00443             // iteration, since size can vary
00444             mesh_ent.getEntAdj( all_verts_ptr[i], TSTTM::EntityType_REGION, dum_hexes, dum_hexes_size );
00445         }
00446     }
00447     catch( TSTT::Error& /* err */ )
00448     {
00449         cerr << "Problem getting connectivity or vertex coords." << endl;
00450         return;
00451     }
00452 }
00453 
00454 void print_time( const bool print_em, double& tot_time, double& utime, double& stime )
00455 {
00456     struct rusage r_usage;
00457     getrusage( RUSAGE_SELF, &r_usage );
00458     utime    = (double)r_usage.ru_utime.tv_sec + ( (double)r_usage.ru_utime.tv_usec / 1.e6 );
00459     stime    = (double)r_usage.ru_stime.tv_sec + ( (double)r_usage.ru_stime.tv_usec / 1.e6 );
00460     tot_time = utime + stime;
00461     if( print_em )
00462         std::cout << "User, system, total time = " << utime << ", " << stime << ", " << tot_time << std::endl;
00463 #ifndef LINUX
00464     std::cout << "Max resident set size = " << r_usage.ru_maxrss * 4096 << " bytes" << std::endl;
00465     std::cout << "Int resident set size = " << r_usage.ru_idrss << std::endl;
00466 #else
00467     system( "ps o args,drs,rss | grep perf | grep -v grep" );  // RedHat 9.0 doesnt fill in actual
00468                                                                // memory data
00469 #endif
00470 }
00471 
00472 void compute_edge( double* start, const int nelem, const double xint, const int stride )
00473 {
00474     for( int i = 1; i < nelem; i++ )
00475     {
00476         start[i * stride]                     = start[0] + i * xint;
00477         start[nelem + 1 + i * stride]         = start[nelem + 1] + i * xint;
00478         start[2 * ( nelem + 1 ) + i * stride] = start[2 * ( nelem + 1 )] + i * xint;
00479     }
00480 }
00481 
00482 void compute_face( double* a, const int nelem, const double xint, const int stride1, const int stride2 )
00483 {
00484     // 2D TFI on a face starting at a, with strides stride1 in ada and stride2 in tse
00485     for( int j = 1; j < nelem; j++ )
00486     {
00487         double tse = j * xint;
00488         for( int i = 1; i < nelem; i++ )
00489         {
00490             double ada = i * xint;
00491 
00492             a[i * stride1 + j * stride2] =
00493                 ( 1.0 - ada ) * a[i * stride1] + ada * a[i * stride1 + nelem * stride2] +
00494                 ( 1.0 - tse ) * a[j * stride2] + tse * a[j * stride2 + nelem * stride1] -
00495                 ( 1.0 - tse ) * ( 1.0 - ada ) * a[0] - ( 1.0 - tse ) * ada * a[nelem * stride1] -
00496                 tse * ( 1.0 - ada ) * a[nelem * stride2] - tse * ada * a[nelem * ( stride1 + stride2 )];
00497             a[nelem + 1 + i * stride1 + j * stride2] =
00498                 ( 1.0 - ada ) * a[nelem + 1 + i * stride1] + ada * a[nelem + 1 + i * stride1 + nelem * stride2] +
00499                 ( 1.0 - tse ) * a[nelem + 1 + j * stride2] + tse * a[nelem + 1 + j * stride2 + nelem * stride1] -
00500                 ( 1.0 - tse ) * ( 1.0 - ada ) * a[nelem + 1 + 0] -
00501                 ( 1.0 - tse ) * ada * a[nelem + 1 + nelem * stride1] -
00502                 tse * ( 1.0 - ada ) * a[nelem + 1 + nelem * stride2] -
00503                 tse * ada * a[nelem + 1 + nelem * ( stride1 + stride2 )];
00504             a[2 * ( nelem + 1 ) + i * stride1 + j * stride2] =
00505                 ( 1.0 - ada ) * a[2 * ( nelem + 1 ) + i * stride1] +
00506                 ada * a[2 * ( nelem + 1 ) + i * stride1 + nelem * stride2] +
00507                 ( 1.0 - tse ) * a[2 * ( nelem + 1 ) + j * stride2] +
00508                 tse * a[2 * ( nelem + 1 ) + j * stride2 + nelem * stride1] -
00509                 ( 1.0 - tse ) * ( 1.0 - ada ) * a[2 * ( nelem + 1 ) + 0] -
00510                 ( 1.0 - tse ) * ada * a[2 * ( nelem + 1 ) + nelem * stride1] -
00511                 tse * ( 1.0 - ada ) * a[2 * ( nelem + 1 ) + nelem * stride2] -
00512                 tse * ada * a[2 * ( nelem + 1 ) + nelem * ( stride1 + stride2 )];
00513         }
00514     }
00515 }
00516 
00517 void build_coords( const int nelem, sidl::array< double >& coords )
00518 {
00519     double ttime0, ttime1, utime1, stime1;
00520     print_time( false, ttime0, utime1, stime1 );
00521 
00522     // allocate the memory
00523     int numv     = nelem + 1;
00524     int numv_sq  = numv * numv;
00525     int tot_numv = numv * numv * numv;
00526     CHECK_SIZE( coords, 3 * tot_numv );
00527     double* coords_ptr = ARRAY_PTR( coords, double );
00528 
00529 // use FORTRAN-like indexing
00530 #define VINDEX( i, j, k ) ( ( i ) + ( (j)*numv ) + ( (k)*numv_sq ) )
00531     int idx;
00532     double scale1, scale2, scale3;
00533     // use these to prevent optimization on 1-scale, etc (real map wouldn't have
00534     // all these equal)
00535     scale1 = LENGTH / nelem;
00536     scale2 = LENGTH / nelem;
00537     scale3 = LENGTH / nelem;
00538 
00539 #ifdef REALTFI
00540     // use a real TFI xform to compute coordinates
00541     // initialize positions of 8 corners
00542     coords_ptr[VINDEX( 0, 0, 0 )]                         = coords_ptr[VINDEX( 0, 0, 0 ) + nelem + 1] =
00543         coords_ptr[VINDEX( 0, 0, 0 ) + 2 * ( nelem + 1 )] = 0.0;
00544     coords_ptr[VINDEX( 0, nelem, 0 )] = coords_ptr[VINDEX( 0, nelem, 0 ) + 2 * ( nelem + 1 )] = 0.0;
00545     coords_ptr[VINDEX( 0, nelem, 0 ) + nelem + 1]                                             = LENGTH;
00546     coords_ptr[VINDEX( 0, 0, nelem )] = coords_ptr[VINDEX( 0, 0, nelem ) + nelem + 1] = 0.0;
00547     coords_ptr[VINDEX( 0, 0, nelem ) + 2 * ( nelem + 1 )]                             = LENGTH;
00548     coords_ptr[VINDEX( 0, nelem, nelem )]                                             = 0.0;
00549     coords_ptr[VINDEX( 0, nelem, nelem ) + nelem + 1] = coords_ptr[VINDEX( 0, nelem, nelem ) + 2 * ( nelem + 1 )] =
00550         LENGTH;
00551     coords_ptr[VINDEX( nelem, 0, 0 )]             = LENGTH;
00552     coords_ptr[VINDEX( nelem, 0, 0 ) + nelem + 1] = coords_ptr[VINDEX( nelem, 0, 0 ) + 2 * ( nelem + 1 )] = 0.0;
00553     coords_ptr[VINDEX( nelem, 0, nelem )] = coords_ptr[VINDEX( nelem, 0, nelem ) + 2 * ( nelem + 1 )] = LENGTH;
00554     coords_ptr[VINDEX( nelem, 0, nelem ) + nelem + 1]                                                 = 0.0;
00555     coords_ptr[VINDEX( nelem, nelem, 0 )] = coords_ptr[VINDEX( nelem, nelem, 0 ) + nelem + 1] = LENGTH;
00556     coords_ptr[VINDEX( nelem, nelem, 0 ) + 2 * ( nelem + 1 )]                                 = 0.0;
00557     coords_ptr[VINDEX( nelem, nelem, nelem )] = coords_ptr[VINDEX( nelem, nelem, nelem ) + nelem + 1] =
00558         coords_ptr[VINDEX( nelem, nelem, nelem ) + 2 * ( nelem + 1 )] = LENGTH;
00559 
00560     // compute edges
00561     // i (stride=1)
00562     compute_edge( &coords_ptr[VINDEX( 0, 0, 0 )], nelem, scale1, 1 );
00563     compute_edge( &coords_ptr[VINDEX( 0, nelem, 0 )], nelem, scale1, 1 );
00564     compute_edge( &coords_ptr[VINDEX( 0, 0, nelem )], nelem, scale1, 1 );
00565     compute_edge( &coords_ptr[VINDEX( 0, nelem, nelem )], nelem, scale1, 1 );
00566     // j (stride=numv)
00567     compute_edge( &coords_ptr[VINDEX( 0, 0, 0 )], nelem, scale1, numv );
00568     compute_edge( &coords_ptr[VINDEX( nelem, 0, 0 )], nelem, scale1, numv );
00569     compute_edge( &coords_ptr[VINDEX( 0, 0, nelem )], nelem, scale1, numv );
00570     compute_edge( &coords_ptr[VINDEX( nelem, 0, nelem )], nelem, scale1, numv );
00571     // k (stride=numv^2)
00572     compute_edge( &coords_ptr[VINDEX( 0, 0, 0 )], nelem, scale1, numv_sq );
00573     compute_edge( &coords_ptr[VINDEX( nelem, 0, 0 )], nelem, scale1, numv_sq );
00574     compute_edge( &coords_ptr[VINDEX( 0, nelem, 0 )], nelem, scale1, numv_sq );
00575     compute_edge( &coords_ptr[VINDEX( nelem, nelem, 0 )], nelem, scale1, numv_sq );
00576 
00577     // compute faces
00578     // i=0, nelem
00579     compute_face( &coords_ptr[VINDEX( 0, 0, 0 )], nelem, scale1, numv, numv_sq );
00580     compute_face( &coords_ptr[VINDEX( nelem, 0, 0 )], nelem, scale1, numv, numv_sq );
00581     // j=0, nelem
00582     compute_face( &coords_ptr[VINDEX( 0, 0, 0 )], nelem, scale1, 1, numv_sq );
00583     compute_face( &coords_ptr[VINDEX( 0, nelem, 0 )], nelem, scale1, 1, numv_sq );
00584     // k=0, nelem
00585     compute_face( &coords_ptr[VINDEX( 0, 0, 0 )], nelem, scale1, 1, numv );
00586     compute_face( &coords_ptr[VINDEX( 0, 0, nelem )], nelem, scale1, 1, numv );
00587 
00588     // initialize corner indices
00589     int i000 = VINDEX( 0, 0, 0 );
00590     int ia00 = VINDEX( nelem, 0, 0 );
00591     int i0t0 = VINDEX( 0, nelem, 0 );
00592     int iat0 = VINDEX( nelem, nelem, 0 );
00593     int i00g = VINDEX( 0, 0, nelem );
00594     int ia0g = VINDEX( nelem, 0, nelem );
00595     int i0tg = VINDEX( 0, nelem, nelem );
00596     int iatg = VINDEX( nelem, nelem, nelem );
00597     double cX, cY, cZ;
00598     int adaInts   = nelem;
00599     int tseInts   = nelem;
00600     int gammaInts = nelem;
00601 
00602     for( int i = 1; i < nelem; i++ )
00603     {
00604         for( int j = 1; j < nelem; j++ )
00605         {
00606             for( int k = 1; k < nelem; k++ )
00607             {
00608                 // idx = VINDEX(i,j,k);
00609                 double tse   = i * scale1;
00610                 double ada   = j * scale2;
00611                 double gamma = k * scale3;
00612                 double tm1   = 1.0 - tse;
00613                 double am1   = 1.0 - ada;
00614                 double gm1   = 1.0 - gamma;
00615 
00616                 cX = gm1 * ( am1 * ( tm1 * coords_ptr[i000] + tse * coords_ptr[i0t0] ) +
00617                              ada * ( tm1 * coords_ptr[ia00] + tse * coords_ptr[iat0] ) ) +
00618                      gamma * ( am1 * ( tm1 * coords_ptr[i00g] + tse * coords_ptr[i0tg] ) +
00619                                ada * ( tm1 * coords_ptr[ia0g] + tse * coords_ptr[iatg] ) );
00620 
00621                 cY = gm1 * ( am1 * ( tm1 * coords_ptr[i000] + tse * coords_ptr[i0t0] ) +
00622                              ada * ( tm1 * coords_ptr[ia00] + tse * coords_ptr[iat0] ) ) +
00623                      gamma * ( am1 * ( tm1 * coords_ptr[i00g] + tse * coords_ptr[i0tg] ) +
00624                                ada * ( tm1 * coords_ptr[ia0g] + tse * coords_ptr[iatg] ) );
00625 
00626                 cZ = gm1 * ( am1 * ( tm1 * coords_ptr[i000] + tse * coords_ptr[i0t0] ) +
00627                              ada * ( tm1 * coords_ptr[ia00] + tse * coords_ptr[iat0] ) ) +
00628                      gamma * ( am1 * ( tm1 * coords_ptr[i00g] + tse * coords_ptr[i0tg] ) +
00629                                ada * ( tm1 * coords_ptr[ia0g] + tse * coords_ptr[iatg] ) );
00630 
00631                 double* ai0k = &coords_ptr[VINDEX( k, 0, i )];
00632                 double* aiak = &coords_ptr[VINDEX( k, adaInts, i )];
00633                 double* a0jk = &coords_ptr[VINDEX( k, j, 0 )];
00634                 double* atjk = &coords_ptr[VINDEX( k, j, tseInts )];
00635                 double* aij0 = &coords_ptr[VINDEX( 0, j, i )];
00636                 double* aijg = &coords_ptr[VINDEX( gammaInts, j, i )];
00637 
00638                 coords_ptr[VINDEX( i, j, k )] = ( am1 * ai0k[0] + ada * aiak[0] + tm1 * a0jk[0] + tse * atjk[0] +
00639                                                   gm1 * aij0[0] + gamma * aijg[0] ) /
00640                                                     2.0 -
00641                                                 cX / 2.0;
00642 
00643                 coords_ptr[nelem + 1 + VINDEX( i, j, k )] =
00644                     ( am1 * ai0k[nelem + 1] + ada * aiak[nelem + 1] + tm1 * a0jk[nelem + 1] + tse * atjk[nelem + 1] +
00645                       gm1 * aij0[nelem + 1] + gamma * aijg[nelem + 1] ) /
00646                         2.0 -
00647                     cY / 2.0;
00648 
00649                 coords_ptr[2 * ( nelem + 1 ) + VINDEX( i, j, k )] =
00650                     ( am1 * ai0k[2 * ( nelem + 1 )] + ada * aiak[2 * ( nelem + 1 )] + tm1 * a0jk[2 * ( nelem + 1 )] +
00651                       tse * atjk[2 * ( nelem + 1 )] + gm1 * aij0[2 * ( nelem + 1 )] +
00652                       gamma * aijg[2 * ( nelem + 1 )] ) /
00653                         2.0 -
00654                     cZ / 2.0;
00655             }
00656         }
00657     }
00658 
00659 #else
00660     for( int i = 0; i < numv; i++ )
00661     {
00662         for( int j = 0; j < numv; j++ )
00663         {
00664             for( int k = 0; k < numv; k++ )
00665             {
00666                 idx = VINDEX( i, j, k );
00667                 // blocked coordinate ordering
00668                 coords_ptr[idx]                = i * scale1;
00669                 coords_ptr[tot_numv + idx]     = j * scale2;
00670                 coords_ptr[2 * tot_numv + idx] = k * scale3;
00671             }
00672         }
00673     }
00674 #endif
00675     print_time( false, ttime1, utime1, stime1 );
00676     std::cout << "TSTT/MOAB: TFI time = " << ttime1 - ttime0 << " sec" << std::endl;
00677 }
00678 
00679 void build_connect( const int nelem, const int vstart, int*& connect )
00680 {
00681     // allocate the memory
00682     int nume_tot = nelem * nelem * nelem;
00683     connect      = new int[8 * nume_tot];
00684 
00685     int vijk;
00686     int numv    = nelem + 1;
00687     int numv_sq = numv * numv;
00688     int idx     = 0;
00689     for( int i = 0; i < nelem; i++ )
00690     {
00691         for( int j = 0; j < nelem; j++ )
00692         {
00693             for( int k = 0; k < nelem; k++ )
00694             {
00695                 vijk           = vstart + VINDEX( i, j, k );
00696                 connect[idx++] = vijk;
00697                 connect[idx++] = vijk + 1;
00698                 connect[idx++] = vijk + 1 + numv;
00699                 connect[idx++] = vijk + numv;
00700                 connect[idx++] = vijk + numv * numv;
00701                 connect[idx++] = vijk + 1 + numv * numv;
00702                 connect[idx++] = vijk + 1 + numv + numv * numv;
00703                 connect[idx++] = vijk + numv + numv * numv;
00704                 assert( i <= numv * numv * numv );
00705             }
00706         }
00707     }
00708 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines