MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 // A test file for Subset Normalization 00002 #include "moab/ParallelComm.hpp" 00003 #include "MBParallelConventions.h" 00004 #include "moab/Core.hpp" 00005 #include "moab/FileOptions.hpp" 00006 #include "ReadParallel.hpp" 00007 #include "Coupler.hpp" 00008 #include "DebugOutput.hpp" 00009 #include "ElemUtil.hpp" 00010 #include <iostream> 00011 #include <iomanip> 00012 #include <sstream> 00013 #include <cstring> 00014 #include <cstdlib> 00015 extern "C" { 00016 #include "moab/FindPtFuncs.h" 00017 } 00018 #include "moab/TupleList.hpp" 00019 #include "moab/gs.hpp" 00020 #include "moab/Types.hpp" 00021 #ifndef IS_BUILDING_MB 00022 #define IS_BUILDING_MB 00023 #include "Internals.hpp" 00024 #undef IS_BUILDING_MB 00025 #else 00026 #include "Internals.hpp" 00027 #endif 00028 00029 using namespace moab; 00030 00031 bool debug = true; 00032 00033 // Forward declarations 00034 void get_file_options( int argc, char** argv, std::vector< const char* >& filenames, std::string& norm_tag, 00035 std::vector< const char* >& tag_names, std::vector< const char* >& tag_values, 00036 std::string& file_opts, int* err ); 00037 00038 void print_tuples( TupleList* tlp ); 00039 00040 int print_vertex_fields( Interface* mbi, std::vector< std::vector< EntityHandle > >& groups, Tag& norm_hdl, 00041 Coupler::IntegType integ_type ); 00042 00043 double const_field( double x, double y, double z ); 00044 double field_1( double x, double y, double z ); 00045 double field_2( double x, double y, double z ); 00046 double field_3( double x, double y, double z ); 00047 double physField( double x, double y, double z ); 00048 00049 ErrorCode integrate_scalar_field_test(); 00050 int pack_tuples( TupleList* tl, void** ptr ); 00051 void unpack_tuples( void* ptr, TupleList** tlp ); 00052 00053 // 00054 // Start of main test program 00055 // 00056 int main( int argc, char** argv ) 00057 { 00058 // need to init MPI first, to tell how many procs and rank 00059 // Used since Coupler is a parallel code. The Coupler may be run 00060 // in parallel or serial mode but will need MPI either way. 00061 int err = MPI_Init( &argc, &argv ); 00062 00063 // Print usage if not enough arguments 00064 if( argc < 3 ) 00065 { 00066 std::cerr << "Usage: "; 00067 std::cerr << argv[0] << " <nfiles> <fname1> ... <fnamen> <norm_tag> <tag_select_opts> <file_opts>" << std::endl; 00068 std::cerr << "nfiles : number of mesh files" << std::endl; 00069 std::cerr << "fname1...fnamen : mesh files" << std::endl; 00070 std::cerr << "norm_tag : name of tag to normalize across meshes" << std::endl; 00071 std::cerr << "tag_select_opts : quoted string of tags and values for subset selection, " 00072 "e.g. \"TAG1=VAL1;TAG2=VAL2;TAG3;TAG4\"" 00073 << std::endl; 00074 std::cerr << "file_opts : quoted string of parallel file read options, e.g. " 00075 "\"OPTION1=VALUE1;OPTION2;OPTION3=VALUE3\"" 00076 << std::endl; 00077 00078 err = integrate_scalar_field_test();MB_CHK_SET_ERR( (ErrorCode)err, "Integrate scalar field test failed" ); 00079 00080 err = MPI_Finalize(); 00081 00082 return err; 00083 } 00084 00085 int nprocs, rank; 00086 err = MPI_Comm_size( MPI_COMM_WORLD, &nprocs ); 00087 assert( MPI_SUCCESS == err ); 00088 err = MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 00089 assert( MPI_SUCCESS == err ); 00090 00091 // Create an ofstream to write output. One file each for each proc. 00092 std::stringstream fname; 00093 fname << argv[0] << rank << ".out"; 00094 if( !std::freopen( fname.str().c_str(), "a", stdout ) ) return -1; 00095 if( !std::freopen( fname.str().c_str(), "a", stderr ) ) return -1; 00096 00097 // Create the moab instance 00098 Interface* mbi = new Core(); 00099 if( NULL == mbi ) 00100 { 00101 std::cerr << "MOAB constructor failed" << std::endl; 00102 return 1; 00103 } 00104 00105 // Get the input options 00106 std::cout << "Getting options..." << std::endl; 00107 std::vector< const char* > filenames; 00108 std::vector< const char* > tagNames; 00109 std::vector< const char* > tagValues; 00110 std::string normTag, fileOpts; 00111 get_file_options( argc, argv, filenames, normTag, tagNames, tagValues, fileOpts, &err );MB_CHK_SET_ERR( (ErrorCode)err, "get_file_options failed" ); 00112 00113 // Print out the input parameters 00114 std::cout << " Input Parameters - " << std::endl; 00115 std::cout << " Filenames: "; 00116 for( std::vector< const char* >::iterator it = filenames.begin(); it != filenames.end(); ++it ) 00117 std::cout << *it << " "; 00118 std::cout << std::endl; 00119 std::cout << " Norm Tag: " << normTag << std::endl; 00120 std::cout << " Selection Data: NumNames=" << tagNames.size() << " NumValues=" << tagValues.size() << std::endl; 00121 std::cout << " TagNames TagValues " << std::endl; 00122 std::cout << " -------------------- --------------------" << std::endl; 00123 std::vector< const char* >::iterator nameIt = tagNames.begin(); 00124 std::vector< const char* >::iterator valIt = tagValues.begin(); 00125 std::cout << std::setiosflags( std::ios::left ); 00126 for( ; nameIt != tagNames.end(); ++nameIt ) 00127 { 00128 std::cout << " " << std::setw( 20 ) << *nameIt; 00129 if( *valIt != 0 ) 00130 { 00131 std::cout << " " << std::setw( 20 ) << *( (int*)( *valIt ) ) << std::endl; 00132 ++valIt; 00133 } 00134 else 00135 std::cout << " NULL " << std::endl; 00136 } 00137 std::cout << std::resetiosflags( std::ios::left ); 00138 std::cout << " File Options: " << fileOpts << std::endl; 00139 00140 // Read in mesh(es) 00141 std::cout << "Reading mesh file(s)..." << std::endl; 00142 std::vector< ParallelComm* > pcs( filenames.size() ); 00143 std::vector< ReadParallel* > rps( filenames.size() ); 00144 00145 // allocate root sets for each mesh for moab 00146 std::vector< EntityHandle > roots( filenames.size() ); 00147 00148 ErrorCode result; 00149 for( unsigned int i = 0; i < filenames.size(); i++ ) 00150 { 00151 pcs[i] = new ParallelComm( mbi, MPI_COMM_WORLD ); 00152 rps[i] = new ReadParallel( mbi, pcs[i] ); 00153 00154 result = mbi->create_meshset( MESHSET_SET, roots[i] ); 00155 00156 MB_CHK_SET_ERR( result, "Creating root set failed" ); 00157 result = rps[i]->load_file( filenames[i], &roots[i], FileOptions( fileOpts.c_str() ) );MB_CHK_SET_ERR( result, "load_file failed" ); 00158 } 00159 00160 // Initialize the debug object for Range printing 00161 DebugOutput debugOut( "ssn_test-", std::cerr ); 00162 debugOut.set_rank( rank ); 00163 debugOut.set_verbosity( 10 ); 00164 00165 // Output what is in root sets 00166 for( unsigned int k = 0; k < filenames.size(); k++ ) 00167 { 00168 00169 Range rootRg; 00170 result = mbi->get_entities_by_handle( roots[k], rootRg );MB_CHK_SET_ERR( result, "can't get entities" ); 00171 debugOut.print( 2, "Root set entities: ", rootRg ); 00172 rootRg.clear(); 00173 00174 Range partRg; 00175 pcs[k]->get_part_entities( partRg ); 00176 debugOut.print( 2, "Partition entities: ", partRg ); 00177 partRg.clear(); 00178 } 00179 00180 // source is 1st mesh, target is 2nd mesh 00181 Range src_elems, targ_elems; 00182 00183 // ****************************** 00184 std::cout << "********** Create Coupler **********" << std::endl; 00185 // Create a coupler 00186 std::cout << "Creating Coupler..." << std::endl; 00187 Coupler mbc( mbi, pcs[0], src_elems, 0 ); 00188 00189 // Get tag handles for passed in tags 00190 std::cout << "Getting tag handles..." << std::endl; 00191 int numTagNames = tagNames.size(); 00192 00193 std::vector< Tag > tagHandles( numTagNames ); 00194 int iTags = 0; 00195 while( iTags < numTagNames ) 00196 { 00197 std::cout << "Getting handle for " << tagNames[iTags] << std::endl; 00198 result = mbi->tag_get_handle( tagNames[iTags], tagHandles[iTags] );MB_CHK_SET_ERR( result, "Retrieving tag handles failed" ); 00199 iTags++; 00200 } 00201 00202 // ****************************** 00203 std::cout << "********** Test create_tuples **********" << std::endl; 00204 // First get some EntitySets for Mesh 1 and Mesh 2 00205 { 00206 00207 Range entsets1, entsets2; 00208 result = mbi->get_entities_by_type_and_tag( roots[0], MBENTITYSET, &tagHandles[0], 00209 (const void* const*)&tagValues[0], tagHandles.size(), entsets1, 00210 Interface::INTERSECT ); // recursive is false 00211 MB_CHK_SET_ERR( result, "sets: get_entities_by_type_and_tag failed on Mesh 1." ); 00212 00213 // Create tuple_list for each mesh's 00214 std::cout << "Creating tuples for mesh 1..." << std::endl; 00215 TupleList* m1TagTuples = NULL; 00216 err = mbc.create_tuples( entsets1, &tagHandles[0], tagHandles.size(), &m1TagTuples );MB_CHK_SET_ERR( (ErrorCode)err, "create_tuples failed" ); 00217 00218 std::cout << " create_tuples returned" << std::endl; 00219 print_tuples( m1TagTuples ); 00220 00221 result = mbi->get_entities_by_type_and_tag( roots[1], MBENTITYSET, &tagHandles[0], 00222 (const void* const*)&tagValues[0], tagHandles.size(), entsets2, 00223 Interface::INTERSECT ); // recursive is false 00224 MB_CHK_SET_ERR( result, "sets: get_entities_by_type_and_tag failed on Mesh 2." ); 00225 00226 std::cout << "Creating tuples for mesh 2..." << std::endl; 00227 TupleList* m2TagTuples = NULL; 00228 err = mbc.create_tuples( entsets2, (Tag*)( &tagHandles[0] ), tagHandles.size(), &m2TagTuples );MB_CHK_SET_ERR( (ErrorCode)err, "create_tuples failed" ); 00229 00230 std::cout << " create_tuples returned" << std::endl; 00231 print_tuples( m2TagTuples ); 00232 00233 // ****************************** 00234 std::cout << "********** Test consolidate_tuples **********" << std::endl; 00235 // In this serial version we only have the tuples from Mesh 1 and Mesh 2. 00236 // Just consolidate those for the test. 00237 std::cout << "Consolidating tuple_lists for Mesh 1 and Mesh 2..." << std::endl; 00238 TupleList** tplp_arr = (TupleList**)malloc( 2 * sizeof( TupleList* ) ); 00239 TupleList* unique_tpl = NULL; 00240 tplp_arr[0] = m1TagTuples; 00241 tplp_arr[1] = m2TagTuples; 00242 00243 err = mbc.consolidate_tuples( tplp_arr, 2, &unique_tpl );MB_CHK_SET_ERR( (ErrorCode)err, "consolidate_tuples failed" ); 00244 std::cout << " consolidate_tuples returned" << std::endl; 00245 print_tuples( unique_tpl ); 00246 } 00247 00248 // ****************************** 00249 std::cout << "********** Test get_matching_entities **********" << std::endl; 00250 std::vector< std::vector< EntityHandle > > m1EntitySets; 00251 std::vector< std::vector< EntityHandle > > m1EntityGroups; 00252 std::vector< std::vector< EntityHandle > > m2EntitySets; 00253 std::vector< std::vector< EntityHandle > > m2EntityGroups; 00254 00255 // Get matching entities for Mesh 1 00256 std::cout << "Get matching entities for mesh 1..." << std::endl; 00257 err = mbc.get_matching_entities( roots[0], &tagHandles[0], &tagValues[0], tagHandles.size(), &m1EntitySets, 00258 &m1EntityGroups );MB_CHK_SET_ERR( (ErrorCode)err, "get_matching_entities failed" ); 00259 00260 std::cout << " get_matching_entities returned " << m1EntityGroups.size() << " entity groups" << std::endl; 00261 00262 // Print out the data in the vector of vectors 00263 std::vector< std::vector< EntityHandle > >::iterator iter_esi; 00264 std::vector< std::vector< EntityHandle > >::iterator iter_egi; 00265 std::vector< EntityHandle >::iterator iter_esj; 00266 std::vector< EntityHandle >::iterator iter_egj; 00267 Range entSetRg; 00268 int icnt; 00269 for( iter_egi = m1EntityGroups.begin(), iter_esi = m1EntitySets.begin(), icnt = 1; 00270 ( iter_egi != m1EntityGroups.end() ) && ( iter_esi != m1EntitySets.end() ); ++iter_egi, ++iter_esi, icnt++ ) 00271 { 00272 std::cout << " EntityGroup(" << icnt << ") = "; 00273 std::cout.flush(); 00274 entSetRg.clear(); 00275 for( iter_egj = ( *iter_egi ).begin(); iter_egj != ( *iter_egi ).end(); ++iter_egj ) 00276 entSetRg.insert( (EntityHandle)*iter_egj ); 00277 debugOut.print( 2, "Mesh1 matching Entities: ", entSetRg ); 00278 std::cout.flush(); 00279 00280 std::cout << " EntitySet(" << icnt << ") = "; 00281 std::cout.flush(); 00282 entSetRg.clear(); 00283 for( iter_esj = ( *iter_esi ).begin(); iter_esj != ( *iter_esi ).end(); ++iter_esj ) 00284 entSetRg.insert( (EntityHandle)*iter_esj ); 00285 debugOut.print( 2, "Mesh1 matching EntitySets: ", entSetRg ); 00286 std::cout.flush(); 00287 } 00288 00289 // Get matching entities for Mesh 2 00290 std::cout << "Get matching entities for mesh 2..." << std::endl; 00291 err = mbc.get_matching_entities( roots[1], &tagHandles[0], &tagValues[0], tagHandles.size(), &m2EntitySets, 00292 &m2EntityGroups );MB_CHK_SET_ERR( (ErrorCode)err, "get_matching_entities failed" ); 00293 00294 std::cout << " get_matching_entities returned " << m2EntityGroups.size() << " entity groups" << std::endl; 00295 for( iter_egi = m2EntityGroups.begin(), iter_esi = m2EntitySets.begin(), icnt = 1; 00296 ( iter_egi != m2EntityGroups.end() ) && ( iter_esi != m2EntitySets.end() ); ++iter_egi, ++iter_esi, icnt++ ) 00297 { 00298 std::cout << " EntityGroup(" << icnt << ") = "; 00299 std::cout.flush(); 00300 entSetRg.clear(); 00301 for( iter_egj = ( *iter_egi ).begin(); iter_egj != ( *iter_egi ).end(); ++iter_egj ) 00302 entSetRg.insert( (EntityHandle)*iter_egj ); 00303 debugOut.print( 2, "Mesh2 matching Entities: ", entSetRg ); 00304 std::cout.flush(); 00305 00306 std::cout << " EntitySet(" << icnt << ") = "; 00307 std::cout.flush(); 00308 entSetRg.clear(); 00309 for( iter_esj = ( *iter_esi ).begin(); iter_esj != ( *iter_esi ).end(); ++iter_esj ) 00310 entSetRg.insert( (EntityHandle)*iter_esj ); 00311 debugOut.print( 2, "Mesh2 matching EntitySets: ", entSetRg ); 00312 std::cout.flush(); 00313 } 00314 00315 if( debug ) 00316 { 00317 // ****************************** 00318 std::cout << "********** Test print_tuples **********" << std::endl; 00319 // temporary test function 00320 std::cout << "Testing print_tuples..." << std::endl; 00321 00322 TupleList test_tuple; 00323 int num_ints = 3, num_longs = 2, num_ulongs = 4, num_reals = 6, num_rows = 10; 00324 00325 std::cout << " print of test_tuples zero init..." << std::endl; 00326 test_tuple.initialize( 0, 0, 0, 0, 0 ); 00327 00328 test_tuple.enableWriteAccess(); 00329 00330 print_tuples( &test_tuple ); 00331 00332 std::cout << " print of test_tuples after setting n to 10..." << std::endl; 00333 test_tuple.set_n( 10 ); 00334 print_tuples( &test_tuple ); 00335 00336 test_tuple.initialize( num_ints, num_longs, num_ulongs, num_reals, num_rows ); 00337 std::cout << " print of test_tuples after init..." << std::endl; 00338 print_tuples( &test_tuple ); 00339 00340 std::cout << " print of test_tuples after setting n to 10..." << std::endl; 00341 test_tuple.set_n( 10 ); 00342 print_tuples( &test_tuple ); 00343 00344 for( int i = 0; i < num_rows; i++ ) 00345 { 00346 int j; 00347 for( j = 0; j < num_ints; j++ ) 00348 test_tuple.vi_wr[i * num_ints + j] = (int)( ( j + 1 ) * ( i + 1 ) ); 00349 00350 for( j = 0; j < num_longs; j++ ) 00351 test_tuple.vl_wr[i * num_longs + j] = (int)( ( j + 1 ) * ( i + 1 ) ); 00352 00353 for( j = 0; j < num_ulongs; j++ ) 00354 test_tuple.vul_wr[i * num_ulongs + j] = (int)( ( j + 1 ) * ( i + 1 ) ); 00355 00356 for( j = 0; j < num_reals; j++ ) 00357 test_tuple.vr_wr[i * num_reals + j] = (int)( ( j + 1 ) * ( i + 1 ) + ( j * 0.01 ) ); 00358 } 00359 std::cout << " print of test_tuples after filling with data..." << std::endl; 00360 print_tuples( &test_tuple ); 00361 00362 // ****************************** 00363 std::cout << "********** Test pack_tuples and unpack_tuples **********" << std::endl; 00364 void* mp_buf; 00365 int buf_sz; 00366 if( rank == 0 ) { buf_sz = pack_tuples( &test_tuple, &mp_buf ); } 00367 00368 // Send buffer size 00369 err = MPI_Bcast( &buf_sz, 1, MPI_INT, 0, MPI_COMM_WORLD ); 00370 00371 if( err != MPI_SUCCESS ) 00372 { 00373 std::cerr << "MPI_Bcast of buffer size failed" << std::endl; 00374 return -1; 00375 } 00376 00377 // Allocate a buffer in the other procs 00378 if( rank != 0 ) { mp_buf = malloc( buf_sz * sizeof( uint ) ); } 00379 00380 err = MPI_Bcast( mp_buf, buf_sz * sizeof( uint ), MPI_UNSIGNED_CHAR, 0, MPI_COMM_WORLD ); 00381 if( err != MPI_SUCCESS ) 00382 { 00383 std::cerr << "MPI_Bcast of buffer failed" << std::endl; 00384 return -1; 00385 } 00386 00387 TupleList* rcv_tuples; 00388 unpack_tuples( mp_buf, &rcv_tuples ); 00389 00390 std::cout << " print of rcv_tuples after unpacking from MPI_Bcast..." << std::endl; 00391 print_tuples( rcv_tuples ); 00392 } 00393 00394 err = integrate_scalar_field_test();MB_CHK_SET_ERR( (ErrorCode)err, "Failure in integrating a scalar_field" ); 00395 00396 // ****************************** 00397 std::cout << "********** Test get_group_integ_vals **********" << std::endl; 00398 std::cout << "Get group integrated field values..." << std::endl; 00399 00400 // print the field values at the vertices before change. 00401 std::cout << " print vertex field values first:" << std::endl; 00402 Tag norm_hdl; 00403 result = mbi->tag_get_handle( normTag.c_str(), norm_hdl );MB_CHK_SET_ERR( (ErrorCode)err, "Failed to get tag handle." ); 00404 00405 Coupler::IntegType integ_type = Coupler::VOLUME; 00406 // Mesh 1 field values 00407 std::cout << " Original entity vertex field values (mesh 1): " << std::endl; 00408 print_vertex_fields( mbi, m1EntityGroups, norm_hdl, integ_type ); 00409 00410 // Mesh 2 field values 00411 std::cout << " Original entity vertex field values (mesh 2): " << std::endl; 00412 print_vertex_fields( mbi, m2EntityGroups, norm_hdl, integ_type ); 00413 00414 // Get the field values 00415 std::vector< double >::iterator iter_ivals; 00416 00417 std::cout << "Get group integrated field values for mesh 1..." << std::endl; 00418 std::vector< double > m1IntegVals( m1EntityGroups.size() ); 00419 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." ); 00420 std::cout << "Mesh 1 integrated field values(" << m1IntegVals.size() << "): "; 00421 for( iter_ivals = m1IntegVals.begin(); iter_ivals != m1IntegVals.end(); ++iter_ivals ) 00422 { 00423 std::cout << ( *iter_ivals ) << " "; 00424 } 00425 std::cout << std::endl; 00426 00427 std::cout << "Get group integrated field values for mesh 2..." << std::endl; 00428 std::vector< double > m2IntegVals( m2EntityGroups.size() ); 00429 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." ); 00430 std::cout << "Mesh 2 integrated field values(" << m2IntegVals.size() << "): "; 00431 for( iter_ivals = m2IntegVals.begin(); iter_ivals != m2IntegVals.end(); ++iter_ivals ) 00432 { 00433 std::cout << ( *iter_ivals ) << " "; 00434 } 00435 std::cout << std::endl; 00436 00437 // ****************************** 00438 std::cout << "********** Test apply_group_norm_factors **********" << std::endl; 00439 // Make the norm factors by inverting the integration values. 00440 double val; 00441 for( unsigned int i = 0; i < m1IntegVals.size(); i++ ) 00442 { 00443 val = m1IntegVals[i]; 00444 m1IntegVals[i] = 1 / val; 00445 } 00446 00447 for( unsigned int i = 0; i < m2IntegVals.size(); i++ ) 00448 { 00449 val = m2IntegVals[i]; 00450 m2IntegVals[i] = 1 / val; 00451 } 00452 00453 std::cout << "Mesh 1 norm factors(" << m1IntegVals.size() << "): "; 00454 for( iter_ivals = m1IntegVals.begin(); iter_ivals != m1IntegVals.end(); ++iter_ivals ) 00455 { 00456 std::cout << ( *iter_ivals ) << " "; 00457 } 00458 std::cout << std::endl; 00459 00460 std::cout << "Mesh 2 norm factors(" << m2IntegVals.size() << "): "; 00461 for( iter_ivals = m2IntegVals.begin(); iter_ivals != m2IntegVals.end(); ++iter_ivals ) 00462 { 00463 std::cout << ( *iter_ivals ) << " "; 00464 } 00465 std::cout << std::endl; 00466 00467 // Apply the factors and reprint the vertices 00468 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." ); 00469 00470 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." ); 00471 00472 // Get the norm_tag_factor on the EntitySets 00473 // Get the handle for the norm factor tag 00474 Tag norm_factor_hdl; 00475 std::string normFactor = normTag + "_normf"; 00476 result = mbi->tag_get_handle( normFactor.c_str(), norm_factor_hdl );MB_CHK_SET_ERR( result, "Failed to get norm factor tag handle." ); 00477 00478 // Mesh 1 values 00479 std::cout << "Mesh 1 norm factors per EntitySet..."; 00480 for( iter_esi = m1EntitySets.begin(); iter_esi != m1EntitySets.end(); ++iter_esi ) 00481 { 00482 for( iter_esj = ( *iter_esi ).begin(); iter_esj != ( *iter_esi ).end(); ++iter_esj ) 00483 { 00484 double data = 0; 00485 EntityHandle eh = *iter_esj; 00486 result = mbi->tag_get_data( norm_factor_hdl, &eh, 1, &data );MB_CHK_SET_ERR( result, "Failed to get tag data." ); 00487 std::cout << data << ", "; 00488 } 00489 } 00490 std::cout << std::endl; 00491 00492 // Mesh 2 values 00493 std::cout << "Mesh 2 norm factors per EntitySet..."; 00494 for( iter_esi = m2EntitySets.begin(); iter_esi != m2EntitySets.end(); ++iter_esi ) 00495 { 00496 for( iter_esj = ( *iter_esi ).begin(); iter_esj != ( *iter_esi ).end(); ++iter_esj ) 00497 { 00498 double data = 0; 00499 EntityHandle eh = *iter_esj; 00500 result = mbi->tag_get_data( norm_factor_hdl, &eh, 1, &data );MB_CHK_SET_ERR( result, "Failed to get tag data." ); 00501 std::cout << data << ", "; 00502 } 00503 } 00504 std::cout << std::endl; 00505 00506 // ****************************** 00507 std::cout << "********** Test normalize_subset **********" << std::endl; 00508 // Now call the Coupler::normalize_subset routine and see if we get an error. 00509 std::cout << "Running Coupler::normalize_subset() on mesh 1" << std::endl; 00510 err = mbc.normalize_subset( (EntityHandle)roots[0], normTag.c_str(), &tagNames[0], numTagNames, &tagValues[0], 00511 Coupler::VOLUME, 4 );MB_CHK_SET_ERR( (ErrorCode)err, "Failure in call to Coupler::normalize_subset() on mesh 1" ); 00512 00513 // Print the normFactor on each EntitySet after the above call. 00514 // Mesh 1 values 00515 std::cout << "Mesh 1 norm factors per EntitySet..."; 00516 for( iter_esi = m1EntitySets.begin(); iter_esi != m1EntitySets.end(); ++iter_esi ) 00517 { 00518 for( iter_esj = ( *iter_esi ).begin(); iter_esj != ( *iter_esi ).end(); ++iter_esj ) 00519 { 00520 double data = 0; 00521 EntityHandle eh = *iter_esj; 00522 result = mbi->tag_get_data( norm_factor_hdl, &eh, 1, &data );MB_CHK_SET_ERR( result, "Failed to get tag data." ); 00523 std::cout << data << ", "; 00524 } 00525 } 00526 std::cout << std::endl; 00527 00528 std::cout << "Running Coupler::normalize_subset() on mesh 2" << std::endl; 00529 err = mbc.normalize_subset( (EntityHandle)roots[1], normTag.c_str(), &tagNames[0], numTagNames, &tagValues[0], 00530 Coupler::VOLUME, 4 );MB_CHK_SET_ERR( (ErrorCode)err, "Failure in call to Coupler::normalize_subset() on mesh 2" ); 00531 00532 // Mesh 2 values 00533 std::cout << "Mesh 2 norm factors per EntitySet..."; 00534 for( iter_esi = m2EntitySets.begin(); iter_esi != m2EntitySets.end(); ++iter_esi ) 00535 { 00536 for( iter_esj = ( *iter_esi ).begin(); iter_esj != ( *iter_esi ).end(); ++iter_esj ) 00537 { 00538 double data = 0; 00539 EntityHandle eh = *iter_esj; 00540 result = mbi->tag_get_data( norm_factor_hdl, &eh, 1, &data );MB_CHK_SET_ERR( result, "Failed to get tag data." ); 00541 00542 std::cout << data << ", "; 00543 } 00544 } 00545 std::cout << std::endl; 00546 00547 // Done, cleanup 00548 std::cout << "********** ssn_test DONE! **********" << std::endl; 00549 MPI_Finalize(); 00550 return 0; 00551 } 00552 00553 ErrorCode integrate_scalar_field_test() 00554 { 00555 // ****************************** 00556 std::cout << "********** Test moab::Element::Map::integrate_scalar_field **********" << std::endl; 00557 // Create a simple hex centered at 0,0,0 with sides of length 2. 00558 std::vector< CartVect > biunit_cube( 8 ); 00559 biunit_cube[0] = CartVect( -1, -1, -1 ); 00560 biunit_cube[1] = CartVect( 1, -1, -1 ); 00561 biunit_cube[2] = CartVect( 1, 1, -1 ); 00562 biunit_cube[3] = CartVect( -1, 1, -1 ); 00563 biunit_cube[4] = CartVect( -1, -1, 1 ); 00564 biunit_cube[5] = CartVect( 1, -1, 1 ); 00565 biunit_cube[6] = CartVect( 1, 1, 1 ); 00566 biunit_cube[7] = CartVect( -1, 1, 1 ); 00567 00568 std::vector< CartVect > zerobase_cube( 8 ); 00569 zerobase_cube[0] = CartVect( 0, 0, 0 ); 00570 zerobase_cube[1] = CartVect( 2, 0, 0 ); 00571 zerobase_cube[2] = CartVect( 2, 2, 0 ); 00572 zerobase_cube[3] = CartVect( 0, 2, 0 ); 00573 zerobase_cube[4] = CartVect( 0, 0, 2 ); 00574 zerobase_cube[5] = CartVect( 2, 0, 2 ); 00575 zerobase_cube[6] = CartVect( 2, 2, 2 ); 00576 zerobase_cube[7] = CartVect( 0, 2, 2 ); 00577 00578 // Calculate field values at the corners of both cubes 00579 double bcf[8], bf1[8], bf2[8], bf3[8], zcf[8], zf1[8], zf2[8], zf3[8]; 00580 for( int i = 0; i < 8; i++ ) 00581 { 00582 bcf[i] = const_field( biunit_cube[i][0], biunit_cube[i][1], biunit_cube[i][2] ); 00583 bf1[i] = field_1( biunit_cube[i][0], biunit_cube[i][1], biunit_cube[i][2] ); 00584 bf2[i] = field_2( biunit_cube[i][0], biunit_cube[i][1], biunit_cube[i][2] ); 00585 bf3[i] = field_3( biunit_cube[i][0], biunit_cube[i][1], biunit_cube[i][2] ); 00586 00587 zcf[i] = const_field( zerobase_cube[i][0], zerobase_cube[i][1], zerobase_cube[i][2] ); 00588 zf1[i] = field_1( zerobase_cube[i][0], zerobase_cube[i][1], zerobase_cube[i][2] ); 00589 zf2[i] = field_2( zerobase_cube[i][0], zerobase_cube[i][1], zerobase_cube[i][2] ); 00590 zf3[i] = field_3( zerobase_cube[i][0], zerobase_cube[i][1], zerobase_cube[i][2] ); 00591 } 00592 00593 std::cout << "Integrated values:" << std::endl; 00594 00595 try 00596 { 00597 double field_const1, field_const2; 00598 double field_linear1, field_linear2; 00599 double field_quad1, field_quad2; 00600 double field_cubic1, field_cubic2; 00601 00602 int ipoints = 0; 00603 Element::LinearHex biunit_hexMap( biunit_cube ); 00604 Element::LinearHex zerobase_hexMap( zerobase_cube ); 00605 00606 field_const1 = biunit_hexMap.integrate_scalar_field( bcf ); 00607 field_const2 = zerobase_hexMap.integrate_scalar_field( zcf ); 00608 std::cout << " binunit_cube, const_field(num_pts=" << ipoints << "): field_val=" << field_const1 00609 << std::endl; 00610 std::cout << " zerobase_cube, const_field(num_pts=" << ipoints << "): field_val=" << field_const2 00611 << std::endl; 00612 00613 field_linear1 = biunit_hexMap.integrate_scalar_field( bf1 ); 00614 field_linear2 = zerobase_hexMap.integrate_scalar_field( zf1 ); 00615 std::cout << " binunit_cube, field_1(num_pts=" << ipoints << "): field_val=" << field_linear1 << std::endl; 00616 std::cout << " zerobase_cube, field_1(num_pts=" << ipoints << "): field_val=" << field_linear2 << std::endl; 00617 00618 field_quad1 = biunit_hexMap.integrate_scalar_field( bf2 ); 00619 field_quad2 = zerobase_hexMap.integrate_scalar_field( zf2 ); 00620 std::cout << " binunit_cube, field_2(num_pts=" << ipoints << "): field_val=" << field_quad1 << std::endl; 00621 std::cout << " zerobase_cube, field_2(num_pts=" << ipoints << "): field_val=" << field_quad2 << std::endl; 00622 00623 field_cubic1 = biunit_hexMap.integrate_scalar_field( bf3 ); 00624 field_cubic2 = zerobase_hexMap.integrate_scalar_field( zf3 ); 00625 std::cout << " binunit_cube, field_3(num_pts=" << ipoints << "): field_val=" << field_cubic1 << std::endl; 00626 std::cout << " zerobase_cube, field_3(num_pts=" << ipoints << "): field_val=" << field_cubic2 << std::endl; 00627 } 00628 catch( moab::Element::Map::ArgError ) 00629 { 00630 MB_CHK_SET_ERR( MB_FAILURE, "Failed to set vertices on Element::Map." ); 00631 } 00632 catch( moab::Element::Map::EvaluationError ) 00633 { 00634 MB_CHK_SET_ERR( MB_FAILURE, "Failed to get inverse evaluation of coordinate on Element::Map." ); 00635 } 00636 return MB_SUCCESS; 00637 } 00638 00639 // Function to parse input parameters 00640 void get_file_options( int argc, char** argv, std::vector< const char* >& filenames, std::string& normTag, 00641 std::vector< const char* >& tagNames, std::vector< const char* >& tagValues, 00642 std::string& fileOpts, int* err ) 00643 { 00644 int npos = 1; 00645 00646 // get number of files 00647 int nfiles = atoi( argv[npos++] ); 00648 00649 // get mesh filenames 00650 filenames.resize( nfiles ); 00651 for( int i = 0; i < nfiles; i++ ) 00652 filenames[i] = argv[npos++]; 00653 00654 // get normTag 00655 if( npos < argc ) 00656 normTag = argv[npos++]; 00657 else 00658 { 00659 std::cerr << "Insufficient parameters: norm_tag missing" << std::endl; 00660 *err = 1; 00661 return; 00662 } 00663 00664 // get tag selection options 00665 if( npos < argc ) 00666 { 00667 char* opts = argv[npos++]; 00668 // char sep1[1] = {';'}; 00669 // char sep2[1] = {'='}; 00670 bool end_vals_seen = false; 00671 std::vector< char* > tmpTagOpts; 00672 00673 // first get the options 00674 for( char* i = strtok( opts, ";" ); i; i = strtok( 0, ";" ) ) 00675 { 00676 if( debug ) std::cout << "get_file_options: i=" << i << std::endl; 00677 tmpTagOpts.push_back( i ); 00678 } 00679 00680 // parse out the name and val or just name. 00681 for( unsigned int j = 0; j < tmpTagOpts.size(); j++ ) 00682 { 00683 char* e = strtok( tmpTagOpts[j], "=" ); 00684 if( debug ) std::cout << "get_file_options: name=" << e << std::endl; 00685 tagNames.push_back( e ); 00686 e = strtok( 0, "=" ); 00687 if( e != NULL ) 00688 { 00689 if( debug ) std::cout << "get_file_options: val=" << e << std::endl; 00690 // We have a value 00691 if( end_vals_seen ) 00692 { 00693 // ERROR we should not have a value after none are seen 00694 std::cerr << "Incorrect parameters: new value seen after end of values" << std::endl; 00695 *err = 1; 00696 return; 00697 } 00698 // Otherwise get the value string from e and convert it to an int 00699 int* valp = new int; 00700 *valp = atoi( e ); 00701 tagValues.push_back( (const char*)valp ); 00702 } 00703 else 00704 { 00705 // Otherwise there is no '=' so push a null on the list 00706 end_vals_seen = true; 00707 tagValues.push_back( (const char*)0 ); 00708 } 00709 } 00710 } 00711 else 00712 { 00713 std::cerr << "Insufficient parameters: tag_select_opts missing" << std::endl; 00714 *err = 1; 00715 return; 00716 } 00717 00718 // get fileOpts 00719 if( npos < argc ) 00720 fileOpts = argv[npos++]; 00721 else 00722 { 00723 std::cerr << "Insufficient parameters: file_opts missing" << std::endl; 00724 *err = 1; 00725 return; 00726 } 00727 } 00728 00729 // Function to print out a tuple_list. 00730 void print_tuples( TupleList* tlp ) 00731 { 00732 uint mi, ml, mul, mr; 00733 tlp->getTupleSize( mi, ml, mul, mr ); 00734 std::cout << " tuple data: (n=" << tlp->get_n() << ")" << std::endl; 00735 std::cout << " mi:" << mi << " ml:" << ml << " mul:" << mul << " mr:" << mr << std::endl; 00736 std::cout << " [" << std::setw( 11 * mi ) << " int data" 00737 << " |" << std::setw( 11 * ml ) << " long data" 00738 << " |" << std::setw( 11 * mul ) << " Ulong data" 00739 << " |" << std::setw( 11 * mr ) << " real data" 00740 << " " << std::endl 00741 << " "; 00742 for( unsigned int i = 0; i < tlp->get_n(); i++ ) 00743 { 00744 if( mi > 0 ) 00745 { 00746 for( unsigned int j = 0; j < mi; j++ ) 00747 { 00748 std::cout << std::setw( 10 ) << tlp->vi_rd[i * mi + j] << " "; 00749 } 00750 } 00751 else 00752 { 00753 std::cout << " "; 00754 } 00755 std::cout << "| "; 00756 00757 if( ml > 0 ) 00758 { 00759 for( unsigned int j = 0; j < ml; j++ ) 00760 { 00761 std::cout << std::setw( 10 ) << tlp->vl_rd[i * ml + j] << " "; 00762 } 00763 } 00764 else 00765 { 00766 std::cout << " "; 00767 } 00768 std::cout << "| "; 00769 00770 if( mul > 0 ) 00771 { 00772 for( unsigned int j = 0; j < mul; j++ ) 00773 { 00774 std::cout << std::setw( 10 ) << tlp->vul_rd[i * mul + j] << " "; 00775 } 00776 } 00777 else 00778 { 00779 std::cout << " "; 00780 } 00781 std::cout << "| "; 00782 00783 if( mr > 0 ) 00784 { 00785 for( unsigned int j = 0; j < mr; j++ ) 00786 { 00787 std::cout << std::setw( 10 ) << tlp->vr_rd[i * mr + j] << " "; 00788 } 00789 } 00790 else 00791 { 00792 std::cout << " "; 00793 } 00794 00795 if( i + 1 < tlp->get_n() ) std::cout << std::endl << " "; 00796 } 00797 std::cout << "]" << std::endl; 00798 } 00799 00800 // Function to print vertex field values 00801 int print_vertex_fields( Interface* mbi, std::vector< std::vector< EntityHandle > >& groups, Tag& norm_hdl, 00802 Coupler::IntegType integ_type ) 00803 { 00804 int err = 0; 00805 ErrorCode result; 00806 std::vector< EntityHandle >::iterator iter_j; 00807 00808 for( unsigned int i = 0; i < groups.size(); i++ ) 00809 { 00810 std::cout << " Group - " << std::endl << " "; 00811 for( iter_j = groups[i].begin(); iter_j != groups[i].end(); ++iter_j ) 00812 { 00813 EntityHandle ehandle = ( *iter_j ); 00814 // Check that the entity in iter_j is of the same dimension as the 00815 // integ_type we are performing 00816 int j_type = mbi->dimension_from_handle( ehandle ); 00817 00818 if( ( integ_type == Coupler::VOLUME ) && ( j_type != 3 ) ) continue; 00819 00820 // Retrieve the vertices from the element 00821 const EntityHandle* conn = NULL; 00822 int num_verts = 0; 00823 result = mbi->get_connectivity( ehandle, conn, num_verts ); 00824 if( MB_SUCCESS != result ) return 1; 00825 std::cout << std::fixed; 00826 for( int iv = 0; iv < num_verts; iv++ ) 00827 { 00828 double data = 0; 00829 result = mbi->tag_get_data( norm_hdl, &conn[iv], 1, &data ); 00830 if( MB_SUCCESS != result ) return 1; 00831 std::cout << std::setprecision( 8 ) << data << ", "; 00832 } 00833 std::cout << std::endl << " "; 00834 } 00835 std::cout << std::endl; 00836 std::cout.unsetf( std::ios_base::floatfield ); // turn off fixed notation 00837 } 00838 00839 return err; 00840 } 00841 00842 // Function for a constant field value 00843 double const_field( double /*x*/, double /*y*/, double /*z*/ ) 00844 { 00845 // return 5.0/40.0; 00846 return 5.0; 00847 } 00848 00849 // Functions for a some field values 00850 double field_1( double x, double y, double z ) 00851 { 00852 double f = fabs( x ) + fabs( y ) + fabs( z ); 00853 // return f/24.0; 00854 return f; 00855 } 00856 00857 double field_2( double x, double y, double z ) 00858 { 00859 double f = x * x + y * y + z * z; 00860 // return f/32.0; 00861 return f; 00862 } 00863 00864 double field_3( double x, double y, double z ) 00865 { 00866 double f = 2 * x + 2 * y + 2 * z; 00867 // return f/48.0; 00868 return f; 00869 } 00870 00871 // Function used to create field on mesh for testing. 00872 double physField( double x, double y, double z ) 00873 { 00874 double out; 00875 00876 // 1/r^2 decay from {0,0,0} 00877 00878 out = x * x + y * y + z * z; 00879 out += 1e-1; // clamp 00880 out = 1 / out; 00881 00882 return out; 00883 } 00884 00885 #define UINT_PER_X( X ) ( ( sizeof( X ) + sizeof( uint ) - 1 ) / sizeof( uint ) ) 00886 #define UINT_PER_REAL UINT_PER_X( realType ) 00887 #define UINT_PER_LONG UINT_PER_X( slong ) 00888 #define UINT_PER_ULONG UINT_PER_X( Ulong ) 00889 #define UINT_PER_UNSIGNED UINT_PER_X( unsigned ) 00890 00891 // Function for packing tuple_list 00892 int pack_tuples( TupleList* tl, void** ptr ) 00893 { 00894 uint mi, ml, mul, mr; 00895 tl->getTupleSize( mi, ml, mul, mr ); 00896 00897 uint n = tl->get_n(); 00898 00899 int sz_buf = 1 + 4 * UINT_PER_UNSIGNED + 00900 tl->get_n() * ( mi + ml * UINT_PER_LONG + mul * UINT_PER_LONG + mr * UINT_PER_REAL ); 00901 00902 uint* buf = (uint*)malloc( sz_buf * sizeof( uint ) ); 00903 *ptr = (void*)buf; 00904 00905 // copy n 00906 memcpy( buf, &n, sizeof( uint ) ), buf += 1; 00907 // copy mi 00908 memcpy( buf, &mi, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED; 00909 // copy ml 00910 memcpy( buf, &ml, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED; 00911 // copy mul 00912 memcpy( buf, &mul, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED; 00913 // copy mr 00914 memcpy( buf, &mr, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED; 00915 // copy vi_rd 00916 memcpy( buf, tl->vi_rd, tl->get_n() * mi * sizeof( sint ) ), buf += tl->get_n() * mi; 00917 // copy vl_rd 00918 memcpy( buf, tl->vl_rd, tl->get_n() * ml * sizeof( slong ) ), buf += tl->get_n() * ml * UINT_PER_LONG; 00919 // copy vul_rd 00920 memcpy( buf, tl->vul_rd, tl->get_n() * mul * sizeof( Ulong ) ), buf += tl->get_n() * mul * UINT_PER_ULONG; 00921 // copy vr_rd 00922 memcpy( buf, tl->vr_rd, tl->get_n() * mr * sizeof( realType ) ), buf += tl->get_n() * mr * UINT_PER_REAL; 00923 00924 return sz_buf; 00925 } 00926 00927 // Function for packing tuple_list 00928 void unpack_tuples( void* ptr, TupleList** tlp ) 00929 { 00930 TupleList* tl = new TupleList(); 00931 *tlp = tl; 00932 00933 uint nt; 00934 unsigned mit, mlt, mult, mrt; 00935 uint* buf = (uint*)ptr; 00936 00937 // get n 00938 memcpy( &nt, buf, sizeof( uint ) ), buf += 1; 00939 // get mi 00940 memcpy( &mit, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED; 00941 // get ml 00942 memcpy( &mlt, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED; 00943 // get mul 00944 memcpy( &mult, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED; 00945 // get mr 00946 memcpy( &mrt, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED; 00947 00948 // initialize tl 00949 tl->initialize( mit, mlt, mult, mrt, nt ); 00950 tl->enableWriteAccess(); 00951 tl->set_n( nt ); 00952 00953 uint mi, ml, mul, mr; 00954 tl->getTupleSize( mi, ml, mul, mr ); 00955 00956 // get vi_wr 00957 memcpy( tl->vi_wr, buf, tl->get_n() * mi * sizeof( sint ) ), buf += tl->get_n() * mi; 00958 // get vl_wr 00959 memcpy( tl->vl_wr, buf, tl->get_n() * ml * sizeof( slong ) ), buf += tl->get_n() * ml * UINT_PER_LONG; 00960 // get vul_wr 00961 memcpy( tl->vul_wr, buf, tl->get_n() * mul * sizeof( Ulong ) ), buf += tl->get_n() * mul * UINT_PER_ULONG; 00962 // get vr_wr 00963 memcpy( tl->vr_wr, buf, tl->get_n() * mr * sizeof( realType ) ), buf += tl->get_n() * mr * UINT_PER_REAL; 00964 00965 tl->disableWriteAccess(); 00966 00967 return; 00968 }