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. More...

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

◆ IntegType

enum moab::Coupler::IntegType

Enumerator
VOLUME

Definition at line 55 of file Coupler.hpp.

  {
  VOLUME
  };

◆ Method

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

◆ Coupler()

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.

  : 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;
 }

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

◆ ~Coupler()

moab::Coupler::~Coupler ( )

virtual

Definition at line 73 of file Coupler.cpp.

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

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

Member Function Documentation

◆ all_boxes()

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

inline

Definition at line 441 of file Coupler.hpp.

  {
  return allBoxes;
  }

References allBoxes.

◆ apply_group_norm_factor()

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 1785 of file Coupler.cpp.

 {
  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;
 }

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().

◆ consolidate_tuples()

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

Definition at line 1561 of file Coupler.cpp.

 {
  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;
 }

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().

◆ constant_interp()

ErrorCode moab::Coupler::constant_interp	(	EntityHandle	elem,
		Tag	tag,
		double &	field
	)

private

Definition at line 1109 of file Coupler.cpp.

 {
  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;
 }

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

Referenced by interpolate().

◆ create_tuples() [1/2]

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

Definition at line 1504 of file Coupler.cpp.

 {
  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 );
 }

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

Referenced by get_matching_entities(), and main().

◆ create_tuples() [2/2]

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

Definition at line 1526 of file Coupler.cpp.

 {
  // 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;
 }

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.

◆ do_normalization()

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 1201 of file Coupler.cpp.

 {
  // 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;
 }

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().

◆ get_gl_points_on_elements()

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

Definition at line 1928 of file Coupler.cpp.

 {
  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;
 }

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().

◆ get_group_integ_vals()

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 1650 of file Coupler.cpp.

 {
  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;
 }

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().

◆ get_matching_entities() [1/2]

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 1291 of file Coupler.cpp.

 {
  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 );
 }

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

Referenced by main(), and normalize_subset().

◆ get_matching_entities() [2/2]

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 1313 of file Coupler.cpp.

 {
  // 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;
 }

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.

◆ initialize_spectral_elements()

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

Definition at line 196 of file Coupler.cpp.

 {
  /*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;
 }

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().

◆ initialize_tree()

ErrorCode moab::Coupler::initialize_tree ( )

Initialize the kdtree, locally and across communicator.

Definition at line 83 of file Coupler.cpp.

 {
  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;
 }

References allBoxes, moab::BoundBox::bMax, moab::BoundBox::bMin, moab::AdaptiveKDTree::build_tree(), 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(), and spherical.

Referenced by Coupler().

◆ interp_field()

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

private

Definition at line 1019 of file Coupler.cpp.

 {
  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;
 }

References _spectralSource, ErrorCode, moab::Element::SpectralHex::evaluate_scalar_field(), moab::Element::Map::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().

◆ interpolate() [1/4]

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`
	)

◆ interpolate() [2/4]

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 663 of file Coupler.cpp.

 {
  // 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;
 }

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, targetPts, moab::TupleList::vi_rd, moab::TupleList::vi_wr, moab::TupleList::vr_rd, moab::TupleList::vr_wr, and moab::TupleList::vul_rd.

◆ interpolate() [3/4]

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

Definition at line 643 of file Coupler.cpp.

 {
  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 );
 }

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

◆ interpolate() [4/4]

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.

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

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

Referenced by interpolate(), and main().

◆ local_root()

EntityHandle moab::Coupler::local_root ( ) const

inline

Definition at line 437 of file Coupler.hpp.

  {
  return localRoot;
  }

References localRoot.

◆ locate_points() [1/2]

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.

 {
  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; ..
  
 #ifdef VERBOSE
  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++;
  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 );
  
 #ifdef VERBOSE
  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++;
  }
  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++;
  }
  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 );
  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;
 }

Referenced by locate_points(), and main().

◆ locate_points() [2/2]

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.

 {
  // 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 );
 }

References moab::Range::begin(), moab::Range::end(), 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.

◆ mapped_pts()

TupleList* moab::Coupler::mapped_pts ( ) const

inline

Definition at line 461 of file Coupler.hpp.

  {
  return mappedPts;
  }

References mappedPts.

◆ mb_impl()

Interface* moab::Coupler::mb_impl ( ) const

inline

Definition at line 429 of file Coupler.hpp.

  {
  return mbImpl;
  }

References mbImpl.

◆ my_id()

int moab::Coupler::my_id ( ) const

inline

Definition at line 453 of file Coupler.hpp.

  {
  return myId;
  }

References myId.

◆ my_pc()

ParallelComm* moab::Coupler::my_pc ( ) const

inline

Definition at line 445 of file Coupler.hpp.

  {
  return myPc;
  }

References myPc.

◆ my_range()

const Range& moab::Coupler::my_range ( ) const

inline

Definition at line 457 of file Coupler.hpp.

  {
  return myRange;
  }

References myRange.

◆ my_tree()

AdaptiveKDTree* moab::Coupler::my_tree ( ) const

inline

Definition at line 433 of file Coupler.hpp.

  {
  return myTree;
  }

References myTree.

◆ nat_param()

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 747 of file Coupler.cpp.

 {
  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;
 }

Referenced by test_local_box().

◆ normalize_mesh()

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

Definition at line 1123 of file Coupler.cpp.

 {
  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;
 }

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

◆ normalize_subset() [1/2]

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 1154 of file Coupler.cpp.

 {
  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 );
 }

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

Referenced by main().

◆ normalize_subset() [2/2]

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 1177 of file Coupler.cpp.

 {
  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;
 }

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

◆ num_its()

int moab::Coupler::num_its ( ) const

inline

Definition at line 465 of file Coupler.hpp.

  {
  return numIts;
  }

References numIts.

◆ set_spherical()

void moab::Coupler::set_spherical ( bool arg1 = true )

inline

Definition at line 470 of file Coupler.hpp.

  {
  spherical = arg1;
  }

References spherical.

◆ target_ents()

const Range& moab::Coupler::target_ents ( ) const

inline

Definition at line 449 of file Coupler.hpp.

  {
  return targetEnts;
  }

References targetEnts.

◆ test_local_box()

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 562 of file Coupler.cpp.

 {
  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;
 }

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().

Member Data Documentation

◆ _ntot

int moab::Coupler::_ntot

private

Definition at line 559 of file Coupler.hpp.

Referenced by get_gl_points_on_elements(), initialize_spectral_elements(), and interpolate().

◆ _spectralSource

void* moab::Coupler::_spectralSource

private

Definition at line 556 of file Coupler.hpp.

Referenced by Coupler(), initialize_spectral_elements(), interp_field(), interpolate(), nat_param(), and ~Coupler().

◆ _spectralTarget

void* moab::Coupler::_spectralTarget

private

Definition at line 557 of file Coupler.hpp.

Referenced by Coupler(), initialize_spectral_elements(), and ~Coupler().

◆ _xm1Tag

moab::Tag moab::Coupler::_xm1Tag

private

Definition at line 558 of file Coupler.hpp.

Referenced by get_gl_points_on_elements(), initialize_spectral_elements(), and nat_param().

◆ _ym1Tag

moab::Tag moab::Coupler::_ym1Tag

private

Definition at line 558 of file Coupler.hpp.

Referenced by get_gl_points_on_elements(), initialize_spectral_elements(), and nat_param().

◆ _zm1Tag

moab::Tag moab::Coupler::_zm1Tag

private

Definition at line 558 of file Coupler.hpp.

Referenced by get_gl_points_on_elements(), initialize_spectral_elements(), and nat_param().

◆ allBoxes

std::vector< double > moab::Coupler::allBoxes

private

Definition at line 510 of file Coupler.hpp.

Referenced by all_boxes(), initialize_tree(), and locate_points().

◆ localRoot

EntityHandle moab::Coupler::localRoot

private

Definition at line 506 of file Coupler.hpp.

Referenced by initialize_tree(), local_root(), nat_param(), and test_local_box().

◆ mappedPts

TupleList* moab::Coupler::mappedPts

private

Definition at line 534 of file Coupler.hpp.

Referenced by Coupler(), interpolate(), locate_points(), mapped_pts(), test_local_box(), and ~Coupler().

◆ max_dim

int moab::Coupler::max_dim

private

Definition at line 551 of file Coupler.hpp.

Referenced by initialize_tree(), and nat_param().

◆ mbImpl

Interface* moab::Coupler::mbImpl

private

Definition at line 498 of file Coupler.hpp.

Referenced by apply_group_norm_factor(), constant_interp(), create_tuples(), get_gl_points_on_elements(), get_group_integ_vals(), get_matching_entities(), initialize_spectral_elements(), initialize_tree(), interp_field(), interpolate(), locate_points(), mb_impl(), nat_param(), normalize_mesh(), and normalize_subset().

◆ myId

int moab::Coupler::myId

private

Definition at line 518 of file Coupler.hpp.

Referenced by my_id().

◆ myPc

ParallelComm* moab::Coupler::myPc

private

Definition at line 514 of file Coupler.hpp.

Referenced by do_normalization(), get_matching_entities(), initialize_tree(), interpolate(), locate_points(), and my_pc().

◆ myRange

Range moab::Coupler::myRange

private

Definition at line 522 of file Coupler.hpp.

Referenced by Coupler(), initialize_tree(), and my_range().

◆ myTree

AdaptiveKDTree* moab::Coupler::myTree

private

Definition at line 502 of file Coupler.hpp.

Referenced by Coupler(), initialize_tree(), my_tree(), nat_param(), test_local_box(), and ~Coupler().

◆ numIts

int moab::Coupler::numIts

private

Definition at line 548 of file Coupler.hpp.

Referenced by initialize_tree(), and num_its().

◆ spherical

bool moab::Coupler::spherical

private

Definition at line 562 of file Coupler.hpp.

Referenced by initialize_tree(), nat_param(), and set_spherical().

◆ targetEnts

Range moab::Coupler::targetEnts

private

Definition at line 526 of file Coupler.hpp.

Referenced by locate_points(), and target_ents().

◆ targetPts

TupleList* moab::Coupler::targetPts

private

Definition at line 543 of file Coupler.hpp.

Referenced by Coupler(), interpolate(), locate_points(), and ~Coupler().

The documentation for this class was generated from the following files:

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Enumeration Documentation

◆ IntegType

◆ Method

Constructor & Destructor Documentation

◆ Coupler()

◆ ~Coupler()

Member Function Documentation

◆ all_boxes()

◆ apply_group_norm_factor()

◆ consolidate_tuples()

◆ constant_interp()

◆ create_tuples() [1/2]

◆ create_tuples() [2/2]

◆ do_normalization()

◆ get_gl_points_on_elements()

◆ get_group_integ_vals()

◆ get_matching_entities() [1/2]

◆ get_matching_entities() [2/2]

◆ initialize_spectral_elements()

◆ initialize_tree()

◆ interp_field()

◆ interpolate() [1/4]

◆ interpolate() [2/4]

◆ interpolate() [3/4]

◆ interpolate() [4/4]

◆ local_root()

◆ locate_points() [1/2]

◆ locate_points() [2/2]

◆ mapped_pts()

◆ mb_impl()

◆ my_id()

◆ my_pc()

◆ my_range()

◆ my_tree()

◆ nat_param()

◆ normalize_mesh()

◆ normalize_subset() [1/2]

◆ normalize_subset() [2/2]

◆ num_its()

◆ set_spherical()

◆ target_ents()

◆ test_local_box()

Member Data Documentation

◆ _ntot

◆ _spectralSource

◆ _spectralTarget

◆ _xm1Tag

◆ _ym1Tag

◆ _zm1Tag

◆ allBoxes

◆ localRoot

◆ mappedPts

◆ max_dim

◆ mbImpl

◆ myId

◆ myPc

◆ myRange

◆ myTree

◆ numIts

◆ spherical

◆ targetEnts

◆ targetPts