MOAB: Mesh Oriented datABase  (version 5.5.0)
perf.cpp
Go to the documentation of this file.
1 /**
2  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3  * storing and accessing finite element mesh data.
4  *
5  * Copyright 2004 Sandia Corporation. Under the terms of Contract
6  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7  * retains certain rights in this software.
8  *
9  * This library is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public
11  * License as published by the Free Software Foundation; either
12  * version 2.1 of the License, or (at your option) any later version.
13  *
14  */
15 
16 // MOAB performance tests building mapped mesh with nodes and
17 // hexes created one at a time. This also creates the node to hex adjacencies.
18 
19 // Different platforms follow different conventions for usage
20 #if !defined( _MSC_VER ) && !defined( __MINGW32__ )
21 #include <sys/resource.h>
22 #else
23 #include <time.h>
24 #endif
25 #ifdef SOLARIS
26 extern "C" int getrusage( int, struct rusage* );
27 #ifndef RUSAGE_SELF
28 #include </usr/ucbinclude/sys/rusage.h>
29 #endif
30 #endif
31 
32 #ifndef IS_BUILDING_MB
33 #define IS_BUILDING_MB
34 #endif
35 
36 #include <cstdlib>
37 #include <cstdio>
38 #include <cassert>
39 #include <iostream>
40 #include "moab/Core.hpp"
41 #include "moab/ReadUtilIface.hpp"
42 #include "VertexSequence.hpp"
43 #include "StructuredElementSeq.hpp"
44 #include "EntitySequence.hpp"
45 #include "SequenceManager.hpp"
46 #include "moab/HomXform.hpp"
47 #include "moab/SetIterator.hpp"
48 
49 using namespace moab;
50 
51 double LENGTH = 1.0;
52 
53 void testA( const int nelem, const double* coords );
54 void testB( const int nelem, const double* coords, int* connect );
55 void testC( const int nelem, const double* coords );
56 void testD( const int nelem, const double* coords, int ver );
57 void testE( const int nelem, const double* coords, int* connect );
58 void print_time( const bool print_em, double& tot_time, double& utime, double& stime, long& imem, long& rmem );
59 void query_vert_to_elem();
60 void query_elem_to_vert();
61 void query_struct_elem_to_vert();
62 void query_elem_to_vert_direct();
63 void query_vert_to_elem_direct();
64 ErrorCode normalize_elems( double* coords );
65 void check_answers( const char* );
66 
67 #define RC( msg ) \
68  if( MB_SUCCESS != result ) do \
69  { \
70  std::cout << "FAIL in " << ( msg ) << std::endl; \
71  return; \
72  } while( true )
73 #define RR( msg ) \
74  if( MB_SUCCESS != result ) do \
75  { \
76  std::cout << "FAIL in " << ( msg ) << std::endl; \
77  return result; \
78  } while( true )
79 
80 void compute_edge( double* start, const int nelem, const double xint, const int stride )
81 {
82  for( int i = 1; i < nelem; i++ )
83  {
84  start[i * stride] = start[0] + i * xint;
85  start[nelem + 1 + i * stride] = start[nelem + 1] + i * xint;
86  start[2 * ( nelem + 1 ) + i * stride] = start[2 * ( nelem + 1 )] + i * xint;
87  }
88 }
89 
90 void compute_face( double* a, const int nelem, const double xint, const int stride1, const int stride2 )
91 {
92  // 2D TFI on a face starting at a, with strides stride1 in ada and stride2 in tse
93  for( int j = 1; j < nelem; j++ )
94  {
95  double tse = j * xint;
96  for( int i = 1; i < nelem; i++ )
97  {
98  double ada = i * xint;
99 
100  a[i * stride1 + j * stride2] =
101  ( 1.0 - ada ) * a[i * stride1] + ada * a[i * stride1 + nelem * stride2] +
102  ( 1.0 - tse ) * a[j * stride2] + tse * a[j * stride2 + nelem * stride1] -
103  ( 1.0 - tse ) * ( 1.0 - ada ) * a[0] - ( 1.0 - tse ) * ada * a[nelem * stride1] -
104  tse * ( 1.0 - ada ) * a[nelem * stride2] - tse * ada * a[nelem * ( stride1 + stride2 )];
105  a[nelem + 1 + i * stride1 + j * stride2] =
106  ( 1.0 - ada ) * a[nelem + 1 + i * stride1] + ada * a[nelem + 1 + i * stride1 + nelem * stride2] +
107  ( 1.0 - tse ) * a[nelem + 1 + j * stride2] + tse * a[nelem + 1 + j * stride2 + nelem * stride1] -
108  ( 1.0 - tse ) * ( 1.0 - ada ) * a[nelem + 1 + 0] -
109  ( 1.0 - tse ) * ada * a[nelem + 1 + nelem * stride1] -
110  tse * ( 1.0 - ada ) * a[nelem + 1 + nelem * stride2] -
111  tse * ada * a[nelem + 1 + nelem * ( stride1 + stride2 )];
112  a[2 * ( nelem + 1 ) + i * stride1 + j * stride2] =
113  ( 1.0 - ada ) * a[2 * ( nelem + 1 ) + i * stride1] +
114  ada * a[2 * ( nelem + 1 ) + i * stride1 + nelem * stride2] +
115  ( 1.0 - tse ) * a[2 * ( nelem + 1 ) + j * stride2] +
116  tse * a[2 * ( nelem + 1 ) + j * stride2 + nelem * stride1] -
117  ( 1.0 - tse ) * ( 1.0 - ada ) * a[2 * ( nelem + 1 ) + 0] -
118  ( 1.0 - tse ) * ada * a[2 * ( nelem + 1 ) + nelem * stride1] -
119  tse * ( 1.0 - ada ) * a[2 * ( nelem + 1 ) + nelem * stride2] -
120  tse * ada * a[2 * ( nelem + 1 ) + nelem * ( stride1 + stride2 )];
121  }
122  }
123 }
124 
125 void build_coords( const int nelem, double*& coords )
126 {
127  double ttime0 = 0.0, ttime1 = 0.0, utime1 = 0.0, stime1 = 0.0;
128  long imem = 0, rmem = 0;
129  print_time( false, ttime0, utime1, stime1, imem, rmem );
130  // allocate the memory
131  int numv = nelem + 1;
132  int numv_sq = numv * numv;
133  int tot_numv = numv * numv * numv;
134  coords = new double[3 * tot_numv];
135 
136 // use FORTRAN-like indexing
137 #define VINDEX( i, j, k ) ( ( i ) + ( (j)*numv ) + ( (k)*numv_sq ) )
138  int idx;
139  double scale1, scale2, scale3;
140  // use these to prevent optimization on 1-scale, etc (real map wouldn't have
141  // all these equal)
142  scale1 = LENGTH / nelem;
143  scale2 = LENGTH / nelem;
144  scale3 = LENGTH / nelem;
145 
146 #ifdef REALTFI
147  // use a real TFI xform to compute coordinates
148  // compute edges
149  // i (stride=1)
150  compute_edge( &coords[VINDEX( 0, 0, 0 )], nelem, scale1, 1 );
151  compute_edge( &coords[VINDEX( 0, nelem, 0 )], nelem, scale1, 1 );
152  compute_edge( &coords[VINDEX( 0, 0, nelem )], nelem, scale1, 1 );
153  compute_edge( &coords[VINDEX( 0, nelem, nelem )], nelem, scale1, 1 );
154  // j (stride=numv)
155  compute_edge( &coords[VINDEX( 0, 0, 0 )], nelem, scale1, numv );
156  compute_edge( &coords[VINDEX( nelem, 0, 0 )], nelem, scale1, numv );
157  compute_edge( &coords[VINDEX( 0, 0, nelem )], nelem, scale1, numv );
158  compute_edge( &coords[VINDEX( nelem, 0, nelem )], nelem, scale1, numv );
159  // k (stride=numv^2)
160  compute_edge( &coords[VINDEX( 0, 0, 0 )], nelem, scale1, numv_sq );
161  compute_edge( &coords[VINDEX( nelem, 0, 0 )], nelem, scale1, numv_sq );
162  compute_edge( &coords[VINDEX( 0, nelem, 0 )], nelem, scale1, numv_sq );
163  compute_edge( &coords[VINDEX( nelem, nelem, 0 )], nelem, scale1, numv_sq );
164 
165  // compute faces
166  // i=0, nelem
167  compute_face( &coords[VINDEX( 0, 0, 0 )], nelem, scale1, numv, numv_sq );
168  compute_face( &coords[VINDEX( nelem, 0, 0 )], nelem, scale1, numv, numv_sq );
169  // j=0, nelem
170  compute_face( &coords[VINDEX( 0, 0, 0 )], nelem, scale1, 1, numv_sq );
171  compute_face( &coords[VINDEX( 0, nelem, 0 )], nelem, scale1, 1, numv_sq );
172  // k=0, nelem
173  compute_face( &coords[VINDEX( 0, 0, 0 )], nelem, scale1, 1, numv );
174  compute_face( &coords[VINDEX( 0, 0, nelem )], nelem, scale1, 1, numv );
175 
176  // initialize corner indices
177  int i000 = VINDEX( 0, 0, 0 );
178  int ia00 = VINDEX( nelem, 0, 0 );
179  int i0t0 = VINDEX( 0, nelem, 0 );
180  int iat0 = VINDEX( nelem, nelem, 0 );
181  int i00g = VINDEX( 0, 0, nelem );
182  int ia0g = VINDEX( nelem, 0, nelem );
183  int i0tg = VINDEX( 0, nelem, nelem );
184  int iatg = VINDEX( nelem, nelem, nelem );
185  double cX, cY, cZ;
186  int adaInts = nelem;
187  int tseInts = nelem;
188  int gammaInts = nelem;
189 
190  for( int i = 1; i < nelem; i++ )
191  {
192  for( int j = 1; j < nelem; j++ )
193  {
194  for( int k = 1; k < nelem; k++ )
195  {
196  // idx = VINDEX(i,j,k);
197  double tse = i * scale1;
198  double ada = j * scale2;
199  double gamma = k * scale3;
200  double tm1 = 1.0 - tse;
201  double am1 = 1.0 - ada;
202  double gm1 = 1.0 - gamma;
203 
204  cX = gm1 * ( am1 * ( tm1 * coords[i000] + tse * coords[i0t0] ) +
205  ada * ( tm1 * coords[ia00] + tse * coords[iat0] ) ) +
206  gamma * ( am1 * ( tm1 * coords[i00g] + tse * coords[i0tg] ) +
207  ada * ( tm1 * coords[ia0g] + tse * coords[iatg] ) );
208 
209  cY = gm1 * ( am1 * ( tm1 * coords[i000] + tse * coords[i0t0] ) +
210  ada * ( tm1 * coords[ia00] + tse * coords[iat0] ) ) +
211  gamma * ( am1 * ( tm1 * coords[i00g] + tse * coords[i0tg] ) +
212  ada * ( tm1 * coords[ia0g] + tse * coords[iatg] ) );
213 
214  cZ = gm1 * ( am1 * ( tm1 * coords[i000] + tse * coords[i0t0] ) +
215  ada * ( tm1 * coords[ia00] + tse * coords[iat0] ) ) +
216  gamma * ( am1 * ( tm1 * coords[i00g] + tse * coords[i0tg] ) +
217  ada * ( tm1 * coords[ia0g] + tse * coords[iatg] ) );
218 
219  double* ai0k = &coords[VINDEX( k, 0, i )];
220  double* aiak = &coords[VINDEX( k, adaInts, i )];
221  double* a0jk = &coords[VINDEX( k, j, 0 )];
222  double* atjk = &coords[VINDEX( k, j, tseInts )];
223  double* aij0 = &coords[VINDEX( 0, j, i )];
224  double* aijg = &coords[VINDEX( gammaInts, j, i )];
225 
226  coords[VINDEX( i, j, k )] = ( am1 * ai0k[0] + ada * aiak[0] + tm1 * a0jk[0] + tse * atjk[0] +
227  gm1 * aij0[0] + gamma * aijg[0] ) /
228  2.0 -
229  cX / 2.0;
230 
231  coords[nelem + 1 + VINDEX( i, j, k )] =
232  ( am1 * ai0k[nelem + 1] + ada * aiak[nelem + 1] + tm1 * a0jk[nelem + 1] + tse * atjk[nelem + 1] +
233  gm1 * aij0[nelem + 1] + gamma * aijg[nelem + 1] ) /
234  2.0 -
235  cY / 2.0;
236 
237  coords[2 * ( nelem + 1 ) + VINDEX( i, j, k )] =
238  ( am1 * ai0k[2 * ( nelem + 1 )] + ada * aiak[2 * ( nelem + 1 )] + tm1 * a0jk[2 * ( nelem + 1 )] +
239  tse * atjk[2 * ( nelem + 1 )] + gm1 * aij0[2 * ( nelem + 1 )] +
240  gamma * aijg[2 * ( nelem + 1 )] ) /
241  2.0 -
242  cZ / 2.0;
243  }
244  }
245  }
246 
247 #else
248  for( int i = 0; i < numv; i++ )
249  {
250  for( int j = 0; j < numv; j++ )
251  {
252  for( int k = 0; k < numv; k++ )
253  {
254  idx = VINDEX( i, j, k );
255  // blocked coordinate ordering
256  coords[idx] = i * scale1;
257  coords[tot_numv + idx] = j * scale2;
258  coords[2 * tot_numv + idx] = k * scale3;
259  }
260  }
261  }
262 #endif
263  print_time( false, ttime1, utime1, stime1, imem, rmem );
264  // std::cout << "MOAB: TFI time = " << ttime1-ttime0 << " sec"
265  // << std::endl;
266 }
267 
268 void build_connect( const int nelem, const EntityHandle /*vstart*/, int*& connect )
269 {
270  // allocate the memory
271  int nume_tot = nelem * nelem * nelem;
272  connect = new int[8 * nume_tot];
273 
274  EntityHandle vijk;
275  int numv = nelem + 1;
276  int numv_sq = numv * numv;
277  int idx = 0;
278  for( int i = 0; i < nelem; i++ )
279  {
280  for( int j = 0; j < nelem; j++ )
281  {
282  for( int k = 0; k < nelem; k++ )
283  {
284  vijk = VINDEX( i, j, k );
285  connect[idx++] = vijk;
286  connect[idx++] = vijk + 1;
287  connect[idx++] = vijk + 1 + numv;
288  connect[idx++] = vijk + numv;
289  connect[idx++] = vijk + numv * numv;
290  connect[idx++] = vijk + 1 + numv * numv;
291  connect[idx++] = vijk + 1 + numv + numv * numv;
292  connect[idx++] = vijk + numv + numv * numv;
293  assert( i <= numv * numv * numv );
294  }
295  }
296  }
297 }
298 
299 Interface* gMB;
300 Tag pos_tag, pos2_tag;
301 
302 void init()
303 {
304  gMB = new Core();
305  double def_val[3] = { 0.0, 0.0, 0.0 };
306  ErrorCode rval =
307  gMB->tag_get_handle( "position_tag", 3, MB_TYPE_DOUBLE, pos_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
308  assert( MB_SUCCESS == rval );
309  if( rval )
310  {
311  } // empty line to remove compiler warning
312  rval = gMB->tag_get_handle( "position2_tag", 3, MB_TYPE_DOUBLE, pos2_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
313  assert( MB_SUCCESS == rval );
314 }
315 
316 int main( int argc, char* argv[] )
317 {
318  int nelem = 20;
319  if( argc < 3 )
320  {
321  std::cout << "Usage: " << argv[0] << " <ints_per_side> [A|B|C|D [1|2|3|4]|E]" << std::endl;
322  return 1;
323  }
324 
325  char which_test = '\0';
326  int ver = 0;
327 
328  sscanf( argv[1], "%d", &nelem );
329  if( argc >= 3 ) sscanf( argv[2], "%c", &which_test );
330  if( argc >= 4 ) sscanf( argv[3], "%d", &ver );
331 
332  if( 3 <= argc && which_test != 'A' && which_test != 'B' && which_test != 'C' && which_test != 'D' &&
333  which_test != 'E' )
334  {
335  std::cout << "Must indicate A or B, C, D or E for test." << std::endl;
336  return 1;
337  }
338 
339  if( 4 <= argc && which_test == 'D' && ( ver < 1 || ver > 4 ) )
340  {
341  std::cout << "Must indicate version 1, 2, 3, or 4 for test D." << std::endl;
342  return 1;
343  }
344 
345  // std::cout << "number of elements: " << nelem << "; test "
346  // << which_test << std::endl;
347 
348  // pre-build the coords array
349  double* coords = NULL;
350  build_coords( nelem, coords );
351  assert( NULL != coords );
352 
353  int* connect = NULL;
354 
355  // test A: create structured mesh
356  if( '\0' == which_test || 'A' == which_test ) testA( nelem, coords );
357 
358  build_connect( nelem, 1, connect );
359 
360  // test B: create mesh using bulk interface
361  if( '\0' == which_test || 'B' == which_test ) testB( nelem, coords, connect );
362 
363  // test C: create mesh using individual interface
364  if( '\0' == which_test || 'C' == which_test ) testC( nelem, coords );
365 
366  // test D: query mesh using iterators
367  if( '\0' == which_test || 'D' == which_test ) testD( nelem, coords, ver );
368 
369  // test E: query mesh using direct access
370  if( '\0' == which_test || 'E' == which_test ) testE( nelem, coords, connect );
371 
372  return 0;
373 }
374 
375 void query_elem_to_vert()
376 {
377  Range all_hexes;
378  ErrorCode result = gMB->get_entities_by_type( 0, MBHEX, all_hexes );
379  RC( "query_elem_to_vert" );
380  const EntityHandle* connect;
381  int num_connect;
382  double dum_coords[24];
383  for( Range::iterator eit = all_hexes.begin(); eit != all_hexes.end(); ++eit )
384  {
385  result = gMB->get_connectivity( *eit, connect, num_connect );
386  RC( "query_elem_to_vert" );
387  result = gMB->get_coords( connect, num_connect, dum_coords );
388  RC( "query_elem_to_vert" );
389 
390  // compute the centroid
391  double centroid[3] = { 0.0, 0.0, 0.0 };
392  for( int j = 0; j < 8; j++ )
393  {
394  centroid[0] += dum_coords[3 * j + 0];
395  centroid[1] += dum_coords[3 * j + 1];
396  centroid[2] += dum_coords[3 * j + 2];
397  }
398  centroid[0] *= 0.125;
399  centroid[1] *= 0.125;
400  centroid[2] *= 0.125;
401  result = gMB->tag_set_data( pos_tag, &( *eit ), 1, centroid );
402  RC( "query_elem_to_vert" );
403  }
404 }
405 
406 void query_vert_to_elem()
407 {
408  Range all_verts;
409  std::vector< EntityHandle > neighbor_hexes;
410  std::vector< double > neighbor_pos;
411  double coords[3];
412  neighbor_pos.resize( 3 * 8 ); // average vertex will have 8 adjacent hexes
413  ErrorCode result = gMB->get_entities_by_type( 0, MBVERTEX, all_verts );
414  RC( "query_vert_to_elem" );
415  for( Range::iterator vit = all_verts.begin(); vit != all_verts.end(); ++vit )
416  {
417  neighbor_hexes.clear();
418  result = gMB->get_coords( &( *vit ), 1, coords );
419  RC( "query_vert_to_elem" );
420  result = gMB->get_adjacencies( &( *vit ), 1, 3, false, neighbor_hexes );
421  RC( "query_vert_to_elem" );
422  assert( neighbor_pos.size() >= 3 * neighbor_hexes.size() );
423  result = gMB->tag_get_data( pos2_tag, &neighbor_hexes[0], neighbor_hexes.size(), &neighbor_pos[0] );
424  RC( "query_vert_to_elem" );
425  for( unsigned int i = 0; i < neighbor_hexes.size(); i++ )
426  {
427  neighbor_pos[3 * i] += coords[0];
428  neighbor_pos[3 * i + 1] += coords[1];
429  neighbor_pos[3 * i + 2] += coords[2];
430  }
431 
432  result = gMB->tag_set_data( pos2_tag, &neighbor_hexes[0], neighbor_hexes.size(), &neighbor_pos[0] );
433  RC( "query_vert_to_elem" );
434  }
435 
436  // get all hexes and divide pos_tag by 8; reuse all_verts
437  result = normalize_elems( coords );
438  RC( "query_vert_to_elem" );
439 }
440 
441 void query_elem_to_vert_iters( int chunk_size,
442  bool check_valid,
443  std::vector< EntityHandle >& connect,
444  double* dum_coords,
445  double* dum_pos )
446 {
447  std::vector< EntityHandle > hexes;
448  SetIterator* iter;
449  ErrorCode result = gMB->create_set_iterator( 0, MBHEX, -1, chunk_size, check_valid, iter );
450  RC( "query_elem_to_vert_iters" );
451  bool atend = false;
452  while( !atend )
453  {
454  hexes.clear();
455  result = iter->get_next_arr( hexes, atend );
456  RC( "query_elem_to_vert_iters" );
457  result = gMB->get_connectivity( &hexes[0], hexes.size(), connect );
458  RC( "query_elem_to_vert_iters" );
459  result = gMB->get_coords( &connect[0], connect.size(), dum_coords );
460  RC( "query_elem_to_vert_iters" );
461  result = gMB->tag_get_data( pos_tag, &hexes[0], hexes.size(), dum_pos );
462  RC( "query_elem_to_vert_iters" );
463  for( unsigned int i = 0; i < hexes.size(); i++ )
464  {
465  // compute the centroid
466  for( int j = 0; j < 8; j++ )
467  {
468  dum_pos[3 * i + 0] += dum_coords[24 * i + 3 * j];
469  dum_pos[3 * i + 1] += dum_coords[24 * i + 3 * j + 1];
470  dum_pos[3 * i + 2] += dum_coords[24 * i + 3 * j + 2];
471  }
472  dum_pos[3 * i + 0] *= 0.125;
473  dum_pos[3 * i + 1] *= 0.125;
474  dum_pos[3 * i + 2] *= 0.125;
475  }
476  result = gMB->tag_set_data( pos_tag, &hexes[0], hexes.size(), dum_pos );
477  RC( "query_elem_to_vert_iters" );
478  }
479 
480  delete iter;
481 }
482 
483 void query_vert_to_elem_iters( int chunk_size,
484  bool check_valid,
485  std::vector< EntityHandle >& /*connect*/,
486  double* dum_coords,
487  double* dum_pos )
488 {
489  std::vector< EntityHandle > verts, neighbor_hexes;
490  SetIterator* iter;
491  ErrorCode result = gMB->create_set_iterator( 0, MBVERTEX, -1, chunk_size, check_valid, iter );
492  RC( "query_vert_to_elem_iters" );
493  assert( MB_SUCCESS == result );
494  bool atend = false;
495  while( !atend )
496  {
497  verts.clear();
498  result = iter->get_next_arr( verts, atend );
499  RC( "query_vert_to_elem_iters" );
500  result = gMB->get_coords( &verts[0], verts.size(), dum_coords );
501  RC( "query_vert_to_elem_iters" );
502  chunk_size = std::min( (int)verts.size(), chunk_size );
503  for( int i = 0; i < chunk_size; i++ )
504  {
505  neighbor_hexes.clear();
506  result = gMB->get_adjacencies( &verts[i], 1, 3, false, neighbor_hexes );
507  RC( "query_vert_to_elem_iters" );
508  result = gMB->tag_get_data( pos2_tag, &neighbor_hexes[0], neighbor_hexes.size(), dum_pos );
509  RC( "query_vert_to_elem_iters" );
510  for( unsigned int j = 0; j < neighbor_hexes.size(); j++ )
511  {
512  dum_pos[3 * j + 0] += dum_coords[3 * i + 0];
513  dum_pos[3 * j + 1] += dum_coords[3 * i + 1];
514  dum_pos[3 * j + 2] += dum_coords[3 * i + 2];
515  }
516  result = gMB->tag_set_data( pos2_tag, &neighbor_hexes[0], neighbor_hexes.size(), dum_pos );
517  RC( "query_vert_to_elem_iters" );
518  }
519  }
520 
521  result = normalize_elems( dum_coords );
522  RC( "query_vert_to_elem_iters" );
523 }
524 
525 // normalize pos_tag on all elems by 1/8; coords assumed large enough to hold 3 doubles
526 ErrorCode normalize_elems( double* coords )
527 {
528  // get all hexes and divide pos_tag by 8
529  Range elems;
530  ErrorCode result = gMB->get_entities_by_type( 0, MBHEX, elems );
531  RR( "normalize" );
532 
533  for( Range::iterator vit = elems.begin(); vit != elems.end(); ++vit )
534  {
535  result = gMB->tag_get_data( pos2_tag, &( *vit ), 1, coords );
536  RR( "normalize" );
537  coords[0] *= 0.125;
538  coords[1] *= 0.125;
539  coords[2] *= 0.125;
540  result = gMB->tag_set_data( pos2_tag, &( *vit ), 1, coords );
541  RR( "normalize" );
542  }
543 
544  return result;
545 }
546 
547 void query_struct_elem_to_vert()
548 {
549  // assumes brick mapped mesh with handles starting at zero
550  Range all_hexes;
551  ErrorCode result = gMB->get_entities_by_type( 0, MBHEX, all_hexes );
552  RC( "query_struct_elem_to_vert" );
553  double dum_coords[24];
554  std::vector< EntityHandle > connect;
555  for( Range::iterator eit = all_hexes.begin(); eit != all_hexes.end(); ++eit )
556  {
557  result = gMB->get_connectivity( &( *eit ), 1, connect );
558  RC( "query_struct_elem_to_vert" );
559  result = gMB->get_coords( &connect[0], connect.size(), dum_coords );
560  RC( "query_struct_elem_to_vert" );
561 
562  double centroid[3] = { 0.0, 0.0, 0.0 };
563  for( int j = 0; j < 8; j++ )
564  {
565  centroid[0] += dum_coords[3 * j + 0];
566  centroid[1] += dum_coords[3 * j + 1];
567  centroid[2] += dum_coords[3 * j + 2];
568  }
569  centroid[0] *= 0.125;
570  centroid[1] *= 0.125;
571  centroid[2] *= 0.125;
572  result = gMB->tag_set_data( pos_tag, &( *eit ), 1, centroid );
573  RC( "query_struct_elem_to_vert" );
574  }
575 }
576 
577 #if defined( _MSC_VER ) || defined( __MINGW32__ )
578 void print_time( const bool print_em, double& tot_time, double& utime, double& stime, long& imem, long& rmem )
579 {
580  utime = (double)clock() / CLOCKS_PER_SEC;
581  if( print_em ) std::cout << "Total wall time = " << utime << std::endl;
582  tot_time = stime = 0;
583  imem = rmem = 0;
584 }
585 #else
586 void print_time( const bool print_em, double& tot_time, double& utime, double& stime, long& imem, long& rmem )
587 {
588  struct rusage r_usage;
589  getrusage( RUSAGE_SELF, &r_usage );
590  utime = (double)r_usage.ru_utime.tv_sec + ( (double)r_usage.ru_utime.tv_usec / 1.e6 );
591  stime = (double)r_usage.ru_stime.tv_sec + ( (double)r_usage.ru_stime.tv_usec / 1.e6 );
592  tot_time = utime + stime;
593  if( print_em )
594  std::cout << "User, system, total time = " << utime << ", " << stime << ", " << tot_time << std::endl;
595 #ifndef LINUX
596  if( print_em )
597  {
598  std::cout << "Max resident set size = " << r_usage.ru_maxrss << " kbytes" << std::endl;
599  std::cout << "Int resident set size = " << r_usage.ru_idrss << std::endl;
600  }
601  imem = r_usage.ru_idrss;
602  rmem = r_usage.ru_maxrss;
603 #else
604  system( "ps o args,drs,rss | grep perf | grep -v grep" ); // RedHat 9.0 doesnt fill in actual
605  // memory data
606  imem = rmem = 0;
607 #endif
608 }
609 #endif
610 
611 void testA( const int nelem, const double* coords )
612 {
613  double ttime0 = 0.0, ttime1 = 0.0, ttime2 = 0.0, ttime3 = 0.0, ttime4 = 0.0, utime = 0.0, stime = 0.0;
614  long imem0 = 0, rmem0 = 0, imem1 = 0, rmem1 = 0, imem2 = 0, rmem2 = 0, imem3 = 0, rmem3 = 0, imem4 = 0, rmem4 = 0;
615 
616  print_time( false, ttime0, utime, stime, imem0, rmem0 );
617 
618  // make a 3d block of vertices
619  EntitySequence* dum_seq = NULL;
620  ScdVertexData* vseq = NULL;
621  StructuredElementSeq* eseq = NULL;
622  init();
623  Core* mbcore = dynamic_cast< Core* >( gMB );
624  assert( mbcore != NULL );
625  SequenceManager* seq_mgr = mbcore->sequence_manager();
626  HomCoord vseq_minmax[2] = { HomCoord( 0, 0, 0 ), HomCoord( nelem, nelem, nelem ) };
627  EntityHandle vstart, estart;
628 
629  ErrorCode result = seq_mgr->create_scd_sequence( vseq_minmax[0], vseq_minmax[1], MBVERTEX, 1, vstart, dum_seq );
630  RC( "testA" );
631  if( NULL != dum_seq ) vseq = dynamic_cast< ScdVertexData* >( dum_seq->data() );
632  assert( MB_FAILURE != result && vstart != 0 && dum_seq != NULL && vseq != NULL );
633  // now the element sequence
634  result = seq_mgr->create_scd_sequence( vseq_minmax[0], vseq_minmax[1], MBHEX, 1, estart, dum_seq );
635  if( NULL != dum_seq ) eseq = dynamic_cast< StructuredElementSeq* >( dum_seq );
636  assert( MB_FAILURE != result && estart != 0 && dum_seq != NULL && eseq != NULL );
637 
638  // only need to add one vseq to this, unity transform
639  // trick: if I know it's going to be unity, just input 3 sets of equivalent points
640  result = eseq->sdata()->add_vsequence( vseq, vseq_minmax[0], vseq_minmax[0], vseq_minmax[0], vseq_minmax[0],
641  vseq_minmax[0], vseq_minmax[0] );
642  assert( MB_SUCCESS == result );
643 
644  // set the coordinates of the vertices
645  EntityHandle handle;
646  int i;
647  double dumv[3];
648  int num_verts = ( nelem + 1 ) * ( nelem + 1 ) * ( nelem + 1 );
649  for( i = 0, handle = vstart; i < num_verts; i++, handle++ )
650  {
651  dumv[0] = coords[i];
652  dumv[1] = coords[num_verts + i];
653  dumv[2] = coords[2 * num_verts + i];
654  result = gMB->set_coords( &handle, 1, dumv );
655  assert( MB_SUCCESS == result );
656  }
657 
658  print_time( false, ttime1, utime, stime, imem1, rmem1 );
659 
660  // query the mesh 2 ways
661  query_struct_elem_to_vert();
662 
663  print_time( false, ttime2, utime, stime, imem2, rmem2 );
664 
665  query_vert_to_elem();
666 
667  print_time( false, ttime3, utime, stime, imem3, rmem3 );
668 
669 #ifndef NDEBUG
670  check_answers( "A" );
671 #endif
672 
673  delete gMB;
674 
675  print_time( false, ttime4, utime, stime, imem4, rmem4 );
676 
677  std::cout << "MOAB_scd:nelem,construct,e_to_v,v_to_e,after_dtor,total= " << nelem << " " << ttime1 - ttime0 << " "
678  << ttime2 - ttime1 << " " << ttime3 - ttime2 << " " << ttime4 - ttime3 << " " << ttime4 - ttime0
679  << " seconds" << std::endl;
680  std::cout << "MOAB_scd_memory(rss):initial,after_construction,e-v,v-e,after_dtor= " << rmem0 << " " << rmem1 << " "
681  << rmem2 << " " << rmem3 << " " << rmem4 << " kb" << std::endl;
682 }
683 
684 void testB( const int nelem, const double* coords, int* connect )
685 {
686  double ttime0 = 0.0, ttime1 = 0.0, ttime2 = 0.0, ttime3 = 0.0, ttime4 = 0.0, utime = 0.0, stime = 0.0;
687  long imem0 = 0, rmem0 = 0, imem1 = 0, rmem1 = 0, imem2 = 0, rmem2 = 0, imem3 = 0, rmem3 = 0, imem4 = 0, rmem4 = 0;
688 
689  print_time( false, ttime0, utime, stime, imem0, rmem0 );
690 
691  int num_verts = ( nelem + 1 ) * ( nelem + 1 ) * ( nelem + 1 );
692  int num_elems = nelem * nelem * nelem;
693  EntityHandle vstart, estart;
694 
695  // get the read interface
696  ReadUtilIface* readMeshIface;
697  init();
698  gMB->query_interface( readMeshIface );
699 
700  // create a sequence to hold the node coordinates
701  // get the current number of entities and start at the next slot
702  std::vector< double* > coord_arrays;
703  ErrorCode result = readMeshIface->get_node_coords( 3, num_verts, 1, vstart, coord_arrays );
704  RC( "testB" );
705  assert( MB_SUCCESS == result && 1 == vstart && coord_arrays[0] && coord_arrays[1] && coord_arrays[2] );
706  // memcpy the coordinate data into place
707  memcpy( coord_arrays[0], coords, sizeof( double ) * num_verts );
708  memcpy( coord_arrays[1], &coords[num_verts], sizeof( double ) * num_verts );
709  memcpy( coord_arrays[2], &coords[2 * num_verts], sizeof( double ) * num_verts );
710 
711  EntityHandle* conn = 0;
712  result = readMeshIface->get_element_connect( num_elems, 8, MBHEX, 1, estart, conn );
713  assert( MB_SUCCESS == result );
714  for( int i = 0; i < num_elems * 8; i++ )
715  conn[i] = vstart + connect[i];
716 
717  delete[] connect;
718 
719  result = readMeshIface->update_adjacencies( estart, num_elems, 8, conn );
720  assert( MB_SUCCESS == result );
721 
722  print_time( false, ttime1, utime, stime, imem1, rmem1 );
723 
724  // query the mesh 2 ways
725  query_elem_to_vert();
726 
727  print_time( false, ttime2, utime, stime, imem2, rmem2 );
728 
729  query_vert_to_elem();
730 
731  print_time( false, ttime3, utime, stime, imem3, rmem3 );
732 
733 #ifndef NDEBUG
734  check_answers( "B" );
735 #endif
736 
737  delete gMB;
738 
739  print_time( false, ttime4, utime, stime, imem4, rmem4 );
740 
741  std::cout << "MOAB_ucd_blocked:nelem,construct,e_to_v,v_to_e,after_dtor,total= " << nelem << " " << ttime1 - ttime0
742  << " " << ttime2 - ttime1 << " " << ttime3 - ttime2 << " " << ttime4 - ttime3 << " " << ttime4 - ttime0
743  << " seconds" << std::endl;
744  std::cout << "MOAB_ucdblocked_memory_(rss):initial,after_construction,e-v,v-e,after_dtor= " << rmem0 << " " << rmem1
745  << " " << rmem2 << " " << rmem3 << " " << rmem4 << " kb" << std::endl;
746 }
747 
748 void testC( const int nelem, const double* coords )
749 {
750  double ttime0 = 0.0, ttime1 = 0.0, ttime2 = 0.0, ttime3 = 0.0, ttime4 = 0.0, utime = 0.0, stime = 0.0;
751  long imem0 = 0, rmem0 = 0, imem1 = 0, rmem1 = 0, imem2 = 0, rmem2 = 0, imem3 = 0, rmem3 = 0, imem4 = 0, rmem4 = 0;
752 
753  print_time( false, ttime0, utime, stime, imem0, rmem0 );
754 
755  // create the vertices; assume we don't need to keep a list of vertex handles, since they'll
756  // be created in sequence
757  int numv = nelem + 1;
758  int numv_sq = numv * numv;
759  int num_verts = numv * numv * numv;
760  double dum_coords[3] = { coords[0], coords[num_verts], coords[2 * num_verts] };
761  EntityHandle vstart;
762 
763  init();
764  ErrorCode result = gMB->create_vertex( dum_coords, vstart );
765  RC( "testC" );
766  assert( MB_SUCCESS == result && 1 == vstart );
767 
768  EntityHandle dum_vert, vijk;
769  int i;
770  for( i = 1; i < num_verts; i++ )
771  {
772  dum_coords[0] = coords[i];
773  dum_coords[1] = coords[num_verts + i];
774  dum_coords[2] = coords[2 * num_verts + i];
775  result = gMB->create_vertex( dum_coords, dum_vert );
776  assert( MB_SUCCESS == result );
777  }
778 
779  EntityHandle dum_conn[8];
780  for( i = 0; i < nelem; i++ )
781  {
782  for( int j = 0; j < nelem; j++ )
783  {
784  for( int k = 0; k < nelem; k++ )
785  {
786  vijk = vstart + VINDEX( i, j, k );
787  dum_conn[0] = vijk;
788  dum_conn[1] = vijk + 1;
789  dum_conn[2] = vijk + 1 + numv;
790  dum_conn[3] = vijk + numv;
791  dum_conn[4] = vijk + numv * numv;
792  dum_conn[5] = vijk + 1 + numv * numv;
793  dum_conn[6] = vijk + 1 + numv + numv * numv;
794  dum_conn[7] = vijk + numv + numv * numv;
795  result = gMB->create_element( MBHEX, dum_conn, 8, dum_vert );
796  assert( MB_SUCCESS == result );
797  }
798  }
799  }
800 
801  print_time( false, ttime1, utime, stime, imem1, rmem1 );
802 
803  // query the mesh 2 ways
804  query_elem_to_vert();
805 
806  print_time( false, ttime2, utime, stime, imem2, rmem2 );
807 
808  query_vert_to_elem();
809 
810  print_time( false, ttime3, utime, stime, imem3, rmem3 );
811 
812 #ifndef NDEBUG
813  check_answers( "C" );
814 #endif
815 
816  delete gMB;
817 
818  print_time( false, ttime4, utime, stime, imem4, rmem4 );
819 
820  std::cout << "MOAB_ucd_indiv:nelem,construct,e_to_v,v_to_e,after_dtor,total= " << nelem << " " << ttime1 - ttime0
821  << " " << ttime2 - ttime1 << " " << ttime3 - ttime2 << " " << ttime4 - ttime3 << " " << ttime4 - ttime0
822  << " seconds" << std::endl;
823  std::cout << "MOAB_ucd_indiv_memory_(rss):initial,after_construction,e-v,v-e,after_dtor= " << rmem0 << " " << rmem1
824  << " " << rmem2 << " " << rmem3 << " " << rmem4 << " kb" << std::endl;
825 }
826 
827 void testD( const int nelem, const double* coords, int ver )
828 {
829  double ttime0 = 0.0, ttime1 = 0.0, ttime2 = 0.0, ttime3 = 0.0, ttime4 = 0.0, ttime5 = 0.0, ttime6 = 0.0,
830  ttime7 = 0.0, ttime8 = 0.0, ttime9 = 0.0, ttime10 = 0.0, utime = 0.0, stime = 0.0;
831  long imem0 = 0, rmem0 = 0, imem1 = 0, rmem1 = 0, imem2 = 0, rmem2 = 0, imem3 = 0, rmem3 = 0, imem4 = 0, rmem4 = 0,
832  imem5 = 0, rmem5 = 0, imem6 = 0, rmem6 = 0, imem7 = 0, rmem7 = 0, imem8 = 0, rmem8 = 0, imem9 = 0, rmem9 = 0,
833  imem10 = 0, rmem10 = 0;
834 
835  print_time( false, ttime0, utime, stime, imem0, rmem0 );
836 
837  // create the vertices; assume we don't need to keep a list of vertex handles, since they'll
838  // be created in sequence
839  int numv = nelem + 1;
840  int numv_sq = numv * numv;
841  int num_verts = numv * numv * numv;
842  std::vector< double > dum_coords( 24 ), dum_pos( 24 );
843  dum_coords[0] = coords[0];
844  dum_coords[1] = coords[num_verts];
845  dum_coords[2] = coords[2 * num_verts];
846  EntityHandle vstart;
847 
848  init();
849  ErrorCode result = gMB->create_vertex( &dum_coords[0], vstart );
850  RC( "testD" );
851  assert( MB_SUCCESS == result && 1 == vstart );
852 
853  EntityHandle dum_vert, vijk;
854  int i;
855  for( i = 1; i < num_verts; i++ )
856  {
857  dum_coords[0] = coords[i];
858  dum_coords[1] = coords[num_verts + i];
859  dum_coords[2] = coords[2 * num_verts + i];
860  result = gMB->create_vertex( &dum_coords[0], dum_vert );
861  assert( MB_SUCCESS == result );
862  }
863 
864  EntityHandle dum_conn[8];
865  for( i = 0; i < nelem; i++ )
866  {
867  for( int j = 0; j < nelem; j++ )
868  {
869  for( int k = 0; k < nelem; k++ )
870  {
871  vijk = vstart + VINDEX( i, j, k );
872  dum_conn[0] = vijk;
873  dum_conn[1] = vijk + 1;
874  dum_conn[2] = vijk + 1 + numv;
875  dum_conn[3] = vijk + numv;
876  dum_conn[4] = vijk + numv * numv;
877  dum_conn[5] = vijk + 1 + numv * numv;
878  dum_conn[6] = vijk + 1 + numv + numv * numv;
879  dum_conn[7] = vijk + numv + numv * numv;
880  result = gMB->create_element( MBHEX, dum_conn, 8, dum_vert );
881  assert( MB_SUCCESS == result );
882  }
883  }
884  }
885 
886  print_time( false, ttime1, utime, stime, imem1, rmem1 );
887 
888  // query the mesh 2 ways with !check_valid
889  std::vector< EntityHandle > connect( 8 );
890 #ifndef NDEBUG
891  // used only in debug mode
892  double def_val[3] = { 0.0, 0.0, 0.0 };
893 #endif
894  if( ver == 0 || ver == 1 )
895  {
896  query_elem_to_vert_iters( 1, false, connect, &dum_coords[0], &dum_pos[0] );
897  print_time( false, ttime2, utime, stime, imem2, rmem2 );
898  query_vert_to_elem_iters( 1, false, connect, &dum_coords[0], &dum_pos[0] );
899  print_time( false, ttime3, utime, stime, imem3, rmem3 );
900 #ifndef NDEBUG
901  check_answers( "D" );
902  result = gMB->tag_delete( pos_tag );
903  assert( MB_SUCCESS == result );
904  result =
905  gMB->tag_get_handle( "position_tag", 3, MB_TYPE_DOUBLE, pos_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
906  assert( MB_SUCCESS == result );
907  result = gMB->tag_delete( pos2_tag );
908  assert( MB_SUCCESS == result );
909  result =
910  gMB->tag_get_handle( "position2_tag", 3, MB_TYPE_DOUBLE, pos2_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
911  assert( MB_SUCCESS == result );
912 #endif
913  }
914 
915  if( ver == 0 || ver == 2 )
916  {
917  if( ver != 0 ) print_time( false, ttime3, utime, stime, imem3, rmem3 );
918  query_elem_to_vert_iters( 1, true, connect, &dum_coords[0], &dum_pos[0] );
919  print_time( false, ttime4, utime, stime, imem4, rmem4 );
920  query_vert_to_elem_iters( 1, true, connect, &dum_coords[0], &dum_pos[0] );
921  print_time( false, ttime5, utime, stime, imem5, rmem5 );
922 #ifndef NDEBUG
923  check_answers( "D" );
924  result = gMB->tag_delete( pos_tag );
925  assert( MB_SUCCESS == result );
926  result =
927  gMB->tag_get_handle( "position_tag", 3, MB_TYPE_DOUBLE, pos_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
928  assert( MB_SUCCESS == result );
929  result = gMB->tag_delete( pos2_tag );
930  assert( MB_SUCCESS == result );
931  result =
932  gMB->tag_get_handle( "position2_tag", 3, MB_TYPE_DOUBLE, pos2_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
933  assert( MB_SUCCESS == result );
934 #endif
935  }
936 
937  if( ver == 0 || ver >= 3 )
938  {
939  dum_coords.resize( 2400 );
940  dum_pos.resize( 300 );
941  }
942  if( ver == 0 || ver == 3 )
943  {
944  if( ver != 0 ) print_time( false, ttime5, utime, stime, imem3, rmem3 );
945  query_elem_to_vert_iters( 100, false, connect, &dum_coords[0], &dum_pos[0] );
946  print_time( false, ttime6, utime, stime, imem6, rmem6 );
947  query_vert_to_elem_iters( 100, false, connect, &dum_coords[0], &dum_pos[0] );
948  print_time( false, ttime7, utime, stime, imem7, rmem7 );
949 #ifndef NDEBUG
950  check_answers( "D" );
951  result = gMB->tag_delete( pos_tag );
952  assert( MB_SUCCESS == result );
953  result =
954  gMB->tag_get_handle( "position_tag", 3, MB_TYPE_DOUBLE, pos_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
955  assert( MB_SUCCESS == result );
956  result = gMB->tag_delete( pos2_tag );
957  assert( MB_SUCCESS == result );
958  result =
959  gMB->tag_get_handle( "position2_tag", 3, MB_TYPE_DOUBLE, pos2_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
960  assert( MB_SUCCESS == result );
961 #endif
962  }
963 
964  if( ver == 0 || ver == 4 )
965  {
966  if( ver != 0 ) print_time( false, ttime7, utime, stime, imem3, rmem3 );
967  query_elem_to_vert_iters( 100, true, connect, &dum_coords[0], &dum_pos[0] );
968  print_time( false, ttime8, utime, stime, imem8, rmem8 );
969  query_vert_to_elem_iters( 100, true, connect, &dum_coords[0], &dum_pos[0] );
970  print_time( false, ttime9, utime, stime, imem9, rmem9 );
971 #ifndef NDEBUG
972  check_answers( "D" );
973  result = gMB->tag_delete( pos_tag );
974  assert( MB_SUCCESS == result );
975  result =
976  gMB->tag_get_handle( "position_tag", 3, MB_TYPE_DOUBLE, pos_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
977  assert( MB_SUCCESS == result );
978  result = gMB->tag_delete( pos2_tag );
979  assert( MB_SUCCESS == result );
980  result =
981  gMB->tag_get_handle( "position2_tag", 3, MB_TYPE_DOUBLE, pos2_tag, MB_TAG_DENSE | MB_TAG_CREAT, def_val );
982  assert( MB_SUCCESS == result );
983 #endif
984  }
985 
986  if( ver > 0 && ver < 4 ) print_time( false, ttime9, utime, stime, imem9, rmem9 );
987  delete gMB;
988 
989  print_time( false, ttime10, utime, stime, imem10, rmem10 );
990 
991  if( ver == 0 || ver == 1 )
992  {
993  std::cout << "MOAB_ucd_iters_!check_valid_1:nelem,construct,e-v,v-e,after_dtor,total= " << nelem << " "
994  << ttime1 - ttime0 << " " << ttime2 - ttime1 << " " << ttime3 - ttime2 << " " << ttime10 - ttime9
995  << " " << ttime3 - ttime0 + ttime10 - ttime9 << std::endl;
996  std::cout << "MOAB_ucd_iters_memory_(rss)_!check_valid_1:initial,after_construction,e-v,v-"
997  "e,after_dtor= "
998  << rmem0 << " " << rmem1 << " " << rmem2 << " " << rmem3 << " " << rmem10 << " kb" << std::endl;
999  }
1000  if( ver == 0 || ver == 2 )
1001  {
1002  std::cout << "MOAB_ucd_iters_check_valid_1:nelem,construct,e-v,v-e,after_dtor,total= " << nelem << " "
1003  << ttime1 - ttime0 << " " << ttime4 - ttime3 << " " << ttime5 - ttime4 << " " << ttime10 - ttime9
1004  << " " << ttime1 - ttime0 + ttime5 - ttime3 + ttime10 - ttime9 << std::endl;
1005  std::cout << "MOAB_ucd_iters_memory_(rss)_check_valid_1:initial,after_construction,e-v,v-e,"
1006  "after_dtor= "
1007  << rmem0 << " " << rmem1 << " " << rmem2 << " " << rmem3 << " " << rmem10 << " kb" << std::endl;
1008  }
1009  if( ver == 0 || ver == 3 )
1010  {
1011  std::cout << "MOAB_ucd_iters_!check_valid_100:nelem,construct,e-v,v-e,after_dtor,total= " << nelem << " "
1012  << ttime1 - ttime0 << " " << ttime6 - ttime5 << " " << ttime7 - ttime6 << " " << ttime10 - ttime9
1013  << " " << ttime1 - ttime0 + ttime7 - ttime5 + ttime10 - ttime9 << std::endl;
1014  std::cout << "MOAB_ucd_iters_memory_(rss)_!check_valid_100:initial,after_construction,e-v,"
1015  "v-e,after_dtor= "
1016  << rmem0 << " " << rmem1 << " " << rmem6 << " " << rmem7 << " " << rmem10 << " kb" << std::endl;
1017  }
1018  if( ver == 0 || ver == 4 )
1019  {
1020  std::cout << "MOAB_ucd_iters_check_valid_100:nelem,construct,e-v,v-e,after_dtor,total= " << nelem << " "
1021  << ttime1 - ttime0 << " " << ttime8 - ttime7 << " " << ttime9 - ttime8 << " " << ttime10 - ttime9
1022  << " " << ttime1 - ttime0 + ttime10 - ttime7 << std::endl;
1023  std::cout << "MOAB_ucd_iters_memory_(rss)_check_valid_100:initial,after_construction,e-v,v-"
1024  "e,after_dtor= "
1025  << rmem0 << " " << rmem1 << " " << rmem8 << " " << rmem9 << " " << rmem10 << " kb" << std::endl;
1026  }
1027 }
1028 
1029 void testE( const int nelem, const double* coords, int* connect )
1030 {
1031  double ttime0 = 0.0, ttime1 = 0.0, ttime2 = 0.0, ttime3 = 0.0, ttime4 = 0.0, ttime5 = 0.0, ttime6 = 0.0,
1032  utime = 0.0, stime = 0.0;
1033  long imem0 = 0, rmem0 = 0, imem1 = 0, rmem1 = 0, imem2 = 0, rmem2 = 0, imem3 = 0, rmem3 = 0, imem4 = 0, rmem4 = 0,
1034  imem5 = 0, rmem5 = 0, imem6 = 0, rmem6 = 0;
1035 
1036  print_time( false, ttime0, utime, stime, imem0, rmem0 );
1037 
1038  int num_verts = ( nelem + 1 ) * ( nelem + 1 ) * ( nelem + 1 );
1039  int num_elems = nelem * nelem * nelem;
1040  EntityHandle vstart, estart;
1041 
1042  // get the read interface
1043  ReadUtilIface* readMeshIface;
1044  init();
1045  gMB->query_interface( readMeshIface );
1046 
1047  // create a sequence to hold the node coordinates
1048  // get the current number of entities and start at the next slot
1049  std::vector< double* > coord_arrays;
1050  ErrorCode result = readMeshIface->get_node_coords( 3, num_verts, 1, vstart, coord_arrays );
1051  RC( "testE" );
1052  assert( MB_SUCCESS == result && 1 == vstart && coord_arrays[0] && coord_arrays[1] && coord_arrays[2] );
1053  // memcpy the coordinate data into place
1054  memcpy( coord_arrays[0], coords, sizeof( double ) * num_verts );
1055  memcpy( coord_arrays[1], &coords[num_verts], sizeof( double ) * num_verts );
1056  memcpy( coord_arrays[2], &coords[2 * num_verts], sizeof( double ) * num_verts );
1057 
1058  EntityHandle* conn = 0;
1059  result = readMeshIface->get_element_connect( num_elems, 8, MBHEX, 1, estart, conn );
1060  assert( MB_SUCCESS == result );
1061  for( int i = 0; i < num_elems * 8; i++ )
1062  conn[i] = vstart + connect[i];
1063 
1064  delete[] connect;
1065 
1066  Range verts( vstart, vstart + num_verts - 1 ), elems( estart, estart + num_elems - 1 );
1067 
1068  result = readMeshIface->update_adjacencies( estart, num_elems, 8, conn );
1069  assert( MB_SUCCESS == result );
1070 
1071  print_time( false, ttime1, utime, stime, imem1, rmem1 );
1072 
1073  // query the mesh 2 ways
1074  query_elem_to_vert_direct();
1075 
1076  print_time( false, ttime2, utime, stime, imem2, rmem2 );
1077 
1078  query_vert_to_elem_direct();
1079 
1080  print_time( false, ttime3, utime, stime, imem3, rmem3 );
1081 
1082 #ifndef NDEBUG
1083  check_answers( "E" );
1084 #endif
1085 
1086  query_elem_to_vert_direct();
1087 
1088  print_time( false, ttime4, utime, stime, imem4, rmem4 );
1089 
1090  query_vert_to_elem_direct();
1091 
1092  print_time( false, ttime5, utime, stime, imem5, rmem5 );
1093 
1094 #ifndef NDEBUG
1095  check_answers( "E" );
1096 #endif
1097 
1098  delete gMB;
1099 
1100  print_time( false, ttime6, utime, stime, imem6, rmem6 );
1101 
1102  std::cout << "MOAB_ucd_direct:nelem,construct,e_to_v,v_to_e,after_dtor,total= " << nelem << " " << ttime1 - ttime0
1103  << " " << ttime2 - ttime1 << " " << ttime3 - ttime2 << " " << ttime6 - ttime5 << " "
1104  << ttime3 - ttime0 + ttime6 - ttime5 << " seconds" << std::endl;
1105  std::cout << "MOAB_ucd_direct_memory_(rss):initial,after_construction,e-v,v-e,after_dtor= " << rmem0 << " " << rmem1
1106  << " " << rmem2 << " " << rmem3 << " " << rmem6 << " kb" << std::endl;
1107 
1108  std::cout << "MOAB_ucd_direct2:nelem,construct,e_to_v,v_to_e,after_dtor,total= " << nelem << " " << ttime1 - ttime0
1109  << " " << ttime4 - ttime3 << " " << ttime5 - ttime4 << " " << ttime6 - ttime5 << " "
1110  << ttime1 - ttime0 + ttime6 - ttime3 << " seconds" << std::endl;
1111  std::cout << "MOAB_ucd_direct2_memory_(rss):initial,after_construction,e-v,v-e,after_dtor= " << rmem0 << " "
1112  << rmem1 << " " << rmem4 << " " << rmem5 << " " << rmem6 << " kb" << std::endl;
1113 }
1114 
1115 void query_elem_to_vert_direct()
1116 {
1117  Range all_hexes, all_verts;
1118  ErrorCode result = gMB->get_entities_by_type( 0, MBHEX, all_hexes );
1119  assert( MB_SUCCESS == result );
1120  result = gMB->get_entities_by_type( 0, MBVERTEX, all_verts );
1121  RC( "query_elem_to_vert_direct" );
1122  EntityHandle* connect;
1123  int ecount, vcount, vpere;
1124  double* coords[3];
1125 
1126  result = gMB->connect_iterate( all_hexes.begin(), all_hexes.end(), connect, vpere, ecount );
1127  if( MB_SUCCESS != result || ecount != (int)all_hexes.size() )
1128  {
1129  std::cout << "FAILED in connect_iterate!" << std::endl;
1130  return;
1131  }
1132  result = gMB->coords_iterate( all_verts.begin(), all_verts.end(), coords[0], coords[1], coords[2], vcount );
1133  if( MB_SUCCESS != result || vcount != (int)all_verts.size() )
1134  {
1135  std::cout << "FAILED in coords_iterate!" << std::endl;
1136  return;
1137  }
1138  double* centroid;
1139  result = gMB->tag_iterate( pos_tag, all_hexes.begin(), all_hexes.end(), ecount, (void*&)centroid );
1140  if( MB_SUCCESS != result || ecount != (int)all_hexes.size() )
1141  {
1142  std::cout << "FAILED in connect_iterate!" << std::endl;
1143  return;
1144  }
1145 
1146  EntityHandle vstart = *all_verts.begin();
1147  for( int i = 0; i < ecount; i++ )
1148  {
1149  // compute the centroid
1150  for( int j = 0; j < vpere; j++ )
1151  {
1152  int vind = *connect - vstart;
1153  connect++;
1154  centroid[3 * i + 0] += coords[0][vind];
1155  centroid[3 * i + 1] += coords[1][vind];
1156  centroid[3 * i + 2] += coords[2][vind];
1157  }
1158 
1159  // now normalize
1160  centroid[3 * i + 0] *= 0.125;
1161  centroid[3 * i + 1] *= 0.125;
1162  centroid[3 * i + 2] *= 0.125;
1163  }
1164 }
1165 
1166 void query_vert_to_elem_direct()
1167 {
1168  Range all_verts, tmp_ents;
1169  ErrorCode result = gMB->get_entities_by_type( 0, MBVERTEX, all_verts );
1170  assert( MB_SUCCESS == result );
1171 
1172  // make sure vertex-element adjacencies are created
1173  result = gMB->get_adjacencies( &( *all_verts.begin() ), 1, 3, false, tmp_ents );
1174  assert( MB_SUCCESS == result );
1175 
1176  const std::vector< EntityHandle >** adjs;
1177  int count;
1178  result = gMB->adjacencies_iterate( all_verts.begin(), all_verts.end(), adjs, count );
1179  if( MB_SUCCESS != result && count != (int)all_verts.size() )
1180  {
1181  std::cout << "FAILED:adjacencies_iterate." << std::endl;
1182  return;
1183  }
1184 
1185  double* coords[3];
1186  result = gMB->coords_iterate( all_verts.begin(), all_verts.end(), coords[0], coords[1], coords[2], count );
1187  if( MB_SUCCESS != result || count != (int)all_verts.size() )
1188  {
1189  std::cout << "FAILED in coords_iterate!" << std::endl;
1190  return;
1191  }
1192  // get all hexes, then iterator over pos2_tag
1193  double* centroid;
1194  int ecount;
1195  tmp_ents.clear();
1196  result = gMB->get_entities_by_type( 0, MBHEX, tmp_ents );
1197  assert( MB_SUCCESS == result );
1198  result = gMB->tag_iterate( pos2_tag, tmp_ents.begin(), tmp_ents.end(), ecount, (void*&)centroid );
1199  if( MB_SUCCESS != result || ecount != (int)tmp_ents.size() )
1200  {
1201  std::cout << "FAILED in tag_iterate!" << std::endl;
1202  return;
1203  }
1204 
1205  int i;
1206  Range::iterator vit;
1207  EntityHandle estart = *tmp_ents.begin();
1208  for( vit = all_verts.begin(), i = 0; vit != all_verts.end(); ++vit, i++ )
1209  {
1210  assert( adjs[i] );
1211  for( std::vector< EntityHandle >::const_iterator vit2 = adjs[i]->begin(); vit2 != adjs[i]->end(); ++vit2 )
1212  if( *vit >= estart )
1213  {
1214  int eind = *vit2 - estart;
1215  centroid[3 * eind + 0] += coords[0][i];
1216  centroid[3 * eind + 1] += coords[1][i];
1217  centroid[3 * eind + 2] += coords[2][i];
1218  }
1219  }
1220 
1221  // now normalize
1222  for( i = 0; i < (int)tmp_ents.size(); i++ )
1223  {
1224  centroid[3 * i + 0] *= 0.125;
1225  centroid[3 * i + 1] *= 0.125;
1226  centroid[3 * i + 2] *= 0.125;
1227  }
1228 }
1229 
1230 void check_answers( const char* /*test_name*/ )
1231 {
1232  Range elems;
1233  ErrorCode result = gMB->get_entities_by_type( 0, MBHEX, elems );
1234  RC( "check_answers" );
1235 
1236  double coords1[3], coords2[3], del[3];
1237  double diff = 0.0;
1238  for( Range::iterator vit = elems.begin(); vit != elems.end(); ++vit )
1239  {
1240  result = gMB->tag_get_data( pos_tag, &( *vit ), 1, coords1 );
1241  RC( "check_answers" );
1242  result = gMB->tag_get_data( pos2_tag, &( *vit ), 1, coords2 );
1243  RC( "check_answers" );
1244  for( int i = 0; i < 3; i++ )
1245  del[i] = fabs( coords1[i] - coords2[i] );
1246  if( del[0] || del[1] || del[2] )
1247  {
1248  double len2 = std::max( coords1[0] * coords1[0] + coords1[1] * coords1[1] + coords1[2] * coords1[2],
1249  coords2[0] * coords2[0] + coords2[1] * coords2[1] + coords2[2] * coords2[2] ),
1250  num = del[0] * del[0] + del[1] * del[1] + del[2] * del[2];
1251  if( len2 > 0.0 )
1252  diff = std::max( diff, num / sqrt( len2 ) );
1253  else if( num > 0.0 )
1254  diff = sqrt( num );
1255  }
1256  }
1257  if( diff > 0.0 ) std::cout << "Max relative difference = " << diff << std::endl;
1258 }