MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 }