ssn_test.cpp

Go to the documentation of this file.

// A test file for Subset Normalization
#include "moab/ParallelComm.hpp"
#include "MBParallelConventions.h"
#include "moab/Core.hpp"
#include "moab/FileOptions.hpp"
#include "ReadParallel.hpp"
#include "Coupler.hpp"
#include "DebugOutput.hpp"
#include "ElemUtil.hpp"
#include <iostream>
#include <iomanip>
#include <sstream>
#include <cstring>
#include <cstdlib>
extern "C" {
#include "moab/FindPtFuncs.h"
}
#include "moab/TupleList.hpp"
#include "moab/gs.hpp"
#include "moab/Types.hpp"
#ifndef IS_BUILDING_MB
#define IS_BUILDING_MB
#include "Internals.hpp"
#undef IS_BUILDING_MB
#else
#include "Internals.hpp"
#endif
 
using namespace moab;
 
bool debug = true;
 
// Forward declarations
void get_file_options( int argc,
 char** argv,
 std::vector< const char* >& filenames,
 std::string& norm_tag,
 std::vector< const char* >& tag_names,
 std::vector< const char* >& tag_values,
 std::string& file_opts,
 int* err );
 
void print_tuples( TupleList* tlp );
 
int print_vertex_fields( Interface* mbi,
 std::vector< std::vector< EntityHandle > >& groups,
 Tag& norm_hdl,
 Coupler::IntegType integ_type );
 
double const_field( double x, double y, double z );
double field_1( double x, double y, double z );
double field_2( double x, double y, double z );
double field_3( double x, double y, double z );
double physField( double x, double y, double z );
 
ErrorCode integrate_scalar_field_test();
int pack_tuples( TupleList* tl, void** ptr );
void unpack_tuples( void* ptr, TupleList** tlp );
 
//
// Start of main test program
//
int main( int argc, char** argv )
{
 // need to init MPI first, to tell how many procs and rank
 // Used since Coupler is a parallel code. The Coupler may be run
 // in parallel or serial mode but will need MPI either way.
 int err = MPI_Init( &argc, &argv );
 
 // Print usage if not enough arguments
 if( argc < 3 )
 {
 std::cerr << "Usage: ";
 std::cerr << argv[0] << " <nfiles> <fname1> ... <fnamen> <norm_tag> <tag_select_opts> <file_opts>" << std::endl;
 std::cerr << "nfiles : number of mesh files" << std::endl;
 std::cerr << "fname1...fnamen : mesh files" << std::endl;
 std::cerr << "norm_tag : name of tag to normalize across meshes" << std::endl;
 std::cerr << "tag_select_opts : quoted string of tags and values for subset selection, "
 "e.g. \"TAG1=VAL1;TAG2=VAL2;TAG3;TAG4\""
 << std::endl;
 std::cerr << "file_opts : quoted string of parallel file read options, e.g. "
 "\"OPTION1=VALUE1;OPTION2;OPTION3=VALUE3\""
 << std::endl;
 
 err = integrate_scalar_field_test();MB_CHK_SET_ERR( (ErrorCode)err, "Integrate scalar field test failed" );
 
 err = MPI_Finalize();
 
 return err;
 }
 
 int nprocs, rank;
 err = MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
 assert( MPI_SUCCESS == err );
 err = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
 assert( MPI_SUCCESS == err );
 
 // Create an ofstream to write output. One file each for each proc.
 std::stringstream fname;
 fname << argv[0] << rank << ".out";
 if( !std::freopen( fname.str().c_str(), "a", stdout ) ) return -1;
 if( !std::freopen( fname.str().c_str(), "a", stderr ) ) return -1;
 
 // Create the moab instance
 Interface* mbi = new Core();
 if( NULL == mbi )
 {
 std::cerr << "MOAB constructor failed" << std::endl;
 return 1;
 }
 
 // Get the input options
 std::cout << "Getting options..." << std::endl;
 std::vector< const char* > filenames;
 std::vector< const char* > tagNames;
 std::vector< const char* > tagValues;
 std::string normTag, fileOpts;
 get_file_options( argc, argv, filenames, normTag, tagNames, tagValues, fileOpts, &err );MB_CHK_SET_ERR( (ErrorCode)err, "get_file_options failed" );
 
 // Print out the input parameters
 std::cout << " Input Parameters - " << std::endl;
 std::cout << " Filenames: ";
 for( std::vector< const char* >::iterator it = filenames.begin(); it != filenames.end(); ++it )
 std::cout << *it << " ";
 std::cout << std::endl;
 std::cout << " Norm Tag: " << normTag << std::endl;
 std::cout << " Selection Data: NumNames=" << tagNames.size() << " NumValues=" << tagValues.size() << std::endl;
 std::cout << " TagNames TagValues " << std::endl;
 std::cout << " -------------------- --------------------" << std::endl;
 std::vector< const char* >::iterator nameIt = tagNames.begin();
 std::vector< const char* >::iterator valIt = tagValues.begin();
 std::cout << std::setiosflags( std::ios::left );
 for( ; nameIt != tagNames.end(); ++nameIt )
 {
 std::cout << " " << std::setw( 20 ) << *nameIt;
 if( *valIt != 0 )
 {
 std::cout << " " << std::setw( 20 ) << *( (int*)( *valIt ) ) << std::endl;
 ++valIt;
 }
 else
 std::cout << " NULL " << std::endl;
 }
 std::cout << std::resetiosflags( std::ios::left );
 std::cout << " File Options: " << fileOpts << std::endl;
 
 // Read in mesh(es)
 std::cout << "Reading mesh file(s)..." << std::endl;
 std::vector< ParallelComm* > pcs( filenames.size() );
 std::vector< ReadParallel* > rps( filenames.size() );
 
 // allocate root sets for each mesh for moab
 std::vector< EntityHandle > roots( filenames.size() );
 
 ErrorCode result;
 for( unsigned int i = 0; i < filenames.size(); i++ )
 {
 pcs[i] = new ParallelComm( mbi, MPI_COMM_WORLD );
 rps[i] = new ReadParallel( mbi, pcs[i] );
 
 result = mbi->create_meshset( MESHSET_SET, roots[i] );
 
 MB_CHK_SET_ERR( result, "Creating root set failed" );
 result = rps[i]->load_file( filenames[i], &roots[i], FileOptions( fileOpts.c_str() ) );MB_CHK_SET_ERR( result, "load_file failed" );
 }
 
 // Initialize the debug object for Range printing
 DebugOutput debugOut( "ssn_test-", std::cerr );
 debugOut.set_rank( rank );
 debugOut.set_verbosity( 10 );
 
 // Output what is in root sets
 for( unsigned int k = 0; k < filenames.size(); k++ )
 {
 
 Range rootRg;
 result = mbi->get_entities_by_handle( roots[k], rootRg );MB_CHK_SET_ERR( result, "can't get entities" );
 debugOut.print( 2, "Root set entities: ", rootRg );
 rootRg.clear();
 
 Range partRg;
 pcs[k]->get_part_entities( partRg );
 debugOut.print( 2, "Partition entities: ", partRg );
 partRg.clear();
 }
 
 // source is 1st mesh, target is 2nd mesh
 Range src_elems, targ_elems;
 
 // ******************************
 std::cout << "********** Create Coupler **********" << std::endl;
 // Create a coupler
 std::cout << "Creating Coupler..." << std::endl;
 Coupler mbc( mbi, pcs[0], src_elems, 0 );
 
 // Get tag handles for passed in tags
 std::cout << "Getting tag handles..." << std::endl;
 int numTagNames = tagNames.size();
 
 std::vector< Tag > tagHandles( numTagNames );
 int iTags = 0;
 while( iTags < numTagNames )
 {
 std::cout << "Getting handle for " << tagNames[iTags] << std::endl;
 result = mbi->tag_get_handle( tagNames[iTags], tagHandles[iTags] );MB_CHK_SET_ERR( result, "Retrieving tag handles failed" );
 iTags++;
 }
 
 // ******************************
 std::cout << "********** Test create_tuples **********" << std::endl;
 // First get some EntitySets for Mesh 1 and Mesh 2
 {
 
 Range entsets1, entsets2;
 result = mbi->get_entities_by_type_and_tag( roots[0], MBENTITYSET, &tagHandles[0],
 (const void* const*)&tagValues[0], tagHandles.size(), entsets1,
 Interface::INTERSECT ); // recursive is false
 MB_CHK_SET_ERR( result, "sets: get_entities_by_type_and_tag failed on Mesh 1." );
 
 // Create tuple_list for each mesh's
 std::cout << "Creating tuples for mesh 1..." << std::endl;
 TupleList* m1TagTuples = NULL;
 err = mbc.create_tuples( entsets1, &tagHandles[0], tagHandles.size(), &m1TagTuples );MB_CHK_SET_ERR( (ErrorCode)err, "create_tuples failed" );
 
 std::cout << " create_tuples returned" << std::endl;
 print_tuples( m1TagTuples );
 
 result = mbi->get_entities_by_type_and_tag( roots[1], MBENTITYSET, &tagHandles[0],
 (const void* const*)&tagValues[0], tagHandles.size(), entsets2,
 Interface::INTERSECT ); // recursive is false
 MB_CHK_SET_ERR( result, "sets: get_entities_by_type_and_tag failed on Mesh 2." );
 
 std::cout << "Creating tuples for mesh 2..." << std::endl;
 TupleList* m2TagTuples = NULL;
 err = mbc.create_tuples( entsets2, (Tag*)( &tagHandles[0] ), tagHandles.size(), &m2TagTuples );MB_CHK_SET_ERR( (ErrorCode)err, "create_tuples failed" );
 
 std::cout << " create_tuples returned" << std::endl;
 print_tuples( m2TagTuples );
 
 // ******************************
 std::cout << "********** Test consolidate_tuples **********" << std::endl;
 // In this serial version we only have the tuples from Mesh 1 and Mesh 2.
 // Just consolidate those for the test.
 std::cout << "Consolidating tuple_lists for Mesh 1 and Mesh 2..." << std::endl;
 TupleList** tplp_arr = (TupleList**)malloc( 2 * sizeof( TupleList* ) );
 TupleList* unique_tpl = NULL;
 tplp_arr[0] = m1TagTuples;
 tplp_arr[1] = m2TagTuples;
 
 err = mbc.consolidate_tuples( tplp_arr, 2, &unique_tpl );MB_CHK_SET_ERR( (ErrorCode)err, "consolidate_tuples failed" );
 std::cout << " consolidate_tuples returned" << std::endl;
 print_tuples( unique_tpl );
 }
 
 // ******************************
 std::cout << "********** Test get_matching_entities **********" << std::endl;
 std::vector< std::vector< EntityHandle > > m1EntitySets;
 std::vector< std::vector< EntityHandle > > m1EntityGroups;
 std::vector< std::vector< EntityHandle > > m2EntitySets;
 std::vector< std::vector< EntityHandle > > m2EntityGroups;
 
 // Get matching entities for Mesh 1
 std::cout << "Get matching entities for mesh 1..." << std::endl;
 err = mbc.get_matching_entities( roots[0], &tagHandles[0], &tagValues[0], tagHandles.size(), &m1EntitySets,
 &m1EntityGroups );MB_CHK_SET_ERR( (ErrorCode)err, "get_matching_entities failed" );
 
 std::cout << " get_matching_entities returned " << m1EntityGroups.size() << " entity groups" << std::endl;
 
 // Print out the data in the vector of vectors
 std::vector< std::vector< EntityHandle > >::iterator iter_esi;
 std::vector< std::vector< EntityHandle > >::iterator iter_egi;
 std::vector< EntityHandle >::iterator iter_esj;
 std::vector< EntityHandle >::iterator iter_egj;
 Range entSetRg;
 int icnt;
 for( iter_egi = m1EntityGroups.begin(), iter_esi = m1EntitySets.begin(), icnt = 1;
 ( iter_egi != m1EntityGroups.end() ) && ( iter_esi != m1EntitySets.end() ); ++iter_egi, ++iter_esi, icnt++ )
 {
 std::cout << " EntityGroup(" << icnt << ") = ";
 std::cout.flush();
 entSetRg.clear();
 for( iter_egj = ( *iter_egi ).begin(); iter_egj != ( *iter_egi ).end(); ++iter_egj )
 entSetRg.insert( (EntityHandle)*iter_egj );
 debugOut.print( 2, "Mesh1 matching Entities: ", entSetRg );
 std::cout.flush();
 
 std::cout << " EntitySet(" << icnt << ") = ";
 std::cout.flush();
 entSetRg.clear();
 for( iter_esj = ( *iter_esi ).begin(); iter_esj != ( *iter_esi ).end(); ++iter_esj )
 entSetRg.insert( (EntityHandle)*iter_esj );
 debugOut.print( 2, "Mesh1 matching EntitySets: ", entSetRg );
 std::cout.flush();
 }
 
 // Get matching entities for Mesh 2
 std::cout << "Get matching entities for mesh 2..." << std::endl;
 err = mbc.get_matching_entities( roots[1], &tagHandles[0], &tagValues[0], tagHandles.size(), &m2EntitySets,
 &m2EntityGroups );MB_CHK_SET_ERR( (ErrorCode)err, "get_matching_entities failed" );
 
 std::cout << " get_matching_entities returned " << m2EntityGroups.size() << " entity groups" << std::endl;
 for( iter_egi = m2EntityGroups.begin(), iter_esi = m2EntitySets.begin(), icnt = 1;
 ( iter_egi != m2EntityGroups.end() ) && ( iter_esi != m2EntitySets.end() ); ++iter_egi, ++iter_esi, icnt++ )
 {
 std::cout << " EntityGroup(" << icnt << ") = ";
 std::cout.flush();
 entSetRg.clear();
 for( iter_egj = ( *iter_egi ).begin(); iter_egj != ( *iter_egi ).end(); ++iter_egj )
 entSetRg.insert( (EntityHandle)*iter_egj );
 debugOut.print( 2, "Mesh2 matching Entities: ", entSetRg );
 std::cout.flush();
 
 std::cout << " EntitySet(" << icnt << ") = ";
 std::cout.flush();
 entSetRg.clear();
 for( iter_esj = ( *iter_esi ).begin(); iter_esj != ( *iter_esi ).end(); ++iter_esj )
 entSetRg.insert( (EntityHandle)*iter_esj );
 debugOut.print( 2, "Mesh2 matching EntitySets: ", entSetRg );
 std::cout.flush();
 }
 
 if( debug )
 {
 // ******************************
 std::cout << "********** Test print_tuples **********" << std::endl;
 // temporary test function
 std::cout << "Testing print_tuples..." << std::endl;
 
 TupleList test_tuple;
 int num_ints = 3, num_longs = 2, num_ulongs = 4, num_reals = 6, num_rows = 10;
 
 std::cout << " print of test_tuples zero init..." << std::endl;
 test_tuple.initialize( 0, 0, 0, 0, 0 );
 
 test_tuple.enableWriteAccess();
 
 print_tuples( &test_tuple );
 
 std::cout << " print of test_tuples after setting n to 10..." << std::endl;
 test_tuple.set_n( 10 );
 print_tuples( &test_tuple );
 
 test_tuple.initialize( num_ints, num_longs, num_ulongs, num_reals, num_rows );
 std::cout << " print of test_tuples after init..." << std::endl;
 print_tuples( &test_tuple );
 
 std::cout << " print of test_tuples after setting n to 10..." << std::endl;
 test_tuple.set_n( 10 );
 print_tuples( &test_tuple );
 
 for( int i = 0; i < num_rows; i++ )
 {
 int j;
 for( j = 0; j < num_ints; j++ )
 test_tuple.vi_wr[i * num_ints + j] = (int)( ( j + 1 ) * ( i + 1 ) );
 
 for( j = 0; j < num_longs; j++ )
 test_tuple.vl_wr[i * num_longs + j] = (int)( ( j + 1 ) * ( i + 1 ) );
 
 for( j = 0; j < num_ulongs; j++ )
 test_tuple.vul_wr[i * num_ulongs + j] = (int)( ( j + 1 ) * ( i + 1 ) );
 
 for( j = 0; j < num_reals; j++ )
 test_tuple.vr_wr[i * num_reals + j] = (int)( ( j + 1 ) * ( i + 1 ) + ( j * 0.01 ) );
 }
 std::cout << " print of test_tuples after filling with data..." << std::endl;
 print_tuples( &test_tuple );
 
 // ******************************
 std::cout << "********** Test pack_tuples and unpack_tuples **********" << std::endl;
 void* mp_buf;
 int buf_sz;
 if( rank == 0 )
 {
 buf_sz = pack_tuples( &test_tuple, &mp_buf );
 }
 
 // Send buffer size
 err = MPI_Bcast( &buf_sz, 1, MPI_INT, 0, MPI_COMM_WORLD );
 
 if( err != MPI_SUCCESS )
 {
 std::cerr << "MPI_Bcast of buffer size failed" << std::endl;
 return -1;
 }
 
 // Allocate a buffer in the other procs
 if( rank != 0 )
 {
 mp_buf = malloc( buf_sz * sizeof( uint ) );
 }
 
 err = MPI_Bcast( mp_buf, buf_sz * sizeof( uint ), MPI_UNSIGNED_CHAR, 0, MPI_COMM_WORLD );
 if( err != MPI_SUCCESS )
 {
 std::cerr << "MPI_Bcast of buffer failed" << std::endl;
 return -1;
 }
 
 TupleList* rcv_tuples;
 unpack_tuples( mp_buf, &rcv_tuples );
 
 std::cout << " print of rcv_tuples after unpacking from MPI_Bcast..." << std::endl;
 print_tuples( rcv_tuples );
 }
 
 err = integrate_scalar_field_test();MB_CHK_SET_ERR( (ErrorCode)err, "Failure in integrating a scalar_field" );
 
 // ******************************
 std::cout << "********** Test get_group_integ_vals **********" << std::endl;
 std::cout << "Get group integrated field values..." << std::endl;
 
 // print the field values at the vertices before change.
 std::cout << " print vertex field values first:" << std::endl;
 Tag norm_hdl;
 result = mbi->tag_get_handle( normTag.c_str(), norm_hdl );MB_CHK_SET_ERR( (ErrorCode)err, "Failed to get tag handle." );
 
 Coupler::IntegType integ_type = Coupler::VOLUME;
 // Mesh 1 field values
 std::cout << " Original entity vertex field values (mesh 1): " << std::endl;
 print_vertex_fields( mbi, m1EntityGroups, norm_hdl, integ_type );
 
 // Mesh 2 field values
 std::cout << " Original entity vertex field values (mesh 2): " << std::endl;
 print_vertex_fields( mbi, m2EntityGroups, norm_hdl, integ_type );
 
 // Get the field values
 std::vector< double >::iterator iter_ivals;
 
 std::cout << "Get group integrated field values for mesh 1..." << std::endl;
 std::vector< double > m1IntegVals( m1EntityGroups.size() );
 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." );
 std::cout << "Mesh 1 integrated field values(" << m1IntegVals.size() << "): ";
 for( iter_ivals = m1IntegVals.begin(); iter_ivals != m1IntegVals.end(); ++iter_ivals )
 {
 std::cout << ( *iter_ivals ) << " ";
 }
 std::cout << std::endl;
 
 std::cout << "Get group integrated field values for mesh 2..." << std::endl;
 std::vector< double > m2IntegVals( m2EntityGroups.size() );
 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." );
 std::cout << "Mesh 2 integrated field values(" << m2IntegVals.size() << "): ";
 for( iter_ivals = m2IntegVals.begin(); iter_ivals != m2IntegVals.end(); ++iter_ivals )
 {
 std::cout << ( *iter_ivals ) << " ";
 }
 std::cout << std::endl;
 
 // ******************************
 std::cout << "********** Test apply_group_norm_factors **********" << std::endl;
 // Make the norm factors by inverting the integration values.
 double val;
 for( unsigned int i = 0; i < m1IntegVals.size(); i++ )
 {
 val = m1IntegVals[i];
 m1IntegVals[i] = 1 / val;
 }
 
 for( unsigned int i = 0; i < m2IntegVals.size(); i++ )
 {
 val = m2IntegVals[i];
 m2IntegVals[i] = 1 / val;
 }
 
 std::cout << "Mesh 1 norm factors(" << m1IntegVals.size() << "): ";
 for( iter_ivals = m1IntegVals.begin(); iter_ivals != m1IntegVals.end(); ++iter_ivals )
 {
 std::cout << ( *iter_ivals ) << " ";
 }
 std::cout << std::endl;
 
 std::cout << "Mesh 2 norm factors(" << m2IntegVals.size() << "): ";
 for( iter_ivals = m2IntegVals.begin(); iter_ivals != m2IntegVals.end(); ++iter_ivals )
 {
 std::cout << ( *iter_ivals ) << " ";
 }
 std::cout << std::endl;
 
 // Apply the factors and reprint the vertices
 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." );
 
 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." );
 
 // Get the norm_tag_factor on the EntitySets
 // Get the handle for the norm factor tag
 Tag norm_factor_hdl;
 std::string normFactor = normTag + "_normf";
 result = mbi->tag_get_handle( normFactor.c_str(), norm_factor_hdl );MB_CHK_SET_ERR( result, "Failed to get norm factor tag handle." );
 
 // Mesh 1 values
 std::cout << "Mesh 1 norm factors per EntitySet...";
 for( iter_esi = m1EntitySets.begin(); iter_esi != m1EntitySets.end(); ++iter_esi )
 {
 for( iter_esj = ( *iter_esi ).begin(); iter_esj != ( *iter_esi ).end(); ++iter_esj )
 {
 double data = 0;
 EntityHandle eh = *iter_esj;
 result = mbi->tag_get_data( norm_factor_hdl, &eh, 1, &data );MB_CHK_SET_ERR( result, "Failed to get tag data." );
 std::cout << data << ", ";
 }
 }
 std::cout << std::endl;
 
 // Mesh 2 values
 std::cout << "Mesh 2 norm factors per EntitySet...";
 for( iter_esi = m2EntitySets.begin(); iter_esi != m2EntitySets.end(); ++iter_esi )
 {
 for( iter_esj = ( *iter_esi ).begin(); iter_esj != ( *iter_esi ).end(); ++iter_esj )
 {
 double data = 0;
 EntityHandle eh = *iter_esj;
 result = mbi->tag_get_data( norm_factor_hdl, &eh, 1, &data );MB_CHK_SET_ERR( result, "Failed to get tag data." );
 std::cout << data << ", ";
 }
 }
 std::cout << std::endl;
 
 // ******************************
 std::cout << "********** Test normalize_subset **********" << std::endl;
 // Now call the Coupler::normalize_subset routine and see if we get an error.
 std::cout << "Running Coupler::normalize_subset() on mesh 1" << std::endl;
 err = mbc.normalize_subset( (EntityHandle)roots[0], normTag.c_str(), &tagNames[0], numTagNames, &tagValues[0],
 Coupler::VOLUME, 4 );MB_CHK_SET_ERR( (ErrorCode)err, "Failure in call to Coupler::normalize_subset() on mesh 1" );
 
 // Print the normFactor on each EntitySet after the above call.
 // Mesh 1 values
 std::cout << "Mesh 1 norm factors per EntitySet...";
 for( iter_esi = m1EntitySets.begin(); iter_esi != m1EntitySets.end(); ++iter_esi )
 {
 for( iter_esj = ( *iter_esi ).begin(); iter_esj != ( *iter_esi ).end(); ++iter_esj )
 {
 double data = 0;
 EntityHandle eh = *iter_esj;
 result = mbi->tag_get_data( norm_factor_hdl, &eh, 1, &data );MB_CHK_SET_ERR( result, "Failed to get tag data." );
 std::cout << data << ", ";
 }
 }
 std::cout << std::endl;
 
 std::cout << "Running Coupler::normalize_subset() on mesh 2" << std::endl;
 err = mbc.normalize_subset( (EntityHandle)roots[1], normTag.c_str(), &tagNames[0], numTagNames, &tagValues[0],
 Coupler::VOLUME, 4 );MB_CHK_SET_ERR( (ErrorCode)err, "Failure in call to Coupler::normalize_subset() on mesh 2" );
 
 // Mesh 2 values
 std::cout << "Mesh 2 norm factors per EntitySet...";
 for( iter_esi = m2EntitySets.begin(); iter_esi != m2EntitySets.end(); ++iter_esi )
 {
 for( iter_esj = ( *iter_esi ).begin(); iter_esj != ( *iter_esi ).end(); ++iter_esj )
 {
 double data = 0;
 EntityHandle eh = *iter_esj;
 result = mbi->tag_get_data( norm_factor_hdl, &eh, 1, &data );MB_CHK_SET_ERR( result, "Failed to get tag data." );
 
 std::cout << data << ", ";
 }
 }
 std::cout << std::endl;
 
 // Done, cleanup
 std::cout << "********** ssn_test DONE! **********" << std::endl;
 MPI_Finalize();
 return 0;
}
 
ErrorCode integrate_scalar_field_test()
{
 // ******************************
 std::cout << "********** Test moab::Element::Map::integrate_scalar_field **********" << std::endl;
 // Create a simple hex centered at 0,0,0 with sides of length 2.
 std::vector< CartVect > biunit_cube( 8 );
 biunit_cube[0] = CartVect( -1, -1, -1 );
 biunit_cube[1] = CartVect( 1, -1, -1 );
 biunit_cube[2] = CartVect( 1, 1, -1 );
 biunit_cube[3] = CartVect( -1, 1, -1 );
 biunit_cube[4] = CartVect( -1, -1, 1 );
 biunit_cube[5] = CartVect( 1, -1, 1 );
 biunit_cube[6] = CartVect( 1, 1, 1 );
 biunit_cube[7] = CartVect( -1, 1, 1 );
 
 std::vector< CartVect > zerobase_cube( 8 );
 zerobase_cube[0] = CartVect( 0, 0, 0 );
 zerobase_cube[1] = CartVect( 2, 0, 0 );
 zerobase_cube[2] = CartVect( 2, 2, 0 );
 zerobase_cube[3] = CartVect( 0, 2, 0 );
 zerobase_cube[4] = CartVect( 0, 0, 2 );
 zerobase_cube[5] = CartVect( 2, 0, 2 );
 zerobase_cube[6] = CartVect( 2, 2, 2 );
 zerobase_cube[7] = CartVect( 0, 2, 2 );
 
 // Calculate field values at the corners of both cubes
 double bcf[8], bf1[8], bf2[8], bf3[8], zcf[8], zf1[8], zf2[8], zf3[8];
 for( int i = 0; i < 8; i++ )
 {
 bcf[i] = const_field( biunit_cube[i][0], biunit_cube[i][1], biunit_cube[i][2] );
 bf1[i] = field_1( biunit_cube[i][0], biunit_cube[i][1], biunit_cube[i][2] );
 bf2[i] = field_2( biunit_cube[i][0], biunit_cube[i][1], biunit_cube[i][2] );
 bf3[i] = field_3( biunit_cube[i][0], biunit_cube[i][1], biunit_cube[i][2] );
 
 zcf[i] = const_field( zerobase_cube[i][0], zerobase_cube[i][1], zerobase_cube[i][2] );
 zf1[i] = field_1( zerobase_cube[i][0], zerobase_cube[i][1], zerobase_cube[i][2] );
 zf2[i] = field_2( zerobase_cube[i][0], zerobase_cube[i][1], zerobase_cube[i][2] );
 zf3[i] = field_3( zerobase_cube[i][0], zerobase_cube[i][1], zerobase_cube[i][2] );
 }
 
 std::cout << "Integrated values:" << std::endl;
 
 try
 {
 double field_const1, field_const2;
 double field_linear1, field_linear2;
 double field_quad1, field_quad2;
 double field_cubic1, field_cubic2;
 
 int ipoints = 0;
 Element::LinearHex biunit_hexMap( biunit_cube );
 Element::LinearHex zerobase_hexMap( zerobase_cube );
 
 field_const1 = biunit_hexMap.integrate_scalar_field( bcf );
 field_const2 = zerobase_hexMap.integrate_scalar_field( zcf );
 std::cout << " binunit_cube, const_field(num_pts=" << ipoints << "): field_val=" << field_const1
 << std::endl;
 std::cout << " zerobase_cube, const_field(num_pts=" << ipoints << "): field_val=" << field_const2
 << std::endl;
 
 field_linear1 = biunit_hexMap.integrate_scalar_field( bf1 );
 field_linear2 = zerobase_hexMap.integrate_scalar_field( zf1 );
 std::cout << " binunit_cube, field_1(num_pts=" << ipoints << "): field_val=" << field_linear1 << std::endl;
 std::cout << " zerobase_cube, field_1(num_pts=" << ipoints << "): field_val=" << field_linear2 << std::endl;
 
 field_quad1 = biunit_hexMap.integrate_scalar_field( bf2 );
 field_quad2 = zerobase_hexMap.integrate_scalar_field( zf2 );
 std::cout << " binunit_cube, field_2(num_pts=" << ipoints << "): field_val=" << field_quad1 << std::endl;
 std::cout << " zerobase_cube, field_2(num_pts=" << ipoints << "): field_val=" << field_quad2 << std::endl;
 
 field_cubic1 = biunit_hexMap.integrate_scalar_field( bf3 );
 field_cubic2 = zerobase_hexMap.integrate_scalar_field( zf3 );
 std::cout << " binunit_cube, field_3(num_pts=" << ipoints << "): field_val=" << field_cubic1 << std::endl;
 std::cout << " zerobase_cube, field_3(num_pts=" << ipoints << "): field_val=" << field_cubic2 << std::endl;
 }
 catch( moab::Element::Map::ArgError )
 {
 MB_CHK_SET_ERR( MB_FAILURE, "Failed to set vertices on Element::Map." );
 }
 catch( moab::Element::Map::EvaluationError )
 {
 MB_CHK_SET_ERR( MB_FAILURE, "Failed to get inverse evaluation of coordinate on Element::Map." );
 }
 return MB_SUCCESS;
}
 
// Function to parse input parameters
void get_file_options( int argc,
 char** argv,
 std::vector< const char* >& filenames,
 std::string& normTag,
 std::vector< const char* >& tagNames,
 std::vector< const char* >& tagValues,
 std::string& fileOpts,
 int* err )
{
 int npos = 1;
 
 // get number of files
 int nfiles = atoi( argv[npos++] );
 
 // get mesh filenames
 filenames.resize( nfiles );
 for( int i = 0; i < nfiles; i++ )
 filenames[i] = argv[npos++];
 
 // get normTag
 if( npos < argc )
 normTag = argv[npos++];
 else
 {
 std::cerr << "Insufficient parameters: norm_tag missing" << std::endl;
 *err = 1;
 return;
 }
 
 // get tag selection options
 if( npos < argc )
 {
 char* opts = argv[npos++];
 // char sep1[1] = {';'};
 // char sep2[1] = {'='};
 bool end_vals_seen = false;
 std::vector< char* > tmpTagOpts;
 
 // first get the options
 for( char* i = strtok( opts, ";" ); i; i = strtok( 0, ";" ) )
 {
 if( debug ) std::cout << "get_file_options: i=" << i << std::endl;
 tmpTagOpts.push_back( i );
 }
 
 // parse out the name and val or just name.
 for( unsigned int j = 0; j < tmpTagOpts.size(); j++ )
 {
 char* e = strtok( tmpTagOpts[j], "=" );
 if( debug ) std::cout << "get_file_options: name=" << e << std::endl;
 tagNames.push_back( e );
 e = strtok( 0, "=" );
 if( e != NULL )
 {
 if( debug ) std::cout << "get_file_options: val=" << e << std::endl;
 // We have a value
 if( end_vals_seen )
 {
 // ERROR we should not have a value after none are seen
 std::cerr << "Incorrect parameters: new value seen after end of values" << std::endl;
 *err = 1;
 return;
 }
 // Otherwise get the value string from e and convert it to an int
 int* valp = new int;
 *valp = atoi( e );
 tagValues.push_back( (const char*)valp );
 }
 else
 {
 // Otherwise there is no '=' so push a null on the list
 end_vals_seen = true;
 tagValues.push_back( (const char*)0 );
 }
 }
 }
 else
 {
 std::cerr << "Insufficient parameters: tag_select_opts missing" << std::endl;
 *err = 1;
 return;
 }
 
 // get fileOpts
 if( npos < argc )
 fileOpts = argv[npos++];
 else
 {
 std::cerr << "Insufficient parameters: file_opts missing" << std::endl;
 *err = 1;
 return;
 }
}
 
// Function to print out a tuple_list.
void print_tuples( TupleList* tlp )
{
 uint mi, ml, mul, mr;
 tlp->getTupleSize( mi, ml, mul, mr );
 std::cout << " tuple data: (n=" << tlp->get_n() << ")" << std::endl;
 std::cout << " mi:" << mi << " ml:" << ml << " mul:" << mul << " mr:" << mr << std::endl;
 std::cout << " [" << std::setw( 11 * mi ) << " int data"
 << " |" << std::setw( 11 * ml ) << " long data"
 << " |" << std::setw( 11 * mul ) << " Ulong data"
 << " |" << std::setw( 11 * mr ) << " real data"
 << " " << std::endl
 << " ";
 for( unsigned int i = 0; i < tlp->get_n(); i++ )
 {
 if( mi > 0 )
 {
 for( unsigned int j = 0; j < mi; j++ )
 {
 std::cout << std::setw( 10 ) << tlp->vi_rd[i * mi + j] << " ";
 }
 }
 else
 {
 std::cout << " ";
 }
 std::cout << "| ";
 
 if( ml > 0 )
 {
 for( unsigned int j = 0; j < ml; j++ )
 {
 std::cout << std::setw( 10 ) << tlp->vl_rd[i * ml + j] << " ";
 }
 }
 else
 {
 std::cout << " ";
 }
 std::cout << "| ";
 
 if( mul > 0 )
 {
 for( unsigned int j = 0; j < mul; j++ )
 {
 std::cout << std::setw( 10 ) << tlp->vul_rd[i * mul + j] << " ";
 }
 }
 else
 {
 std::cout << " ";
 }
 std::cout << "| ";
 
 if( mr > 0 )
 {
 for( unsigned int j = 0; j < mr; j++ )
 {
 std::cout << std::setw( 10 ) << tlp->vr_rd[i * mr + j] << " ";
 }
 }
 else
 {
 std::cout << " ";
 }
 
 if( i + 1 < tlp->get_n() ) std::cout << std::endl << " ";
 }
 std::cout << "]" << std::endl;
}
 
// Function to print vertex field values
int print_vertex_fields( Interface* mbi,
 std::vector< std::vector< EntityHandle > >& groups,
 Tag& norm_hdl,
 Coupler::IntegType integ_type )
{
 int err = 0;
 ErrorCode result;
 std::vector< EntityHandle >::iterator iter_j;
 
 for( unsigned int i = 0; i < groups.size(); i++ )
 {
 std::cout << " Group - " << std::endl << " ";
 for( iter_j = groups[i].begin(); iter_j != groups[i].end(); ++iter_j )
 {
 EntityHandle ehandle = ( *iter_j );
 // Check that the entity in iter_j is of the same dimension as the
 // integ_type we are performing
 int j_type = mbi->dimension_from_handle( ehandle );
 
 if( ( integ_type == Coupler::VOLUME ) && ( j_type != 3 ) ) continue;
 
 // Retrieve the vertices from the element
 const EntityHandle* conn = NULL;
 int num_verts = 0;
 result = mbi->get_connectivity( ehandle, conn, num_verts );
 if( MB_SUCCESS != result ) return 1;
 std::cout << std::fixed;
 for( int iv = 0; iv < num_verts; iv++ )
 {
 double data = 0;
 result = mbi->tag_get_data( norm_hdl, &conn[iv], 1, &data );
 if( MB_SUCCESS != result ) return 1;
 std::cout << std::setprecision( 8 ) << data << ", ";
 }
 std::cout << std::endl << " ";
 }
 std::cout << std::endl;
 std::cout.unsetf( std::ios_base::floatfield ); // turn off fixed notation
 }
 
 return err;
}
 
// Function for a constant field value
double const_field( double /*x*/, double /*y*/, double /*z*/ )
{
 // return 5.0/40.0;
 return 5.0;
}
 
// Functions for a some field values
double field_1( double x, double y, double z )
{
 double f = fabs( x ) + fabs( y ) + fabs( z );
 // return f/24.0;
 return f;
}
 
double field_2( double x, double y, double z )
{
 double f = x * x + y * y + z * z;
 // return f/32.0;
 return f;
}
 
double field_3( double x, double y, double z )
{
 double f = 2 * x + 2 * y + 2 * z;
 // return f/48.0;
 return f;
}
 
// Function used to create field on mesh for testing.
double physField( double x, double y, double z )
{
 double out;
 
 // 1/r^2 decay from {0,0,0}
 
 out = x * x + y * y + z * z;
 out += 1e-1; // clamp
 out = 1 / out;
 
 return out;
}
 
#define UINT_PER_X( X ) ( ( sizeof( X ) + sizeof( uint ) - 1 ) / sizeof( uint ) )
#define UINT_PER_REAL UINT_PER_X( realType )
#define UINT_PER_LONG UINT_PER_X( slong )
#define UINT_PER_ULONG UINT_PER_X( Ulong )
#define UINT_PER_UNSIGNED UINT_PER_X( unsigned )
 
// Function for packing tuple_list
int pack_tuples( TupleList* tl, void** ptr )
{
 uint mi, ml, mul, mr;
 tl->getTupleSize( mi, ml, mul, mr );
 
 uint n = tl->get_n();
 
 int sz_buf = 1 + 4 * UINT_PER_UNSIGNED +
 tl->get_n() * ( mi + ml * UINT_PER_LONG + mul * UINT_PER_LONG + mr * UINT_PER_REAL );
 
 uint* buf = (uint*)malloc( sz_buf * sizeof( uint ) );
 *ptr = (void*)buf;
 
 // copy n
 memcpy( buf, &n, sizeof( uint ) ), buf += 1;
 // copy mi
 memcpy( buf, &mi, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
 // copy ml
 memcpy( buf, &ml, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
 // copy mul
 memcpy( buf, &mul, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
 // copy mr
 memcpy( buf, &mr, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
 // copy vi_rd
 memcpy( buf, tl->vi_rd, tl->get_n() * mi * sizeof( sint ) ), buf += tl->get_n() * mi;
 // copy vl_rd
 memcpy( buf, tl->vl_rd, tl->get_n() * ml * sizeof( slong ) ), buf += tl->get_n() * ml * UINT_PER_LONG;
 // copy vul_rd
 memcpy( buf, tl->vul_rd, tl->get_n() * mul * sizeof( Ulong ) ), buf += tl->get_n() * mul * UINT_PER_ULONG;
 // copy vr_rd
 memcpy( buf, tl->vr_rd, tl->get_n() * mr * sizeof( realType ) ), buf += tl->get_n() * mr * UINT_PER_REAL;
 
 return sz_buf;
}
 
// Function for packing tuple_list
void unpack_tuples( void* ptr, TupleList** tlp )
{
 TupleList* tl = new TupleList();
 *tlp = tl;
 
 uint nt;
 unsigned mit, mlt, mult, mrt;
 uint* buf = (uint*)ptr;
 
 // get n
 memcpy( &nt, buf, sizeof( uint ) ), buf += 1;
 // get mi
 memcpy( &mit, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
 // get ml
 memcpy( &mlt, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
 // get mul
 memcpy( &mult, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
 // get mr
 memcpy( &mrt, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
 
 // initialize tl
 tl->initialize( mit, mlt, mult, mrt, nt );
 tl->enableWriteAccess();
 tl->set_n( nt );
 
 uint mi, ml, mul, mr;
 tl->getTupleSize( mi, ml, mul, mr );
 
 // get vi_wr
 memcpy( tl->vi_wr, buf, tl->get_n() * mi * sizeof( sint ) ), buf += tl->get_n() * mi;
 // get vl_wr
 memcpy( tl->vl_wr, buf, tl->get_n() * ml * sizeof( slong ) ), buf += tl->get_n() * ml * UINT_PER_LONG;
 // get vul_wr
 memcpy( tl->vul_wr, buf, tl->get_n() * mul * sizeof( Ulong ) ), buf += tl->get_n() * mul * UINT_PER_ULONG;
 // get vr_wr
 memcpy( tl->vr_wr, buf, tl->get_n() * mr * sizeof( realType ) ), buf += tl->get_n() * mr * UINT_PER_REAL;
 
 tl->disableWriteAccess();
 
 return;
}