This class couples data between meshes. More...

#include <Coupler.hpp>

Collaboration diagram for moab::Coupler:

Public Types
enum	Method { CONSTANT, LINEAR_FE, QUADRATIC_FE, SPECTRAL, SPHERICAL }
enum	IntegType { VOLUME }
Public Member Functions
	Coupler (Interface impl, ParallelComm pc, Range &local_elems, int coupler_id, bool init_tree=true, int max_ent_dim=3)
virtual	~Coupler ()
ErrorCode	initialize_tree ()
	Initialize the kdtree, locally and across communicator.
ErrorCode	locate_points (double xyz, unsigned int num_points, double rel_eps=0.0, double abs_eps=0.0, TupleList tl=NULL, bool store_local=true)
ErrorCode	locate_points (Range &ents, double rel_eps=0.0, double abs_eps=0.0, TupleList *tl=NULL, bool store_local=true)
ErrorCode	interpolate (Coupler::Method method, Tag tag, double interp_vals, TupleList tl=NULL, bool normalize=true)
ErrorCode	interpolate (Coupler::Method method, const std::string &tag_name, double interp_vals, TupleList tl=NULL, bool normalize=true)
ErrorCode	interpolate (Coupler::Method methods, const std::string tag_names, int points_per_method, int num_methods, double interp_vals, TupleList *tl=NULL, bool normalize=true)
ErrorCode	interpolate (Coupler::Method methods, Tag tag_names, int points_per_method, int num_methods, double interp_vals, TupleList *tl=NULL, bool normalize=true)
ErrorCode	normalize_mesh (EntityHandle root_set, const char *norm_tag, Coupler::IntegType integ_type, int num_integ_pts)
ErrorCode	normalize_subset (EntityHandle root_set, const char norm_tag, const char tag_names, int num_tags, const char *tag_values, Coupler::IntegType integ_type, int num_integ_pts)
ErrorCode	normalize_subset (EntityHandle root_set, const char norm_tag, Tag tag_handles, int num_tags, const char **tag_values, Coupler::IntegType integ_type, int num_integ_pts)
ErrorCode	do_normalization (const char *norm_tag, std::vector< std::vector< EntityHandle > > &entity_sets, std::vector< std::vector< EntityHandle > > &entity_groups, Coupler::IntegType integ_type, int num_integ_pts)
ErrorCode	get_matching_entities (EntityHandle root_set, const char tag_names, const char tag_values, int num_tags, std::vector< std::vector< EntityHandle > > entity_sets, std::vector< std::vector< EntityHandle > > entity_groups)
ErrorCode	get_matching_entities (EntityHandle root_set, Tag tag_handles, const char tag_values, int num_tags, std::vector< std::vector< EntityHandle > > entity_sets, std::vector< std::vector< EntityHandle > > *entity_groups)
ErrorCode	create_tuples (Range &ent_sets, const char tag_names, unsigned int num_tags, TupleList tuples)
ErrorCode	create_tuples (Range &ent_sets, Tag tag_handles, unsigned int num_tags, TupleList *tuples)
ErrorCode	consolidate_tuples (TupleList all_tuples, unsigned int num_tuples, TupleList unique_tuples)
ErrorCode	get_group_integ_vals (std::vector< std::vector< EntityHandle > > &groups, std::vector< double > &integ_vals, const char *norm_tag, int num_integ_pts, Coupler::IntegType integ_type)
ErrorCode	apply_group_norm_factor (std::vector< std::vector< EntityHandle > > &entity_sets, std::vector< double > &norm_factors, const char *norm_tag, Coupler::IntegType integ_type)
ErrorCode	initialize_spectral_elements (EntityHandle rootSource, EntityHandle rootTarget, bool &specSou, bool &specTar)
ErrorCode	get_gl_points_on_elements (Range &targ_elems, std::vector< double > &vpos, int &numPointsOfInterest)
Interface *	mb_impl () const
AdaptiveKDTree *	my_tree () const
EntityHandle	local_root () const
const std::vector< double > &	all_boxes () const
ParallelComm *	my_pc () const
const Range &	target_ents () const
int	my_id () const
const Range &	my_range () const
TupleList *	mapped_pts () const
int	num_its () const
void	set_spherical (bool arg1=true)
Private Member Functions
ErrorCode	nat_param (double xyz[3], std::vector< EntityHandle > &entities, std::vector< CartVect > &nat_coords, double epsilon=0.0)
ErrorCode	interp_field (EntityHandle elem, CartVect nat_coord, Tag tag, double &field)
ErrorCode	constant_interp (EntityHandle elem, Tag tag, double &field)
ErrorCode	test_local_box (double xyz, int from_proc, int remote_index, int index, bool &point_located, double rel_eps=0.0, double abs_eps=0.0, TupleList tl=NULL)
Private Attributes
Interface *	mbImpl
AdaptiveKDTree *	myTree
EntityHandle	localRoot
std::vector< double >	allBoxes
ParallelComm *	myPc
int	myId
Range	myRange
Range	targetEnts
TupleList *	mappedPts
TupleList *	targetPts
int	numIts
int	max_dim
void *	_spectralSource
void *	_spectralTarget
moab::Tag	_xm1Tag
moab::Tag	_ym1Tag
moab::Tag	_zm1Tag
int	_ntot
bool	spherical

Detailed Description

This class couples data between meshes.

Author:: Tim Tautges

The coupler interpolates solution data at a set of points. Data being interpolated resides on a source mesh, in a tag. Applications calling this coupler send in entities, usually points or vertices, and receive back the tag value interpolated at those points. Entities in the source mesh containing those points do not have to reside on the same processor.

To use, an application should:

instantiate this coupler by calling the constructor collectively on all processors in the communicator
call locate_points, which locates the points to be interpolated and (optionally) caches the results in this object
call interpolate, which does the interpolation

Multiple interpolations can be done after locating the points.

Definition at line 43 of file Coupler.hpp.

Member Enumeration Documentation

enum moab::Coupler::IntegType

Enumerator:

VOLUME

Definition at line 55 of file Coupler.hpp.

    {
        VOLUME
    };

enum moab::Coupler::Method

Enumerator:

CONSTANT
LINEAR_FE
QUADRATIC_FE
SPECTRAL
SPHERICAL

Definition at line 46 of file Coupler.hpp.

    {
        CONSTANT,
        LINEAR_FE,
        QUADRATIC_FE,
        SPECTRAL,
        SPHERICAL
    };

Constructor & Destructor Documentation

moab::Coupler::Coupler	(	Interface *	impl,
		ParallelComm *	pc,
		Range &	local_elems,
		int	coupler_id,
		bool	init_tree = `true`,
		int	max_ent_dim = `3`
	)

Definition at line 47 of file Coupler.cpp.

References _spectralSource, _spectralTarget, moab::Range::empty(), initialize_tree(), mappedPts, myRange, myTree, and targetPts.

    : mbImpl( impl ), myPc( pc ), myId( coupler_id ), numIts( 3 ), max_dim( max_ent_dim ), _ntot( 0 ),
      spherical( false )
{
    assert( NULL != impl && ( pc || !local_elems.empty() ) );

    // Keep track of the local points, at least for now
    myRange = local_elems;
    myTree  = NULL;

    // Now initialize the tree
    if( init_tree ) initialize_tree();

    // Initialize tuple lists to indicate not initialized
    mappedPts       = NULL;
    targetPts       = NULL;
    _spectralSource = _spectralTarget = NULL;
}

moab::Coupler::~Coupler ( ) [virtual]

Definition at line 73 of file Coupler.cpp.

References _spectralSource, _spectralTarget, mappedPts, myTree, and targetPts.

{
    // This will clear the cache
    delete(moab::Element::SpectralHex*)_spectralSource;
    delete(moab::Element::SpectralHex*)_spectralTarget;
    delete myTree;
    delete targetPts;
    delete mappedPts;
}

Member Function Documentation

const std::vector< double >& moab::Coupler::all_boxes ( ) const [inline]

Definition at line 441 of file Coupler.hpp.

References allBoxes.

    {
        return allBoxes;
    }

ErrorCode moab::Coupler::apply_group_norm_factor	(	std::vector< std::vector< EntityHandle > > &	entity_sets,
		std::vector< double > &	norm_factors,
		const char *	norm_tag,
		Coupler::IntegType	integ_type
	)

Definition at line 1787 of file Coupler.cpp.

References ErrorCode, ERRORR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_SPARSE, MB_TYPE_DOUBLE, mbImpl, moab::Interface::tag_get_handle(), and moab::Interface::tag_set_data().

Referenced by do_normalization(), and main().

{
    ErrorCode err;

    // Construct the new tag for the normalization factor from the norm_tag name
    // and "_normf".
    int norm_tag_len       = strlen( norm_tag );
    const char* normf_appd = "_normf";
    int normf_appd_len     = strlen( normf_appd );

    char* normf_tag = (char*)malloc( norm_tag_len + normf_appd_len + 1 );
    char* tmp_ptr   = normf_tag;

    memcpy( tmp_ptr, norm_tag, norm_tag_len );
    tmp_ptr += norm_tag_len;
    memcpy( tmp_ptr, normf_appd, normf_appd_len );
    tmp_ptr += normf_appd_len;
    *tmp_ptr = '\0';

    Tag normf_hdl;
    // Check to see if the tag exists.  If not then create it and get the handle.
    err = mbImpl->tag_get_handle( normf_tag, 1, moab::MB_TYPE_DOUBLE, normf_hdl,
                                  moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT );ERRORR( "Failed to create normalization factor tag.", err );
    if( normf_hdl == NULL )
    {
        std::string msg( "Failed to create normalization factor tag named '" );
        msg += std::string( normf_tag ) + std::string( "'" );ERRORR( msg.c_str(), MB_FAILURE );
    }
    free( normf_tag );

    std::vector< std::vector< EntityHandle > >::iterator iter_i;
    std::vector< EntityHandle >::iterator iter_j;
    std::vector< double >::iterator iter_f;
    double grp_norm_factor = 0.0;

    // Loop over the entity sets
    for( iter_i = entity_sets.begin(), iter_f = norm_factors.begin();
         ( iter_i != entity_sets.end() ) && ( iter_f != norm_factors.end() ); ++iter_i, ++iter_f )
    {
        grp_norm_factor = *iter_f;

        // Loop over the all the entity sets in iter_i and set the
        // new normf_tag with the norm factor value on each
        for( iter_j = ( *iter_i ).begin(); iter_j != ( *iter_i ).end(); ++iter_j )
        {
            EntityHandle entset = *iter_j;

            std::cout << "Coupler: applying normalization for entity set=" << entset
                      << ",  normalization_factor=" << grp_norm_factor << std::endl;

            err = mbImpl->tag_set_data( normf_hdl, &entset, 1, &grp_norm_factor );ERRORR( "Failed to set normalization factor on entity set.", err );
        }
    }

    return MB_SUCCESS;
}

ErrorCode moab::Coupler::consolidate_tuples	(	TupleList **	all_tuples,
		unsigned int	num_tuples,
		TupleList **	unique_tuples
	)

Definition at line 1563 of file Coupler.cpp.

References moab::TupleList::enableWriteAccess(), moab::TupleList::get_n(), moab::TupleList::getTupleSize(), MB_SUCCESS, moab::TupleList::resize(), moab::TupleList::set_n(), sint, moab::TupleList::sort(), moab::TupleList::vi_rd, and moab::TupleList::vi_wr.

Referenced by get_matching_entities(), and main().

{
    int total_rcv_tuples = 0;
    int offset = 0, copysz = 0;
    unsigned num_tags = 0;

    uint ml, mul, mr;
    uint* mi = (uint*)malloc( sizeof( uint ) * num_tuples );

    for( unsigned int i = 0; i < num_tuples; i++ )
    {
        all_tuples[i]->getTupleSize( mi[i], ml, mul, mr );
    }

    for( unsigned int i = 0; i < num_tuples; i++ )
    {
        if( all_tuples[i] != NULL )
        {
            total_rcv_tuples += all_tuples[i]->get_n();
            num_tags = mi[i];
        }
    }
    const unsigned int_size  = sizeof( sint );
    const unsigned int_width = num_tags * int_size;

    // Get the total size of all of the tuple_lists in all_tuples.
    for( unsigned int i = 0; i < num_tuples; i++ )
    {
        if( all_tuples[i] != NULL ) total_rcv_tuples += all_tuples[i]->get_n();
    }

    // Copy the tuple_lists into a single tuple_list.
    TupleList* all_tuples_list = new TupleList( num_tags, 0, 0, 0, total_rcv_tuples );
    all_tuples_list->enableWriteAccess();
    // all_tuples_list->initialize(num_tags, 0, 0, 0, total_rcv_tuples);
    for( unsigned int i = 0; i < num_tuples; i++ )
    {
        if( all_tuples[i] != NULL )
        {
            copysz = all_tuples[i]->get_n() * int_width;
            memcpy( all_tuples_list->vi_wr + offset, all_tuples[i]->vi_rd, copysz );
            offset = offset + ( all_tuples[i]->get_n() * mi[i] );
            all_tuples_list->set_n( all_tuples_list->get_n() + all_tuples[i]->get_n() );
        }
    }

    // Sort the new tuple_list. Use a radix type sort, starting with the last (or least significant)
    // tag column in the vi array and working towards the first (or most significant) tag column.
    TupleList::buffer sort_buffer;
    sort_buffer.buffer_init( 2 * total_rcv_tuples * int_width );
    for( int i = num_tags - 1; i >= 0; i-- )
    {
        all_tuples_list->sort( i, &sort_buffer );
    }

    // Cycle through the sorted list eliminating duplicates.
    // Keep counters to the current end of the tuple_list (w/out dups) and the last tuple examined.
    unsigned int end_idx = 0, last_idx = 1;
    while( last_idx < all_tuples_list->get_n() )
    {
        if( memcmp( all_tuples_list->vi_rd + ( end_idx * num_tags ), all_tuples_list->vi_rd + ( last_idx * num_tags ),
                    int_width ) == 0 )
        {
            // Values equal - skip
            last_idx += 1;
        }
        else
        {
            // Values different - copy
            // Move up the end index
            end_idx += 1;
            memcpy( all_tuples_list->vi_wr + ( end_idx * num_tags ), all_tuples_list->vi_rd + ( last_idx * num_tags ),
                    int_width );
            last_idx += 1;
        }
    }
    // Update the count in all_tuples_list
    all_tuples_list->set_n( end_idx + 1 );

    // Resize the tuple_list
    all_tuples_list->resize( all_tuples_list->get_n() );

    // Set the output parameter
    *unique_tuples = all_tuples_list;

    return MB_SUCCESS;
}

ErrorCode moab::Coupler::constant_interp	(	EntityHandle	elem,
		Tag	tag,
		double &	field
	)		`[private]`

Definition at line 1111 of file Coupler.cpp.

References ErrorCode, MB_SUCCESS, mbImpl, and moab::Interface::tag_get_data().

Referenced by interpolate().

{
    double tempField;

    // Get the tag values at the vertices
    ErrorCode result = mbImpl->tag_get_data( tag, &elem, 1, &tempField );
    if( MB_SUCCESS != result ) return result;

    field = tempField;

    return MB_SUCCESS;
}

ErrorCode moab::Coupler::create_tuples	(	Range &	ent_sets,
		const char **	tag_names,
		unsigned int	num_tags,
		TupleList **	tuples
	)

Definition at line 1506 of file Coupler.cpp.

References ErrorCode, ERRORR, MB_TAG_ANY, MB_TYPE_DOUBLE, mbImpl, t, and moab::Interface::tag_get_handle().

Referenced by get_matching_entities(), and main().

{
    ErrorCode err;
    std::vector< Tag > tag_handles;

    for( unsigned int t = 0; t < num_tags; t++ )
    {
        // Get tag handle & size
        Tag th;
        err = mbImpl->tag_get_handle( tag_names[t], 1, moab::MB_TYPE_DOUBLE, th, moab::MB_TAG_ANY );ERRORR( "Failed to get tag handle.", err );
        tag_handles.push_back( th );
    }

    return create_tuples( ent_sets, &tag_handles[0], num_tags, tuple_list );
}

ErrorCode moab::Coupler::create_tuples	(	Range &	ent_sets,
		Tag *	tag_handles,
		unsigned int	num_tags,
		TupleList **	tuples
	)

Definition at line 1528 of file Coupler.cpp.

References moab::TupleList::disableWriteAccess(), moab::TupleList::enableWriteAccess(), ErrorCode, ERRORR, moab::TupleList::getTupleSize(), moab::TupleList::inc_n(), MB_SUCCESS, mbImpl, moab::Range::size(), moab::Interface::tag_get_data(), and moab::TupleList::vi_wr.

{
    // ASSUMPTION: All tags are of type integer.  This may need to be expanded in future.
    ErrorCode err;

    // Allocate a tuple_list for the number of entity sets passed in
    TupleList* tag_tuples = new TupleList( num_tags, 0, 0, 0, (int)ent_sets.size() );
    // tag_tuples->initialize(num_tags, 0, 0, 0, num_sets);
    uint mi, ml, mul, mr;
    tag_tuples->getTupleSize( mi, ml, mul, mr );
    tag_tuples->enableWriteAccess();

    if( mi == 0 ) ERRORR( "Failed to initialize tuple_list.", MB_FAILURE );

    // Loop over the filtered entity sets retrieving each matching tag value one by one.
    int val;
    for( unsigned int i = 0; i < ent_sets.size(); i++ )
    {
        for( unsigned int j = 0; j < num_tags; j++ )
        {
            EntityHandle set_handle = ent_sets[i];
            err                     = mbImpl->tag_get_data( tag_handles[j], &set_handle, 1, &val );ERRORR( "Failed to get integer tag data.", err );
            tag_tuples->vi_wr[i * mi + j] = val;
        }

        // If we get here there was no error so increment n in the tuple_list
        tag_tuples->inc_n();
    }
    tag_tuples->disableWriteAccess();
    *tuples = tag_tuples;

    return MB_SUCCESS;
}

ErrorCode moab::Coupler::do_normalization	(	const char *	norm_tag,
		std::vector< std::vector< EntityHandle > > &	entity_sets,
		std::vector< std::vector< EntityHandle > > &	entity_groups,
		Coupler::IntegType	integ_type,
		int	num_integ_pts
	)

Definition at line 1203 of file Coupler.cpp.

References apply_group_norm_factor(), ErrorCode, ERRORMPI, ERRORR, get_group_integ_vals(), MASTER_PROC, myPc, moab::ProcConfig::proc_comm(), and moab::ParallelComm::proc_config().

Referenced by normalize_mesh(), and normalize_subset().

{
    // SLAVE START ****************************************************************
    ErrorCode err;
    int ierr = 0;

    // Setup data for parallel computing
    int nprocs, rank;
    ierr = MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
    ERRORMPI( "Getting number of procs failed.", ierr );
    ierr = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
    ERRORMPI( "Getting rank failed.", ierr );

    // Get the integrated field value for each group(vector) of entities.
    // If no entities are in a group then a zero will be put in the list
    // of return values.
    unsigned int num_ent_grps = entity_groups.size();
    std::vector< double > integ_vals( num_ent_grps );

    err = get_group_integ_vals( entity_groups, integ_vals, norm_tag, num_integ_pts, integ_type );ERRORR( "Failed to get integrated field values for groups in mesh.", err );
    // SLAVE END   ****************************************************************

    // SLAVE/MASTER START #########################################################
    // Send list of integrated values back to master proc. The ordering of the
    // values will match the ordering of the entity groups (i.e. vector of vectors)
    // sent from master to slaves earlier.  The values for each entity group will
    // be summed during the transfer.
    std::vector< double > sum_integ_vals( num_ent_grps );

    if( nprocs > 1 )
    {
        // If parallel then send the values back to the master.
        ierr = MPI_Reduce( &integ_vals[0], &sum_integ_vals[0], num_ent_grps, MPI_DOUBLE, MPI_SUM, MASTER_PROC,
                           myPc->proc_config().proc_comm() );
        ERRORMPI( "Transfer and reduction of integrated values failed.", ierr );
    }
    else
    {
        // Otherwise just copy the vector
        sum_integ_vals = integ_vals;
    }
    // SLAVE/MASTER END   #########################################################

    // MASTER START ***************************************************************
    // Calculate the normalization factor for each group by taking the
    // inverse of each integrated field value. Put the normalization factor
    // for each group back into the list in the same order.
    for( unsigned int i = 0; i < num_ent_grps; i++ )
    {
        double val = sum_integ_vals[i];
        if( fabs( val ) > 1e-8 )
            sum_integ_vals[i] = 1.0 / val;
        else
        {
            sum_integ_vals[i] = 0.0; /* VSM: not sure what we should do here ? */
            /* commenting out error below since if integral(value)=0.0, then normalization
               is probably unnecessary to start with ? */
            /* ERRORR("Integrating an invalid field -- integral("<<norm_tag<<") = "<<val<<".", err);
             */
        }
    }
    // MASTER END   ***************************************************************

    // MASTER/SLAVE START #########################################################
    if( nprocs > 1 )
    {
        // If parallel then broadcast the normalization factors to the procs.
        ierr = MPI_Bcast( &sum_integ_vals[0], num_ent_grps, MPI_DOUBLE, MASTER_PROC, myPc->proc_config().proc_comm() );
        ERRORMPI( "Broadcast of normalization factors failed.", ierr );
    }
    // MASTER/SLAVE END   #########################################################

    // SLAVE START ****************************************************************
    // Save the normalization factors to a new tag with name of norm_tag's value
    // and the string "_normF" appended.  This new tag will be created on the entity
    // set that contains all of the entities from a group.

    err = apply_group_norm_factor( entity_sets, sum_integ_vals, norm_tag, integ_type );ERRORR( "Failed to set the normalization factor for groups in mesh.", err );
    // SLAVE END   ****************************************************************

    return err;
}

ErrorCode moab::Coupler::get_gl_points_on_elements	(	Range &	targ_elems,
		std::vector< double > &	vpos,
		int &	numPointsOfInterest
	)

Definition at line 1930 of file Coupler.cpp.

References _ntot, _xm1Tag, _ym1Tag, _zm1Tag, moab::Range::begin(), moab::Range::end(), ErrorCode, MB_SUCCESS, mbImpl, moab::Range::size(), and moab::Interface::tag_get_by_ptr().

{
    numPointsOfInterest = targ_elems.size() * _ntot;
    vpos.resize( 3 * numPointsOfInterest );
    int ielem = 0;
    for( Range::iterator eit = targ_elems.begin(); eit != targ_elems.end(); ++eit, ielem += _ntot * 3 )
    {
        EntityHandle eh = *eit;
        const double* xval;
        const double* yval;
        const double* zval;
        ErrorCode rval = mbImpl->tag_get_by_ptr( _xm1Tag, &eh, 1, (const void**)&xval );
        if( moab::MB_SUCCESS != rval )
        {
            std::cout << "Can't get xm1 values \n";
            return MB_FAILURE;
        }
        rval = mbImpl->tag_get_by_ptr( _ym1Tag, &eh, 1, (const void**)&yval );
        if( moab::MB_SUCCESS != rval )
        {
            std::cout << "Can't get ym1 values \n";
            return MB_FAILURE;
        }
        rval = mbImpl->tag_get_by_ptr( _zm1Tag, &eh, 1, (const void**)&zval );
        if( moab::MB_SUCCESS != rval )
        {
            std::cout << "Can't get zm1 values \n";
            return MB_FAILURE;
        }
        // Now, in a stride, populate vpos
        for( int i = 0; i < _ntot; i++ )
        {
            vpos[ielem + 3 * i]     = xval[i];
            vpos[ielem + 3 * i + 1] = yval[i];
            vpos[ielem + 3 * i + 2] = zval[i];
        }
    }

    return MB_SUCCESS;
}

ErrorCode moab::Coupler::get_group_integ_vals	(	std::vector< std::vector< EntityHandle > > &	groups,
		std::vector< double > &	integ_vals,
		const char *	norm_tag,
		int	num_integ_pts,
		Coupler::IntegType	integ_type
	)

Definition at line 1652 of file Coupler.cpp.

References ErrorCode, ERRORR, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Element::Map::integrate_scalar_field(), MB_SUCCESS, MB_TAG_CREAT, MB_TAG_SPARSE, MB_TYPE_DOUBLE, MB_UNSUPPORTED_OPERATION, MBEDGE, MBENTITYSET, MBHEX, mbImpl, MBQUAD, MBTET, moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), moab::Interface::type_from_handle(), and VOLUME.

Referenced by do_normalization(), and main().

{
    ErrorCode err;

    std::vector< std::vector< EntityHandle > >::iterator iter_i;
    std::vector< EntityHandle >::iterator iter_j;
    double grp_intrgr_val, intgr_val;

    // Get the tag handle for norm_tag
    Tag norm_hdl;
    err =
        mbImpl->tag_get_handle( norm_tag, 1, moab::MB_TYPE_DOUBLE, norm_hdl, moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT );ERRORR( "Failed to get norm_tag handle.", err );

    // Check size of integ_vals vector
    if( integ_vals.size() != groups.size() ) integ_vals.resize( groups.size() );

    // Loop over the groups(vectors) of entities
    unsigned int i;
    for( i = 0, iter_i = groups.begin(); iter_i != groups.end(); i++, ++iter_i )
    {
        grp_intrgr_val = 0;

        // Loop over the all the entities in the group, integrating
        // the field_fn over the entity in iter_j
        for( iter_j = ( *iter_i ).begin(); iter_j != ( *iter_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
            EntityType j_type;
            j_type = mbImpl->type_from_handle( ehandle );ERRORR( "Failed to get entity type.", err );
            // Skip any entities in the group that are not of the type being considered
            if( ( integ_type == VOLUME ) && ( j_type < MBTET || j_type >= MBENTITYSET ) ) continue;

            intgr_val = 0;

            // Retrieve the vertices from the element
            const EntityHandle* verts = NULL;
            int connectivity_size     = 0;

            err = mbImpl->get_connectivity( ehandle, verts, connectivity_size, false );ERRORR( "Failed to get vertices from entity.", err );

            // Get the vertex coordinates and the field values at the vertices.
            double* coords = (double*)malloc( sizeof( double ) * ( 3 * connectivity_size ) );
            /* TODO: VSM: check if this works for lower dimensions also without problems */
            /* if (3 == geom_dim) */
            err = mbImpl->get_coords( verts, connectivity_size, coords );ERRORR( "Failed to get vertex coordinates.", err );

            /* allocate the field data array */
            double* vfield = (double*)malloc( sizeof( double ) * ( connectivity_size ) );
            err            = mbImpl->tag_get_data( norm_hdl, verts, connectivity_size, vfield );
            if( MB_SUCCESS != err )
            {
                free( coords );
            }
            ERRORR( "Failed to get vertex coordinates.", err );

            // Get coordinates of all corner vertices (in normal order) and
            // put in array of CartVec.
            std::vector< CartVect > vertices( connectivity_size );

            // Put the vertices into a CartVect vector
            double* x = coords;
            for( int j = 0; j < connectivity_size; j++, x += 3 )
            {
                vertices[j] = CartVect( x );
            }
            free( coords );

            moab::Element::Map* elemMap;
            if( j_type == MBHEX )
            {
                if( connectivity_size == 8 )
                    elemMap = new moab::Element::LinearHex( vertices );
                else
                    elemMap = new moab::Element::QuadraticHex( vertices );
            }
            else if( j_type == MBTET )
            {
                elemMap = new moab::Element::LinearTet( vertices );
            }
            else if( j_type == MBQUAD )
            {
                elemMap = new moab::Element::LinearQuad( vertices );
            }
            /*
            else if (j_type == MBTRI) {
              elemMap = new moab::Element::LinearTri(vertices);
            }
            */
            else if( j_type == MBEDGE )
            {
                elemMap = new moab::Element::LinearEdge( vertices );
            }
            else
                ERRORR( "Unknown topology type.", MB_UNSUPPORTED_OPERATION );

            // Set the vertices in the Map and perform the integration
            try
            {
                /* VSM: Do we need this call ?? */
                // elemMap->set_vertices(vertices);

                // Perform the actual integration over the element
                intgr_val = elemMap->integrate_scalar_field( vfield );

                // Combine the result with those of the group
                grp_intrgr_val += intgr_val;
            }
            catch( moab::Element::Map::ArgError& )
            {
                std::cerr << "Failed to set vertices on Element::Map." << std::endl;
            }
            catch( moab::Element::Map::EvaluationError& )
            {
                std::cerr << "Failed to get inverse evaluation of coordinate on Element::Map." << std::endl;
            }

            delete( elemMap );
            free( vfield );
        }

        // Set the group integrated value in the vector
        integ_vals[i] = grp_intrgr_val;
    }

    return err;
}

ErrorCode moab::Coupler::get_matching_entities	(	EntityHandle	root_set,
		const char **	tag_names,
		const char **	tag_values,
		int	num_tags,
		std::vector< std::vector< EntityHandle > > *	entity_sets,
		std::vector< std::vector< EntityHandle > > *	entity_groups
	)

Definition at line 1293 of file Coupler.cpp.

References ErrorCode, ERRORR, MB_TAG_ANY, MB_TYPE_DOUBLE, mbImpl, t, and moab::Interface::tag_get_handle().

Referenced by main(), and normalize_subset().

{
    ErrorCode err;
    std::vector< Tag > tag_handles;

    for( int t = 0; t < num_tags; t++ )
    {
        // Get tag handle & size
        Tag th;
        err = mbImpl->tag_get_handle( tag_names[t], 1, moab::MB_TYPE_DOUBLE, th, moab::MB_TAG_ANY );ERRORR( "Failed to get tag handle.", err );
        tag_handles.push_back( th );
    }

    return get_matching_entities( root_set, &tag_handles[0], tag_values, num_tags, entity_sets, entity_groups );
}

ErrorCode moab::Coupler::get_matching_entities	(	EntityHandle	root_set,
		Tag *	tag_handles,
		const char **	tag_values,
		int	num_tags,
		std::vector< std::vector< EntityHandle > > *	entity_sets,
		std::vector< std::vector< EntityHandle > > *	entity_groups
	)

Definition at line 1315 of file Coupler.cpp.

References moab::Range::clear(), consolidate_tuples(), create_tuples(), moab::debug, ErrorCode, ERRORMPI, ERRORR, moab::Interface::get_entities_by_handle(), moab::Interface::get_entities_by_type_and_tag(), moab::TupleList::get_n(), moab::TupleList::getTupleSize(), moab::Interface::INTERSECT, MASTER_PROC, MBENTITYSET, mbImpl, myPc, moab::pack_tuples(), moab::ProcConfig::proc_comm(), moab::ParallelComm::proc_config(), moab::TupleList::reset(), moab::Range::size(), moab::unpack_tuples(), and moab::TupleList::vi_rd.

{
    // SLAVE START ****************************************************************

    // Setup data for parallel computing
    ErrorCode err;
    int ierr = 0;
    int nprocs, rank;
    ierr = MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
    ERRORMPI( "Getting number of procs failed.", ierr );
    ierr = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
    ERRORMPI( "Getting rank failed.", ierr );

    Range ent_sets;
    err =
        mbImpl->get_entities_by_type_and_tag( root_set, moab::MBENTITYSET, tag_handles, (const void* const*)tag_values,
                                              num_tags, ent_sets, Interface::INTERSECT, false );ERRORR( "Core::get_entities_by_type_and_tag failed.", err );

    TupleList* tag_list = NULL;
    err                 = create_tuples( ent_sets, tag_handles, num_tags, &tag_list );ERRORR( "Failed to create tuples from entity sets.", err );

    // Free up range
    ent_sets.clear();
    // SLAVE END   ****************************************************************

    // If we are running in a multi-proc session then send tuple list back to master
    // proc for consolidation. Otherwise just copy the pointer to the tuple_list.
    TupleList* cons_tuples;
    if( nprocs > 1 )
    {
        // SLAVE/MASTER START #########################################################

        // Pack the tuple_list in a buffer.
        uint* tuple_buf;
        int tuple_buf_len;
        tuple_buf_len = pack_tuples( tag_list, (void**)&tuple_buf );

        // Free tag_list here as its not used again if nprocs > 1
        tag_list->reset();

        // Send back the buffer sizes to the master proc
        int* recv_cnts       = (int*)malloc( nprocs * sizeof( int ) );
        int* offsets         = (int*)malloc( nprocs * sizeof( int ) );
        uint* all_tuples_buf = NULL;

        MPI_Gather( &tuple_buf_len, 1, MPI_INT, recv_cnts, 1, MPI_INT, MASTER_PROC, myPc->proc_config().proc_comm() );
        ERRORMPI( "Gathering buffer sizes failed.", err );

        // Allocate a buffer large enough for all the data
        if( rank == MASTER_PROC )
        {
            int all_tuples_len = recv_cnts[0];
            offsets[0]         = 0;
            for( int i = 1; i < nprocs; i++ )
            {
                offsets[i] = offsets[i - 1] + recv_cnts[i - 1];
                all_tuples_len += recv_cnts[i];
            }

            all_tuples_buf = (uint*)malloc( all_tuples_len * sizeof( uint ) );
        }

        // Send all buffers to the master proc for consolidation
        MPI_Gatherv( (void*)tuple_buf, tuple_buf_len, MPI_INT, (void*)all_tuples_buf, recv_cnts, offsets, MPI_INT,
                     MASTER_PROC, myPc->proc_config().proc_comm() );
        ERRORMPI( "Gathering tuple_lists failed.", err );
        free( tuple_buf );  // malloc'd in pack_tuples

        if( rank == MASTER_PROC )
        {
            // Unpack the tuple_list from the buffer.
            TupleList** tl_array = (TupleList**)malloc( nprocs * sizeof( TupleList* ) );
            for( int i = 0; i < nprocs; i++ )
                unpack_tuples( (void*)&all_tuples_buf[offsets[i]], &tl_array[i] );

            // Free all_tuples_buf here as it is only allocated on the MASTER_PROC
            free( all_tuples_buf );
            // SLAVE/MASTER END   #########################################################

            // MASTER START ***************************************************************
            // Consolidate all tuple_lists into one tuple_list with no duplicates.
            err = consolidate_tuples( tl_array, nprocs, &cons_tuples );ERRORR( "Failed to consolidate tuples.", err );

            for( int i = 0; i < nprocs; i++ )
                tl_array[i]->reset();
            free( tl_array );
            // MASTER END   ***************************************************************
        }

        // Free offsets and recv_cnts as they are allocated on all procs
        free( offsets );
        free( recv_cnts );

        // MASTER/SLAVE START #########################################################
        // Broadcast condensed tuple list back to all procs.
        uint* ctl_buf;
        int ctl_buf_sz;
        if( rank == MASTER_PROC ) ctl_buf_sz = pack_tuples( cons_tuples, (void**)&ctl_buf );

        // Send buffer size
        ierr = MPI_Bcast( &ctl_buf_sz, 1, MPI_INT, MASTER_PROC, myPc->proc_config().proc_comm() );
        ERRORMPI( "Broadcasting tuple_list size failed.", ierr );

        // Allocate a buffer in the other procs
        if( rank != MASTER_PROC ) ctl_buf = (uint*)malloc( ctl_buf_sz * sizeof( uint ) );

        ierr = MPI_Bcast( (void*)ctl_buf, ctl_buf_sz, MPI_INT, MASTER_PROC, myPc->proc_config().proc_comm() );
        ERRORMPI( "Broadcasting tuple_list failed.", ierr );

        if( rank != MASTER_PROC ) unpack_tuples( ctl_buf, &cons_tuples );
        free( ctl_buf );
        // MASTER/SLAVE END   #########################################################
    }
    else
        cons_tuples = tag_list;

    // SLAVE START ****************************************************************
    // Loop over the tuple list getting the entities with the tags in the tuple_list entry
    uint mi, ml, mul, mr;
    cons_tuples->getTupleSize( mi, ml, mul, mr );

    for( unsigned int i = 0; i < cons_tuples->get_n(); i++ )
    {
        // Get Entity Sets that match the tags and values.

        // Convert the data in the tuple_list to an array of pointers to the data
        // in the tuple_list as that is what the iMesh API call is expecting.
        int** vals = (int**)malloc( mi * sizeof( int* ) );
        for( unsigned int j = 0; j < mi; j++ )
            vals[j] = (int*)&( cons_tuples->vi_rd[( i * mi ) + j] );

        // Get entities recursively based on type and tag data
        err = mbImpl->get_entities_by_type_and_tag( root_set, moab::MBENTITYSET, tag_handles, (const void* const*)vals,
                                                    mi, ent_sets, Interface::INTERSECT, false );ERRORR( "Core::get_entities_by_type_and_tag failed.", err );
        if( debug ) std::cout << "ent_sets_size=" << ent_sets.size() << std::endl;

        // Free up the array of pointers
        free( vals );

        // Loop over the entity sets and then free the memory for ent_sets.
        std::vector< EntityHandle > ent_set_hdls;
        std::vector< EntityHandle > ent_hdls;
        for( unsigned int j = 0; j < ent_sets.size(); j++ )
        {
            // Save the entity set
            ent_set_hdls.push_back( ent_sets[j] );

            // Get all entities for the entity set
            Range ents;

            /* VSM: do we need to filter out entity sets ? */
            err = mbImpl->get_entities_by_handle( ent_sets[j], ents, false );ERRORR( "Core::get_entities_by_handle failed.", err );
            if( debug ) std::cout << "ents_size=" << ents.size() << std::endl;

            // Save all of the entities from the entity set and free the memory for ents.
            for( unsigned int k = 0; k < ents.size(); k++ )
            {
                ent_hdls.push_back( ents[k] );
            }
            ents.clear();
            if( debug ) std::cout << "ent_hdls.size=" << ent_hdls.size() << std::endl;
        }

        // Free the entity set list for next tuple iteration.
        ent_sets.clear();

        // Push ent_set_hdls onto entity_sets, ent_hdls onto entity_groups
        // and clear both ent_set_hdls and ent_hdls.
        entity_sets->push_back( ent_set_hdls );
        ent_set_hdls.clear();
        entity_groups->push_back( ent_hdls );
        ent_hdls.clear();
        if( debug )
            std::cout << "entity_sets->size=" << entity_sets->size()
                      << ", entity_groups->size=" << entity_groups->size() << std::endl;
    }

    cons_tuples->reset();
    // SLAVE END   ****************************************************************

    return err;
}

ErrorCode moab::Coupler::initialize_spectral_elements	(	EntityHandle	rootSource,
		EntityHandle	rootTarget,
		bool &	specSou,
		bool &	specTar
	)

Definition at line 196 of file Coupler.cpp.

References _ntot, _spectralSource, _spectralTarget, _xm1Tag, _ym1Tag, _zm1Tag, moab::Range::empty(), ErrorCode, moab::Interface::get_entities_by_type_and_tag(), MB_SUCCESS, MB_TYPE_DOUBLE, MB_TYPE_INTEGER, MBENTITYSET, mbImpl, moab::Interface::tag_get_data(), and moab::Interface::tag_get_handle().

{
    /*void * _spectralSource;
      void * _spectralTarget;*/

    moab::Range spectral_sets;
    moab::Tag sem_tag;
    int sem_dims[3];
    ErrorCode rval = mbImpl->tag_get_handle( "SEM_DIMS", 3, moab::MB_TYPE_INTEGER, sem_tag );
    if( moab::MB_SUCCESS != rval )
    {
        std::cout << "Can't find tag, no spectral set\n";
        return MB_SUCCESS;  // Nothing to do, no spectral elements
    }
    rval = mbImpl->get_entities_by_type_and_tag( rootSource, moab::MBENTITYSET, &sem_tag, NULL, 1, spectral_sets );
    if( moab::MB_SUCCESS != rval || spectral_sets.empty() )
        std::cout << "Can't get sem set on source\n";
    else
    {
        moab::EntityHandle firstSemSet = spectral_sets[0];
        rval                           = mbImpl->tag_get_data( sem_tag, &firstSemSet, 1, (void*)sem_dims );
        if( moab::MB_SUCCESS != rval ) return MB_FAILURE;

        if( sem_dims[0] != sem_dims[1] || sem_dims[0] != sem_dims[2] )
        {
            std::cout << " dimensions are different. bail out\n";
            return MB_FAILURE;
        }
        // Repeat for target sets
        spectral_sets.empty();
        // Now initialize a source spectral element !
        _spectralSource = new moab::Element::SpectralHex( sem_dims[0] );
        specSou         = true;
    }
    // Repeat for target source
    rval = mbImpl->get_entities_by_type_and_tag( rootTarget, moab::MBENTITYSET, &sem_tag, NULL, 1, spectral_sets );
    if( moab::MB_SUCCESS != rval || spectral_sets.empty() )
        std::cout << "Can't get sem set on target\n";
    else
    {
        moab::EntityHandle firstSemSet = spectral_sets[0];
        rval                           = mbImpl->tag_get_data( sem_tag, &firstSemSet, 1, (void*)sem_dims );
        if( moab::MB_SUCCESS != rval ) return MB_FAILURE;

        if( sem_dims[0] != sem_dims[1] || sem_dims[0] != sem_dims[2] )
        {
            std::cout << " dimensions are different. bail out\n";
            return MB_FAILURE;
        }
        // Repeat for target sets
        spectral_sets.empty();
        // Now initialize a target spectral element !
        _spectralTarget = new moab::Element::SpectralHex( sem_dims[0] );
        specTar         = true;
    }

    _ntot = sem_dims[0] * sem_dims[1] * sem_dims[2];
    rval  = mbImpl->tag_get_handle( "SEM_X", _ntot, moab::MB_TYPE_DOUBLE, _xm1Tag );
    if( moab::MB_SUCCESS != rval )
    {
        std::cout << "Can't get xm1tag \n";
        return MB_FAILURE;
    }
    rval = mbImpl->tag_get_handle( "SEM_Y", _ntot, moab::MB_TYPE_DOUBLE, _ym1Tag );
    if( moab::MB_SUCCESS != rval )
    {
        std::cout << "Can't get ym1tag \n";
        return MB_FAILURE;
    }
    rval = mbImpl->tag_get_handle( "SEM_Z", _ntot, moab::MB_TYPE_DOUBLE, _zm1Tag );
    if( moab::MB_SUCCESS != rval )
    {
        std::cout << "Can't get zm1tag \n";
        return MB_FAILURE;
    }

    return MB_SUCCESS;
}

ErrorCode moab::Coupler::initialize_tree ( )

Initialize the kdtree, locally and across communicator.

Definition at line 83 of file Coupler.cpp.

References allBoxes, moab::BoundBox::bMax, moab::BoundBox::bMin, moab::Range::empty(), ErrorCode, moab::CartVect::get(), moab::Tree::get_bounding_box(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::AdaptiveKDTree::get_info(), moab::ParallelComm::get_part_entities(), moab::CartVect::length(), localRoot, max_dim, MB_SUCCESS, mbImpl, myPc, myRange, myTree, numIts, moab::ProcConfig::proc_comm(), moab::ParallelComm::proc_config(), moab::ProcConfig::proc_rank(), moab::ProcConfig::proc_size(), radius, and spherical.

Referenced by Coupler().

{
    Range local_ents;

    // Get entities on the local part
    ErrorCode result = MB_SUCCESS;
    if( myPc )
    {
        result = myPc->get_part_entities( local_ents, max_dim );
        if( local_ents.empty() )
        {
            max_dim--;
            result = myPc->get_part_entities( local_ents, max_dim );  // go one dimension lower
            // right now, this is used for spherical meshes only
        }
    }
    else
        local_ents = myRange;
    if( MB_SUCCESS != result || local_ents.empty() )
    {
        std::cout << "Problems getting source entities" << std::endl;
        return result;
    }

    // Build the tree for local processor
    int max_per_leaf = 6;
    for( int i = 0; i < numIts; i++ )
    {
        std::ostringstream str;
        str << "PLANE_SET=0;"
            << "MAX_PER_LEAF=" << max_per_leaf << ";";
        if( spherical && !local_ents.empty() )
        {
            // get coordinates of one vertex, and use for radius computation
            EntityHandle elem = local_ents[0];
            const EntityHandle* conn;
            int numn = 0;
            mbImpl->get_connectivity( elem, conn, numn );
            CartVect pos0;
            mbImpl->get_coords( conn, 1, &( pos0[0] ) );
            double radius = pos0.length();
            str << "SPHERICAL=true;RADIUS=" << radius << ";";
        }
        FileOptions opts( str.str().c_str() );
        myTree = new AdaptiveKDTree( mbImpl );
        result = myTree->build_tree( local_ents, &localRoot, &opts );
        if( MB_SUCCESS != result )
        {
            std::cout << "Problems building tree";
            if( numIts != i )
            {
                delete myTree;
                max_per_leaf *= 2;
                std::cout << "; increasing elements/leaf to " << max_per_leaf << std::endl;
            }
            else
            {
                std::cout << "; exiting" << std::endl;
                return result;
            }
        }
        else
            break;  // Get out of tree building
    }

    // Get the bounding box for local tree
    if( myPc )
        allBoxes.resize( 6 * myPc->proc_config().proc_size() );
    else
        allBoxes.resize( 6 );

    unsigned int my_rank = ( myPc ? myPc->proc_config().proc_rank() : 0 );
    BoundBox box;
    result = myTree->get_bounding_box( box, &localRoot );
    if( MB_SUCCESS != result ) return result;
    box.bMin.get( &allBoxes[6 * my_rank] );
    box.bMax.get( &allBoxes[6 * my_rank + 3] );

    // Now communicate to get all boxes
    if( myPc )
    {
        int mpi_err;
#if( MPI_VERSION >= 2 )
        // Use "in place" option
        mpi_err = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &allBoxes[0], 6, MPI_DOUBLE,
                                 myPc->proc_config().proc_comm() );
#else
        {
            std::vector< double > allBoxes_tmp( 6 * myPc->proc_config().proc_size() );
            mpi_err  = MPI_Allgather( &allBoxes[6 * my_rank], 6, MPI_DOUBLE, &allBoxes_tmp[0], 6, MPI_DOUBLE,
                                      myPc->proc_config().proc_comm() );
            allBoxes = allBoxes_tmp;
        }
#endif
        if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
    }

    /*  std::ostringstream blah;
    for(int i = 0; i < allBoxes.size(); i++)
    blah << allBoxes[i] << " ";
    std::cout << blah.str() << "\n";*/

#ifdef VERBOSE
    double min[3] = { 0, 0, 0 }, max[3] = { 0, 0, 0 };
    unsigned int dep;
    myTree->get_info( localRoot, min, max, dep );
    std::cout << "Proc " << my_rank << ": box min/max, tree depth = (" << min[0] << "," << min[1] << "," << min[2]
              << "), (" << max[0] << "," << max[1] << "," << max[2] << "), " << dep << std::endl;
#endif

    return result;
}

ErrorCode moab::Coupler::interp_field	(	EntityHandle	elem,
		CartVect	nat_coord,
		Tag	tag,
		double &	field
	)		`[private]`

Definition at line 1021 of file Coupler.cpp.

References _spectralSource, ErrorCode, moab::Element::Map::evaluate_scalar_field(), moab::Element::SpectralHex::evaluate_scalar_field(), moab::Interface::get_connectivity(), MB_SUCCESS, MBHEX, mbImpl, MBQUAD, MBTET, MBTRI, moab::Interface::tag_get_by_ptr(), moab::Interface::tag_get_data(), and moab::Interface::type_from_handle().

Referenced by interpolate().

{
    if( _spectralSource )
    {
        // Get tag values at the GL points for some field (Tag)
        const double* vx;
        ErrorCode rval = mbImpl->tag_get_by_ptr( tag, &elem, 1, (const void**)&vx );
        if( moab::MB_SUCCESS != rval )
        {
            std::cout << "Can't get field values for the tag \n";
            return MB_FAILURE;
        }
        Element::SpectralHex* spcHex = (Element::SpectralHex*)_spectralSource;
        field                        = spcHex->evaluate_scalar_field( nat_coord, vx );
    }
    else
    {
        double vfields[27];  // Will work for linear hex, quadratic hex or Tets
        moab::Element::Map* elemMap = NULL;
        int num_verts               = 0;
        // Get the EntityType
        // Get the tag values at the vertices
        const EntityHandle* connect;
        int num_connect;
        ErrorCode result = mbImpl->get_connectivity( elem, connect, num_connect );
        if( MB_SUCCESS != result ) return result;
        EntityType etype = mbImpl->type_from_handle( elem );
        if( MBHEX == etype )
        {
            if( 8 == num_connect )
            {
                elemMap   = new moab::Element::LinearHex();
                num_verts = 8;
            }
            else
            { /* (MBHEX == etype && 27 == num_connect) */
                elemMap   = new moab::Element::QuadraticHex();
                num_verts = 27;
            }
        }
        else if( MBTET == etype )
        {
            elemMap   = new moab::Element::LinearTet();
            num_verts = 4;
        }
        else if( MBQUAD == etype )
        {
            elemMap   = new moab::Element::LinearQuad();
            num_verts = 4;
        }
        else if( MBTRI == etype )
        {
            elemMap   = new moab::Element::LinearTri();
            num_verts = 3;
        }
        else
            return MB_FAILURE;

        result = mbImpl->tag_get_data( tag, connect, std::min( num_verts, num_connect ), vfields );
        if( MB_SUCCESS != result )
        {
            delete elemMap;
            return result;
        }

        // Function for the interpolation
        field = 0;

        // Check the number of vertices
        assert( num_connect >= num_verts );

        // Calculate the field
        try
        {
            field = elemMap->evaluate_scalar_field( nat_coord, vfields );
        }
        catch( moab::Element::Map::EvaluationError& )
        {
            delete elemMap;
            return MB_FAILURE;
        }

        delete elemMap;
    }

    return MB_SUCCESS;
}

ErrorCode moab::Coupler::interpolate	(	Coupler::Method	method,
		Tag	tag,
		double *	interp_vals,
		TupleList *	tl = `NULL`,
		bool	normalize = `true`
	)		`[inline]`

Definition at line 565 of file Coupler.hpp.

References moab::TupleList::get_n(), and targetPts.

Referenced by interpolate(), and main().

{
    int num_pts = ( tl ? tl->get_n() : targetPts->get_n() );
    return interpolate( &method, &tag, &num_pts, 1, interp_vals, tl, normalize );
}

ErrorCode moab::Coupler::interpolate	(	Coupler::Method	method,
		const std::string &	tag_name,
		double *	interp_vals,
		TupleList *	tl = `NULL`,
		bool	normalize = `true`
	)

Definition at line 645 of file Coupler.cpp.

References _ntot, _spectralSource, ErrorCode, interpolate(), MB_CHK_SET_ERR, MB_TYPE_DOUBLE, mbImpl, and moab::Interface::tag_get_handle().

{
    Tag tag;
    ErrorCode result;
    if( _spectralSource )
    {
        result = mbImpl->tag_get_handle( interp_tag.c_str(), _ntot, MB_TYPE_DOUBLE, tag );MB_CHK_SET_ERR( result, "Failed to get handle for interpolation tag \"" << interp_tag << "\"" );
    }
    else
    {
        result = mbImpl->tag_get_handle( interp_tag.c_str(), 1, MB_TYPE_DOUBLE, tag );MB_CHK_SET_ERR( result, "Failed to get handle for interpolation tag \"" << interp_tag << "\"" );
    }

    return interpolate( method, tag, interp_vals, tl, normalize );
}

ErrorCode moab::Coupler::interpolate	(	Coupler::Method *	methods,
		const std::string *	tag_names,
		int *	points_per_method,
		int	num_methods,
		double *	interp_vals,
		TupleList *	tl = `NULL`,
		bool	normalize = `true`
	)

ErrorCode moab::Coupler::interpolate	(	Coupler::Method *	methods,
		Tag *	tag_names,
		int *	points_per_method,
		int	num_methods,
		double *	interp_vals,
		TupleList *	tl = `NULL`,
		bool	normalize = `true`
	)

Definition at line 665 of file Coupler.cpp.

References CONSTANT, constant_interp(), moab::ProcConfig::crystal_router(), moab::TupleList::enableWriteAccess(), ErrorCode, moab::TupleList::get_n(), moab::TupleList::inc_n(), moab::TupleList::initialize(), interp_field(), LINEAR_FE, mappedPts, MB_SUCCESS, myPc, moab::ParallelComm::proc_config(), QUADRATIC_FE, SPHERICAL, t, targetPts, moab::TupleList::vi_rd, moab::TupleList::vi_wr, moab::TupleList::vr_rd, moab::TupleList::vr_wr, and moab::TupleList::vul_rd.

{
    // if (!((LINEAR_FE == method) || (CONSTANT == method)))
    // return MB_FAILURE;

    // remote pts first
    TupleList* tl_tmp = ( tl ? tl : targetPts );

    ErrorCode result = MB_SUCCESS;

    unsigned int pts_total = 0;
    for( int i = 0; i < num_methods; i++ )
        pts_total += points_per_method[i];

    // If tl was passed in non-NULL, just have those points, otherwise have targetPts plus
    // locally mapped pts
    if( pts_total != tl_tmp->get_n() ) return MB_FAILURE;

    TupleList tinterp;
    tinterp.initialize( 5, 0, 0, 1, tl_tmp->get_n() );
    int t = 0;
    tinterp.enableWriteAccess();
    for( int i = 0; i < num_methods; i++ )
    {
        for( int j = 0; j < points_per_method[i]; j++ )
        {
            tinterp.vi_wr[5 * t]     = tl_tmp->vi_rd[3 * t];
            tinterp.vi_wr[5 * t + 1] = tl_tmp->vi_rd[3 * t + 1];
            tinterp.vi_wr[5 * t + 2] = tl_tmp->vi_rd[3 * t + 2];
            tinterp.vi_wr[5 * t + 3] = methods[i];
            tinterp.vi_wr[5 * t + 4] = i;
            tinterp.vr_wr[t]         = 0.0;
            tinterp.inc_n();
            t++;
        }
    }

    // Scatter/gather interpolation points
    if( myPc )
    {
        ( myPc->proc_config().crystal_router() )->gs_transfer( 1, tinterp, 0 );

        // Perform interpolation on local source mesh; put results into
        // tinterp.vr_wr[i]

        mappedPts->enableWriteAccess();
        for( unsigned int i = 0; i < tinterp.get_n(); i++ )
        {
            int mindex    = tinterp.vi_rd[5 * i + 2];
            Method method = (Method)tinterp.vi_rd[5 * i + 3];
            Tag tag       = tags[tinterp.vi_rd[5 * i + 4]];

            result = MB_FAILURE;
            if( LINEAR_FE == method || QUADRATIC_FE == method || SPHERICAL == method )
            {
                result = interp_field( mappedPts->vul_rd[mindex], CartVect( mappedPts->vr_wr + 3 * mindex ), tag,
                                       tinterp.vr_wr[i] );
            }
            else if( CONSTANT == method )
            {
                result = constant_interp( mappedPts->vul_rd[mindex], tag, tinterp.vr_wr[i] );
            }

            if( MB_SUCCESS != result ) return result;
        }

        // Scatter/gather interpolation data
        myPc->proc_config().crystal_router()->gs_transfer( 1, tinterp, 0 );
    }

    // Copy the interpolated field as a unit
    for( unsigned int i = 0; i < tinterp.get_n(); i++ )
        interp_vals[tinterp.vi_rd[5 * i + 1]] = tinterp.vr_rd[i];

    // Done
    return MB_SUCCESS;
}

EntityHandle moab::Coupler::local_root ( ) const [inline]

Definition at line 437 of file Coupler.hpp.

References localRoot.

    {
        return localRoot;
    }

ErrorCode moab::Coupler::locate_points	(	double *	xyz,
		unsigned int	num_points,
		double	rel_eps = `0.0`,
		double	abs_eps = `0.0`,
		TupleList *	tl = `NULL`,
		bool	store_local = `true`
	)

Definition at line 319 of file Coupler.cpp.

Referenced by locate_points(), and main().

{
    assert( tl || store_local );

    // target_pts: TL(to_proc, tgt_index, x, y, z): tuples sent to source mesh procs representing
    // pts to be located source_pts: TL(from_proc, tgt_index, src_index): results of source mesh
    // proc point location, ready to send
    //             back to tgt procs; src_index of -1 indicates point not located (arguably not
    //             useful...)
    TupleList target_pts;
    target_pts.initialize( 2, 0, 0, 3, num_points );
    target_pts.enableWriteAccess();

    TupleList source_pts;
    mappedPts = new TupleList( 0, 0, 1, 3, target_pts.get_max() );
    mappedPts->enableWriteAccess();

    source_pts.initialize( 3, 0, 0, 0, target_pts.get_max() );
    source_pts.enableWriteAccess();

    mappedPts->set_n( 0 );
    source_pts.set_n( 0 );
    ErrorCode result;

    unsigned int my_rank = ( myPc ? myPc->proc_config().proc_rank() : 0 );
    bool point_located;

    // For each point, find box(es) containing the point,
    // appending results to tuple_list;
    // keep local points separately, in local_pts, which has pairs
    // of <local_index, mapped_index>, where mapped_index is the index
    // into the mappedPts tuple list
    for( unsigned int i = 0; i < 3 * num_points; i += 3 )
    {

        std::vector< int > procs_to_send_to;
        for( unsigned int j = 0; j < ( myPc ? myPc->proc_config().proc_size() : 0 ); j++ )
        {
            // Test if point is in proc's box
            if( ( allBoxes[6 * j] <= xyz[i] + abs_eps ) && ( xyz[i] <= allBoxes[6 * j + 3] + abs_eps ) &&
                ( allBoxes[6 * j + 1] <= xyz[i + 1] + abs_eps ) && ( xyz[i + 1] <= allBoxes[6 * j + 4] + abs_eps ) &&
                ( allBoxes[6 * j + 2] <= xyz[i + 2] + abs_eps ) && ( xyz[i + 2] <= allBoxes[6 * j + 5] + abs_eps ) )
            {
                // If in this proc's box, will send to proc to test further
                procs_to_send_to.push_back( j );
            }
        }
        if( procs_to_send_to.empty() )
        {
#ifdef VERBOSE
            std::cout << " point index " << i / 3 << ": " << xyz[i] << " " << xyz[i + 1] << " " << xyz[i + 2]
                      << " not found in any box\n";
#endif
            // try to find the closest box, and put it in that box, anyway
            double min_dist = 1.e+20;
            int index       = -1;
            for( unsigned int j = 0; j < ( myPc ? myPc->proc_config().proc_size() : 0 ); j++ )
            {
                BoundBox box( &allBoxes[6 * j] );           // form back the box
                double distance = box.distance( &xyz[i] );  // will compute the distance in 3d, from the box
                if( distance < min_dist )
                {
                    index    = j;
                    min_dist = distance;
                }
            }
            if( index == -1 )
            {
                // need to abort early, nothing we can do
                assert( "cannot locate any box for some points" );
                // need a better exit strategy
            }
#ifdef VERBOSE
            std::cout << " point index " << i / 3 << " added to box for proc j:" << index << "\n";
#endif
            procs_to_send_to.push_back( index );  // will send to just one proc, that has the closest box
        }
        // we finally decided to populate the tuple list for a list of processors
        for( size_t k = 0; k < procs_to_send_to.size(); k++ )
        {
            unsigned int j = procs_to_send_to[k];
            // Check size of tuple list, grow if we're at max
            if( target_pts.get_n() == target_pts.get_max() )
                target_pts.resize( std::max( 10.0, 1.5 * target_pts.get_max() ) );

            target_pts.vi_wr[2 * target_pts.get_n()]     = j;
            target_pts.vi_wr[2 * target_pts.get_n() + 1] = i / 3;

            target_pts.vr_wr[3 * target_pts.get_n()]     = xyz[i];
            target_pts.vr_wr[3 * target_pts.get_n() + 1] = xyz[i + 1];
            target_pts.vr_wr[3 * target_pts.get_n() + 2] = xyz[i + 2];
            target_pts.inc_n();
        }

    }  // end for (unsigned int i = 0; ..

    int num_to_me = 0;
    for( unsigned int i = 0; i < target_pts.get_n(); i++ )
        if( target_pts.vi_rd[2 * i] == (int)my_rank ) num_to_me++;
#ifdef VERBOSE
    printf( "rank: %u local points: %u, nb sent target pts: %u mappedPts: %u num to me: %d \n", my_rank, num_points,
            target_pts.get_n(), mappedPts->get_n(), num_to_me );
#endif
    // Perform scatter/gather, to gather points to source mesh procs
    if( myPc )
    {
        ( myPc->proc_config().crystal_router() )->gs_transfer( 1, target_pts, 0 );

        num_to_me = 0;
        for( unsigned int i = 0; i < target_pts.get_n(); i++ )
        {
            if( target_pts.vi_rd[2 * i] == (int)my_rank ) num_to_me++;
        }
#ifdef VERBOSE
        printf( "rank: %u after first gs nb received_pts: %u; num_from_me = %d\n", my_rank, target_pts.get_n(),
                num_to_me );
#endif
        // After scatter/gather:
        // target_pts.set_n(# points local proc has to map);
        // target_pts.vi_wr[2*i] = proc sending point i
        // target_pts.vi_wr[2*i + 1] = index of point i on sending proc
        // target_pts.vr_wr[3*i..3*i + 2] = xyz of point i
        //
        // Mapping builds the tuple list:
        // source_pts.set_n(target_pts.get_n())
        // source_pts.vi_wr[3*i] = target_pts.vi_wr[2*i] = sending proc
        // source_pts.vi_wr[3*i + 1] = index of point i on sending proc
        // source_pts.vi_wr[3*i + 2] = index of mapped point (-1 if not mapped)
        //
        // Also, mapping builds local tuple_list mappedPts:
        // mappedPts->set_n( # mapped points );
        // mappedPts->vul_wr[i] = local handle of mapped entity
        // mappedPts->vr_wr[3*i..3*i + 2] = natural coordinates in mapped entity

        // Test target points against my elements
        for( unsigned i = 0; i < target_pts.get_n(); i++ )
        {
            result = test_local_box( target_pts.vr_wr + 3 * i, target_pts.vi_rd[2 * i], target_pts.vi_rd[2 * i + 1], i,
                                     point_located, rel_eps, abs_eps, &source_pts );
            if( MB_SUCCESS != result ) return result;
        }

        // No longer need target_pts
        target_pts.reset();
#ifdef VERBOSE
        printf( "rank: %u nb sent source pts: %u, mappedPts now: %u\n", my_rank, source_pts.get_n(),
                mappedPts->get_n() );
#endif
        // Send target points back to target procs
        ( myPc->proc_config().crystal_router() )->gs_transfer( 1, source_pts, 0 );

#ifdef VERBOSE
        printf( "rank: %u nb received source pts: %u\n", my_rank, source_pts.get_n() );
#endif
    }

    // Store proc/index tuples in targetPts, and/or pass back to application;
    // the tuple this gets stored to looks like:
    // tl.set_n(# mapped points);
    // tl.vi_wr[3*i] = remote proc mapping point
    // tl.vi_wr[3*i + 1] = local index of mapped point
    // tl.vi_wr[3*i + 2] = remote index of mapped point
    //
    // Local index is mapped into either myRange, holding the handles of
    // local mapped entities, or myXyz, holding locations of mapped pts

    // Store information about located points
    TupleList* tl_tmp;
    if( !store_local )
        tl_tmp = tl;
    else
    {
        targetPts = new TupleList();
        tl_tmp    = targetPts;
    }

    tl_tmp->initialize( 3, 0, 0, 0, num_points );
    tl_tmp->set_n( num_points );  // Automatically sets tl to write_enabled
    // Initialize so we know afterwards how many pts weren't located
    std::fill( tl_tmp->vi_wr, tl_tmp->vi_wr + 3 * num_points, -1 );

    unsigned int local_pts = 0;
    for( unsigned int i = 0; i < source_pts.get_n(); i++ )
    {
        if( -1 != source_pts.vi_rd[3 * i + 2] )
        {  // Why bother sending message saying "i don't have the point" if it gets discarded?
            int tgt_index = 3 * source_pts.vi_rd[3 * i + 1];
            // Prefer always entities that are local, from the source_pts
            // if a local entity was already found to contain the target point, skip
            // tl_tmp->vi_wr[tgt_index] is -1 initially, but it could already be set with
            // a remote processor
            if( tl_tmp->vi_wr[tgt_index] != (int)my_rank )
            {
                tl_tmp->vi_wr[tgt_index]     = source_pts.vi_rd[3 * i];
                tl_tmp->vi_wr[tgt_index + 1] = source_pts.vi_rd[3 * i + 1];
                tl_tmp->vi_wr[tgt_index + 2] = source_pts.vi_rd[3 * i + 2];
            }
        }
    }

    // Count missing points
    unsigned int missing_pts = 0;
    for( unsigned int i = 0; i < num_points; i++ )
    {
        if( tl_tmp->vi_rd[3 * i + 1] == -1 )
        {
            missing_pts++;
#ifdef VERBOSE
            printf( "missing point at index i:  %d -> %15.10f %15.10f %15.10f\n", i, xyz[3 * i], xyz[3 * i + 1],
                    xyz[3 * i + 2] );
#endif
        }
        else if( tl_tmp->vi_rd[3 * i] == (int)my_rank )
            local_pts++;
    }
#ifdef VERBOSE
    printf( "rank: %u point location: wanted %u got %u locally, %u remote, missing %u\n", my_rank, num_points,
            local_pts, num_points - missing_pts - local_pts, missing_pts );
#endif
    assert( 0 == missing_pts );  // Will likely break on curved geometries

    // No longer need source_pts
    source_pts.reset();

    // Copy into tl if passed in and storing locally
    if( tl && store_local )
    {
        tl->initialize( 3, 0, 0, 0, num_points );
        tl->enableWriteAccess();
        memcpy( tl->vi_wr, tl_tmp->vi_rd, 3 * tl_tmp->get_n() * sizeof( int ) );
        tl->set_n( tl_tmp->get_n() );
        tl->disableWriteAccess();
    }

    tl_tmp->disableWriteAccess();

    // Done
    return MB_SUCCESS;
}

ErrorCode moab::Coupler::locate_points	(	Range &	ents,
		double	rel_eps = `0.0`,
		double	abs_eps = `0.0`,
		TupleList *	tl = `NULL`,
		bool	store_local = `true`
	)

Definition at line 278 of file Coupler.cpp.

References ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), locate_points(), moab::CN::MAX_NODES_PER_ELEMENT, MB_SUCCESS, mbImpl, MBVERTEX, moab::Range::size(), moab::Range::subset_by_type(), moab::subtract(), and targetEnts.

{
    // Get locations
    std::vector< double > locs( 3 * targ_ents.size() );
    Range verts    = targ_ents.subset_by_type( MBVERTEX );
    ErrorCode rval = mbImpl->get_coords( verts, &locs[0] );
    if( MB_SUCCESS != rval ) return rval;
    // Now get other ents; reuse verts
    unsigned int num_verts = verts.size();
    verts                  = subtract( targ_ents, verts );
    // Compute centroids
    std::vector< EntityHandle > dum_conn( CN::MAX_NODES_PER_ELEMENT );
    std::vector< double > dum_pos( CN::MAX_NODES_PER_ELEMENT );
    const EntityHandle* conn;
    int num_conn;
    double* coords = &locs[num_verts];
    // Do this here instead of a function to allow reuse of dum_pos and dum_conn
    for( Range::const_iterator rit = verts.begin(); rit != verts.end(); ++rit )
    {
        rval = mbImpl->get_connectivity( *rit, conn, num_conn, false, &dum_conn );
        if( MB_SUCCESS != rval ) return rval;
        rval = mbImpl->get_coords( conn, num_conn, &dum_pos[0] );
        if( MB_SUCCESS != rval ) return rval;
        coords[0] = coords[1] = coords[2] = 0.0;
        for( int i = 0; i < num_conn; i++ )
        {
            coords[0] += dum_pos[3 * i];
            coords[1] += dum_pos[3 * i + 1];
            coords[2] += dum_pos[3 * i + 2];
        }
        coords[0] /= num_conn;
        coords[1] /= num_conn;
        coords[2] /= num_conn;
        coords += 3;
    }

    if( store_local ) targetEnts = targ_ents;

    return locate_points( &locs[0], targ_ents.size(), rel_eps, abs_eps, tl, store_local );
}

TupleList* moab::Coupler::mapped_pts ( ) const [inline]

Definition at line 461 of file Coupler.hpp.

References mappedPts.

    {
        return mappedPts;
    }

Interface* moab::Coupler::mb_impl ( ) const [inline]

Definition at line 429 of file Coupler.hpp.

References mbImpl.

    {
        return mbImpl;
    }

int moab::Coupler::my_id ( ) const [inline]

Definition at line 453 of file Coupler.hpp.

References myId.

    {
        return myId;
    }

ParallelComm* moab::Coupler::my_pc ( ) const [inline]

Definition at line 445 of file Coupler.hpp.

References myPc.

    {
        return myPc;
    }

const Range& moab::Coupler::my_range ( ) const [inline]

Definition at line 457 of file Coupler.hpp.

References myRange.

    {
        return myRange;
    }

AdaptiveKDTree* moab::Coupler::my_tree ( ) const [inline]

Definition at line 433 of file Coupler.hpp.

References myTree.

    {
        return myTree;
    }

ErrorCode moab::Coupler::nat_param	(	double	xyz[3],
		std::vector< EntityHandle > &	entities,
		std::vector< CartVect > &	nat_coords,
		double	epsilon = `0.0`
	)		`[private]`

Definition at line 749 of file Coupler.cpp.

Referenced by test_local_box().

{
    if( !myTree ) return MB_FAILURE;

    AdaptiveKDTreeIter treeiter;
    ErrorCode result = myTree->get_tree_iterator( localRoot, treeiter );
    if( MB_SUCCESS != result )
    {
        std::cout << "Problems getting iterator" << std::endl;
        return result;
    }

    EntityHandle closest_leaf;
    if( epsilon )
    {
        std::vector< double > dists;
        std::vector< EntityHandle > leaves;

        // Two tolerances
        result = myTree->distance_search( xyz, epsilon, leaves,
                                          /*iter_tol*/ epsilon,
                                          /*inside_tol*/ 10 * epsilon, &dists, NULL, &localRoot );
        if( leaves.empty() )
            // Not found returns success here, with empty list, just like case with no epsilon
            return MB_SUCCESS;
        // Get closest leaf
        double min_dist                           = *dists.begin();
        closest_leaf                              = *leaves.begin();
        std::vector< EntityHandle >::iterator vit = leaves.begin() + 1;
        std::vector< double >::iterator dit       = dists.begin() + 1;
        for( ; vit != leaves.end() && min_dist; ++vit, ++dit )
        {
            if( *dit < min_dist )
            {
                min_dist     = *dit;
                closest_leaf = *vit;
            }
        }
    }
    else
    {
        result = myTree->point_search( xyz, treeiter, 1.0e-10, 1.0e-6, NULL, &localRoot );
        if( MB_ENTITY_NOT_FOUND == result )  // Point is outside of myTree's bounding box
            return MB_SUCCESS;
        else if( MB_SUCCESS != result )
        {
            std::cout << "Problems getting leaf \n";
            return result;
        }
        closest_leaf = treeiter.handle();
    }

    // Find natural coordinates of point in element(s) in that leaf
    CartVect tmp_nat_coords;
    Range range_leaf;
    result = mbImpl->get_entities_by_dimension( closest_leaf, max_dim, range_leaf, false );
    if( MB_SUCCESS != result ) std::cout << "Problem getting leaf in a range" << std::endl;

    // Loop over the range_leaf
    for( Range::iterator iter = range_leaf.begin(); iter != range_leaf.end(); ++iter )
    {
        // Test to find out in which entity the point is
        // Get the EntityType and create the appropriate Element::Map subtype
        // If spectral, do not need coordinates, just the GL points
        EntityType etype = mbImpl->type_from_handle( *iter );
        if( NULL != this->_spectralSource && MBHEX == etype )
        {
            EntityHandle eh = *iter;
            const double* xval;
            const double* yval;
            const double* zval;
            ErrorCode rval = mbImpl->tag_get_by_ptr( _xm1Tag, &eh, 1, (const void**)&xval );
            if( moab::MB_SUCCESS != rval )
            {
                std::cout << "Can't get xm1 values \n";
                return MB_FAILURE;
            }
            rval = mbImpl->tag_get_by_ptr( _ym1Tag, &eh, 1, (const void**)&yval );
            if( moab::MB_SUCCESS != rval )
            {
                std::cout << "Can't get ym1 values \n";
                return MB_FAILURE;
            }
            rval = mbImpl->tag_get_by_ptr( _zm1Tag, &eh, 1, (const void**)&zval );
            if( moab::MB_SUCCESS != rval )
            {
                std::cout << "Can't get zm1 values \n";
                return MB_FAILURE;
            }
            Element::SpectralHex* spcHex = (Element::SpectralHex*)_spectralSource;

            spcHex->set_gl_points( (double*)xval, (double*)yval, (double*)zval );
            try
            {
                tmp_nat_coords = spcHex->ievaluate( CartVect( xyz ), epsilon );  // introduce
                bool inside    = spcHex->inside_nat_space( CartVect( tmp_nat_coords ), epsilon );
                if( !inside )
                {
#ifdef VERBOSE
                    std::cout << "point " << xyz[0] << " " << xyz[1] << " " << xyz[2]
                              << " is not converging inside hex " << mbImpl->id_from_handle( eh ) << "\n";
#endif
                    continue;  // It is possible that the point is outside, so it will not converge
                }
            }
            catch( Element::Map::EvaluationError& )
            {
                continue;
            }
        }
        else
        {
            const EntityHandle* connect;
            int num_connect;

            // Get connectivity
            result = mbImpl->get_connectivity( *iter, connect, num_connect, true );
            if( MB_SUCCESS != result ) return result;

            // Get coordinates of the vertices
            std::vector< CartVect > coords_vert( num_connect );
            result = mbImpl->get_coords( connect, num_connect, &( coords_vert[0][0] ) );
            if( MB_SUCCESS != result )
            {
                std::cout << "Problems getting coordinates of vertices\n";
                return result;
            }
            CartVect pos( xyz );
            if( MBHEX == etype )
            {
                if( 8 == num_connect )
                {
                    Element::LinearHex hexmap( coords_vert );
                    if( !hexmap.inside_box( pos, epsilon ) ) continue;
                    try
                    {
                        tmp_nat_coords = hexmap.ievaluate( pos, epsilon );
                        bool inside    = hexmap.inside_nat_space( tmp_nat_coords, epsilon );
                        if( !inside ) continue;
                    }
                    catch( Element::Map::EvaluationError& )
                    {
                        continue;
                    }
                }
                else if( 27 == num_connect )
                {
                    Element::QuadraticHex hexmap( coords_vert );
                    if( !hexmap.inside_box( pos, epsilon ) ) continue;
                    try
                    {
                        tmp_nat_coords = hexmap.ievaluate( pos, epsilon );
                        bool inside    = hexmap.inside_nat_space( tmp_nat_coords, epsilon );
                        if( !inside ) continue;
                    }
                    catch( Element::Map::EvaluationError& )
                    {
                        continue;
                    }
                }
                else  // TODO this case not treated yet, no interpolation
                    continue;
            }
            else if( MBTET == etype )
            {
                Element::LinearTet tetmap( coords_vert );
                // This is just a linear solve; unless degenerate, will not except
                tmp_nat_coords = tetmap.ievaluate( pos );
                bool inside    = tetmap.inside_nat_space( tmp_nat_coords, epsilon );
                if( !inside ) continue;
            }
            else if( MBQUAD == etype && spherical )
            {
                Element::SphericalQuad sphermap( coords_vert );
                /* skip box test, because it can filter out good elements with high curvature
                 * if (!sphermap.inside_box(pos, epsilon))
                  continue;*/
                try
                {
                    tmp_nat_coords = sphermap.ievaluate( pos, epsilon );
                    bool inside    = sphermap.inside_nat_space( tmp_nat_coords, epsilon );
                    if( !inside ) continue;
                }
                catch( Element::Map::EvaluationError& )
                {
                    continue;
                }
            }
            else if( MBTRI == etype && spherical )
            {
                Element::SphericalTri sphermap( coords_vert );
                /* skip box test, because it can filter out good elements with high curvature
                 * if (!sphermap.inside_box(pos, epsilon))
                    continue;*/
                try
                {
                    tmp_nat_coords = sphermap.ievaluate( pos, epsilon );
                    bool inside    = sphermap.inside_nat_space( tmp_nat_coords, epsilon );
                    if( !inside ) continue;
                }
                catch( Element::Map::EvaluationError& )
                {
                    continue;
                }
            }

            else if( MBQUAD == etype )
            {
                Element::LinearQuad quadmap( coords_vert );
                if( !quadmap.inside_box( pos, epsilon ) ) continue;
                try
                {
                    tmp_nat_coords = quadmap.ievaluate( pos, epsilon );
                    bool inside    = quadmap.inside_nat_space( tmp_nat_coords, epsilon );
                    if( !inside ) continue;
                }
                catch( Element::Map::EvaluationError& )
                {
                    continue;
                }
                if( !quadmap.inside_nat_space( tmp_nat_coords, epsilon ) ) continue;
            }
            /*
            else if (etype == MBTRI){
              Element::LinearTri trimap(coords_vert);
              if (!trimap.inside_box( pos, epsilon))
                continue;
              try {
                tmp_nat_coords = trimap.ievaluate(pos, epsilon);
                bool inside = trimap.inside_nat_space(tmp_nat_coords, epsilon);
                if (!inside) continue;
              }
              catch (Element::Map::EvaluationError) {
                continue;
              }
              if (!trimap.inside_nat_space(tmp_nat_coords, epsilon))
                continue;
            }
            */
            else if( etype == MBEDGE )
            {
                Element::LinearEdge edgemap( coords_vert );
                try
                {
                    tmp_nat_coords = edgemap.ievaluate( CartVect( xyz ), epsilon );
                }
                catch( Element::Map::EvaluationError )
                {
                    continue;
                }
                if( !edgemap.inside_nat_space( tmp_nat_coords, epsilon ) ) continue;
            }
            else
            {
                std::cout << "Entity not Hex/Tet/Quad/Tri/Edge. Please verify." << std::endl;
                continue;
            }
        }
        // If we get here then we've found the coordinates.
        // Save them and the entity and return success.
        entities.push_back( *iter );
        nat_coords.push_back( tmp_nat_coords );
        return MB_SUCCESS;
    }

    // Didn't find any elements containing the point
    return MB_SUCCESS;
}

ErrorCode moab::Coupler::normalize_mesh	(	EntityHandle	root_set,
		const char *	norm_tag,
		Coupler::IntegType	integ_type,
		int	num_integ_pts
	)

Definition at line 1125 of file Coupler.cpp.

References do_normalization(), entities, ErrorCode, ERRORR, moab::Interface::get_entities_by_handle(), and mbImpl.

{
    ErrorCode err;

    // SLAVE START ****************************************************************
    // Search for entities based on tag_handles and tag_values
    std::vector< std::vector< EntityHandle > > entity_sets;
    std::vector< std::vector< EntityHandle > > entity_groups;

    // put the root_set into entity_sets
    std::vector< EntityHandle > ent_set;
    ent_set.push_back( root_set );
    entity_sets.push_back( ent_set );

    // get all entities from root_set and put into entity_groups
    std::vector< EntityHandle > entities;
    err = mbImpl->get_entities_by_handle( root_set, entities, true );ERRORR( "Failed to get entities in root_set.", err );

    entity_groups.push_back( entities );

    // Call do_normalization() to continue common normalization processing
    err = do_normalization( norm_tag, entity_sets, entity_groups, integ_type, num_integ_pts );ERRORR( "Failure in do_normalization().", err );
    // SLAVE END   ****************************************************************

    return err;
}

ErrorCode moab::Coupler::normalize_subset	(	EntityHandle	root_set,
		const char *	norm_tag,
		const char **	tag_names,
		int	num_tags,
		const char **	tag_values,
		Coupler::IntegType	integ_type,
		int	num_integ_pts
	)

Definition at line 1156 of file Coupler.cpp.

References ErrorCode, ERRORR, MB_TAG_ANY, MB_TYPE_DOUBLE, mbImpl, t, and moab::Interface::tag_get_handle().

Referenced by main().

{
    moab::ErrorCode err;
    std::vector< Tag > tag_handles;

    // Lookup tag handles from tag names
    for( int t = 0; t < num_tags; t++ )
    {
        // get tag handle & size
        Tag th;
        err = mbImpl->tag_get_handle( tag_names[t], 1, moab::MB_TYPE_DOUBLE, th, moab::MB_TAG_ANY );ERRORR( "Failed to get tag handle.", err );
        tag_handles.push_back( th );
    }

    return normalize_subset( root_set, norm_tag, &tag_handles[0], num_tags, tag_values, integ_type, num_integ_pts );
}

ErrorCode moab::Coupler::normalize_subset	(	EntityHandle	root_set,
		const char *	norm_tag,
		Tag *	tag_handles,
		int	num_tags,
		const char **	tag_values,
		Coupler::IntegType	integ_type,
		int	num_integ_pts
	)

Definition at line 1179 of file Coupler.cpp.

References do_normalization(), ErrorCode, ERRORR, and get_matching_entities().

{
    ErrorCode err;

    // SLAVE START ****************************************************************
    // Search for entities based on tag_handles and tag_values
    std::vector< std::vector< EntityHandle > > entity_sets;
    std::vector< std::vector< EntityHandle > > entity_groups;

    err = get_matching_entities( root_set, tag_handles, tag_values, num_tags, &entity_sets, &entity_groups );ERRORR( "Failed to get matching entities.", err );

    // Call do_normalization() to continue common normalization processing
    err = do_normalization( norm_tag, entity_sets, entity_groups, integ_type, num_integ_pts );ERRORR( "Failure in do_normalization().", err );
    // SLAVE END   ****************************************************************

    return err;
}

int moab::Coupler::num_its ( ) const [inline]

Definition at line 465 of file Coupler.hpp.

References numIts.

    {
        return numIts;
    }

void moab::Coupler::set_spherical ( bool arg1 = true ) [inline]

Definition at line 470 of file Coupler.hpp.

References spherical.

    {
        spherical = arg1;
    }

const Range& moab::Coupler::target_ents ( ) const [inline]

Definition at line 449 of file Coupler.hpp.

References targetEnts.

    {
        return targetEnts;
    }

ErrorCode moab::Coupler::test_local_box	(	double *	xyz,
		int	from_proc,
		int	remote_index,
		int	index,
		bool &	point_located,
		double	rel_eps = `0.0`,
		double	abs_eps = `0.0`,
		TupleList *	tl = `NULL`
	)		`[private]`

Definition at line 564 of file Coupler.cpp.

References moab::BoundBox::diagonal_length(), moab::TupleList::disableWriteAccess(), moab::TupleList::enableWriteAccess(), entities, ErrorCode, moab::Tree::get_bounding_box(), moab::TupleList::get_max(), moab::TupleList::get_n(), moab::TupleList::get_writeEnabled(), moab::TupleList::inc_n(), localRoot, mappedPts, MB_SUCCESS, myTree, nat_param(), moab::TupleList::resize(), moab::TupleList::vi_wr, moab::TupleList::vr_wr, and moab::TupleList::vul_wr.

Referenced by locate_points().

{
    std::vector< EntityHandle > entities;
    std::vector< CartVect > nat_coords;
    bool canWrite = false;
    if( tl )
    {
        canWrite = tl->get_writeEnabled();
        if( !canWrite )
        {
            tl->enableWriteAccess();
            canWrite = true;
        }
    }

    if( rel_eps && !abs_eps )
    {
        // Relative epsilon given, translate to absolute epsilon using box dimensions
        BoundBox box;
        myTree->get_bounding_box( box, &localRoot );
        abs_eps = rel_eps * box.diagonal_length();
    }

    ErrorCode result = nat_param( xyz, entities, nat_coords, abs_eps );
    if( MB_SUCCESS != result ) return result;

    // If we didn't find any ents and we're looking locally, nothing more to do
    if( entities.empty() )
    {
        if( tl->get_n() == tl->get_max() ) tl->resize( std::max( 10.0, 1.5 * tl->get_max() ) );

        tl->vi_wr[3 * tl->get_n()]     = from_proc;
        tl->vi_wr[3 * tl->get_n() + 1] = remote_index;
        tl->vi_wr[3 * tl->get_n() + 2] = -1;
        tl->inc_n();

        point_located = false;
        return MB_SUCCESS;
    }

    // Grow if we know we'll exceed size
    if( mappedPts->get_n() + entities.size() >= mappedPts->get_max() )
        mappedPts->resize( std::max( 10.0, 1.5 * mappedPts->get_max() ) );

    std::vector< EntityHandle >::iterator eit = entities.begin();
    std::vector< CartVect >::iterator ncit    = nat_coords.begin();

    mappedPts->enableWriteAccess();
    for( ; eit != entities.end(); ++eit, ++ncit )
    {
        // Store in tuple mappedPts
        mappedPts->vr_wr[3 * mappedPts->get_n()]     = ( *ncit )[0];
        mappedPts->vr_wr[3 * mappedPts->get_n() + 1] = ( *ncit )[1];
        mappedPts->vr_wr[3 * mappedPts->get_n() + 2] = ( *ncit )[2];
        mappedPts->vul_wr[mappedPts->get_n()]        = *eit;
        mappedPts->inc_n();

        // Also store local point, mapped point indices
        if( tl->get_n() == tl->get_max() ) tl->resize( std::max( 10.0, 1.5 * tl->get_max() ) );

        // Store in tuple source_pts
        tl->vi_wr[3 * tl->get_n()]     = from_proc;
        tl->vi_wr[3 * tl->get_n() + 1] = remote_index;
        tl->vi_wr[3 * tl->get_n() + 2] = mappedPts->get_n() - 1;
        tl->inc_n();
    }

    point_located = true;

    if( tl && !canWrite ) tl->disableWriteAccess();

    return MB_SUCCESS;
}