Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
moab::Coupler Class Reference

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)
 
Interfacemb_impl () const
 
AdaptiveKDTreemy_tree () const
 
EntityHandle local_root () const
 
const std::vector< double > & all_boxes () const
 
ParallelCommmy_pc () const
 
const Rangetarget_ents () const
 
int my_id () const
 
const Rangemy_range () const
 
TupleListmapped_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

InterfacembImpl
 
AdaptiveKDTreemyTree
 
EntityHandle localRoot
 
std::vector< double > allBoxes
 
ParallelCommmyPc
 
int myId
 
Range myRange
 
Range targetEnts
 
TupleListmappedPts
 
TupleListtargetPts
 
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

Enumerator
VOLUME 

Definition at line 55 of file Coupler.hpp.

56  {
57  VOLUME
58  };

◆ Method

Enumerator
CONSTANT 
LINEAR_FE 
QUADRATIC_FE 
SPECTRAL 
SPHERICAL 

Definition at line 46 of file Coupler.hpp.

47  {
48  CONSTANT,
49  LINEAR_FE,
51  SPECTRAL,
52  SPHERICAL
53  };

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.

53  : mbImpl( impl ), myPc( pc ), myId( coupler_id ), numIts( 3 ), max_dim( max_ent_dim ), _ntot( 0 ),
54  spherical( false )
55 {
56  assert( NULL != impl && ( pc || !local_elems.empty() ) );
57 
58  // Keep track of the local points, at least for now
59  myRange = local_elems;
60  myTree = NULL;
61 
62  // Now initialize the tree
63  if( init_tree ) initialize_tree();
64 
65  // Initialize tuple lists to indicate not initialized
66  mappedPts = NULL;
67  targetPts = NULL;
69 }

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.

74 {
75  // This will clear the cache
78  delete myTree;
79  delete targetPts;
80  delete mappedPts;
81 }

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.

442  {
443  return allBoxes;
444  }

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

1791 {
1792  ErrorCode err;
1793 
1794  // Construct the new tag for the normalization factor from the norm_tag name
1795  // and "_normf".
1796  int norm_tag_len = strlen( norm_tag );
1797  const char* normf_appd = "_normf";
1798  int normf_appd_len = strlen( normf_appd );
1799 
1800  char* normf_tag = (char*)malloc( norm_tag_len + normf_appd_len + 1 );
1801  char* tmp_ptr = normf_tag;
1802 
1803  memcpy( tmp_ptr, norm_tag, norm_tag_len );
1804  tmp_ptr += norm_tag_len;
1805  memcpy( tmp_ptr, normf_appd, normf_appd_len );
1806  tmp_ptr += normf_appd_len;
1807  *tmp_ptr = '\0';
1808 
1809  Tag normf_hdl;
1810  // Check to see if the tag exists. If not then create it and get the handle.
1811  err = mbImpl->tag_get_handle( normf_tag, 1, moab::MB_TYPE_DOUBLE, normf_hdl,
1812  moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT );ERRORR( "Failed to create normalization factor tag.", err );
1813  if( normf_hdl == NULL )
1814  {
1815  std::string msg( "Failed to create normalization factor tag named '" );
1816  msg += std::string( normf_tag ) + std::string( "'" );ERRORR( msg.c_str(), MB_FAILURE );
1817  }
1818  free( normf_tag );
1819 
1820  std::vector< std::vector< EntityHandle > >::iterator iter_i;
1821  std::vector< EntityHandle >::iterator iter_j;
1822  std::vector< double >::iterator iter_f;
1823  double grp_norm_factor = 0.0;
1824 
1825  // Loop over the entity sets
1826  for( iter_i = entity_sets.begin(), iter_f = norm_factors.begin();
1827  ( iter_i != entity_sets.end() ) && ( iter_f != norm_factors.end() ); ++iter_i, ++iter_f )
1828  {
1829  grp_norm_factor = *iter_f;
1830 
1831  // Loop over the all the entity sets in iter_i and set the
1832  // new normf_tag with the norm factor value on each
1833  for( iter_j = ( *iter_i ).begin(); iter_j != ( *iter_i ).end(); ++iter_j )
1834  {
1835  EntityHandle entset = *iter_j;
1836 
1837  std::cout << "Coupler: applying normalization for entity set=" << entset
1838  << ", normalization_factor=" << grp_norm_factor << std::endl;
1839 
1840  err = mbImpl->tag_set_data( normf_hdl, &entset, 1, &grp_norm_factor );ERRORR( "Failed to set normalization factor on entity set.", err );
1841  }
1842  }
1843 
1844  return MB_SUCCESS;
1845 }

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

1564 {
1565  int total_rcv_tuples = 0;
1566  int offset = 0, copysz = 0;
1567  unsigned num_tags = 0;
1568 
1569  uint ml, mul, mr;
1570  uint* mi = (uint*)malloc( sizeof( uint ) * num_tuples );
1571 
1572  for( unsigned int i = 0; i < num_tuples; i++ )
1573  {
1574  all_tuples[i]->getTupleSize( mi[i], ml, mul, mr );
1575  }
1576 
1577  for( unsigned int i = 0; i < num_tuples; i++ )
1578  {
1579  if( all_tuples[i] != NULL )
1580  {
1581  total_rcv_tuples += all_tuples[i]->get_n();
1582  num_tags = mi[i];
1583  }
1584  }
1585  const unsigned int_size = sizeof( sint );
1586  const unsigned int_width = num_tags * int_size;
1587 
1588  // Get the total size of all of the tuple_lists in all_tuples.
1589  for( unsigned int i = 0; i < num_tuples; i++ )
1590  {
1591  if( all_tuples[i] != NULL ) total_rcv_tuples += all_tuples[i]->get_n();
1592  }
1593 
1594  // Copy the tuple_lists into a single tuple_list.
1595  TupleList* all_tuples_list = new TupleList( num_tags, 0, 0, 0, total_rcv_tuples );
1596  all_tuples_list->enableWriteAccess();
1597  // all_tuples_list->initialize(num_tags, 0, 0, 0, total_rcv_tuples);
1598  for( unsigned int i = 0; i < num_tuples; i++ )
1599  {
1600  if( all_tuples[i] != NULL )
1601  {
1602  copysz = all_tuples[i]->get_n() * int_width;
1603  memcpy( all_tuples_list->vi_wr + offset, all_tuples[i]->vi_rd, copysz );
1604  offset = offset + ( all_tuples[i]->get_n() * mi[i] );
1605  all_tuples_list->set_n( all_tuples_list->get_n() + all_tuples[i]->get_n() );
1606  }
1607  }
1608 
1609  // Sort the new tuple_list. Use a radix type sort, starting with the last (or least significant)
1610  // tag column in the vi array and working towards the first (or most significant) tag column.
1611  TupleList::buffer sort_buffer;
1612  sort_buffer.buffer_init( 2 * total_rcv_tuples * int_width );
1613  for( int i = num_tags - 1; i >= 0; i-- )
1614  {
1615  all_tuples_list->sort( i, &sort_buffer );
1616  }
1617 
1618  // Cycle through the sorted list eliminating duplicates.
1619  // Keep counters to the current end of the tuple_list (w/out dups) and the last tuple examined.
1620  unsigned int end_idx = 0, last_idx = 1;
1621  while( last_idx < all_tuples_list->get_n() )
1622  {
1623  if( memcmp( all_tuples_list->vi_rd + ( end_idx * num_tags ), all_tuples_list->vi_rd + ( last_idx * num_tags ),
1624  int_width ) == 0 )
1625  {
1626  // Values equal - skip
1627  last_idx += 1;
1628  }
1629  else
1630  {
1631  // Values different - copy
1632  // Move up the end index
1633  end_idx += 1;
1634  memcpy( all_tuples_list->vi_wr + ( end_idx * num_tags ), all_tuples_list->vi_rd + ( last_idx * num_tags ),
1635  int_width );
1636  last_idx += 1;
1637  }
1638  }
1639  // Update the count in all_tuples_list
1640  all_tuples_list->set_n( end_idx + 1 );
1641 
1642  // Resize the tuple_list
1643  all_tuples_list->resize( all_tuples_list->get_n() );
1644 
1645  // Set the output parameter
1646  *unique_tuples = all_tuples_list;
1647 
1648  return MB_SUCCESS;
1649 }

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

1112 {
1113  double tempField;
1114 
1115  // Get the tag values at the vertices
1116  ErrorCode result = mbImpl->tag_get_data( tag, &elem, 1, &tempField );
1117  if( MB_SUCCESS != result ) return result;
1118 
1119  field = tempField;
1120 
1121  return MB_SUCCESS;
1122 }

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

1510 {
1511  ErrorCode err;
1512  std::vector< Tag > tag_handles;
1513 
1514  for( unsigned int t = 0; t < num_tags; t++ )
1515  {
1516  // Get tag handle & size
1517  Tag th;
1518  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 );
1519  tag_handles.push_back( th );
1520  }
1521 
1522  return create_tuples( ent_sets, &tag_handles[0], num_tags, tuple_list );
1523 }

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

1529 {
1530  // ASSUMPTION: All tags are of type integer. This may need to be expanded in future.
1531  ErrorCode err;
1532 
1533  // Allocate a tuple_list for the number of entity sets passed in
1534  TupleList* tag_tuples = new TupleList( num_tags, 0, 0, 0, (int)ent_sets.size() );
1535  // tag_tuples->initialize(num_tags, 0, 0, 0, num_sets);
1536  uint mi, ml, mul, mr;
1537  tag_tuples->getTupleSize( mi, ml, mul, mr );
1538  tag_tuples->enableWriteAccess();
1539 
1540  if( mi == 0 ) ERRORR( "Failed to initialize tuple_list.", MB_FAILURE );
1541 
1542  // Loop over the filtered entity sets retrieving each matching tag value one by one.
1543  int val;
1544  for( unsigned int i = 0; i < ent_sets.size(); i++ )
1545  {
1546  for( unsigned int j = 0; j < num_tags; j++ )
1547  {
1548  EntityHandle set_handle = ent_sets[i];
1549  err = mbImpl->tag_get_data( tag_handles[j], &set_handle, 1, &val );ERRORR( "Failed to get integer tag data.", err );
1550  tag_tuples->vi_wr[i * mi + j] = val;
1551  }
1552 
1553  // If we get here there was no error so increment n in the tuple_list
1554  tag_tuples->inc_n();
1555  }
1556  tag_tuples->disableWriteAccess();
1557  *tuples = tag_tuples;
1558 
1559  return MB_SUCCESS;
1560 }

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

1208 {
1209  // SLAVE START ****************************************************************
1210  ErrorCode err;
1211  int ierr = 0;
1212 
1213  // Setup data for parallel computing
1214  int nprocs, rank;
1215  ierr = MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
1216  ERRORMPI( "Getting number of procs failed.", ierr );
1217  ierr = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
1218  ERRORMPI( "Getting rank failed.", ierr );
1219 
1220  // Get the integrated field value for each group(vector) of entities.
1221  // If no entities are in a group then a zero will be put in the list
1222  // of return values.
1223  unsigned int num_ent_grps = entity_groups.size();
1224  std::vector< double > integ_vals( num_ent_grps );
1225 
1226  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 );
1227  // SLAVE END ****************************************************************
1228 
1229  // SLAVE/MASTER START #########################################################
1230  // Send list of integrated values back to master proc. The ordering of the
1231  // values will match the ordering of the entity groups (i.e. vector of vectors)
1232  // sent from master to slaves earlier. The values for each entity group will
1233  // be summed during the transfer.
1234  std::vector< double > sum_integ_vals( num_ent_grps );
1235 
1236  if( nprocs > 1 )
1237  {
1238  // If parallel then send the values back to the master.
1239  ierr = MPI_Reduce( &integ_vals[0], &sum_integ_vals[0], num_ent_grps, MPI_DOUBLE, MPI_SUM, MASTER_PROC,
1240  myPc->proc_config().proc_comm() );
1241  ERRORMPI( "Transfer and reduction of integrated values failed.", ierr );
1242  }
1243  else
1244  {
1245  // Otherwise just copy the vector
1246  sum_integ_vals = integ_vals;
1247  }
1248  // SLAVE/MASTER END #########################################################
1249 
1250  // MASTER START ***************************************************************
1251  // Calculate the normalization factor for each group by taking the
1252  // inverse of each integrated field value. Put the normalization factor
1253  // for each group back into the list in the same order.
1254  for( unsigned int i = 0; i < num_ent_grps; i++ )
1255  {
1256  double val = sum_integ_vals[i];
1257  if( fabs( val ) > 1e-8 )
1258  sum_integ_vals[i] = 1.0 / val;
1259  else
1260  {
1261  sum_integ_vals[i] = 0.0; /* VSM: not sure what we should do here ? */
1262  /* commenting out error below since if integral(value)=0.0, then normalization
1263  is probably unnecessary to start with ? */
1264  /* ERRORR("Integrating an invalid field -- integral("<<norm_tag<<") = "<<val<<".", err);
1265  */
1266  }
1267  }
1268  // MASTER END ***************************************************************
1269 
1270  // MASTER/SLAVE START #########################################################
1271  if( nprocs > 1 )
1272  {
1273  // If parallel then broadcast the normalization factors to the procs.
1274  ierr = MPI_Bcast( &sum_integ_vals[0], num_ent_grps, MPI_DOUBLE, MASTER_PROC, myPc->proc_config().proc_comm() );
1275  ERRORMPI( "Broadcast of normalization factors failed.", ierr );
1276  }
1277  // MASTER/SLAVE END #########################################################
1278 
1279  // SLAVE START ****************************************************************
1280  // Save the normalization factors to a new tag with name of norm_tag's value
1281  // and the string "_normF" appended. This new tag will be created on the entity
1282  // set that contains all of the entities from a group.
1283 
1284  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 );
1285  // SLAVE END ****************************************************************
1286 
1287  return err;
1288 }

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

1931 {
1932  numPointsOfInterest = targ_elems.size() * _ntot;
1933  vpos.resize( 3 * numPointsOfInterest );
1934  int ielem = 0;
1935  for( Range::iterator eit = targ_elems.begin(); eit != targ_elems.end(); ++eit, ielem += _ntot * 3 )
1936  {
1937  EntityHandle eh = *eit;
1938  const double* xval;
1939  const double* yval;
1940  const double* zval;
1941  ErrorCode rval = mbImpl->tag_get_by_ptr( _xm1Tag, &eh, 1, (const void**)&xval );
1942  if( moab::MB_SUCCESS != rval )
1943  {
1944  std::cout << "Can't get xm1 values \n";
1945  return MB_FAILURE;
1946  }
1947  rval = mbImpl->tag_get_by_ptr( _ym1Tag, &eh, 1, (const void**)&yval );
1948  if( moab::MB_SUCCESS != rval )
1949  {
1950  std::cout << "Can't get ym1 values \n";
1951  return MB_FAILURE;
1952  }
1953  rval = mbImpl->tag_get_by_ptr( _zm1Tag, &eh, 1, (const void**)&zval );
1954  if( moab::MB_SUCCESS != rval )
1955  {
1956  std::cout << "Can't get zm1 values \n";
1957  return MB_FAILURE;
1958  }
1959  // Now, in a stride, populate vpos
1960  for( int i = 0; i < _ntot; i++ )
1961  {
1962  vpos[ielem + 3 * i] = xval[i];
1963  vpos[ielem + 3 * i + 1] = yval[i];
1964  vpos[ielem + 3 * i + 2] = zval[i];
1965  }
1966  }
1967 
1968  return MB_SUCCESS;
1969 }

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

1657 {
1658  ErrorCode err;
1659 
1660  std::vector< std::vector< EntityHandle > >::iterator iter_i;
1661  std::vector< EntityHandle >::iterator iter_j;
1662  double grp_intrgr_val, intgr_val;
1663 
1664  // Get the tag handle for norm_tag
1665  Tag norm_hdl;
1666  err =
1667  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 );
1668 
1669  // Check size of integ_vals vector
1670  if( integ_vals.size() != groups.size() ) integ_vals.resize( groups.size() );
1671 
1672  // Loop over the groups(vectors) of entities
1673  unsigned int i;
1674  for( i = 0, iter_i = groups.begin(); iter_i != groups.end(); i++, ++iter_i )
1675  {
1676  grp_intrgr_val = 0;
1677 
1678  // Loop over the all the entities in the group, integrating
1679  // the field_fn over the entity in iter_j
1680  for( iter_j = ( *iter_i ).begin(); iter_j != ( *iter_i ).end(); ++iter_j )
1681  {
1682  EntityHandle ehandle = ( *iter_j );
1683 
1684  // Check that the entity in iter_j is of the same dimension as the
1685  // integ_type we are performing
1686  EntityType j_type;
1687  j_type = mbImpl->type_from_handle( ehandle );ERRORR( "Failed to get entity type.", err );
1688  // Skip any entities in the group that are not of the type being considered
1689  if( ( integ_type == VOLUME ) && ( j_type < MBTET || j_type >= MBENTITYSET ) ) continue;
1690 
1691  intgr_val = 0;
1692 
1693  // Retrieve the vertices from the element
1694  const EntityHandle* verts = NULL;
1695  int connectivity_size = 0;
1696 
1697  err = mbImpl->get_connectivity( ehandle, verts, connectivity_size, false );ERRORR( "Failed to get vertices from entity.", err );
1698 
1699  // Get the vertex coordinates and the field values at the vertices.
1700  double* coords = (double*)malloc( sizeof( double ) * ( 3 * connectivity_size ) );
1701  /* TODO: VSM: check if this works for lower dimensions also without problems */
1702  /* if (3 == geom_dim) */
1703  err = mbImpl->get_coords( verts, connectivity_size, coords );ERRORR( "Failed to get vertex coordinates.", err );
1704 
1705  /* allocate the field data array */
1706  double* vfield = (double*)malloc( sizeof( double ) * ( connectivity_size ) );
1707  err = mbImpl->tag_get_data( norm_hdl, verts, connectivity_size, vfield );
1708  if( MB_SUCCESS != err )
1709  {
1710  free( coords );
1711  }
1712  ERRORR( "Failed to get vertex coordinates.", err );
1713 
1714  // Get coordinates of all corner vertices (in normal order) and
1715  // put in array of CartVec.
1716  std::vector< CartVect > vertices( connectivity_size );
1717 
1718  // Put the vertices into a CartVect vector
1719  double* x = coords;
1720  for( int j = 0; j < connectivity_size; j++, x += 3 )
1721  {
1722  vertices[j] = CartVect( x );
1723  }
1724  free( coords );
1725 
1726  moab::Element::Map* elemMap;
1727  if( j_type == MBHEX )
1728  {
1729  if( connectivity_size == 8 )
1730  elemMap = new moab::Element::LinearHex( vertices );
1731  else
1732  elemMap = new moab::Element::QuadraticHex( vertices );
1733  }
1734  else if( j_type == MBTET )
1735  {
1736  elemMap = new moab::Element::LinearTet( vertices );
1737  }
1738  else if( j_type == MBQUAD )
1739  {
1740  elemMap = new moab::Element::LinearQuad( vertices );
1741  }
1742  /*
1743  else if (j_type == MBTRI) {
1744  elemMap = new moab::Element::LinearTri(vertices);
1745  }
1746  */
1747  else if( j_type == MBEDGE )
1748  {
1749  elemMap = new moab::Element::LinearEdge( vertices );
1750  }
1751  else
1752  ERRORR( "Unknown topology type.", MB_UNSUPPORTED_OPERATION );
1753 
1754  // Set the vertices in the Map and perform the integration
1755  try
1756  {
1757  /* VSM: Do we need this call ?? */
1758  // elemMap->set_vertices(vertices);
1759 
1760  // Perform the actual integration over the element
1761  intgr_val = elemMap->integrate_scalar_field( vfield );
1762 
1763  // Combine the result with those of the group
1764  grp_intrgr_val += intgr_val;
1765  }
1767  {
1768  std::cerr << "Failed to set vertices on Element::Map." << std::endl;
1769  }
1771  {
1772  std::cerr << "Failed to get inverse evaluation of coordinate on Element::Map." << std::endl;
1773  }
1774 
1775  delete( elemMap );
1776  free( vfield );
1777  }
1778 
1779  // Set the group integrated value in the vector
1780  integ_vals[i] = grp_intrgr_val;
1781  }
1782 
1783  return err;
1784 }

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

1299 {
1300  ErrorCode err;
1301  std::vector< Tag > tag_handles;
1302 
1303  for( int t = 0; t < num_tags; t++ )
1304  {
1305  // Get tag handle & size
1306  Tag th;
1307  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 );
1308  tag_handles.push_back( th );
1309  }
1310 
1311  return get_matching_entities( root_set, &tag_handles[0], tag_values, num_tags, entity_sets, entity_groups );
1312 }

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

1321 {
1322  // SLAVE START ****************************************************************
1323 
1324  // Setup data for parallel computing
1325  ErrorCode err;
1326  int ierr = 0;
1327  int nprocs, rank;
1328  ierr = MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
1329  ERRORMPI( "Getting number of procs failed.", ierr );
1330  ierr = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
1331  ERRORMPI( "Getting rank failed.", ierr );
1332 
1333  Range ent_sets;
1334  err =
1335  mbImpl->get_entities_by_type_and_tag( root_set, moab::MBENTITYSET, tag_handles, (const void* const*)tag_values,
1336  num_tags, ent_sets, Interface::INTERSECT, false );ERRORR( "Core::get_entities_by_type_and_tag failed.", err );
1337 
1338  TupleList* tag_list = NULL;
1339  err = create_tuples( ent_sets, tag_handles, num_tags, &tag_list );ERRORR( "Failed to create tuples from entity sets.", err );
1340 
1341  // Free up range
1342  ent_sets.clear();
1343  // SLAVE END ****************************************************************
1344 
1345  // If we are running in a multi-proc session then send tuple list back to master
1346  // proc for consolidation. Otherwise just copy the pointer to the tuple_list.
1347  TupleList* cons_tuples;
1348  if( nprocs > 1 )
1349  {
1350  // SLAVE/MASTER START #########################################################
1351 
1352  // Pack the tuple_list in a buffer.
1353  uint* tuple_buf;
1354  int tuple_buf_len;
1355  tuple_buf_len = pack_tuples( tag_list, (void**)&tuple_buf );
1356 
1357  // Free tag_list here as its not used again if nprocs > 1
1358  tag_list->reset();
1359 
1360  // Send back the buffer sizes to the master proc
1361  int* recv_cnts = (int*)malloc( nprocs * sizeof( int ) );
1362  int* offsets = (int*)malloc( nprocs * sizeof( int ) );
1363  uint* all_tuples_buf = NULL;
1364 
1365  MPI_Gather( &tuple_buf_len, 1, MPI_INT, recv_cnts, 1, MPI_INT, MASTER_PROC, myPc->proc_config().proc_comm() );
1366  ERRORMPI( "Gathering buffer sizes failed.", err );
1367 
1368  // Allocate a buffer large enough for all the data
1369  if( rank == MASTER_PROC )
1370  {
1371  int all_tuples_len = recv_cnts[0];
1372  offsets[0] = 0;
1373  for( int i = 1; i < nprocs; i++ )
1374  {
1375  offsets[i] = offsets[i - 1] + recv_cnts[i - 1];
1376  all_tuples_len += recv_cnts[i];
1377  }
1378 
1379  all_tuples_buf = (uint*)malloc( all_tuples_len * sizeof( uint ) );
1380  }
1381 
1382  // Send all buffers to the master proc for consolidation
1383  MPI_Gatherv( (void*)tuple_buf, tuple_buf_len, MPI_INT, (void*)all_tuples_buf, recv_cnts, offsets, MPI_INT,
1385  ERRORMPI( "Gathering tuple_lists failed.", err );
1386  free( tuple_buf ); // malloc'd in pack_tuples
1387 
1388  if( rank == MASTER_PROC )
1389  {
1390  // Unpack the tuple_list from the buffer.
1391  TupleList** tl_array = (TupleList**)malloc( nprocs * sizeof( TupleList* ) );
1392  for( int i = 0; i < nprocs; i++ )
1393  unpack_tuples( (void*)&all_tuples_buf[offsets[i]], &tl_array[i] );
1394 
1395  // Free all_tuples_buf here as it is only allocated on the MASTER_PROC
1396  free( all_tuples_buf );
1397  // SLAVE/MASTER END #########################################################
1398 
1399  // MASTER START ***************************************************************
1400  // Consolidate all tuple_lists into one tuple_list with no duplicates.
1401  err = consolidate_tuples( tl_array, nprocs, &cons_tuples );ERRORR( "Failed to consolidate tuples.", err );
1402 
1403  for( int i = 0; i < nprocs; i++ )
1404  tl_array[i]->reset();
1405  free( tl_array );
1406  // MASTER END ***************************************************************
1407  }
1408 
1409  // Free offsets and recv_cnts as they are allocated on all procs
1410  free( offsets );
1411  free( recv_cnts );
1412 
1413  // MASTER/SLAVE START #########################################################
1414  // Broadcast condensed tuple list back to all procs.
1415  uint* ctl_buf;
1416  int ctl_buf_sz;
1417  if( rank == MASTER_PROC ) ctl_buf_sz = pack_tuples( cons_tuples, (void**)&ctl_buf );
1418 
1419  // Send buffer size
1420  ierr = MPI_Bcast( &ctl_buf_sz, 1, MPI_INT, MASTER_PROC, myPc->proc_config().proc_comm() );
1421  ERRORMPI( "Broadcasting tuple_list size failed.", ierr );
1422 
1423  // Allocate a buffer in the other procs
1424  if( rank != MASTER_PROC ) ctl_buf = (uint*)malloc( ctl_buf_sz * sizeof( uint ) );
1425 
1426  ierr = MPI_Bcast( (void*)ctl_buf, ctl_buf_sz, MPI_INT, MASTER_PROC, myPc->proc_config().proc_comm() );
1427  ERRORMPI( "Broadcasting tuple_list failed.", ierr );
1428 
1429  if( rank != MASTER_PROC ) unpack_tuples( ctl_buf, &cons_tuples );
1430  free( ctl_buf );
1431  // MASTER/SLAVE END #########################################################
1432  }
1433  else
1434  cons_tuples = tag_list;
1435 
1436  // SLAVE START ****************************************************************
1437  // Loop over the tuple list getting the entities with the tags in the tuple_list entry
1438  uint mi, ml, mul, mr;
1439  cons_tuples->getTupleSize( mi, ml, mul, mr );
1440 
1441  for( unsigned int i = 0; i < cons_tuples->get_n(); i++ )
1442  {
1443  // Get Entity Sets that match the tags and values.
1444 
1445  // Convert the data in the tuple_list to an array of pointers to the data
1446  // in the tuple_list as that is what the iMesh API call is expecting.
1447  int** vals = (int**)malloc( mi * sizeof( int* ) );
1448  for( unsigned int j = 0; j < mi; j++ )
1449  vals[j] = (int*)&( cons_tuples->vi_rd[( i * mi ) + j] );
1450 
1451  // Get entities recursively based on type and tag data
1452  err = mbImpl->get_entities_by_type_and_tag( root_set, moab::MBENTITYSET, tag_handles, (const void* const*)vals,
1453  mi, ent_sets, Interface::INTERSECT, false );ERRORR( "Core::get_entities_by_type_and_tag failed.", err );
1454  if( debug ) std::cout << "ent_sets_size=" << ent_sets.size() << std::endl;
1455 
1456  // Free up the array of pointers
1457  free( vals );
1458 
1459  // Loop over the entity sets and then free the memory for ent_sets.
1460  std::vector< EntityHandle > ent_set_hdls;
1461  std::vector< EntityHandle > ent_hdls;
1462  for( unsigned int j = 0; j < ent_sets.size(); j++ )
1463  {
1464  // Save the entity set
1465  ent_set_hdls.push_back( ent_sets[j] );
1466 
1467  // Get all entities for the entity set
1468  Range ents;
1469 
1470  /* VSM: do we need to filter out entity sets ? */
1471  err = mbImpl->get_entities_by_handle( ent_sets[j], ents, false );ERRORR( "Core::get_entities_by_handle failed.", err );
1472  if( debug ) std::cout << "ents_size=" << ents.size() << std::endl;
1473 
1474  // Save all of the entities from the entity set and free the memory for ents.
1475  for( unsigned int k = 0; k < ents.size(); k++ )
1476  {
1477  ent_hdls.push_back( ents[k] );
1478  }
1479  ents.clear();
1480  if( debug ) std::cout << "ent_hdls.size=" << ent_hdls.size() << std::endl;
1481  }
1482 
1483  // Free the entity set list for next tuple iteration.
1484  ent_sets.clear();
1485 
1486  // Push ent_set_hdls onto entity_sets, ent_hdls onto entity_groups
1487  // and clear both ent_set_hdls and ent_hdls.
1488  entity_sets->push_back( ent_set_hdls );
1489  ent_set_hdls.clear();
1490  entity_groups->push_back( ent_hdls );
1491  ent_hdls.clear();
1492  if( debug )
1493  std::cout << "entity_sets->size=" << entity_sets->size()
1494  << ", entity_groups->size=" << entity_groups->size() << std::endl;
1495  }
1496 
1497  cons_tuples->reset();
1498  // SLAVE END ****************************************************************
1499 
1500  return err;
1501 }

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.

200 {
201  /*void * _spectralSource;
202  void * _spectralTarget;*/
203 
204  moab::Range spectral_sets;
205  moab::Tag sem_tag;
206  int sem_dims[3];
207  ErrorCode rval = mbImpl->tag_get_handle( "SEM_DIMS", 3, moab::MB_TYPE_INTEGER, sem_tag );
208  if( moab::MB_SUCCESS != rval )
209  {
210  std::cout << "Can't find tag, no spectral set\n";
211  return MB_SUCCESS; // Nothing to do, no spectral elements
212  }
213  rval = mbImpl->get_entities_by_type_and_tag( rootSource, moab::MBENTITYSET, &sem_tag, NULL, 1, spectral_sets );
214  if( moab::MB_SUCCESS != rval || spectral_sets.empty() )
215  std::cout << "Can't get sem set on source\n";
216  else
217  {
218  moab::EntityHandle firstSemSet = spectral_sets[0];
219  rval = mbImpl->tag_get_data( sem_tag, &firstSemSet, 1, (void*)sem_dims );
220  if( moab::MB_SUCCESS != rval ) return MB_FAILURE;
221 
222  if( sem_dims[0] != sem_dims[1] || sem_dims[0] != sem_dims[2] )
223  {
224  std::cout << " dimensions are different. bail out\n";
225  return MB_FAILURE;
226  }
227  // Repeat for target sets
228  spectral_sets.empty();
229  // Now initialize a source spectral element !
230  _spectralSource = new moab::Element::SpectralHex( sem_dims[0] );
231  specSou = true;
232  }
233  // Repeat for target source
234  rval = mbImpl->get_entities_by_type_and_tag( rootTarget, moab::MBENTITYSET, &sem_tag, NULL, 1, spectral_sets );
235  if( moab::MB_SUCCESS != rval || spectral_sets.empty() )
236  std::cout << "Can't get sem set on target\n";
237  else
238  {
239  moab::EntityHandle firstSemSet = spectral_sets[0];
240  rval = mbImpl->tag_get_data( sem_tag, &firstSemSet, 1, (void*)sem_dims );
241  if( moab::MB_SUCCESS != rval ) return MB_FAILURE;
242 
243  if( sem_dims[0] != sem_dims[1] || sem_dims[0] != sem_dims[2] )
244  {
245  std::cout << " dimensions are different. bail out\n";
246  return MB_FAILURE;
247  }
248  // Repeat for target sets
249  spectral_sets.empty();
250  // Now initialize a target spectral element !
251  _spectralTarget = new moab::Element::SpectralHex( sem_dims[0] );
252  specTar = true;
253  }
254 
255  _ntot = sem_dims[0] * sem_dims[1] * sem_dims[2];
257  if( moab::MB_SUCCESS != rval )
258  {
259  std::cout << "Can't get xm1tag \n";
260  return MB_FAILURE;
261  }
263  if( moab::MB_SUCCESS != rval )
264  {
265  std::cout << "Can't get ym1tag \n";
266  return MB_FAILURE;
267  }
269  if( moab::MB_SUCCESS != rval )
270  {
271  std::cout << "Can't get zm1tag \n";
272  return MB_FAILURE;
273  }
274 
275  return MB_SUCCESS;
276 }

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.

84 {
85  Range local_ents;
86 
87  // Get entities on the local part
88  ErrorCode result = MB_SUCCESS;
89  if( myPc )
90  {
91  result = myPc->get_part_entities( local_ents, max_dim );
92  if( local_ents.empty() )
93  {
94  max_dim--;
95  result = myPc->get_part_entities( local_ents, max_dim ); // go one dimension lower
96  // right now, this is used for spherical meshes only
97  }
98  }
99  else
100  local_ents = myRange;
101  if( MB_SUCCESS != result || local_ents.empty() )
102  {
103  std::cout << "Problems getting source entities" << std::endl;
104  return result;
105  }
106 
107  // Build the tree for local processor
108  int max_per_leaf = 6;
109  for( int i = 0; i < numIts; i++ )
110  {
111  std::ostringstream str;
112  str << "PLANE_SET=0;"
113  << "MAX_PER_LEAF=" << max_per_leaf << ";";
114  if( spherical && !local_ents.empty() )
115  {
116  // get coordinates of one vertex, and use for radius computation
117  EntityHandle elem = local_ents[0];
118  const EntityHandle* conn;
119  int numn = 0;
120  mbImpl->get_connectivity( elem, conn, numn );
121  CartVect pos0;
122  mbImpl->get_coords( conn, 1, &( pos0[0] ) );
123  double radius = pos0.length();
124  str << "SPHERICAL=true;RADIUS=" << radius << ";";
125  }
126  FileOptions opts( str.str().c_str() );
127  myTree = new AdaptiveKDTree( mbImpl );
128  result = myTree->build_tree( local_ents, &localRoot, &opts );
129  if( MB_SUCCESS != result )
130  {
131  std::cout << "Problems building tree";
132  if( numIts != i )
133  {
134  delete myTree;
135  max_per_leaf *= 2;
136  std::cout << "; increasing elements/leaf to " << max_per_leaf << std::endl;
137  }
138  else
139  {
140  std::cout << "; exiting" << std::endl;
141  return result;
142  }
143  }
144  else
145  break; // Get out of tree building
146  }
147 
148  // Get the bounding box for local tree
149  if( myPc )
150  allBoxes.resize( 6 * myPc->proc_config().proc_size() );
151  else
152  allBoxes.resize( 6 );
153 
154  unsigned int my_rank = ( myPc ? myPc->proc_config().proc_rank() : 0 );
155  BoundBox box;
156  result = myTree->get_bounding_box( box, &localRoot );
157  if( MB_SUCCESS != result ) return result;
158  box.bMin.get( &allBoxes[6 * my_rank] );
159  box.bMax.get( &allBoxes[6 * my_rank + 3] );
160 
161  // Now communicate to get all boxes
162  if( myPc )
163  {
164  int mpi_err;
165 #if( MPI_VERSION >= 2 )
166  // Use "in place" option
167  mpi_err = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &allBoxes[0], 6, MPI_DOUBLE,
168  myPc->proc_config().proc_comm() );
169 #else
170  {
171  std::vector< double > allBoxes_tmp( 6 * myPc->proc_config().proc_size() );
172  mpi_err = MPI_Allgather( &allBoxes[6 * my_rank], 6, MPI_DOUBLE, &allBoxes_tmp[0], 6, MPI_DOUBLE,
173  myPc->proc_config().proc_comm() );
174  allBoxes = allBoxes_tmp;
175  }
176 #endif
177  if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
178  }
179 
180  /* std::ostringstream blah;
181  for(int i = 0; i < allBoxes.size(); i++)
182  blah << allBoxes[i] << " ";
183  std::cout << blah.str() << "\n";*/
184 
185 #ifdef VERBOSE
186  double min[3] = { 0, 0, 0 }, max[3] = { 0, 0, 0 };
187  unsigned int dep;
188  myTree->get_info( localRoot, min, max, dep );
189  std::cout << "Proc " << my_rank << ": box min/max, tree depth = (" << min[0] << "," << min[1] << "," << min[2]
190  << "), (" << max[0] << "," << max[1] << "," << max[2] << "), " << dep << std::endl;
191 #endif
192 
193  return result;
194 }

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

1022 {
1023  if( _spectralSource )
1024  {
1025  // Get tag values at the GL points for some field (Tag)
1026  const double* vx;
1027  ErrorCode rval = mbImpl->tag_get_by_ptr( tag, &elem, 1, (const void**)&vx );
1028  if( moab::MB_SUCCESS != rval )
1029  {
1030  std::cout << "Can't get field values for the tag \n";
1031  return MB_FAILURE;
1032  }
1033  Element::SpectralHex* spcHex = (Element::SpectralHex*)_spectralSource;
1034  field = spcHex->evaluate_scalar_field( nat_coord, vx );
1035  }
1036  else
1037  {
1038  double vfields[27]; // Will work for linear hex, quadratic hex or Tets
1039  moab::Element::Map* elemMap = NULL;
1040  int num_verts = 0;
1041  // Get the EntityType
1042  // Get the tag values at the vertices
1043  const EntityHandle* connect;
1044  int num_connect;
1045  ErrorCode result = mbImpl->get_connectivity( elem, connect, num_connect );
1046  if( MB_SUCCESS != result ) return result;
1047  EntityType etype = mbImpl->type_from_handle( elem );
1048  if( MBHEX == etype )
1049  {
1050  if( 8 == num_connect )
1051  {
1052  elemMap = new moab::Element::LinearHex();
1053  num_verts = 8;
1054  }
1055  else
1056  { /* (MBHEX == etype && 27 == num_connect) */
1057  elemMap = new moab::Element::QuadraticHex();
1058  num_verts = 27;
1059  }
1060  }
1061  else if( MBTET == etype )
1062  {
1063  elemMap = new moab::Element::LinearTet();
1064  num_verts = 4;
1065  }
1066  else if( MBQUAD == etype )
1067  {
1068  elemMap = new moab::Element::LinearQuad();
1069  num_verts = 4;
1070  }
1071  else if( MBTRI == etype )
1072  {
1073  elemMap = new moab::Element::LinearTri();
1074  num_verts = 3;
1075  }
1076  else
1077  return MB_FAILURE;
1078 
1079  result = mbImpl->tag_get_data( tag, connect, std::min( num_verts, num_connect ), vfields );
1080  if( MB_SUCCESS != result )
1081  {
1082  delete elemMap;
1083  return result;
1084  }
1085 
1086  // Function for the interpolation
1087  field = 0;
1088 
1089  // Check the number of vertices
1090  assert( num_connect >= num_verts );
1091 
1092  // Calculate the field
1093  try
1094  {
1095  field = elemMap->evaluate_scalar_field( nat_coord, vfields );
1096  }
1098  {
1099  delete elemMap;
1100  return MB_FAILURE;
1101  }
1102 
1103  delete elemMap;
1104  }
1105 
1106  return MB_SUCCESS;
1107 }

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

672 {
673  // if (!((LINEAR_FE == method) || (CONSTANT == method)))
674  // return MB_FAILURE;
675 
676  // remote pts first
677  TupleList* tl_tmp = ( tl ? tl : targetPts );
678 
679  ErrorCode result = MB_SUCCESS;
680 
681  unsigned int pts_total = 0;
682  for( int i = 0; i < num_methods; i++ )
683  pts_total += points_per_method[i];
684 
685  // If tl was passed in non-NULL, just have those points, otherwise have targetPts plus
686  // locally mapped pts
687  if( pts_total != tl_tmp->get_n() ) return MB_FAILURE;
688 
689  TupleList tinterp;
690  tinterp.initialize( 5, 0, 0, 1, tl_tmp->get_n() );
691  int t = 0;
692  tinterp.enableWriteAccess();
693  for( int i = 0; i < num_methods; i++ )
694  {
695  for( int j = 0; j < points_per_method[i]; j++ )
696  {
697  tinterp.vi_wr[5 * t] = tl_tmp->vi_rd[3 * t];
698  tinterp.vi_wr[5 * t + 1] = tl_tmp->vi_rd[3 * t + 1];
699  tinterp.vi_wr[5 * t + 2] = tl_tmp->vi_rd[3 * t + 2];
700  tinterp.vi_wr[5 * t + 3] = methods[i];
701  tinterp.vi_wr[5 * t + 4] = i;
702  tinterp.vr_wr[t] = 0.0;
703  tinterp.inc_n();
704  t++;
705  }
706  }
707 
708  // Scatter/gather interpolation points
709  if( myPc )
710  {
711  ( myPc->proc_config().crystal_router() )->gs_transfer( 1, tinterp, 0 );
712 
713  // Perform interpolation on local source mesh; put results into
714  // tinterp.vr_wr[i]
715 
717  for( unsigned int i = 0; i < tinterp.get_n(); i++ )
718  {
719  int mindex = tinterp.vi_rd[5 * i + 2];
720  Method method = (Method)tinterp.vi_rd[5 * i + 3];
721  Tag tag = tags[tinterp.vi_rd[5 * i + 4]];
722 
723  result = MB_FAILURE;
724  if( LINEAR_FE == method || QUADRATIC_FE == method || SPHERICAL == method )
725  {
726  result = interp_field( mappedPts->vul_rd[mindex], CartVect( mappedPts->vr_wr + 3 * mindex ), tag,
727  tinterp.vr_wr[i] );
728  }
729  else if( CONSTANT == method )
730  {
731  result = constant_interp( mappedPts->vul_rd[mindex], tag, tinterp.vr_wr[i] );
732  }
733 
734  if( MB_SUCCESS != result ) return result;
735  }
736 
737  // Scatter/gather interpolation data
738  myPc->proc_config().crystal_router()->gs_transfer( 1, tinterp, 0 );
739  }
740 
741  // Copy the interpolated field as a unit
742  for( unsigned int i = 0; i < tinterp.get_n(); i++ )
743  interp_vals[tinterp.vi_rd[5 * i + 1]] = tinterp.vr_rd[i];
744 
745  // Done
746  return MB_SUCCESS;
747 }

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

650 {
651  Tag tag;
652  ErrorCode result;
653  if( _spectralSource )
654  {
655  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 << "\"" );
656  }
657  else
658  {
659  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 << "\"" );
660  }
661 
662  return interpolate( method, tag, interp_vals, tl, normalize );
663 }

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.

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

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.

438  {
439  return localRoot;
440  }

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.

325 {
326  assert( tl || store_local );
327 
328  // target_pts: TL(to_proc, tgt_index, x, y, z): tuples sent to source mesh procs representing
329  // pts to be located source_pts: TL(from_proc, tgt_index, src_index): results of source mesh
330  // proc point location, ready to send
331  // back to tgt procs; src_index of -1 indicates point not located (arguably not
332  // useful...)
333  TupleList target_pts;
334  target_pts.initialize( 2, 0, 0, 3, num_points );
335  target_pts.enableWriteAccess();
336 
337  TupleList source_pts;
338  mappedPts = new TupleList( 0, 0, 1, 3, target_pts.get_max() );
340 
341  source_pts.initialize( 3, 0, 0, 0, target_pts.get_max() );
342  source_pts.enableWriteAccess();
343 
344  mappedPts->set_n( 0 );
345  source_pts.set_n( 0 );
346  ErrorCode result;
347 
348  unsigned int my_rank = ( myPc ? myPc->proc_config().proc_rank() : 0 );
349  bool point_located;
350 
351  // For each point, find box(es) containing the point,
352  // appending results to tuple_list;
353  // keep local points separately, in local_pts, which has pairs
354  // of <local_index, mapped_index>, where mapped_index is the index
355  // into the mappedPts tuple list
356  for( unsigned int i = 0; i < 3 * num_points; i += 3 )
357  {
358 
359  std::vector< int > procs_to_send_to;
360  for( unsigned int j = 0; j < ( myPc ? myPc->proc_config().proc_size() : 0 ); j++ )
361  {
362  // Test if point is in proc's box
363  if( ( allBoxes[6 * j] <= xyz[i] + abs_eps ) && ( xyz[i] <= allBoxes[6 * j + 3] + abs_eps ) &&
364  ( allBoxes[6 * j + 1] <= xyz[i + 1] + abs_eps ) && ( xyz[i + 1] <= allBoxes[6 * j + 4] + abs_eps ) &&
365  ( allBoxes[6 * j + 2] <= xyz[i + 2] + abs_eps ) && ( xyz[i + 2] <= allBoxes[6 * j + 5] + abs_eps ) )
366  {
367  // If in this proc's box, will send to proc to test further
368  procs_to_send_to.push_back( j );
369  }
370  }
371  if( procs_to_send_to.empty() )
372  {
373 #ifdef VERBOSE
374  std::cout << " point index " << i / 3 << ": " << xyz[i] << " " << xyz[i + 1] << " " << xyz[i + 2]
375  << " not found in any box\n";
376 #endif
377  // try to find the closest box, and put it in that box, anyway
378  double min_dist = 1.e+20;
379  int index = -1;
380  for( unsigned int j = 0; j < ( myPc ? myPc->proc_config().proc_size() : 0 ); j++ )
381  {
382  BoundBox box( &allBoxes[6 * j] ); // form back the box
383  double distance = box.distance( &xyz[i] ); // will compute the distance in 3d, from the box
384  if( distance < min_dist )
385  {
386  index = j;
387  min_dist = distance;
388  }
389  }
390  if( index == -1 )
391  {
392  // need to abort early, nothing we can do
393  assert( "cannot locate any box for some points" );
394  // need a better exit strategy
395  }
396 #ifdef VERBOSE
397  std::cout << " point index " << i / 3 << " added to box for proc j:" << index << "\n";
398 #endif
399  procs_to_send_to.push_back( index ); // will send to just one proc, that has the closest box
400  }
401  // we finally decided to populate the tuple list for a list of processors
402  for( size_t k = 0; k < procs_to_send_to.size(); k++ )
403  {
404  unsigned int j = procs_to_send_to[k];
405  // Check size of tuple list, grow if we're at max
406  if( target_pts.get_n() == target_pts.get_max() )
407  target_pts.resize( std::max( 10.0, 1.5 * target_pts.get_max() ) );
408 
409  target_pts.vi_wr[2 * target_pts.get_n()] = j;
410  target_pts.vi_wr[2 * target_pts.get_n() + 1] = i / 3;
411 
412  target_pts.vr_wr[3 * target_pts.get_n()] = xyz[i];
413  target_pts.vr_wr[3 * target_pts.get_n() + 1] = xyz[i + 1];
414  target_pts.vr_wr[3 * target_pts.get_n() + 2] = xyz[i + 2];
415  target_pts.inc_n();
416  }
417 
418  } // end for (unsigned int i = 0; ..
419 
420  int num_to_me = 0;
421  for( unsigned int i = 0; i < target_pts.get_n(); i++ )
422  if( target_pts.vi_rd[2 * i] == (int)my_rank ) num_to_me++;
423 #ifdef VERBOSE
424  printf( "rank: %u local points: %u, nb sent target pts: %u mappedPts: %u num to me: %d \n", my_rank, num_points,
425  target_pts.get_n(), mappedPts->get_n(), num_to_me );
426 #endif
427  // Perform scatter/gather, to gather points to source mesh procs
428  if( myPc )
429  {
430  ( myPc->proc_config().crystal_router() )->gs_transfer( 1, target_pts, 0 );
431 
432  num_to_me = 0;
433  for( unsigned int i = 0; i < target_pts.get_n(); i++ )
434  {
435  if( target_pts.vi_rd[2 * i] == (int)my_rank ) num_to_me++;
436  }
437 #ifdef VERBOSE
438  printf( "rank: %u after first gs nb received_pts: %u; num_from_me = %d\n", my_rank, target_pts.get_n(),
439  num_to_me );
440 #endif
441  // After scatter/gather:
442  // target_pts.set_n(# points local proc has to map);
443  // target_pts.vi_wr[2*i] = proc sending point i
444  // target_pts.vi_wr[2*i + 1] = index of point i on sending proc
445  // target_pts.vr_wr[3*i..3*i + 2] = xyz of point i
446  //
447  // Mapping builds the tuple list:
448  // source_pts.set_n(target_pts.get_n())
449  // source_pts.vi_wr[3*i] = target_pts.vi_wr[2*i] = sending proc
450  // source_pts.vi_wr[3*i + 1] = index of point i on sending proc
451  // source_pts.vi_wr[3*i + 2] = index of mapped point (-1 if not mapped)
452  //
453  // Also, mapping builds local tuple_list mappedPts:
454  // mappedPts->set_n( # mapped points );
455  // mappedPts->vul_wr[i] = local handle of mapped entity
456  // mappedPts->vr_wr[3*i..3*i + 2] = natural coordinates in mapped entity
457 
458  // Test target points against my elements
459  for( unsigned i = 0; i < target_pts.get_n(); i++ )
460  {
461  result = test_local_box( target_pts.vr_wr + 3 * i, target_pts.vi_rd[2 * i], target_pts.vi_rd[2 * i + 1], i,
462  point_located, rel_eps, abs_eps, &source_pts );
463  if( MB_SUCCESS != result ) return result;
464  }
465 
466  // No longer need target_pts
467  target_pts.reset();
468 #ifdef VERBOSE
469  printf( "rank: %u nb sent source pts: %u, mappedPts now: %u\n", my_rank, source_pts.get_n(),
470  mappedPts->get_n() );
471 #endif
472  // Send target points back to target procs
473  ( myPc->proc_config().crystal_router() )->gs_transfer( 1, source_pts, 0 );
474 
475 #ifdef VERBOSE
476  printf( "rank: %u nb received source pts: %u\n", my_rank, source_pts.get_n() );
477 #endif
478  }
479 
480  // Store proc/index tuples in targetPts, and/or pass back to application;
481  // the tuple this gets stored to looks like:
482  // tl.set_n(# mapped points);
483  // tl.vi_wr[3*i] = remote proc mapping point
484  // tl.vi_wr[3*i + 1] = local index of mapped point
485  // tl.vi_wr[3*i + 2] = remote index of mapped point
486  //
487  // Local index is mapped into either myRange, holding the handles of
488  // local mapped entities, or myXyz, holding locations of mapped pts
489 
490  // Store information about located points
491  TupleList* tl_tmp;
492  if( !store_local )
493  tl_tmp = tl;
494  else
495  {
496  targetPts = new TupleList();
497  tl_tmp = targetPts;
498  }
499 
500  tl_tmp->initialize( 3, 0, 0, 0, num_points );
501  tl_tmp->set_n( num_points ); // Automatically sets tl to write_enabled
502  // Initialize so we know afterwards how many pts weren't located
503  std::fill( tl_tmp->vi_wr, tl_tmp->vi_wr + 3 * num_points, -1 );
504 
505  unsigned int local_pts = 0;
506  for( unsigned int i = 0; i < source_pts.get_n(); i++ )
507  {
508  if( -1 != source_pts.vi_rd[3 * i + 2] )
509  { // Why bother sending message saying "i don't have the point" if it gets discarded?
510  int tgt_index = 3 * source_pts.vi_rd[3 * i + 1];
511  // Prefer always entities that are local, from the source_pts
512  // if a local entity was already found to contain the target point, skip
513  // tl_tmp->vi_wr[tgt_index] is -1 initially, but it could already be set with
514  // a remote processor
515  if( tl_tmp->vi_wr[tgt_index] != (int)my_rank )
516  {
517  tl_tmp->vi_wr[tgt_index] = source_pts.vi_rd[3 * i];
518  tl_tmp->vi_wr[tgt_index + 1] = source_pts.vi_rd[3 * i + 1];
519  tl_tmp->vi_wr[tgt_index + 2] = source_pts.vi_rd[3 * i + 2];
520  }
521  }
522  }
523 
524  // Count missing points
525  unsigned int missing_pts = 0;
526  for( unsigned int i = 0; i < num_points; i++ )
527  {
528  if( tl_tmp->vi_rd[3 * i + 1] == -1 )
529  {
530  missing_pts++;
531 #ifdef VERBOSE
532  printf( "missing point at index i: %d -> %15.10f %15.10f %15.10f\n", i, xyz[3 * i], xyz[3 * i + 1],
533  xyz[3 * i + 2] );
534 #endif
535  }
536  else if( tl_tmp->vi_rd[3 * i] == (int)my_rank )
537  local_pts++;
538  }
539 #ifdef VERBOSE
540  printf( "rank: %u point location: wanted %u got %u locally, %u remote, missing %u\n", my_rank, num_points,
541  local_pts, num_points - missing_pts - local_pts, missing_pts );
542 #endif
543  assert( 0 == missing_pts ); // Will likely break on curved geometries
544 
545  // No longer need source_pts
546  source_pts.reset();
547 
548  // Copy into tl if passed in and storing locally
549  if( tl && store_local )
550  {
551  tl->initialize( 3, 0, 0, 0, num_points );
552  tl->enableWriteAccess();
553  memcpy( tl->vi_wr, tl_tmp->vi_rd, 3 * tl_tmp->get_n() * sizeof( int ) );
554  tl->set_n( tl_tmp->get_n() );
555  tl->disableWriteAccess();
556  }
557 
558  tl_tmp->disableWriteAccess();
559 
560  // Done
561  return MB_SUCCESS;
562 }

References allBoxes, moab::ProcConfig::crystal_router(), moab::TupleList::disableWriteAccess(), moab::BoundBox::distance(), moab::TupleList::enableWriteAccess(), ErrorCode, moab::TupleList::get_max(), moab::TupleList::get_n(), moab::TupleList::inc_n(), moab::TupleList::initialize(), mappedPts, MB_SUCCESS, myPc, moab::ParallelComm::proc_config(), moab::ProcConfig::proc_rank(), moab::ProcConfig::proc_size(), moab::TupleList::reset(), moab::TupleList::resize(), moab::TupleList::set_n(), targetPts, test_local_box(), moab::TupleList::vi_rd, moab::TupleList::vi_wr, and moab::TupleList::vr_wr.

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.

279 {
280  // Get locations
281  std::vector< double > locs( 3 * targ_ents.size() );
282  Range verts = targ_ents.subset_by_type( MBVERTEX );
283  ErrorCode rval = mbImpl->get_coords( verts, &locs[0] );
284  if( MB_SUCCESS != rval ) return rval;
285  // Now get other ents; reuse verts
286  unsigned int num_verts = verts.size();
287  verts = subtract( targ_ents, verts );
288  // Compute centroids
289  std::vector< EntityHandle > dum_conn( CN::MAX_NODES_PER_ELEMENT );
290  std::vector< double > dum_pos( CN::MAX_NODES_PER_ELEMENT );
291  const EntityHandle* conn;
292  int num_conn;
293  double* coords = &locs[num_verts];
294  // Do this here instead of a function to allow reuse of dum_pos and dum_conn
295  for( Range::const_iterator rit = verts.begin(); rit != verts.end(); ++rit )
296  {
297  rval = mbImpl->get_connectivity( *rit, conn, num_conn, false, &dum_conn );
298  if( MB_SUCCESS != rval ) return rval;
299  rval = mbImpl->get_coords( conn, num_conn, &dum_pos[0] );
300  if( MB_SUCCESS != rval ) return rval;
301  coords[0] = coords[1] = coords[2] = 0.0;
302  for( int i = 0; i < num_conn; i++ )
303  {
304  coords[0] += dum_pos[3 * i];
305  coords[1] += dum_pos[3 * i + 1];
306  coords[2] += dum_pos[3 * i + 2];
307  }
308  coords[0] /= num_conn;
309  coords[1] /= num_conn;
310  coords[2] /= num_conn;
311  coords += 3;
312  }
313 
314  if( store_local ) targetEnts = targ_ents;
315 
316  return locate_points( &locs[0], targ_ents.size(), rel_eps, abs_eps, tl, store_local );
317 }

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.

462  {
463  return mappedPts;
464  }

References mappedPts.

◆ mb_impl()

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

Definition at line 429 of file Coupler.hpp.

430  {
431  return mbImpl;
432  }

References mbImpl.

◆ my_id()

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

Definition at line 453 of file Coupler.hpp.

454  {
455  return myId;
456  }

References myId.

◆ my_pc()

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

Definition at line 445 of file Coupler.hpp.

446  {
447  return myPc;
448  }

References myPc.

◆ my_range()

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

Definition at line 457 of file Coupler.hpp.

458  {
459  return myRange;
460  }

References myRange.

◆ my_tree()

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

Definition at line 433 of file Coupler.hpp.

434  {
435  return myTree;
436  }

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

753 {
754  if( !myTree ) return MB_FAILURE;
755 
756  AdaptiveKDTreeIter treeiter;
757  ErrorCode result = myTree->get_tree_iterator( localRoot, treeiter );
758  if( MB_SUCCESS != result )
759  {
760  std::cout << "Problems getting iterator" << std::endl;
761  return result;
762  }
763 
764  EntityHandle closest_leaf;
765  if( epsilon )
766  {
767  std::vector< double > dists;
768  std::vector< EntityHandle > leaves;
769 
770  // Two tolerances
771  result = myTree->distance_search( xyz, epsilon, leaves,
772  /*iter_tol*/ epsilon,
773  /*inside_tol*/ 10 * epsilon, &dists, NULL, &localRoot );
774  if( leaves.empty() )
775  // Not found returns success here, with empty list, just like case with no epsilon
776  return MB_SUCCESS;
777  // Get closest leaf
778  double min_dist = *dists.begin();
779  closest_leaf = *leaves.begin();
780  std::vector< EntityHandle >::iterator vit = leaves.begin() + 1;
781  std::vector< double >::iterator dit = dists.begin() + 1;
782  for( ; vit != leaves.end() && min_dist; ++vit, ++dit )
783  {
784  if( *dit < min_dist )
785  {
786  min_dist = *dit;
787  closest_leaf = *vit;
788  }
789  }
790  }
791  else
792  {
793  result = myTree->point_search( xyz, treeiter, 1.0e-10, 1.0e-6, NULL, &localRoot );
794  if( MB_ENTITY_NOT_FOUND == result ) // Point is outside of myTree's bounding box
795  return MB_SUCCESS;
796  else if( MB_SUCCESS != result )
797  {
798  std::cout << "Problems getting leaf \n";
799  return result;
800  }
801  closest_leaf = treeiter.handle();
802  }
803 
804  // Find natural coordinates of point in element(s) in that leaf
805  CartVect tmp_nat_coords;
806  Range range_leaf;
807  result = mbImpl->get_entities_by_dimension( closest_leaf, max_dim, range_leaf, false );
808  if( MB_SUCCESS != result ) std::cout << "Problem getting leaf in a range" << std::endl;
809 
810  // Loop over the range_leaf
811  for( Range::iterator iter = range_leaf.begin(); iter != range_leaf.end(); ++iter )
812  {
813  // Test to find out in which entity the point is
814  // Get the EntityType and create the appropriate Element::Map subtype
815  // If spectral, do not need coordinates, just the GL points
816  EntityType etype = mbImpl->type_from_handle( *iter );
817  if( NULL != this->_spectralSource && MBHEX == etype )
818  {
819  EntityHandle eh = *iter;
820  const double* xval;
821  const double* yval;
822  const double* zval;
823  ErrorCode rval = mbImpl->tag_get_by_ptr( _xm1Tag, &eh, 1, (const void**)&xval );
824  if( moab::MB_SUCCESS != rval )
825  {
826  std::cout << "Can't get xm1 values \n";
827  return MB_FAILURE;
828  }
829  rval = mbImpl->tag_get_by_ptr( _ym1Tag, &eh, 1, (const void**)&yval );
830  if( moab::MB_SUCCESS != rval )
831  {
832  std::cout << "Can't get ym1 values \n";
833  return MB_FAILURE;
834  }
835  rval = mbImpl->tag_get_by_ptr( _zm1Tag, &eh, 1, (const void**)&zval );
836  if( moab::MB_SUCCESS != rval )
837  {
838  std::cout << "Can't get zm1 values \n";
839  return MB_FAILURE;
840  }
841  Element::SpectralHex* spcHex = (Element::SpectralHex*)_spectralSource;
842 
843  spcHex->set_gl_points( (double*)xval, (double*)yval, (double*)zval );
844  try
845  {
846  tmp_nat_coords = spcHex->ievaluate( CartVect( xyz ), epsilon ); // introduce
847  bool inside = spcHex->inside_nat_space( CartVect( tmp_nat_coords ), epsilon );
848  if( !inside )
849  {
850 #ifdef VERBOSE
851  std::cout << "point " << xyz[0] << " " << xyz[1] << " " << xyz[2]
852  << " is not converging inside hex " << mbImpl->id_from_handle( eh ) << "\n";
853 #endif
854  continue; // It is possible that the point is outside, so it will not converge
855  }
856  }
857  catch( Element::Map::EvaluationError& )
858  {
859  continue;
860  }
861  }
862  else
863  {
864  const EntityHandle* connect;
865  int num_connect;
866 
867  // Get connectivity
868  result = mbImpl->get_connectivity( *iter, connect, num_connect, true );
869  if( MB_SUCCESS != result ) return result;
870 
871  // Get coordinates of the vertices
872  std::vector< CartVect > coords_vert( num_connect );
873  result = mbImpl->get_coords( connect, num_connect, &( coords_vert[0][0] ) );
874  if( MB_SUCCESS != result )
875  {
876  std::cout << "Problems getting coordinates of vertices\n";
877  return result;
878  }
879  CartVect pos( xyz );
880  if( MBHEX == etype )
881  {
882  if( 8 == num_connect )
883  {
884  Element::LinearHex hexmap( coords_vert );
885  if( !hexmap.inside_box( pos, epsilon ) ) continue;
886  try
887  {
888  tmp_nat_coords = hexmap.ievaluate( pos, epsilon );
889  bool inside = hexmap.inside_nat_space( tmp_nat_coords, epsilon );
890  if( !inside ) continue;
891  }
892  catch( Element::Map::EvaluationError& )
893  {
894  continue;
895  }
896  }
897  else if( 27 == num_connect )
898  {
899  Element::QuadraticHex hexmap( coords_vert );
900  if( !hexmap.inside_box( pos, epsilon ) ) continue;
901  try
902  {
903  tmp_nat_coords = hexmap.ievaluate( pos, epsilon );
904  bool inside = hexmap.inside_nat_space( tmp_nat_coords, epsilon );
905  if( !inside ) continue;
906  }
907  catch( Element::Map::EvaluationError& )
908  {
909  continue;
910  }
911  }
912  else // TODO this case not treated yet, no interpolation
913  continue;
914  }
915  else if( MBTET == etype )
916  {
917  Element::LinearTet tetmap( coords_vert );
918  // This is just a linear solve; unless degenerate, will not except
919  tmp_nat_coords = tetmap.ievaluate( pos );
920  bool inside = tetmap.inside_nat_space( tmp_nat_coords, epsilon );
921  if( !inside ) continue;
922  }
923  else if( MBQUAD == etype && spherical )
924  {
925  Element::SphericalQuad sphermap( coords_vert );
926  /* skip box test, because it can filter out good elements with high curvature
927  * if (!sphermap.inside_box(pos, epsilon))
928  continue;*/
929  try
930  {
931  tmp_nat_coords = sphermap.ievaluate( pos, epsilon );
932  bool inside = sphermap.inside_nat_space( tmp_nat_coords, epsilon );
933  if( !inside ) continue;
934  }
935  catch( Element::Map::EvaluationError& )
936  {
937  continue;
938  }
939  }
940  else if( MBTRI == etype && spherical )
941  {
942  Element::SphericalTri sphermap( coords_vert );
943  /* skip box test, because it can filter out good elements with high curvature
944  * if (!sphermap.inside_box(pos, epsilon))
945  continue;*/
946  try
947  {
948  tmp_nat_coords = sphermap.ievaluate( pos, epsilon );
949  bool inside = sphermap.inside_nat_space( tmp_nat_coords, epsilon );
950  if( !inside ) continue;
951  }
952  catch( Element::Map::EvaluationError& )
953  {
954  continue;
955  }
956  }
957 
958  else if( MBQUAD == etype )
959  {
960  Element::LinearQuad quadmap( coords_vert );
961  if( !quadmap.inside_box( pos, epsilon ) ) continue;
962  try
963  {
964  tmp_nat_coords = quadmap.ievaluate( pos, epsilon );
965  bool inside = quadmap.inside_nat_space( tmp_nat_coords, epsilon );
966  if( !inside ) continue;
967  }
968  catch( Element::Map::EvaluationError& )
969  {
970  continue;
971  }
972  if( !quadmap.inside_nat_space( tmp_nat_coords, epsilon ) ) continue;
973  }
974  /*
975  else if (etype == MBTRI){
976  Element::LinearTri trimap(coords_vert);
977  if (!trimap.inside_box( pos, epsilon))
978  continue;
979  try {
980  tmp_nat_coords = trimap.ievaluate(pos, epsilon);
981  bool inside = trimap.inside_nat_space(tmp_nat_coords, epsilon);
982  if (!inside) continue;
983  }
984  catch (Element::Map::EvaluationError) {
985  continue;
986  }
987  if (!trimap.inside_nat_space(tmp_nat_coords, epsilon))
988  continue;
989  }
990  */
991  else if( etype == MBEDGE )
992  {
993  Element::LinearEdge edgemap( coords_vert );
994  try
995  {
996  tmp_nat_coords = edgemap.ievaluate( CartVect( xyz ), epsilon );
997  }
998  catch( Element::Map::EvaluationError )
999  {
1000  continue;
1001  }
1002  if( !edgemap.inside_nat_space( tmp_nat_coords, epsilon ) ) continue;
1003  }
1004  else
1005  {
1006  std::cout << "Entity not Hex/Tet/Quad/Tri/Edge. Please verify." << std::endl;
1007  continue;
1008  }
1009  }
1010  // If we get here then we've found the coordinates.
1011  // Save them and the entity and return success.
1012  entities.push_back( *iter );
1013  nat_coords.push_back( tmp_nat_coords );
1014  return MB_SUCCESS;
1015  }
1016 
1017  // Didn't find any elements containing the point
1018  return MB_SUCCESS;
1019 }

References _spectralSource, _xm1Tag, _ym1Tag, _zm1Tag, moab::Range::begin(), moab::AdaptiveKDTree::distance_search(), moab::Range::end(), entities, ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::AdaptiveKDTree::get_tree_iterator(), moab::AdaptiveKDTreeIter::handle(), moab::Interface::id_from_handle(), moab::Element::Map::ievaluate(), moab::Element::LinearTet::ievaluate(), moab::Element::SpectralHex::ievaluate(), moab::Element::SphericalQuad::ievaluate(), moab::Element::SphericalTri::ievaluate(), moab::Element::Map::inside_box(), moab::Element::LinearHex::inside_nat_space(), moab::Element::QuadraticHex::inside_nat_space(), moab::Element::LinearTet::inside_nat_space(), moab::Element::SpectralHex::inside_nat_space(), moab::Element::LinearQuad::inside_nat_space(), moab::Element::LinearTri::inside_nat_space(), moab::Element::LinearEdge::inside_nat_space(), localRoot, max_dim, MB_ENTITY_NOT_FOUND, MB_SUCCESS, MBEDGE, MBHEX, mbImpl, MBQUAD, MBTET, MBTRI, myTree, moab::AdaptiveKDTree::point_search(), moab::Element::SpectralHex::set_gl_points(), spherical, moab::Interface::tag_get_by_ptr(), and moab::Interface::type_from_handle().

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

1129 {
1130  ErrorCode err;
1131 
1132  // SLAVE START ****************************************************************
1133  // Search for entities based on tag_handles and tag_values
1134  std::vector< std::vector< EntityHandle > > entity_sets;
1135  std::vector< std::vector< EntityHandle > > entity_groups;
1136 
1137  // put the root_set into entity_sets
1138  std::vector< EntityHandle > ent_set;
1139  ent_set.push_back( root_set );
1140  entity_sets.push_back( ent_set );
1141 
1142  // get all entities from root_set and put into entity_groups
1143  std::vector< EntityHandle > entities;
1144  err = mbImpl->get_entities_by_handle( root_set, entities, true );ERRORR( "Failed to get entities in root_set.", err );
1145 
1146  entity_groups.push_back( entities );
1147 
1148  // Call do_normalization() to continue common normalization processing
1149  err = do_normalization( norm_tag, entity_sets, entity_groups, integ_type, num_integ_pts );ERRORR( "Failure in do_normalization().", err );
1150  // SLAVE END ****************************************************************
1151 
1152  return err;
1153 }

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

1163 {
1164  moab::ErrorCode err;
1165  std::vector< Tag > tag_handles;
1166 
1167  // Lookup tag handles from tag names
1168  for( int t = 0; t < num_tags; t++ )
1169  {
1170  // get tag handle & size
1171  Tag th;
1172  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 );
1173  tag_handles.push_back( th );
1174  }
1175 
1176  return normalize_subset( root_set, norm_tag, &tag_handles[0], num_tags, tag_values, integ_type, num_integ_pts );
1177 }

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

1186 {
1187  ErrorCode err;
1188 
1189  // SLAVE START ****************************************************************
1190  // Search for entities based on tag_handles and tag_values
1191  std::vector< std::vector< EntityHandle > > entity_sets;
1192  std::vector< std::vector< EntityHandle > > entity_groups;
1193 
1194  err = get_matching_entities( root_set, tag_handles, tag_values, num_tags, &entity_sets, &entity_groups );ERRORR( "Failed to get matching entities.", err );
1195 
1196  // Call do_normalization() to continue common normalization processing
1197  err = do_normalization( norm_tag, entity_sets, entity_groups, integ_type, num_integ_pts );ERRORR( "Failure in do_normalization().", err );
1198  // SLAVE END ****************************************************************
1199 
1200  return err;
1201 }

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.

466  {
467  return numIts;
468  }

References numIts.

◆ set_spherical()

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

Definition at line 470 of file Coupler.hpp.

471  {
472  spherical = arg1;
473  }

References spherical.

◆ target_ents()

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

Definition at line 449 of file Coupler.hpp.

450  {
451  return targetEnts;
452  }

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

572 {
573  std::vector< EntityHandle > entities;
574  std::vector< CartVect > nat_coords;
575  bool canWrite = false;
576  if( tl )
577  {
578  canWrite = tl->get_writeEnabled();
579  if( !canWrite )
580  {
581  tl->enableWriteAccess();
582  canWrite = true;
583  }
584  }
585 
586  if( rel_eps && !abs_eps )
587  {
588  // Relative epsilon given, translate to absolute epsilon using box dimensions
589  BoundBox box;
591  abs_eps = rel_eps * box.diagonal_length();
592  }
593 
594  ErrorCode result = nat_param( xyz, entities, nat_coords, abs_eps );
595  if( MB_SUCCESS != result ) return result;
596 
597  // If we didn't find any ents and we're looking locally, nothing more to do
598  if( entities.empty() )
599  {
600  if( tl->get_n() == tl->get_max() ) tl->resize( std::max( 10.0, 1.5 * tl->get_max() ) );
601 
602  tl->vi_wr[3 * tl->get_n()] = from_proc;
603  tl->vi_wr[3 * tl->get_n() + 1] = remote_index;
604  tl->vi_wr[3 * tl->get_n() + 2] = -1;
605  tl->inc_n();
606 
607  point_located = false;
608  return MB_SUCCESS;
609  }
610 
611  // Grow if we know we'll exceed size
612  if( mappedPts->get_n() + entities.size() >= mappedPts->get_max() )
613  mappedPts->resize( std::max( 10.0, 1.5 * mappedPts->get_max() ) );
614 
615  std::vector< EntityHandle >::iterator eit = entities.begin();
616  std::vector< CartVect >::iterator ncit = nat_coords.begin();
617 
619  for( ; eit != entities.end(); ++eit, ++ncit )
620  {
621  // Store in tuple mappedPts
622  mappedPts->vr_wr[3 * mappedPts->get_n()] = ( *ncit )[0];
623  mappedPts->vr_wr[3 * mappedPts->get_n() + 1] = ( *ncit )[1];
624  mappedPts->vr_wr[3 * mappedPts->get_n() + 2] = ( *ncit )[2];
625  mappedPts->vul_wr[mappedPts->get_n()] = *eit;
626  mappedPts->inc_n();
627 
628  // Also store local point, mapped point indices
629  if( tl->get_n() == tl->get_max() ) tl->resize( std::max( 10.0, 1.5 * tl->get_max() ) );
630 
631  // Store in tuple source_pts
632  tl->vi_wr[3 * tl->get_n()] = from_proc;
633  tl->vi_wr[3 * tl->get_n() + 1] = remote_index;
634  tl->vi_wr[3 * tl->get_n() + 2] = mappedPts->get_n() - 1;
635  tl->inc_n();
636  }
637 
638  point_located = true;
639 
640  if( tl && !canWrite ) tl->disableWriteAccess();
641 
642  return MB_SUCCESS;
643 }

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

◆ _spectralSource

void* moab::Coupler::_spectralSource
private

◆ _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

◆ myId

int moab::Coupler::myId
private

Definition at line 518 of file Coupler.hpp.

Referenced by my_id().

◆ myPc

ParallelComm* moab::Coupler::myPc
private

◆ 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: