MOAB: Mesh Oriented datABase  (version 5.2.1)
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, 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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines