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