MOAB: Mesh Oriented datABase  (version 5.5.0)
elem_eval_time.cpp
Go to the documentation of this file.
1 #include "moab/Core.hpp"
2 #include "moab/Range.hpp"
3 #include "moab/ProgOptions.hpp"
4 #include "moab/CartVect.hpp"
5 #include "moab/ElemEvaluator.hpp"
6 #include "moab/CN.hpp"
7 #include "ElemUtil.hpp"
8 #include <iostream>
9 #include <ctime>
10 #include <cstdlib>
11 #include <cassert>
12 #include <sstream>
13 #include <string>
14 
15 // Different platforms follow different conventions for usage
16 #if !defined( _MSC_VER ) && !defined( __MINGW32__ )
17 #include <sys/resource.h>
18 #endif
19 #ifdef SOLARIS
20 extern "C" int getrusage( int, struct rusage* );
21 #ifndef RUSAGE_SELF
22 #include </usr/ucbinclude/sys/rusage.h>
23 #endif
24 #endif
25 
26 using namespace moab;
27 #define CHK( r, s ) \
28  do \
29  { \
30  if( MB_SUCCESS != ( r ) ) \
31  { \
32  fail( ( r ), ( s ), __FILE__, __LINE__ ); \
33  return ( r ); \
34  } \
35  } while( false )
36 
37 const int ELEMEVAL = 0, ELEMUTIL = 1;
38 
39 static void fail( ErrorCode error_code, const char* str, const char* file_name, int line_number )
40 {
41  std::cerr << str << ", line " << line_number << " of file " << file_name << ", error code " << error_code
42  << std::endl;
43 }
44 
45 double mytime();
46 
47 ErrorCode get_ents( Interface& mbi, std::string& filename, int& dim, Range& elems, EntityType& tp, int& nv )
48 {
49  ErrorCode rval = mbi.load_file( filename.c_str(), 0 );
50  while( elems.empty() && dim >= 1 )
51  {
52  rval = mbi.get_entities_by_dimension( 0, dim, elems );
53  CHK( rval, "get_entities_by_dimension" );
54  if( elems.empty() ) dim--;
55  }
56  if( elems.empty() )
57  {
58  CHK( MB_FAILURE, "No elements in file" );
59  }
60 
61  // check to see they're all the same type & #vertices
62  tp = mbi.type_from_handle( *elems.begin() );
63  EntityType tp2 = mbi.type_from_handle( *elems.rbegin() );
64  if( tp != tp2 ) CHK( MB_FAILURE, "Elements must have same type" );
65 
66  int nv2;
67  const EntityHandle* c1;
68  rval = mbi.get_connectivity( *elems.begin(), c1, nv );
69  CHK( rval, "get_connectivity" );
70  rval = mbi.get_connectivity( *elems.rbegin(), c1, nv2 );
71  CHK( rval, "get_connectivity" );
72  if( nv != nv2 )
73  {
74  CHK( MB_FAILURE, "Elements must have same #vertices" );
75  }
76 
77  return MB_SUCCESS;
78 }
79 
80 void parse_options( ProgOptions& opts, int& dim, std::string& filename )
81 {
82  opts.addOpt< int >( std::string( "dim" ),
83  std::string( "Evaluate 1d/2d/3d elements (default: maximal dimension in mesh)" ), &dim,
84  ProgOptions::int_flag );
85  opts.addOpt< std::string >( std::string( "filename,f" ), std::string( "Filename containing mesh" ), &filename );
86 }
87 
88 ErrorCode get_elem_map( EntityType tp, std::vector< CartVect >& vcoords, int nconn, Element::Map*& elemmap )
89 {
90  switch( tp )
91  {
92  case MBHEX:
93  if( nconn == 8 )
94  {
95  elemmap = new Element::LinearHex( vcoords );
96  break;
97  }
98  else if( nconn == 27 )
99  {
100  elemmap = new Element::QuadraticHex( vcoords );
101  break;
102  }
103  else
104  return MB_FAILURE;
105  case MBTET:
106  if( nconn == 4 )
107  {
108  elemmap = new Element::LinearTet( vcoords );
109  break;
110  }
111  else
112  return MB_FAILURE;
113  case MBQUAD:
114  if( nconn == 4 )
115  {
116  elemmap = new Element::LinearQuad( vcoords );
117  break;
118  }
119  else
120  return MB_FAILURE;
121  default:
122  return MB_FAILURE;
123  }
124 
125  return MB_SUCCESS;
126 }
127 
128 ErrorCode time_forward_eval( Interface* mbi,
129  int method,
130  Range& elems,
131  std::vector< CartVect >& params,
132  std::vector< CartVect >& coords,
133  double& evtime )
134 {
135  evtime = mytime();
136  ErrorCode rval;
137  Range::iterator rit;
138  unsigned int i;
139  if( ELEMEVAL == method )
140  {
141  // construct ElemEvaluator
142  EvalSet eset;
143  ElemEvaluator eeval( mbi );
144  eeval.set_eval_set( *elems.begin() );
145  eeval.set_tag_handle( 0, 0 ); // indicates coordinates as the field to evaluate
146 
147  for( rit = elems.begin(), i = 0; rit != elems.end(); ++rit, i++ )
148  {
149  eeval.set_ent_handle( *rit );
150  rval = eeval.eval( params[i].array(), coords[i].array(), 3 );
151 #ifndef NDEBUG
152  if( MB_SUCCESS != rval ) return rval;
153 #endif
154  }
155  }
156  else if( ELEMUTIL == method )
157  {
158  std::vector< CartVect > vcoords( CN::MAX_NODES_PER_ELEMENT );
159  const EntityHandle* connect;
160  int nconn;
161  Element::Map* elemmap = NULL;
162  for( rit = elems.begin(), i = 0; rit != elems.end(); ++rit, i++ )
163  {
164  rval = mbi->get_connectivity( *rit, connect, nconn );
165  CHK( rval, "get_connectivity" );
166  rval = mbi->get_coords( connect, nconn, vcoords[0].array() );
167  CHK( rval, "get_coords" );
168  rval = get_elem_map( mbi->type_from_handle( *rit ), vcoords, nconn, elemmap );
169  CHK( rval, "get_elem_map" );
170  coords[i] = elemmap->evaluate( params[i] );
171  }
172  }
173 
174  evtime = mytime() - evtime;
175  return MB_SUCCESS;
176 }
177 
178 ErrorCode time_reverse_eval( Interface* mbi,
179  int method,
180  Range& elems,
181  std::vector< CartVect >& coords,
182  std::vector< CartVect >& params,
183  double& retime )
184 {
185  retime = mytime();
186  ErrorCode rval;
187  Range::iterator rit;
188  unsigned int i;
189  if( ELEMEVAL == method )
190  {
191  EvalSet eset;
192  ElemEvaluator eeval( mbi );
193  eeval.set_eval_set( *elems.begin() );
194  eeval.set_tag_handle( 0, 0 ); // indicates coordinates as the field to evaluate
195  int ins;
196  for( rit = elems.begin(), i = 0; rit != elems.end(); ++rit, i++ )
197  {
198  eeval.set_ent_handle( *rit );
199  rval = eeval.reverse_eval( coords[i].array(), 1.0e-10, 1.0e-6, params[i].array(), &ins );
200  assert( ins );
201 #ifndef NDEBUG
202  if( MB_SUCCESS != rval ) return rval;
203 #endif
204  }
205  }
206  else if( ELEMUTIL == method )
207  {
208  std::vector< CartVect > vcoords( CN::MAX_NODES_PER_ELEMENT );
209  const EntityHandle* connect;
210  int nconn;
211  Element::Map* elemmap = NULL;
212  for( rit = elems.begin(), i = 0; rit != elems.end(); ++rit, i++ )
213  {
214  rval = mbi->get_connectivity( *rit, connect, nconn );
215  CHK( rval, "get_connectivity" );
216  rval = mbi->get_coords( connect, nconn, vcoords[0].array() );
217  CHK( rval, "get_coords" );
218  rval = get_elem_map( mbi->type_from_handle( *rit ), vcoords, nconn, elemmap );
219  CHK( rval, "get_elem_map" );
220  coords[i] = elemmap->ievaluate( coords[i], 1.0e-6 );
221  }
222  }
223  retime = mytime() - retime;
224  return MB_SUCCESS;
225 }
226 
227 ErrorCode time_jacobian( Interface* mbi, int method, Range& elems, std::vector< CartVect >& params, double& jactime )
228 {
229  jactime = mytime();
230  ErrorCode rval;
231  Range::iterator rit;
232  unsigned int i;
233  Matrix3 jac;
234  if( ELEMEVAL == method )
235  {
236  EvalSet eset;
237  ElemEvaluator eeval( mbi );
238  eeval.set_eval_set( *elems.begin() );
239  eeval.set_tag_handle( 0, 0 ); // indicates coordinates as the field to evaluate
240  for( rit = elems.begin(), i = 0; rit != elems.end(); ++rit, i++ )
241  {
242  eeval.set_ent_handle( *rit );
243  rval = eeval.jacobian( params[i].array(), jac.array() );
244 #ifndef NDEBUG
245  if( MB_SUCCESS != rval ) return rval;
246 #endif
247  }
248  }
249  else if( ELEMUTIL == method )
250  {
251  std::vector< CartVect > vcoords( CN::MAX_NODES_PER_ELEMENT );
252  const EntityHandle* connect;
253  int nconn;
254  Element::Map* elemmap = NULL;
255  for( rit = elems.begin(), i = 0; rit != elems.end(); ++rit, i++ )
256  {
257  rval = mbi->get_connectivity( *rit, connect, nconn );
258  CHK( rval, "get_connectivity" );
259  rval = mbi->get_coords( connect, nconn, vcoords[0].array() );
260  CHK( rval, "get_coords" );
261  rval = get_elem_map( mbi->type_from_handle( *rit ), vcoords, nconn, elemmap );
262  CHK( rval, "get_elem_map" );
263  jac = elemmap->jacobian( params[i] );
264  }
265  }
266  jactime = mytime() - jactime;
267  return MB_SUCCESS;
268 }
269 
270 ErrorCode time_integrate( Interface* mbi, int method, Tag tag, Range& elems, double& inttime )
271 {
272  inttime = mytime();
273  ErrorCode rval;
274  double integral = 0.0;
275  if( ELEMEVAL == method )
276  {
277  EvalSet eset;
278  ElemEvaluator eeval( mbi );
279  eeval.set_eval_set( *elems.begin() );
280  eeval.set_tag_handle( 0, 0 ); // indicates coordinates as the field to evaluate
281  rval = eeval.set_tag_handle( tag, 0 );
282  CHK( rval, "set_tag_handle" );
283  for( Range::iterator rit = elems.begin(); rit != elems.end(); ++rit )
284  {
285  eeval.set_ent_handle( *rit );
286  rval = eeval.integrate( &integral );
287 #ifndef NDEBUG
288  if( MB_SUCCESS != rval ) return rval;
289 #endif
290  }
291  }
292  else if( ELEMUTIL == method )
293  {
294  std::vector< CartVect > vcoords( CN::MAX_NODES_PER_ELEMENT );
295  std::vector< double > tagval( CN::MAX_NODES_PER_ELEMENT );
296  const EntityHandle* connect;
297  int nconn;
298  Element::Map* elemmap = NULL;
299  Range::iterator rit;
300  unsigned int i;
301  for( rit = elems.begin(), i = 0; rit != elems.end(); ++rit, i++ )
302  {
303  rval = mbi->get_connectivity( *rit, connect, nconn );
304  CHK( rval, "get_connectivity" );
305  rval = mbi->get_coords( connect, nconn, vcoords[0].array() );
306  CHK( rval, "get_coords" );
307  rval = get_elem_map( mbi->type_from_handle( *rit ), vcoords, nconn, elemmap );
308  CHK( rval, "get_elem_map" );
309  rval = mbi->tag_get_data( tag, connect, nconn, &tagval[0] );
310  CHK( rval, "tag_get_data" );
311  integral = elemmap->integrate_scalar_field( &tagval[0] );
312  }
313  }
314  inttime = mytime() - inttime;
315  std::cout << "integral = " << integral << std::endl;
316  return MB_SUCCESS;
317 }
318 
319 ErrorCode put_random_field( Interface& mbi, Tag& tag, Range& elems )
320 {
321  Range verts;
322  ErrorCode rval = mbi.get_adjacencies( elems, 0, false, verts, Interface::UNION );
323  CHK( rval, "get_adjacencies" );
324  rval = mbi.tag_get_handle( "mytag", 1, MB_TYPE_DOUBLE, tag, MB_TAG_CREAT | MB_TAG_DENSE );
325  CHK( rval, "tag_get_handle" );
326  std::vector< double > tag_vals( verts.size() );
327  for( unsigned int i = 0; i < verts.size(); i++ )
328  tag_vals[i] = ( (double)rand() ) / RAND_MAX;
329  rval = mbi.tag_set_data( tag, verts, &tag_vals[0] );
330  return rval;
331 }
332 
333 ErrorCode elem_evals( Interface* mbi,
334  int method,
335  Range& elems,
336  Tag tag,
337  std::vector< CartVect >& params,
338  std::vector< CartVect >& coords,
339  double& evtime,
340  double& retime,
341  double& jactime,
342  double& inttime )
343 {
344  evtime = 0, retime = 0, jactime = 0, inttime = 0; // initializations to avoid compiler warning
345 
346  // time/test forward evaluation, putting results into vector
347  ErrorCode rval = time_forward_eval( mbi, method, elems, params, coords, evtime );
348  CHK( rval, "time_forward_eval" );
349 
350  // time/test reverse evaluation, putting results into vector
351  rval = time_reverse_eval( mbi, method, elems, coords, params, retime );
352  CHK( rval, "time_reverse_eval" );
353 
354  // time/test Jacobian evaluation
355  rval = time_jacobian( mbi, method, elems, params, jactime );
356  CHK( rval, "time_jacobian" );
357 
358  // time/test integration
359  rval = time_integrate( mbi, method, tag, elems, inttime );
360  CHK( rval, "time_integrate" );
361 
362  return rval;
363 }
364 
365 int main( int argc, char* argv[] )
366 {
367  // parse options
368  ProgOptions opts;
369  std::string filename;
370  int dim = 3;
371  parse_options( opts, dim, filename );
372  opts.parseCommandLine( argc, argv );
373  if( filename.empty() )
374  CHK( MB_FAILURE, "No file specified" );
375  else if( dim < 1 || dim > 3 )
376  CHK( MB_FAILURE, "Dimension must be > 0 and <= 3" );
377 
378  // read mesh file & gather element handles
379  Core mbi;
380  Range elems;
381  int nv;
382  EntityType tp;
383  ErrorCode rval = get_ents( mbi, filename, dim, elems, tp, nv );
384  if( MB_SUCCESS != rval ) return 1;
385 
386  // construct (random) parameter values for queries
387  unsigned int ne = elems.size();
388  std::vector< CartVect > params( ne ), coords( ne );
389  srand( time( NULL ) );
390  for( unsigned int i = 0; i < ne; i++ )
391  {
392  params[i][0] = -1 + 2 * ( (double)rand() ) / RAND_MAX;
393  if( dim > 1 ) params[i][1] = -1 + 2 * ( (double)rand() ) / RAND_MAX;
394  if( dim > 2 ) params[i][2] = -1 + 2 * ( (double)rand() ) / RAND_MAX;
395  }
396 
397  // put random field on vertices
398  Tag tag;
399  rval = put_random_field( mbi, tag, elems );
400  CHK( rval, "put_random_field" );
401  double evtime[2], retime[2], jactime[2],
402  inttime[2]; // initializations to avoid compiler warning
403 
404  rval = elem_evals( &mbi, ELEMEVAL, elems, tag, params, coords, evtime[0], retime[0], jactime[0], inttime[0] );
405  CHK( rval, "new elem_evals" );
406 
407  rval = elem_evals( &mbi, ELEMUTIL, elems, tag, params, coords, evtime[1], retime[1], jactime[1], inttime[1] );
408  CHK( rval, "old elem_evals" );
409 
410  std::cout << filename << ": " << elems.size() << " " << CN::EntityTypeName( tp ) << " elements, " << nv
411  << " vertices per element." << std::endl
412  << std::endl;
413  std::cout << "New, old element evaluation code:" << std::endl;
414  std::cout << "Evaluation type, time, time per element:" << std::endl;
415  std::cout << " New Old (New/Old)*100" << std::endl;
416  std::cout << "Forward evaluation " << evtime[0] << ", " << evtime[0] / elems.size() << " " << evtime[1] << ", "
417  << evtime[1] / elems.size() << " " << ( evtime[0] / ( evtime[1] ? evtime[1] : 1 ) ) * 100.0
418  << std::endl;
419  std::cout << "Reverse evaluation " << retime[0] << ", " << retime[0] / elems.size() << " " << retime[1] << ", "
420  << retime[1] / elems.size() << " " << ( retime[0] / ( retime[1] ? retime[1] : 1 ) ) * 100.0
421  << std::endl;
422  std::cout << "Jacobian " << jactime[0] << ", " << jactime[0] / elems.size() << " " << jactime[1]
423  << ", " << jactime[1] / elems.size() << " " << ( jactime[0] / ( jactime[1] ? jactime[1] : 1 ) ) * 100.0
424  << std::endl;
425  std::cout << "Integration " << inttime[0] << ", " << inttime[0] / elems.size() << " " << inttime[1]
426  << ", " << inttime[1] / elems.size() << " " << ( inttime[0] / ( inttime[1] ? inttime[1] : 1 ) ) * 100.0
427  << std::endl;
428 }
429 
430 #if defined( _MSC_VER ) || defined( __MINGW32__ )
431 double mytime2( double& tot_time, double& utime, double& stime, long& imem, long& rmem )
432 {
433  utime = (double)clock() / CLOCKS_PER_SEC;
434  tot_time = stime = 0;
435  imem = rmem = 0;
436  return tot_time;
437 }
438 #else
439 double mytime2( double& tot_time, double& utime, double& stime, long& imem, long& rmem )
440 {
441  struct rusage r_usage;
442  getrusage( RUSAGE_SELF, &r_usage );
443  utime = (double)r_usage.ru_utime.tv_sec + ( (double)r_usage.ru_utime.tv_usec / 1.e6 );
444  stime = (double)r_usage.ru_stime.tv_sec + ( (double)r_usage.ru_stime.tv_usec / 1.e6 );
445  tot_time = utime + stime;
446 #ifndef LINUX
447  imem = r_usage.ru_idrss;
448  rmem = r_usage.ru_maxrss;
449 #else
450  system( "ps o args,drs,rss | grep perf | grep -v grep" ); // RedHat 9.0 doesnt fill in actual
451  // memory data
452  imem = rmem = 0;
453 #endif
454  return tot_time;
455 }
456 #endif
457 double mytime()
458 {
459  double ttime, utime, stime;
460  long imem, rmem;
461  return mytime2( ttime, utime, stime, imem, rmem );
462 }