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

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

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.

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

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.

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

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.

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

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.

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

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.

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

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.

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

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.

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

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.

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

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.

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

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

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

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.

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

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.

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

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 #ifdef VERBOSE
421  int num_to_me = 0;
422  for( unsigned int i = 0; i < target_pts.get_n(); i++ )
423  if( target_pts.vi_rd[2 * i] == (int)my_rank ) num_to_me++;
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 #ifdef VERBOSE
433  num_to_me = 0;
434  for( unsigned int i = 0; i < target_pts.get_n(); i++ )
435  {
436  if( target_pts.vi_rd[2 * i] == (int)my_rank ) num_to_me++;
437  }
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  printf( "rank: %u point location: wanted %u got %u locally, %u remote, missing %u\n", my_rank, num_points,
540  local_pts, num_points - missing_pts - local_pts, missing_pts );
541  assert( 0 == missing_pts ); // Will likely break on curved geometries
542 
543  // No longer need source_pts
544  source_pts.reset();
545 
546  // Copy into tl if passed in and storing locally
547  if( tl && store_local )
548  {
549  tl->initialize( 3, 0, 0, 0, num_points );
550  tl->enableWriteAccess();
551  memcpy( tl->vi_wr, tl_tmp->vi_rd, 3 * tl_tmp->get_n() * sizeof( int ) );
552  tl->set_n( tl_tmp->get_n() );
553  tl->disableWriteAccess();
554  }
555 
556  tl_tmp->disableWriteAccess();
557 
558  // Done
559  return MB_SUCCESS;
560 }

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

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

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

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

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.

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

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.

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

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

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

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: