 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
This class couples data between meshes. More...
#include <Coupler.hpp>
Collaboration diagram for moab::Coupler:Public Types | |
| enum | Method { CONSTANT , LINEAR_FE , QUADRATIC_FE , SPECTRAL , SPHERICAL } |
| enum | IntegType { VOLUME } |
Public Member Functions | |
| Coupler (Interface *impl, ParallelComm *pc, Range &local_elems, int coupler_id, bool init_tree=true, int max_ent_dim=3) | |
| virtual | ~Coupler () |
| ErrorCode | initialize_tree () |
| Initialize the kdtree, locally and across communicator. More... | |
| ErrorCode | locate_points (double *xyz, unsigned int num_points, double rel_eps=0.0, double abs_eps=0.0, TupleList *tl=NULL, bool store_local=true) |
| ErrorCode | locate_points (Range &ents, double rel_eps=0.0, double abs_eps=0.0, TupleList *tl=NULL, bool store_local=true) |
| ErrorCode | interpolate (Coupler::Method method, Tag tag, double *interp_vals, TupleList *tl=NULL, bool normalize=true) |
| ErrorCode | interpolate (Coupler::Method method, const std::string &tag_name, double *interp_vals, TupleList *tl=NULL, bool normalize=true) |
| ErrorCode | interpolate (Coupler::Method *methods, const std::string *tag_names, int *points_per_method, int num_methods, double *interp_vals, TupleList *tl=NULL, bool normalize=true) |
| ErrorCode | interpolate (Coupler::Method *methods, Tag *tag_names, int *points_per_method, int num_methods, double *interp_vals, TupleList *tl=NULL, bool normalize=true) |
| ErrorCode | normalize_mesh (EntityHandle root_set, const char *norm_tag, Coupler::IntegType integ_type, int num_integ_pts) |
| ErrorCode | normalize_subset (EntityHandle root_set, const char *norm_tag, const char **tag_names, int num_tags, const char **tag_values, Coupler::IntegType integ_type, int num_integ_pts) |
| ErrorCode | normalize_subset (EntityHandle root_set, const char *norm_tag, Tag *tag_handles, int num_tags, const char **tag_values, Coupler::IntegType integ_type, int num_integ_pts) |
| ErrorCode | do_normalization (const char *norm_tag, std::vector< std::vector< EntityHandle > > &entity_sets, std::vector< std::vector< EntityHandle > > &entity_groups, Coupler::IntegType integ_type, int num_integ_pts) |
| ErrorCode | get_matching_entities (EntityHandle root_set, const char **tag_names, const char **tag_values, int num_tags, std::vector< std::vector< EntityHandle > > *entity_sets, std::vector< std::vector< EntityHandle > > *entity_groups) |
| ErrorCode | get_matching_entities (EntityHandle root_set, Tag *tag_handles, const char **tag_values, int num_tags, std::vector< std::vector< EntityHandle > > *entity_sets, std::vector< std::vector< EntityHandle > > *entity_groups) |
| ErrorCode | create_tuples (Range &ent_sets, const char **tag_names, unsigned int num_tags, TupleList **tuples) |
| ErrorCode | create_tuples (Range &ent_sets, Tag *tag_handles, unsigned int num_tags, TupleList **tuples) |
| ErrorCode | consolidate_tuples (TupleList **all_tuples, unsigned int num_tuples, TupleList **unique_tuples) |
| ErrorCode | get_group_integ_vals (std::vector< std::vector< EntityHandle > > &groups, std::vector< double > &integ_vals, const char *norm_tag, int num_integ_pts, Coupler::IntegType integ_type) |
| ErrorCode | apply_group_norm_factor (std::vector< std::vector< EntityHandle > > &entity_sets, std::vector< double > &norm_factors, const char *norm_tag, Coupler::IntegType integ_type) |
| ErrorCode | initialize_spectral_elements (EntityHandle rootSource, EntityHandle rootTarget, bool &specSou, bool &specTar) |
| ErrorCode | get_gl_points_on_elements (Range &targ_elems, std::vector< double > &vpos, int &numPointsOfInterest) |
| Interface * | mb_impl () const |
| AdaptiveKDTree * | my_tree () const |
| EntityHandle | local_root () const |
| const std::vector< double > & | all_boxes () const |
| ParallelComm * | my_pc () const |
| const Range & | target_ents () const |
| int | my_id () const |
| const Range & | my_range () const |
| TupleList * | mapped_pts () const |
| int | num_its () const |
| void | set_spherical (bool arg1=true) |
Private Member Functions | |
| ErrorCode | nat_param (double xyz[3], std::vector< EntityHandle > &entities, std::vector< CartVect > &nat_coords, double epsilon=0.0) |
| ErrorCode | interp_field (EntityHandle elem, CartVect nat_coord, Tag tag, double &field) |
| ErrorCode | constant_interp (EntityHandle elem, Tag tag, double &field) |
| ErrorCode | test_local_box (double *xyz, int from_proc, int remote_index, int index, bool &point_located, double rel_eps=0.0, double abs_eps=0.0, TupleList *tl=NULL) |
Private Attributes | |
| Interface * | mbImpl |
| AdaptiveKDTree * | myTree |
| EntityHandle | localRoot |
| std::vector< double > | allBoxes |
| ParallelComm * | myPc |
| int | myId |
| Range | myRange |
| Range | targetEnts |
| TupleList * | mappedPts |
| TupleList * | targetPts |
| int | numIts |
| int | max_dim |
| void * | _spectralSource |
| void * | _spectralTarget |
| moab::Tag | _xm1Tag |
| moab::Tag | _ym1Tag |
| moab::Tag | _zm1Tag |
| int | _ntot |
| bool | spherical |
This class couples data between meshes.
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:
Multiple interpolations can be done after locating the points.
Definition at line 43 of file Coupler.hpp.
| Enumerator | |
|---|---|
| CONSTANT | |
| LINEAR_FE | |
| QUADRATIC_FE | |
| SPECTRAL | |
| SPHERICAL | |
Definition at line 46 of file Coupler.hpp.
47 { 48 CONSTANT, 49 LINEAR_FE, 50 QUADRATIC_FE, 51 SPECTRAL, 52 SPHERICAL 53 };
| 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;
68 _spectralSource = _spectralTarget = NULL;
69 }
References _spectralSource, _spectralTarget, moab::Range::empty(), initialize_tree(), mappedPts, myRange, myTree, and targetPts.
|
virtual |
Definition at line 73 of file Coupler.cpp.
74 {
75 // This will clear the cache
76 delete(moab::Element::SpectralHex*)_spectralSource;
77 delete(moab::Element::SpectralHex*)_spectralTarget;
78 delete myTree;
79 delete targetPts;
80 delete mappedPts;
81 }
References _spectralSource, _spectralTarget, mappedPts, myTree, and targetPts.
|
inline |
| 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().
| 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().
|
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().
| 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().
| 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.
| 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().
| 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().
| 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 }
1764 catch( moab::Element::Map::ArgError& )
1765 {
1766 std::cerr << "Failed to set vertices on Element::Map." << std::endl;
1767 }
1768 catch( moab::Element::Map::EvaluationError& )
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().
| 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().
| 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,
1382 MASTER_PROC, myPc->proc_config().proc_comm() );
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.
| 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];
256 rval = mbImpl->tag_get_handle( "SEM_X", _ntot, moab::MB_TYPE_DOUBLE, _xm1Tag );
257 if( moab::MB_SUCCESS != rval )
258 {
259 std::cout << "Can't get xm1tag \n";
260 return MB_FAILURE;
261 }
262 rval = mbImpl->tag_get_handle( "SEM_Y", _ntot, moab::MB_TYPE_DOUBLE, _ym1Tag );
263 if( moab::MB_SUCCESS != rval )
264 {
265 std::cout << "Can't get ym1tag \n";
266 return MB_FAILURE;
267 }
268 rval = mbImpl->tag_get_handle( "SEM_Z", _ntot, moab::MB_TYPE_DOUBLE, _zm1Tag );
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().
| 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().
|
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 }
1095 catch( moab::Element::Map::EvaluationError& )
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().
| ErrorCode moab::Coupler::interpolate | ( | Coupler::Method * | methods, |
| const std::string * | tag_names, | ||
| int * | points_per_method, | ||
| int | num_methods, | ||
| double * | interp_vals, | ||
| TupleList * | tl = NULL, |
||
| bool | normalize = true |
||
| ) |
| ErrorCode moab::Coupler::interpolate | ( | Coupler::Method * | methods, |
| Tag * | tag_names, | ||
| int * | points_per_method, | ||
| int | num_methods, | ||
| double * | interp_vals, | ||
| TupleList * | tl = NULL, |
||
| bool | normalize = true |
||
| ) |
Definition at line 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
714 mappedPts->enableWriteAccess();
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.
| 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().
|
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().
|
inline |
| 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() );
339 mappedPts->enableWriteAccess();
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().
| 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.
|
inline |
|
inline |
|
inline |
|
inline |
|
inline |
|
inline |
|
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().
| 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.
| 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().
| 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().
|
inline |
|
inline |
|
inline |
Definition at line 449 of file Coupler.hpp.
450 {
451 return targetEnts;
452 }
References targetEnts.
|
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;
588 myTree->get_bounding_box( box, &localRoot );
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
616 mappedPts->enableWriteAccess();
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().
|
private |
Definition at line 559 of file Coupler.hpp.
Referenced by get_gl_points_on_elements(), initialize_spectral_elements(), and interpolate().
|
private |
Definition at line 556 of file Coupler.hpp.
Referenced by Coupler(), initialize_spectral_elements(), interp_field(), interpolate(), nat_param(), and ~Coupler().
|
private |
Definition at line 557 of file Coupler.hpp.
Referenced by Coupler(), initialize_spectral_elements(), and ~Coupler().
|
private |
Definition at line 558 of file Coupler.hpp.
Referenced by get_gl_points_on_elements(), initialize_spectral_elements(), and nat_param().
|
private |
Definition at line 558 of file Coupler.hpp.
Referenced by get_gl_points_on_elements(), initialize_spectral_elements(), and nat_param().
|
private |
Definition at line 558 of file Coupler.hpp.
Referenced by get_gl_points_on_elements(), initialize_spectral_elements(), and nat_param().
|
private |
Definition at line 510 of file Coupler.hpp.
Referenced by all_boxes(), initialize_tree(), and locate_points().
|
private |
Definition at line 506 of file Coupler.hpp.
Referenced by initialize_tree(), local_root(), nat_param(), and test_local_box().
|
private |
Definition at line 534 of file Coupler.hpp.
Referenced by Coupler(), interpolate(), locate_points(), mapped_pts(), test_local_box(), and ~Coupler().
|
private |
Definition at line 551 of file Coupler.hpp.
Referenced by initialize_tree(), and nat_param().
|
private |
Definition at line 498 of file Coupler.hpp.
Referenced by apply_group_norm_factor(), constant_interp(), create_tuples(), get_gl_points_on_elements(), get_group_integ_vals(), get_matching_entities(), initialize_spectral_elements(), initialize_tree(), interp_field(), interpolate(), locate_points(), mb_impl(), nat_param(), normalize_mesh(), and normalize_subset().
|
private |
Definition at line 518 of file Coupler.hpp.
Referenced by my_id().
|
private |
Definition at line 514 of file Coupler.hpp.
Referenced by do_normalization(), get_matching_entities(), initialize_tree(), interpolate(), locate_points(), and my_pc().
|
private |
Definition at line 522 of file Coupler.hpp.
Referenced by Coupler(), initialize_tree(), and my_range().
|
private |
Definition at line 502 of file Coupler.hpp.
Referenced by Coupler(), initialize_tree(), my_tree(), nat_param(), test_local_box(), and ~Coupler().
|
private |
Definition at line 548 of file Coupler.hpp.
Referenced by initialize_tree(), and num_its().
|
private |
Definition at line 562 of file Coupler.hpp.
Referenced by initialize_tree(), nat_param(), and set_spherical().
|
private |
Definition at line 526 of file Coupler.hpp.
Referenced by locate_points(), and target_ents().
|
private |
Definition at line 543 of file Coupler.hpp.
Referenced by Coupler(), interpolate(), locate_points(), and ~Coupler().
1.9.1.