 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 // A test file for Subset Normalization
2 #include "moab/ParallelComm.hpp"
3 #include "MBParallelConventions.h"
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
45 int print_vertex_fields( Interface* mbi,
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
56 ErrorCode integrate_scalar_field_test();
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
566 ErrorCode integrate_scalar_field_test()
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 }
641 catch( moab::Element::Map::ArgError )
642 {
643 MB_CHK_SET_ERR( MB_FAILURE, "Failed to set vertices on Element::Map." );
644 }
645 catch( moab::Element::Map::EvaluationError )
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.
748 void print_tuples( TupleList* tlp )
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
819 int print_vertex_fields( Interface* mbi,
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 }
1.9.1.