Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ssn_test.cpp
Go to the documentation of this file.
1 // A test file for Subset Normalization
2 #include "moab/ParallelComm.hpp"
4 #include "moab/Core.hpp"
5 #include "moab/FileOptions.hpp"
6 #include "ReadParallel.hpp"
7 #include "Coupler.hpp"
8 #include "DebugOutput.hpp"
9 #include "ElemUtil.hpp"
10 #include <iostream>
11 #include <iomanip>
12 #include <sstream>
13 #include <cstring>
14 #include <cstdlib>
15 extern "C" {
16 #include "moab/FindPtFuncs.h"
17 }
18 #include "moab/TupleList.hpp"
19 #include "moab/gs.hpp"
20 #include "moab/Types.hpp"
21 #ifndef IS_BUILDING_MB
22 #define IS_BUILDING_MB
23 #include "Internals.hpp"
24 #undef IS_BUILDING_MB
25 #else
26 #include "Internals.hpp"
27 #endif
28 
29 using namespace moab;
30 
31 bool debug = true;
32 
33 // Forward declarations
34 void get_file_options( int argc,
35  char** argv,
36  std::vector< const char* >& filenames,
37  std::string& norm_tag,
38  std::vector< const char* >& tag_names,
39  std::vector< const char* >& tag_values,
40  std::string& file_opts,
41  int* err );
42 
43 void print_tuples( TupleList* tlp );
44 
46  std::vector< std::vector< EntityHandle > >& groups,
47  Tag& norm_hdl,
48  Coupler::IntegType integ_type );
49 
50 double const_field( double x, double y, double z );
51 double field_1( double x, double y, double z );
52 double field_2( double x, double y, double z );
53 double field_3( double x, double y, double z );
54 double physField( double x, double y, double z );
55 
57 int pack_tuples( TupleList* tl, void** ptr );
58 void unpack_tuples( void* ptr, TupleList** tlp );
59 
60 //
61 // Start of main test program
62 //
63 int main( int argc, char** argv )
64 {
65  // need to init MPI first, to tell how many procs and rank
66  // Used since Coupler is a parallel code. The Coupler may be run
67  // in parallel or serial mode but will need MPI either way.
68  int err = MPI_Init( &argc, &argv );
69 
70  // Print usage if not enough arguments
71  if( argc < 3 )
72  {
73  std::cerr << "Usage: ";
74  std::cerr << argv[0] << " <nfiles> <fname1> ... <fnamen> <norm_tag> <tag_select_opts> <file_opts>" << std::endl;
75  std::cerr << "nfiles : number of mesh files" << std::endl;
76  std::cerr << "fname1...fnamen : mesh files" << std::endl;
77  std::cerr << "norm_tag : name of tag to normalize across meshes" << std::endl;
78  std::cerr << "tag_select_opts : quoted string of tags and values for subset selection, "
79  "e.g. \"TAG1=VAL1;TAG2=VAL2;TAG3;TAG4\""
80  << std::endl;
81  std::cerr << "file_opts : quoted string of parallel file read options, e.g. "
82  "\"OPTION1=VALUE1;OPTION2;OPTION3=VALUE3\""
83  << std::endl;
84 
85  err = integrate_scalar_field_test();MB_CHK_SET_ERR( (ErrorCode)err, "Integrate scalar field test failed" );
86 
87  err = MPI_Finalize();
88 
89  return err;
90  }
91 
92  int nprocs, rank;
93  err = MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
94  assert( MPI_SUCCESS == err );
95  err = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
96  assert( MPI_SUCCESS == err );
97 
98  // Create an ofstream to write output. One file each for each proc.
99  std::stringstream fname;
100  fname << argv[0] << rank << ".out";
101  if( !std::freopen( fname.str().c_str(), "a", stdout ) ) return -1;
102  if( !std::freopen( fname.str().c_str(), "a", stderr ) ) return -1;
103 
104  // Create the moab instance
105  Interface* mbi = new Core();
106  if( NULL == mbi )
107  {
108  std::cerr << "MOAB constructor failed" << std::endl;
109  return 1;
110  }
111 
112  // Get the input options
113  std::cout << "Getting options..." << std::endl;
114  std::vector< const char* > filenames;
115  std::vector< const char* > tagNames;
116  std::vector< const char* > tagValues;
117  std::string normTag, fileOpts;
118  get_file_options( argc, argv, filenames, normTag, tagNames, tagValues, fileOpts, &err );MB_CHK_SET_ERR( (ErrorCode)err, "get_file_options failed" );
119 
120  // Print out the input parameters
121  std::cout << " Input Parameters - " << std::endl;
122  std::cout << " Filenames: ";
123  for( std::vector< const char* >::iterator it = filenames.begin(); it != filenames.end(); ++it )
124  std::cout << *it << " ";
125  std::cout << std::endl;
126  std::cout << " Norm Tag: " << normTag << std::endl;
127  std::cout << " Selection Data: NumNames=" << tagNames.size() << " NumValues=" << tagValues.size() << std::endl;
128  std::cout << " TagNames TagValues " << std::endl;
129  std::cout << " -------------------- --------------------" << std::endl;
130  std::vector< const char* >::iterator nameIt = tagNames.begin();
131  std::vector< const char* >::iterator valIt = tagValues.begin();
132  std::cout << std::setiosflags( std::ios::left );
133  for( ; nameIt != tagNames.end(); ++nameIt )
134  {
135  std::cout << " " << std::setw( 20 ) << *nameIt;
136  if( *valIt != 0 )
137  {
138  std::cout << " " << std::setw( 20 ) << *( (int*)( *valIt ) ) << std::endl;
139  ++valIt;
140  }
141  else
142  std::cout << " NULL " << std::endl;
143  }
144  std::cout << std::resetiosflags( std::ios::left );
145  std::cout << " File Options: " << fileOpts << std::endl;
146 
147  // Read in mesh(es)
148  std::cout << "Reading mesh file(s)..." << std::endl;
149  std::vector< ParallelComm* > pcs( filenames.size() );
150  std::vector< ReadParallel* > rps( filenames.size() );
151 
152  // allocate root sets for each mesh for moab
153  std::vector< EntityHandle > roots( filenames.size() );
154 
155  ErrorCode result;
156  for( unsigned int i = 0; i < filenames.size(); i++ )
157  {
158  pcs[i] = new ParallelComm( mbi, MPI_COMM_WORLD );
159  rps[i] = new ReadParallel( mbi, pcs[i] );
160 
161  result = mbi->create_meshset( MESHSET_SET, roots[i] );
162 
163  MB_CHK_SET_ERR( result, "Creating root set failed" );
164  result = rps[i]->load_file( filenames[i], &roots[i], FileOptions( fileOpts.c_str() ) );MB_CHK_SET_ERR( result, "load_file failed" );
165  }
166 
167  // Initialize the debug object for Range printing
168  DebugOutput debugOut( "ssn_test-", std::cerr );
169  debugOut.set_rank( rank );
170  debugOut.set_verbosity( 10 );
171 
172  // Output what is in root sets
173  for( unsigned int k = 0; k < filenames.size(); k++ )
174  {
175 
176  Range rootRg;
177  result = mbi->get_entities_by_handle( roots[k], rootRg );MB_CHK_SET_ERR( result, "can't get entities" );
178  debugOut.print( 2, "Root set entities: ", rootRg );
179  rootRg.clear();
180 
181  Range partRg;
182  pcs[k]->get_part_entities( partRg );
183  debugOut.print( 2, "Partition entities: ", partRg );
184  partRg.clear();
185  }
186 
187  // source is 1st mesh, target is 2nd mesh
188  Range src_elems, targ_elems;
189 
190  // ******************************
191  std::cout << "********** Create Coupler **********" << std::endl;
192  // Create a coupler
193  std::cout << "Creating Coupler..." << std::endl;
194  Coupler mbc( mbi, pcs[0], src_elems, 0 );
195 
196  // Get tag handles for passed in tags
197  std::cout << "Getting tag handles..." << std::endl;
198  int numTagNames = tagNames.size();
199 
200  std::vector< Tag > tagHandles( numTagNames );
201  int iTags = 0;
202  while( iTags < numTagNames )
203  {
204  std::cout << "Getting handle for " << tagNames[iTags] << std::endl;
205  result = mbi->tag_get_handle( tagNames[iTags], tagHandles[iTags] );MB_CHK_SET_ERR( result, "Retrieving tag handles failed" );
206  iTags++;
207  }
208 
209  // ******************************
210  std::cout << "********** Test create_tuples **********" << std::endl;
211  // First get some EntitySets for Mesh 1 and Mesh 2
212  {
213 
214  Range entsets1, entsets2;
215  result = mbi->get_entities_by_type_and_tag( roots[0], MBENTITYSET, &tagHandles[0],
216  (const void* const*)&tagValues[0], tagHandles.size(), entsets1,
217  Interface::INTERSECT ); // recursive is false
218  MB_CHK_SET_ERR( result, "sets: get_entities_by_type_and_tag failed on Mesh 1." );
219 
220  // Create tuple_list for each mesh's
221  std::cout << "Creating tuples for mesh 1..." << std::endl;
222  TupleList* m1TagTuples = NULL;
223  err = mbc.create_tuples( entsets1, &tagHandles[0], tagHandles.size(), &m1TagTuples );MB_CHK_SET_ERR( (ErrorCode)err, "create_tuples failed" );
224 
225  std::cout << " create_tuples returned" << std::endl;
226  print_tuples( m1TagTuples );
227 
228  result = mbi->get_entities_by_type_and_tag( roots[1], MBENTITYSET, &tagHandles[0],
229  (const void* const*)&tagValues[0], tagHandles.size(), entsets2,
230  Interface::INTERSECT ); // recursive is false
231  MB_CHK_SET_ERR( result, "sets: get_entities_by_type_and_tag failed on Mesh 2." );
232 
233  std::cout << "Creating tuples for mesh 2..." << std::endl;
234  TupleList* m2TagTuples = NULL;
235  err = mbc.create_tuples( entsets2, (Tag*)( &tagHandles[0] ), tagHandles.size(), &m2TagTuples );MB_CHK_SET_ERR( (ErrorCode)err, "create_tuples failed" );
236 
237  std::cout << " create_tuples returned" << std::endl;
238  print_tuples( m2TagTuples );
239 
240  // ******************************
241  std::cout << "********** Test consolidate_tuples **********" << std::endl;
242  // In this serial version we only have the tuples from Mesh 1 and Mesh 2.
243  // Just consolidate those for the test.
244  std::cout << "Consolidating tuple_lists for Mesh 1 and Mesh 2..." << std::endl;
245  TupleList** tplp_arr = (TupleList**)malloc( 2 * sizeof( TupleList* ) );
246  TupleList* unique_tpl = NULL;
247  tplp_arr[0] = m1TagTuples;
248  tplp_arr[1] = m2TagTuples;
249 
250  err = mbc.consolidate_tuples( tplp_arr, 2, &unique_tpl );MB_CHK_SET_ERR( (ErrorCode)err, "consolidate_tuples failed" );
251  std::cout << " consolidate_tuples returned" << std::endl;
252  print_tuples( unique_tpl );
253  }
254 
255  // ******************************
256  std::cout << "********** Test get_matching_entities **********" << std::endl;
257  std::vector< std::vector< EntityHandle > > m1EntitySets;
258  std::vector< std::vector< EntityHandle > > m1EntityGroups;
259  std::vector< std::vector< EntityHandle > > m2EntitySets;
260  std::vector< std::vector< EntityHandle > > m2EntityGroups;
261 
262  // Get matching entities for Mesh 1
263  std::cout << "Get matching entities for mesh 1..." << std::endl;
264  err = mbc.get_matching_entities( roots[0], &tagHandles[0], &tagValues[0], tagHandles.size(), &m1EntitySets,
265  &m1EntityGroups );MB_CHK_SET_ERR( (ErrorCode)err, "get_matching_entities failed" );
266 
267  std::cout << " get_matching_entities returned " << m1EntityGroups.size() << " entity groups" << std::endl;
268 
269  // Print out the data in the vector of vectors
270  std::vector< std::vector< EntityHandle > >::iterator iter_esi;
271  std::vector< std::vector< EntityHandle > >::iterator iter_egi;
272  std::vector< EntityHandle >::iterator iter_esj;
273  std::vector< EntityHandle >::iterator iter_egj;
274  Range entSetRg;
275  int icnt;
276  for( iter_egi = m1EntityGroups.begin(), iter_esi = m1EntitySets.begin(), icnt = 1;
277  ( iter_egi != m1EntityGroups.end() ) && ( iter_esi != m1EntitySets.end() ); ++iter_egi, ++iter_esi, icnt++ )
278  {
279  std::cout << " EntityGroup(" << icnt << ") = ";
280  std::cout.flush();
281  entSetRg.clear();
282  for( iter_egj = ( *iter_egi ).begin(); iter_egj != ( *iter_egi ).end(); ++iter_egj )
283  entSetRg.insert( (EntityHandle)*iter_egj );
284  debugOut.print( 2, "Mesh1 matching Entities: ", entSetRg );
285  std::cout.flush();
286 
287  std::cout << " EntitySet(" << icnt << ") = ";
288  std::cout.flush();
289  entSetRg.clear();
290  for( iter_esj = ( *iter_esi ).begin(); iter_esj != ( *iter_esi ).end(); ++iter_esj )
291  entSetRg.insert( (EntityHandle)*iter_esj );
292  debugOut.print( 2, "Mesh1 matching EntitySets: ", entSetRg );
293  std::cout.flush();
294  }
295 
296  // Get matching entities for Mesh 2
297  std::cout << "Get matching entities for mesh 2..." << std::endl;
298  err = mbc.get_matching_entities( roots[1], &tagHandles[0], &tagValues[0], tagHandles.size(), &m2EntitySets,
299  &m2EntityGroups );MB_CHK_SET_ERR( (ErrorCode)err, "get_matching_entities failed" );
300 
301  std::cout << " get_matching_entities returned " << m2EntityGroups.size() << " entity groups" << std::endl;
302  for( iter_egi = m2EntityGroups.begin(), iter_esi = m2EntitySets.begin(), icnt = 1;
303  ( iter_egi != m2EntityGroups.end() ) && ( iter_esi != m2EntitySets.end() ); ++iter_egi, ++iter_esi, icnt++ )
304  {
305  std::cout << " EntityGroup(" << icnt << ") = ";
306  std::cout.flush();
307  entSetRg.clear();
308  for( iter_egj = ( *iter_egi ).begin(); iter_egj != ( *iter_egi ).end(); ++iter_egj )
309  entSetRg.insert( (EntityHandle)*iter_egj );
310  debugOut.print( 2, "Mesh2 matching Entities: ", entSetRg );
311  std::cout.flush();
312 
313  std::cout << " EntitySet(" << icnt << ") = ";
314  std::cout.flush();
315  entSetRg.clear();
316  for( iter_esj = ( *iter_esi ).begin(); iter_esj != ( *iter_esi ).end(); ++iter_esj )
317  entSetRg.insert( (EntityHandle)*iter_esj );
318  debugOut.print( 2, "Mesh2 matching EntitySets: ", entSetRg );
319  std::cout.flush();
320  }
321 
322  if( debug )
323  {
324  // ******************************
325  std::cout << "********** Test print_tuples **********" << std::endl;
326  // temporary test function
327  std::cout << "Testing print_tuples..." << std::endl;
328 
329  TupleList test_tuple;
330  int num_ints = 3, num_longs = 2, num_ulongs = 4, num_reals = 6, num_rows = 10;
331 
332  std::cout << " print of test_tuples zero init..." << std::endl;
333  test_tuple.initialize( 0, 0, 0, 0, 0 );
334 
335  test_tuple.enableWriteAccess();
336 
337  print_tuples( &test_tuple );
338 
339  std::cout << " print of test_tuples after setting n to 10..." << std::endl;
340  test_tuple.set_n( 10 );
341  print_tuples( &test_tuple );
342 
343  test_tuple.initialize( num_ints, num_longs, num_ulongs, num_reals, num_rows );
344  std::cout << " print of test_tuples after init..." << std::endl;
345  print_tuples( &test_tuple );
346 
347  std::cout << " print of test_tuples after setting n to 10..." << std::endl;
348  test_tuple.set_n( 10 );
349  print_tuples( &test_tuple );
350 
351  for( int i = 0; i < num_rows; i++ )
352  {
353  int j;
354  for( j = 0; j < num_ints; j++ )
355  test_tuple.vi_wr[i * num_ints + j] = (int)( ( j + 1 ) * ( i + 1 ) );
356 
357  for( j = 0; j < num_longs; j++ )
358  test_tuple.vl_wr[i * num_longs + j] = (int)( ( j + 1 ) * ( i + 1 ) );
359 
360  for( j = 0; j < num_ulongs; j++ )
361  test_tuple.vul_wr[i * num_ulongs + j] = (int)( ( j + 1 ) * ( i + 1 ) );
362 
363  for( j = 0; j < num_reals; j++ )
364  test_tuple.vr_wr[i * num_reals + j] = (int)( ( j + 1 ) * ( i + 1 ) + ( j * 0.01 ) );
365  }
366  std::cout << " print of test_tuples after filling with data..." << std::endl;
367  print_tuples( &test_tuple );
368 
369  // ******************************
370  std::cout << "********** Test pack_tuples and unpack_tuples **********" << std::endl;
371  void* mp_buf;
372  int buf_sz;
373  if( rank == 0 )
374  {
375  buf_sz = pack_tuples( &test_tuple, &mp_buf );
376  }
377 
378  // Send buffer size
379  err = MPI_Bcast( &buf_sz, 1, MPI_INT, 0, MPI_COMM_WORLD );
380 
381  if( err != MPI_SUCCESS )
382  {
383  std::cerr << "MPI_Bcast of buffer size failed" << std::endl;
384  return -1;
385  }
386 
387  // Allocate a buffer in the other procs
388  if( rank != 0 )
389  {
390  mp_buf = malloc( buf_sz * sizeof( uint ) );
391  }
392 
393  err = MPI_Bcast( mp_buf, buf_sz * sizeof( uint ), MPI_UNSIGNED_CHAR, 0, MPI_COMM_WORLD );
394  if( err != MPI_SUCCESS )
395  {
396  std::cerr << "MPI_Bcast of buffer failed" << std::endl;
397  return -1;
398  }
399 
400  TupleList* rcv_tuples;
401  unpack_tuples( mp_buf, &rcv_tuples );
402 
403  std::cout << " print of rcv_tuples after unpacking from MPI_Bcast..." << std::endl;
404  print_tuples( rcv_tuples );
405  }
406 
407  err = integrate_scalar_field_test();MB_CHK_SET_ERR( (ErrorCode)err, "Failure in integrating a scalar_field" );
408 
409  // ******************************
410  std::cout << "********** Test get_group_integ_vals **********" << std::endl;
411  std::cout << "Get group integrated field values..." << std::endl;
412 
413  // print the field values at the vertices before change.
414  std::cout << " print vertex field values first:" << std::endl;
415  Tag norm_hdl;
416  result = mbi->tag_get_handle( normTag.c_str(), norm_hdl );MB_CHK_SET_ERR( (ErrorCode)err, "Failed to get tag handle." );
417 
418  Coupler::IntegType integ_type = Coupler::VOLUME;
419  // Mesh 1 field values
420  std::cout << " Original entity vertex field values (mesh 1): " << std::endl;
421  print_vertex_fields( mbi, m1EntityGroups, norm_hdl, integ_type );
422 
423  // Mesh 2 field values
424  std::cout << " Original entity vertex field values (mesh 2): " << std::endl;
425  print_vertex_fields( mbi, m2EntityGroups, norm_hdl, integ_type );
426 
427  // Get the field values
428  std::vector< double >::iterator iter_ivals;
429 
430  std::cout << "Get group integrated field values for mesh 1..." << std::endl;
431  std::vector< double > m1IntegVals( m1EntityGroups.size() );
432  err = mbc.get_group_integ_vals( m1EntityGroups, m1IntegVals, normTag.c_str(), 4, integ_type );MB_CHK_SET_ERR( (ErrorCode)err, "Failed to get the Mesh 1 group integration values." );
433  std::cout << "Mesh 1 integrated field values(" << m1IntegVals.size() << "): ";
434  for( iter_ivals = m1IntegVals.begin(); iter_ivals != m1IntegVals.end(); ++iter_ivals )
435  {
436  std::cout << ( *iter_ivals ) << " ";
437  }
438  std::cout << std::endl;
439 
440  std::cout << "Get group integrated field values for mesh 2..." << std::endl;
441  std::vector< double > m2IntegVals( m2EntityGroups.size() );
442  err = mbc.get_group_integ_vals( m2EntityGroups, m2IntegVals, normTag.c_str(), 4, integ_type );MB_CHK_SET_ERR( (ErrorCode)err, "Failed to get the Mesh 2 group integration values." );
443  std::cout << "Mesh 2 integrated field values(" << m2IntegVals.size() << "): ";
444  for( iter_ivals = m2IntegVals.begin(); iter_ivals != m2IntegVals.end(); ++iter_ivals )
445  {
446  std::cout << ( *iter_ivals ) << " ";
447  }
448  std::cout << std::endl;
449 
450  // ******************************
451  std::cout << "********** Test apply_group_norm_factors **********" << std::endl;
452  // Make the norm factors by inverting the integration values.
453  double val;
454  for( unsigned int i = 0; i < m1IntegVals.size(); i++ )
455  {
456  val = m1IntegVals[i];
457  m1IntegVals[i] = 1 / val;
458  }
459 
460  for( unsigned int i = 0; i < m2IntegVals.size(); i++ )
461  {
462  val = m2IntegVals[i];
463  m2IntegVals[i] = 1 / val;
464  }
465 
466  std::cout << "Mesh 1 norm factors(" << m1IntegVals.size() << "): ";
467  for( iter_ivals = m1IntegVals.begin(); iter_ivals != m1IntegVals.end(); ++iter_ivals )
468  {
469  std::cout << ( *iter_ivals ) << " ";
470  }
471  std::cout << std::endl;
472 
473  std::cout << "Mesh 2 norm factors(" << m2IntegVals.size() << "): ";
474  for( iter_ivals = m2IntegVals.begin(); iter_ivals != m2IntegVals.end(); ++iter_ivals )
475  {
476  std::cout << ( *iter_ivals ) << " ";
477  }
478  std::cout << std::endl;
479 
480  // Apply the factors and reprint the vertices
481  err = mbc.apply_group_norm_factor( m1EntitySets, m1IntegVals, normTag.c_str(), integ_type );MB_CHK_SET_ERR( (ErrorCode)err, "Failed to apply norm factors to Mesh 1." );
482 
483  err = mbc.apply_group_norm_factor( m2EntitySets, m2IntegVals, normTag.c_str(), integ_type );MB_CHK_SET_ERR( (ErrorCode)err, "Failed to apply norm factors to Mesh 2." );
484 
485  // Get the norm_tag_factor on the EntitySets
486  // Get the handle for the norm factor tag
487  Tag norm_factor_hdl;
488  std::string normFactor = normTag + "_normf";
489  result = mbi->tag_get_handle( normFactor.c_str(), norm_factor_hdl );MB_CHK_SET_ERR( result, "Failed to get norm factor tag handle." );
490 
491  // Mesh 1 values
492  std::cout << "Mesh 1 norm factors per EntitySet...";
493  for( iter_esi = m1EntitySets.begin(); iter_esi != m1EntitySets.end(); ++iter_esi )
494  {
495  for( iter_esj = ( *iter_esi ).begin(); iter_esj != ( *iter_esi ).end(); ++iter_esj )
496  {
497  double data = 0;
498  EntityHandle eh = *iter_esj;
499  result = mbi->tag_get_data( norm_factor_hdl, &eh, 1, &data );MB_CHK_SET_ERR( result, "Failed to get tag data." );
500  std::cout << data << ", ";
501  }
502  }
503  std::cout << std::endl;
504 
505  // Mesh 2 values
506  std::cout << "Mesh 2 norm factors per EntitySet...";
507  for( iter_esi = m2EntitySets.begin(); iter_esi != m2EntitySets.end(); ++iter_esi )
508  {
509  for( iter_esj = ( *iter_esi ).begin(); iter_esj != ( *iter_esi ).end(); ++iter_esj )
510  {
511  double data = 0;
512  EntityHandle eh = *iter_esj;
513  result = mbi->tag_get_data( norm_factor_hdl, &eh, 1, &data );MB_CHK_SET_ERR( result, "Failed to get tag data." );
514  std::cout << data << ", ";
515  }
516  }
517  std::cout << std::endl;
518 
519  // ******************************
520  std::cout << "********** Test normalize_subset **********" << std::endl;
521  // Now call the Coupler::normalize_subset routine and see if we get an error.
522  std::cout << "Running Coupler::normalize_subset() on mesh 1" << std::endl;
523  err = mbc.normalize_subset( (EntityHandle)roots[0], normTag.c_str(), &tagNames[0], numTagNames, &tagValues[0],
524  Coupler::VOLUME, 4 );MB_CHK_SET_ERR( (ErrorCode)err, "Failure in call to Coupler::normalize_subset() on mesh 1" );
525 
526  // Print the normFactor on each EntitySet after the above call.
527  // Mesh 1 values
528  std::cout << "Mesh 1 norm factors per EntitySet...";
529  for( iter_esi = m1EntitySets.begin(); iter_esi != m1EntitySets.end(); ++iter_esi )
530  {
531  for( iter_esj = ( *iter_esi ).begin(); iter_esj != ( *iter_esi ).end(); ++iter_esj )
532  {
533  double data = 0;
534  EntityHandle eh = *iter_esj;
535  result = mbi->tag_get_data( norm_factor_hdl, &eh, 1, &data );MB_CHK_SET_ERR( result, "Failed to get tag data." );
536  std::cout << data << ", ";
537  }
538  }
539  std::cout << std::endl;
540 
541  std::cout << "Running Coupler::normalize_subset() on mesh 2" << std::endl;
542  err = mbc.normalize_subset( (EntityHandle)roots[1], normTag.c_str(), &tagNames[0], numTagNames, &tagValues[0],
543  Coupler::VOLUME, 4 );MB_CHK_SET_ERR( (ErrorCode)err, "Failure in call to Coupler::normalize_subset() on mesh 2" );
544 
545  // Mesh 2 values
546  std::cout << "Mesh 2 norm factors per EntitySet...";
547  for( iter_esi = m2EntitySets.begin(); iter_esi != m2EntitySets.end(); ++iter_esi )
548  {
549  for( iter_esj = ( *iter_esi ).begin(); iter_esj != ( *iter_esi ).end(); ++iter_esj )
550  {
551  double data = 0;
552  EntityHandle eh = *iter_esj;
553  result = mbi->tag_get_data( norm_factor_hdl, &eh, 1, &data );MB_CHK_SET_ERR( result, "Failed to get tag data." );
554 
555  std::cout << data << ", ";
556  }
557  }
558  std::cout << std::endl;
559 
560  // Done, cleanup
561  std::cout << "********** ssn_test DONE! **********" << std::endl;
562  MPI_Finalize();
563  return 0;
564 }
565 
567 {
568  // ******************************
569  std::cout << "********** Test moab::Element::Map::integrate_scalar_field **********" << std::endl;
570  // Create a simple hex centered at 0,0,0 with sides of length 2.
571  std::vector< CartVect > biunit_cube( 8 );
572  biunit_cube[0] = CartVect( -1, -1, -1 );
573  biunit_cube[1] = CartVect( 1, -1, -1 );
574  biunit_cube[2] = CartVect( 1, 1, -1 );
575  biunit_cube[3] = CartVect( -1, 1, -1 );
576  biunit_cube[4] = CartVect( -1, -1, 1 );
577  biunit_cube[5] = CartVect( 1, -1, 1 );
578  biunit_cube[6] = CartVect( 1, 1, 1 );
579  biunit_cube[7] = CartVect( -1, 1, 1 );
580 
581  std::vector< CartVect > zerobase_cube( 8 );
582  zerobase_cube[0] = CartVect( 0, 0, 0 );
583  zerobase_cube[1] = CartVect( 2, 0, 0 );
584  zerobase_cube[2] = CartVect( 2, 2, 0 );
585  zerobase_cube[3] = CartVect( 0, 2, 0 );
586  zerobase_cube[4] = CartVect( 0, 0, 2 );
587  zerobase_cube[5] = CartVect( 2, 0, 2 );
588  zerobase_cube[6] = CartVect( 2, 2, 2 );
589  zerobase_cube[7] = CartVect( 0, 2, 2 );
590 
591  // Calculate field values at the corners of both cubes
592  double bcf[8], bf1[8], bf2[8], bf3[8], zcf[8], zf1[8], zf2[8], zf3[8];
593  for( int i = 0; i < 8; i++ )
594  {
595  bcf[i] = const_field( biunit_cube[i][0], biunit_cube[i][1], biunit_cube[i][2] );
596  bf1[i] = field_1( biunit_cube[i][0], biunit_cube[i][1], biunit_cube[i][2] );
597  bf2[i] = field_2( biunit_cube[i][0], biunit_cube[i][1], biunit_cube[i][2] );
598  bf3[i] = field_3( biunit_cube[i][0], biunit_cube[i][1], biunit_cube[i][2] );
599 
600  zcf[i] = const_field( zerobase_cube[i][0], zerobase_cube[i][1], zerobase_cube[i][2] );
601  zf1[i] = field_1( zerobase_cube[i][0], zerobase_cube[i][1], zerobase_cube[i][2] );
602  zf2[i] = field_2( zerobase_cube[i][0], zerobase_cube[i][1], zerobase_cube[i][2] );
603  zf3[i] = field_3( zerobase_cube[i][0], zerobase_cube[i][1], zerobase_cube[i][2] );
604  }
605 
606  std::cout << "Integrated values:" << std::endl;
607 
608  try
609  {
610  double field_const1, field_const2;
611  double field_linear1, field_linear2;
612  double field_quad1, field_quad2;
613  double field_cubic1, field_cubic2;
614 
615  int ipoints = 0;
616  Element::LinearHex biunit_hexMap( biunit_cube );
617  Element::LinearHex zerobase_hexMap( zerobase_cube );
618 
619  field_const1 = biunit_hexMap.integrate_scalar_field( bcf );
620  field_const2 = zerobase_hexMap.integrate_scalar_field( zcf );
621  std::cout << " binunit_cube, const_field(num_pts=" << ipoints << "): field_val=" << field_const1
622  << std::endl;
623  std::cout << " zerobase_cube, const_field(num_pts=" << ipoints << "): field_val=" << field_const2
624  << std::endl;
625 
626  field_linear1 = biunit_hexMap.integrate_scalar_field( bf1 );
627  field_linear2 = zerobase_hexMap.integrate_scalar_field( zf1 );
628  std::cout << " binunit_cube, field_1(num_pts=" << ipoints << "): field_val=" << field_linear1 << std::endl;
629  std::cout << " zerobase_cube, field_1(num_pts=" << ipoints << "): field_val=" << field_linear2 << std::endl;
630 
631  field_quad1 = biunit_hexMap.integrate_scalar_field( bf2 );
632  field_quad2 = zerobase_hexMap.integrate_scalar_field( zf2 );
633  std::cout << " binunit_cube, field_2(num_pts=" << ipoints << "): field_val=" << field_quad1 << std::endl;
634  std::cout << " zerobase_cube, field_2(num_pts=" << ipoints << "): field_val=" << field_quad2 << std::endl;
635 
636  field_cubic1 = biunit_hexMap.integrate_scalar_field( bf3 );
637  field_cubic2 = zerobase_hexMap.integrate_scalar_field( zf3 );
638  std::cout << " binunit_cube, field_3(num_pts=" << ipoints << "): field_val=" << field_cubic1 << std::endl;
639  std::cout << " zerobase_cube, field_3(num_pts=" << ipoints << "): field_val=" << field_cubic2 << std::endl;
640  }
642  {
643  MB_CHK_SET_ERR( MB_FAILURE, "Failed to set vertices on Element::Map." );
644  }
646  {
647  MB_CHK_SET_ERR( MB_FAILURE, "Failed to get inverse evaluation of coordinate on Element::Map." );
648  }
649  return MB_SUCCESS;
650 }
651 
652 // Function to parse input parameters
653 void get_file_options( int argc,
654  char** argv,
655  std::vector< const char* >& filenames,
656  std::string& normTag,
657  std::vector< const char* >& tagNames,
658  std::vector< const char* >& tagValues,
659  std::string& fileOpts,
660  int* err )
661 {
662  int npos = 1;
663 
664  // get number of files
665  int nfiles = atoi( argv[npos++] );
666 
667  // get mesh filenames
668  filenames.resize( nfiles );
669  for( int i = 0; i < nfiles; i++ )
670  filenames[i] = argv[npos++];
671 
672  // get normTag
673  if( npos < argc )
674  normTag = argv[npos++];
675  else
676  {
677  std::cerr << "Insufficient parameters: norm_tag missing" << std::endl;
678  *err = 1;
679  return;
680  }
681 
682  // get tag selection options
683  if( npos < argc )
684  {
685  char* opts = argv[npos++];
686  // char sep1[1] = {';'};
687  // char sep2[1] = {'='};
688  bool end_vals_seen = false;
689  std::vector< char* > tmpTagOpts;
690 
691  // first get the options
692  for( char* i = strtok( opts, ";" ); i; i = strtok( 0, ";" ) )
693  {
694  if( debug ) std::cout << "get_file_options: i=" << i << std::endl;
695  tmpTagOpts.push_back( i );
696  }
697 
698  // parse out the name and val or just name.
699  for( unsigned int j = 0; j < tmpTagOpts.size(); j++ )
700  {
701  char* e = strtok( tmpTagOpts[j], "=" );
702  if( debug ) std::cout << "get_file_options: name=" << e << std::endl;
703  tagNames.push_back( e );
704  e = strtok( 0, "=" );
705  if( e != NULL )
706  {
707  if( debug ) std::cout << "get_file_options: val=" << e << std::endl;
708  // We have a value
709  if( end_vals_seen )
710  {
711  // ERROR we should not have a value after none are seen
712  std::cerr << "Incorrect parameters: new value seen after end of values" << std::endl;
713  *err = 1;
714  return;
715  }
716  // Otherwise get the value string from e and convert it to an int
717  int* valp = new int;
718  *valp = atoi( e );
719  tagValues.push_back( (const char*)valp );
720  }
721  else
722  {
723  // Otherwise there is no '=' so push a null on the list
724  end_vals_seen = true;
725  tagValues.push_back( (const char*)0 );
726  }
727  }
728  }
729  else
730  {
731  std::cerr << "Insufficient parameters: tag_select_opts missing" << std::endl;
732  *err = 1;
733  return;
734  }
735 
736  // get fileOpts
737  if( npos < argc )
738  fileOpts = argv[npos++];
739  else
740  {
741  std::cerr << "Insufficient parameters: file_opts missing" << std::endl;
742  *err = 1;
743  return;
744  }
745 }
746 
747 // Function to print out a tuple_list.
749 {
750  uint mi, ml, mul, mr;
751  tlp->getTupleSize( mi, ml, mul, mr );
752  std::cout << " tuple data: (n=" << tlp->get_n() << ")" << std::endl;
753  std::cout << " mi:" << mi << " ml:" << ml << " mul:" << mul << " mr:" << mr << std::endl;
754  std::cout << " [" << std::setw( 11 * mi ) << " int data"
755  << " |" << std::setw( 11 * ml ) << " long data"
756  << " |" << std::setw( 11 * mul ) << " Ulong data"
757  << " |" << std::setw( 11 * mr ) << " real data"
758  << " " << std::endl
759  << " ";
760  for( unsigned int i = 0; i < tlp->get_n(); i++ )
761  {
762  if( mi > 0 )
763  {
764  for( unsigned int j = 0; j < mi; j++ )
765  {
766  std::cout << std::setw( 10 ) << tlp->vi_rd[i * mi + j] << " ";
767  }
768  }
769  else
770  {
771  std::cout << " ";
772  }
773  std::cout << "| ";
774 
775  if( ml > 0 )
776  {
777  for( unsigned int j = 0; j < ml; j++ )
778  {
779  std::cout << std::setw( 10 ) << tlp->vl_rd[i * ml + j] << " ";
780  }
781  }
782  else
783  {
784  std::cout << " ";
785  }
786  std::cout << "| ";
787 
788  if( mul > 0 )
789  {
790  for( unsigned int j = 0; j < mul; j++ )
791  {
792  std::cout << std::setw( 10 ) << tlp->vul_rd[i * mul + j] << " ";
793  }
794  }
795  else
796  {
797  std::cout << " ";
798  }
799  std::cout << "| ";
800 
801  if( mr > 0 )
802  {
803  for( unsigned int j = 0; j < mr; j++ )
804  {
805  std::cout << std::setw( 10 ) << tlp->vr_rd[i * mr + j] << " ";
806  }
807  }
808  else
809  {
810  std::cout << " ";
811  }
812 
813  if( i + 1 < tlp->get_n() ) std::cout << std::endl << " ";
814  }
815  std::cout << "]" << std::endl;
816 }
817 
818 // Function to print vertex field values
820  std::vector< std::vector< EntityHandle > >& groups,
821  Tag& norm_hdl,
822  Coupler::IntegType integ_type )
823 {
824  int err = 0;
825  ErrorCode result;
826  std::vector< EntityHandle >::iterator iter_j;
827 
828  for( unsigned int i = 0; i < groups.size(); i++ )
829  {
830  std::cout << " Group - " << std::endl << " ";
831  for( iter_j = groups[i].begin(); iter_j != groups[i].end(); ++iter_j )
832  {
833  EntityHandle ehandle = ( *iter_j );
834  // Check that the entity in iter_j is of the same dimension as the
835  // integ_type we are performing
836  int j_type = mbi->dimension_from_handle( ehandle );
837 
838  if( ( integ_type == Coupler::VOLUME ) && ( j_type != 3 ) ) continue;
839 
840  // Retrieve the vertices from the element
841  const EntityHandle* conn = NULL;
842  int num_verts = 0;
843  result = mbi->get_connectivity( ehandle, conn, num_verts );
844  if( MB_SUCCESS != result ) return 1;
845  std::cout << std::fixed;
846  for( int iv = 0; iv < num_verts; iv++ )
847  {
848  double data = 0;
849  result = mbi->tag_get_data( norm_hdl, &conn[iv], 1, &data );
850  if( MB_SUCCESS != result ) return 1;
851  std::cout << std::setprecision( 8 ) << data << ", ";
852  }
853  std::cout << std::endl << " ";
854  }
855  std::cout << std::endl;
856  std::cout.unsetf( std::ios_base::floatfield ); // turn off fixed notation
857  }
858 
859  return err;
860 }
861 
862 // Function for a constant field value
863 double const_field( double /*x*/, double /*y*/, double /*z*/ )
864 {
865  // return 5.0/40.0;
866  return 5.0;
867 }
868 
869 // Functions for a some field values
870 double field_1( double x, double y, double z )
871 {
872  double f = fabs( x ) + fabs( y ) + fabs( z );
873  // return f/24.0;
874  return f;
875 }
876 
877 double field_2( double x, double y, double z )
878 {
879  double f = x * x + y * y + z * z;
880  // return f/32.0;
881  return f;
882 }
883 
884 double field_3( double x, double y, double z )
885 {
886  double f = 2 * x + 2 * y + 2 * z;
887  // return f/48.0;
888  return f;
889 }
890 
891 // Function used to create field on mesh for testing.
892 double physField( double x, double y, double z )
893 {
894  double out;
895 
896  // 1/r^2 decay from {0,0,0}
897 
898  out = x * x + y * y + z * z;
899  out += 1e-1; // clamp
900  out = 1 / out;
901 
902  return out;
903 }
904 
905 #define UINT_PER_X( X ) ( ( sizeof( X ) + sizeof( uint ) - 1 ) / sizeof( uint ) )
906 #define UINT_PER_REAL UINT_PER_X( realType )
907 #define UINT_PER_LONG UINT_PER_X( slong )
908 #define UINT_PER_ULONG UINT_PER_X( Ulong )
909 #define UINT_PER_UNSIGNED UINT_PER_X( unsigned )
910 
911 // Function for packing tuple_list
912 int pack_tuples( TupleList* tl, void** ptr )
913 {
914  uint mi, ml, mul, mr;
915  tl->getTupleSize( mi, ml, mul, mr );
916 
917  uint n = tl->get_n();
918 
919  int sz_buf = 1 + 4 * UINT_PER_UNSIGNED +
920  tl->get_n() * ( mi + ml * UINT_PER_LONG + mul * UINT_PER_LONG + mr * UINT_PER_REAL );
921 
922  uint* buf = (uint*)malloc( sz_buf * sizeof( uint ) );
923  *ptr = (void*)buf;
924 
925  // copy n
926  memcpy( buf, &n, sizeof( uint ) ), buf += 1;
927  // copy mi
928  memcpy( buf, &mi, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
929  // copy ml
930  memcpy( buf, &ml, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
931  // copy mul
932  memcpy( buf, &mul, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
933  // copy mr
934  memcpy( buf, &mr, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
935  // copy vi_rd
936  memcpy( buf, tl->vi_rd, tl->get_n() * mi * sizeof( sint ) ), buf += tl->get_n() * mi;
937  // copy vl_rd
938  memcpy( buf, tl->vl_rd, tl->get_n() * ml * sizeof( slong ) ), buf += tl->get_n() * ml * UINT_PER_LONG;
939  // copy vul_rd
940  memcpy( buf, tl->vul_rd, tl->get_n() * mul * sizeof( Ulong ) ), buf += tl->get_n() * mul * UINT_PER_ULONG;
941  // copy vr_rd
942  memcpy( buf, tl->vr_rd, tl->get_n() * mr * sizeof( realType ) ), buf += tl->get_n() * mr * UINT_PER_REAL;
943 
944  return sz_buf;
945 }
946 
947 // Function for packing tuple_list
948 void unpack_tuples( void* ptr, TupleList** tlp )
949 {
950  TupleList* tl = new TupleList();
951  *tlp = tl;
952 
953  uint nt;
954  unsigned mit, mlt, mult, mrt;
955  uint* buf = (uint*)ptr;
956 
957  // get n
958  memcpy( &nt, buf, sizeof( uint ) ), buf += 1;
959  // get mi
960  memcpy( &mit, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
961  // get ml
962  memcpy( &mlt, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
963  // get mul
964  memcpy( &mult, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
965  // get mr
966  memcpy( &mrt, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
967 
968  // initialize tl
969  tl->initialize( mit, mlt, mult, mrt, nt );
970  tl->enableWriteAccess();
971  tl->set_n( nt );
972 
973  uint mi, ml, mul, mr;
974  tl->getTupleSize( mi, ml, mul, mr );
975 
976  // get vi_wr
977  memcpy( tl->vi_wr, buf, tl->get_n() * mi * sizeof( sint ) ), buf += tl->get_n() * mi;
978  // get vl_wr
979  memcpy( tl->vl_wr, buf, tl->get_n() * ml * sizeof( slong ) ), buf += tl->get_n() * ml * UINT_PER_LONG;
980  // get vul_wr
981  memcpy( tl->vul_wr, buf, tl->get_n() * mul * sizeof( Ulong ) ), buf += tl->get_n() * mul * UINT_PER_ULONG;
982  // get vr_wr
983  memcpy( tl->vr_wr, buf, tl->get_n() * mr * sizeof( realType ) ), buf += tl->get_n() * mr * UINT_PER_REAL;
984 
985  tl->disableWriteAccess();
986 
987  return;
988 }