MOAB: Mesh Oriented datABase
(version 5.3.1)
|
00001 #include "moab/ScdInterface.hpp" 00002 #include "moab/Core.hpp" 00003 #include "SequenceManager.hpp" 00004 #include "EntitySequence.hpp" 00005 #include "StructuredElementSeq.hpp" 00006 #include "VertexSequence.hpp" 00007 #include "ScdVertexData.hpp" 00008 #include "MBTagConventions.hpp" 00009 #ifdef MOAB_HAVE_MPI 00010 #include "moab/ParallelComm.hpp" 00011 #include "moab/TupleList.hpp" 00012 #include "moab/gs.hpp" 00013 #endif 00014 #include <cassert> 00015 #include <iostream> 00016 #include <functional> 00017 00018 #define ERRORR( rval, str ) \ 00019 { \ 00020 if( MB_SUCCESS != ( rval ) ) \ 00021 { \ 00022 std::cerr << ( str ); \ 00023 return rval; \ 00024 } \ 00025 } 00026 00027 namespace moab 00028 { 00029 00030 const char* ScdParData::PartitionMethodNames[] = { "alljorkori", "alljkbal", "sqij", "sqjk", 00031 "sqijk", "trivial", "rcbzoltan", "nopart" }; 00032 00033 ScdInterface::ScdInterface( Interface* imp, bool boxes ) 00034 : mbImpl( imp ), searchedBoxes( false ), boxPeriodicTag( 0 ), boxDimsTag( 0 ), globalBoxDimsTag( 0 ), 00035 partMethodTag( 0 ), boxSetTag( 0 ) 00036 { 00037 if( boxes ) find_boxes( scdBoxes ); 00038 } 00039 00040 // Destructor 00041 ScdInterface::~ScdInterface() 00042 { 00043 std::vector< ScdBox* > tmp_boxes; 00044 tmp_boxes.swap( scdBoxes ); 00045 00046 for( std::vector< ScdBox* >::iterator rit = tmp_boxes.begin(); rit != tmp_boxes.end(); ++rit ) 00047 delete *rit; 00048 00049 if( box_set_tag( false ) ) mbImpl->tag_delete( box_set_tag() ); 00050 } 00051 00052 Interface* ScdInterface::impl() const 00053 { 00054 return mbImpl; 00055 } 00056 00057 ErrorCode ScdInterface::find_boxes( std::vector< ScdBox* >& scd_boxes ) 00058 { 00059 Range tmp_boxes; 00060 ErrorCode rval = find_boxes( tmp_boxes ); 00061 if( MB_SUCCESS != rval ) return rval; 00062 00063 for( Range::iterator rit = tmp_boxes.begin(); rit != tmp_boxes.end(); ++rit ) 00064 { 00065 ScdBox* tmp_box = get_scd_box( *rit ); 00066 if( tmp_box ) 00067 scd_boxes.push_back( tmp_box ); 00068 else 00069 rval = MB_FAILURE; 00070 } 00071 00072 return rval; 00073 } 00074 00075 ErrorCode ScdInterface::find_boxes( Range& scd_boxes ) 00076 { 00077 ErrorCode rval = MB_SUCCESS; 00078 box_dims_tag(); 00079 Range boxes; 00080 if( !searchedBoxes ) 00081 { 00082 rval = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &boxDimsTag, NULL, 1, boxes, Interface::UNION ); 00083 searchedBoxes = true; 00084 if( !boxes.empty() ) 00085 { 00086 scdBoxes.resize( boxes.size() ); 00087 rval = mbImpl->tag_get_data( boxSetTag, boxes, &scdBoxes[0] ); 00088 ScdBox* dum = NULL; 00089 // std::remove_if(scdBoxes.begin(), scdBoxes.end(), 00090 // std::bind2nd(std::equal_to<ScdBox*>(), dum) ) ; 00091 std::remove_if( scdBoxes.begin(), scdBoxes.end(), 00092 std::bind( std::equal_to< ScdBox* >(), std::placeholders::_1, dum ) ); 00093 // https://stackoverflow.com/questions/32739018/a-replacement-for-stdbind2nd 00094 } 00095 } 00096 00097 for( std::vector< ScdBox* >::iterator vit = scdBoxes.begin(); vit != scdBoxes.end(); ++vit ) 00098 scd_boxes.insert( ( *vit )->box_set() ); 00099 00100 return rval; 00101 } 00102 00103 ScdBox* ScdInterface::get_scd_box( EntityHandle eh ) 00104 { 00105 ScdBox* scd_box = NULL; 00106 if( !box_set_tag( false ) ) return scd_box; 00107 00108 mbImpl->tag_get_data( box_set_tag(), &eh, 1, &scd_box ); 00109 return scd_box; 00110 } 00111 00112 ErrorCode ScdInterface::construct_box( HomCoord low, HomCoord high, const double* const coords, unsigned int num_coords, 00113 ScdBox*& new_box, int* const lperiodic, ScdParData* par_data, bool assign_gids, 00114 int tag_shared_ents ) 00115 { 00116 // create a rectangular structured mesh block 00117 ErrorCode rval; 00118 00119 int tmp_lper[3] = { 0, 0, 0 }; 00120 if( lperiodic ) std::copy( lperiodic, lperiodic + 3, tmp_lper ); 00121 00122 #ifndef MOAB_HAVE_MPI 00123 if( -1 != tag_shared_ents ) ERRORR( MB_FAILURE, "Parallel capability requested but MOAB not compiled parallel." ); 00124 if( -1 == tag_shared_ents && !assign_gids ) assign_gids = true; // need to assign gids in order to tag shared verts 00125 #else 00126 if( par_data && low == high && ScdParData::NOPART != par_data->partMethod ) 00127 { 00128 // requesting creation of parallel mesh, so need to compute partition 00129 if( !par_data->pComm ) 00130 { 00131 // this is a really boneheaded way to have to create a PC 00132 par_data->pComm = ParallelComm::get_pcomm( mbImpl, 0 ); 00133 if( NULL == par_data->pComm ) par_data->pComm = new ParallelComm( mbImpl, MPI_COMM_WORLD ); 00134 } 00135 int ldims[6]; 00136 rval = compute_partition( par_data->pComm->size(), par_data->pComm->rank(), *par_data, ldims, tmp_lper, 00137 par_data->pDims );ERRORR( rval, "Error returned from compute_partition." ); 00138 low.set( ldims[0], ldims[1], ldims[2] ); 00139 high.set( ldims[3], ldims[4], ldims[5] ); 00140 if( par_data->pComm->get_debug_verbosity() > 0 ) 00141 { 00142 std::cout << "Proc " << par_data->pComm->rank() << ": " << *par_data; 00143 std::cout << "Proc " << par_data->pComm->rank() << " local dims: " << low << "-" << high << std::endl; 00144 } 00145 } 00146 #endif 00147 00148 HomCoord tmp_size = high - low + HomCoord( 1, 1, 1, 0 ); 00149 if( ( tmp_size[1] && num_coords && (int)num_coords < tmp_size[0] ) || 00150 ( tmp_size[2] && num_coords && (int)num_coords < tmp_size[0] * tmp_size[1] ) ) 00151 return MB_FAILURE; 00152 00153 rval = create_scd_sequence( low, high, MBVERTEX, 0, new_box );ERRORR( rval, "Trouble creating scd vertex sequence." ); 00154 00155 // set the vertex coordinates 00156 double *xc, *yc, *zc; 00157 rval = new_box->get_coordinate_arrays( xc, yc, zc );ERRORR( rval, "Couldn't get vertex coordinate arrays." ); 00158 00159 if( coords && num_coords ) 00160 { 00161 unsigned int i = 0; 00162 for( int kl = low[2]; kl <= high[2]; kl++ ) 00163 { 00164 for( int jl = low[1]; jl <= high[1]; jl++ ) 00165 { 00166 for( int il = low[0]; il <= high[0]; il++ ) 00167 { 00168 xc[i] = coords[3 * i]; 00169 if( new_box->box_size()[1] ) yc[i] = coords[3 * i + 1]; 00170 if( new_box->box_size()[2] ) zc[i] = coords[3 * i + 2]; 00171 i++; 00172 } 00173 } 00174 } 00175 } 00176 else 00177 { 00178 unsigned int i = 0; 00179 for( int kl = low[2]; kl <= high[2]; kl++ ) 00180 { 00181 for( int jl = low[1]; jl <= high[1]; jl++ ) 00182 { 00183 for( int il = low[0]; il <= high[0]; il++ ) 00184 { 00185 xc[i] = (double)il; 00186 if( new_box->box_size()[1] ) 00187 yc[i] = (double)jl; 00188 else 00189 yc[i] = 0.0; 00190 if( new_box->box_size()[2] ) 00191 zc[i] = (double)kl; 00192 else 00193 zc[i] = 0.0; 00194 i++; 00195 } 00196 } 00197 } 00198 } 00199 00200 // create element sequence 00201 Core* mbcore = dynamic_cast< Core* >( mbImpl ); 00202 SequenceManager* seq_mgr = mbcore->sequence_manager(); 00203 00204 EntitySequence* tmp_seq; 00205 EntityHandle start_ent; 00206 00207 // construct the sequence 00208 EntityType this_tp = MBHEX; 00209 if( 1 >= tmp_size[2] ) this_tp = MBQUAD; 00210 if( 1 >= tmp_size[2] && 1 >= tmp_size[1] ) this_tp = MBEDGE; 00211 rval = seq_mgr->create_scd_sequence( low, high, this_tp, 0, start_ent, tmp_seq, tmp_lper );ERRORR( rval, "Trouble creating scd element sequence." ); 00212 00213 new_box->elem_seq( tmp_seq ); 00214 new_box->start_element( start_ent ); 00215 00216 // add vertex seq to element seq, forward orientation, unity transform 00217 rval = new_box->add_vbox( new_box, 00218 // p1: imin,jmin 00219 low, low, 00220 // p2: imax,jmin 00221 low + HomCoord( 1, 0, 0 ), low + HomCoord( 1, 0, 0 ), 00222 // p3: imin,jmax 00223 low + HomCoord( 0, 1, 0 ), low + HomCoord( 0, 1, 0 ) );ERRORR( rval, "Error constructing structured element sequence." ); 00224 00225 // add the new hexes to the scd box set; vertices were added in call to create_scd_sequence 00226 Range tmp_range( new_box->start_element(), new_box->start_element() + new_box->num_elements() - 1 ); 00227 rval = mbImpl->add_entities( new_box->box_set(), tmp_range );ERRORR( rval, "Couldn't add new hexes to box set." ); 00228 00229 if( par_data ) new_box->par_data( *par_data ); 00230 00231 if( assign_gids ) 00232 { 00233 rval = assign_global_ids( new_box );ERRORR( rval, "Trouble assigning global ids" ); 00234 } 00235 00236 #ifdef MOAB_HAVE_MPI 00237 if( par_data && -1 != tag_shared_ents ) { rval = tag_shared_vertices( par_data->pComm, new_box ); } 00238 #endif 00239 00240 return MB_SUCCESS; 00241 } 00242 00243 ErrorCode ScdInterface::assign_global_ids( ScdBox* box ) 00244 { 00245 // Get a ptr to global id memory 00246 void* data; 00247 int count = 0; 00248 Tag gid_tag = mbImpl->globalId_tag(); 00249 Range tmp_range( box->start_vertex(), box->start_vertex() + box->num_vertices() ); 00250 ErrorCode rval = mbImpl->tag_iterate( gid_tag, tmp_range.begin(), tmp_range.end(), count, data );ERRORR( rval, "Failed to get tag iterator." ); 00251 assert( count == box->num_vertices() ); 00252 int* gid_data = (int*)data; 00253 int di = box->par_data().gDims[3] - box->par_data().gDims[0] + 1; 00254 int dj = box->par_data().gDims[4] - box->par_data().gDims[1] + 1; 00255 00256 for( int kl = box->box_dims()[2]; kl <= box->box_dims()[5]; kl++ ) 00257 { 00258 for( int jl = box->box_dims()[1]; jl <= box->box_dims()[4]; jl++ ) 00259 { 00260 for( int il = box->box_dims()[0]; il <= box->box_dims()[3]; il++ ) 00261 { 00262 int itmp = 00263 ( !box->locally_periodic()[0] && box->par_data().gPeriodic[0] && il == box->par_data().gDims[3] 00264 ? box->par_data().gDims[0] 00265 : il ); 00266 *gid_data = ( -1 != kl ? kl * di * dj : 0 ) + jl * di + itmp + 1; 00267 gid_data++; 00268 } 00269 } 00270 } 00271 00272 return MB_SUCCESS; 00273 } 00274 00275 ErrorCode ScdInterface::create_scd_sequence( const HomCoord& low, const HomCoord& high, EntityType tp, int starting_id, 00276 ScdBox*& new_box, int* is_periodic ) 00277 { 00278 HomCoord tmp_size = high - low + HomCoord( 1, 1, 1, 0 ); 00279 if( ( tp == MBHEX && 1 >= tmp_size[2] ) || ( tp == MBQUAD && 1 >= tmp_size[1] ) || 00280 ( tp == MBEDGE && 1 >= tmp_size[0] ) ) 00281 return MB_TYPE_OUT_OF_RANGE; 00282 00283 Core* mbcore = dynamic_cast< Core* >( mbImpl ); 00284 assert( mbcore != NULL ); 00285 SequenceManager* seq_mgr = mbcore->sequence_manager(); 00286 00287 EntitySequence* tmp_seq; 00288 EntityHandle start_ent, scd_set; 00289 00290 // construct the sequence 00291 ErrorCode rval = seq_mgr->create_scd_sequence( low, high, tp, starting_id, start_ent, tmp_seq, is_periodic ); 00292 if( MB_SUCCESS != rval ) return rval; 00293 00294 // create the set for this rectangle 00295 rval = create_box_set( low, high, scd_set ); 00296 if( MB_SUCCESS != rval ) return rval; 00297 00298 // make the ScdBox 00299 new_box = new ScdBox( this, scd_set, tmp_seq ); 00300 if( !new_box ) return MB_FAILURE; 00301 00302 // set the start vertex/element 00303 Range new_range; 00304 if( MBVERTEX == tp ) { new_range.insert( start_ent, start_ent + new_box->num_vertices() - 1 ); } 00305 else 00306 { 00307 new_range.insert( start_ent, start_ent + new_box->num_elements() - 1 ); 00308 } 00309 00310 // put the entities in the box set 00311 rval = mbImpl->add_entities( scd_set, new_range ); 00312 if( MB_SUCCESS != rval ) return rval; 00313 00314 // tag the set with the box 00315 rval = mbImpl->tag_set_data( box_set_tag(), &scd_set, 1, &new_box ); 00316 if( MB_SUCCESS != rval ) return rval; 00317 00318 return MB_SUCCESS; 00319 } 00320 00321 ErrorCode ScdInterface::create_box_set( const HomCoord& low, const HomCoord& high, EntityHandle& scd_set, 00322 int* is_periodic ) 00323 { 00324 // create the set and put the entities in it 00325 ErrorCode rval = mbImpl->create_meshset( MESHSET_SET, scd_set ); 00326 if( MB_SUCCESS != rval ) return rval; 00327 00328 // tag the set with parameter extents 00329 int boxdims[6]; 00330 for( int i = 0; i < 3; i++ ) 00331 boxdims[i] = low[i]; 00332 for( int i = 0; i < 3; i++ ) 00333 boxdims[3 + i] = high[i]; 00334 rval = mbImpl->tag_set_data( box_dims_tag(), &scd_set, 1, boxdims ); 00335 if( MB_SUCCESS != rval ) return rval; 00336 00337 if( is_periodic ) 00338 { 00339 rval = mbImpl->tag_set_data( box_periodic_tag(), &scd_set, 1, is_periodic ); 00340 if( MB_SUCCESS != rval ) return rval; 00341 } 00342 00343 return rval; 00344 } 00345 00346 Tag ScdInterface::box_periodic_tag( bool create_if_missing ) 00347 { 00348 // Reset boxPeriodicTag in case it has been deleted (e.g. by Core::clean_up_failed_read) 00349 if( boxPeriodicTag ) 00350 { 00351 std::string tag_name; 00352 if( MB_TAG_NOT_FOUND == mbImpl->tag_get_name( boxPeriodicTag, tag_name ) ) boxPeriodicTag = NULL; 00353 } 00354 00355 if( boxPeriodicTag || !create_if_missing ) return boxPeriodicTag; 00356 00357 ErrorCode rval = 00358 mbImpl->tag_get_handle( "BOX_PERIODIC", 3, MB_TYPE_INTEGER, boxPeriodicTag, MB_TAG_SPARSE | MB_TAG_CREAT ); 00359 if( MB_SUCCESS != rval ) return 0; 00360 return boxPeriodicTag; 00361 } 00362 00363 Tag ScdInterface::box_dims_tag( bool create_if_missing ) 00364 { 00365 // Reset boxDimsTag in case it has been deleted (e.g. by clean_up_failed_read) 00366 if( boxDimsTag ) 00367 { 00368 std::string tag_name; 00369 if( MB_TAG_NOT_FOUND == mbImpl->tag_get_name( boxDimsTag, tag_name ) ) boxDimsTag = NULL; 00370 } 00371 00372 if( boxDimsTag || !create_if_missing ) return boxDimsTag; 00373 00374 ErrorCode rval = mbImpl->tag_get_handle( "BOX_DIMS", 6, MB_TYPE_INTEGER, boxDimsTag, MB_TAG_SPARSE | MB_TAG_CREAT ); 00375 if( MB_SUCCESS != rval ) return 0; 00376 return boxDimsTag; 00377 } 00378 00379 Tag ScdInterface::global_box_dims_tag( bool create_if_missing ) 00380 { 00381 // Reset globalBoxDimsTag in case it has been deleted (e.g. by Core::clean_up_failed_read) 00382 if( globalBoxDimsTag ) 00383 { 00384 std::string tag_name; 00385 if( MB_TAG_NOT_FOUND == mbImpl->tag_get_name( globalBoxDimsTag, tag_name ) ) globalBoxDimsTag = NULL; 00386 } 00387 00388 if( globalBoxDimsTag || !create_if_missing ) return globalBoxDimsTag; 00389 00390 ErrorCode rval = 00391 mbImpl->tag_get_handle( "GLOBAL_BOX_DIMS", 6, MB_TYPE_INTEGER, globalBoxDimsTag, MB_TAG_SPARSE | MB_TAG_CREAT ); 00392 if( MB_SUCCESS != rval ) return 0; 00393 return globalBoxDimsTag; 00394 } 00395 00396 Tag ScdInterface::part_method_tag( bool create_if_missing ) 00397 { 00398 // Reset partMethodTag in case it has been deleted (e.g. by Core::clean_up_failed_read) 00399 if( partMethodTag ) 00400 { 00401 std::string tag_name; 00402 if( MB_TAG_NOT_FOUND == mbImpl->tag_get_name( partMethodTag, tag_name ) ) partMethodTag = NULL; 00403 } 00404 00405 if( partMethodTag || !create_if_missing ) return partMethodTag; 00406 00407 ErrorCode rval = 00408 mbImpl->tag_get_handle( "PARTITION_METHOD", 1, MB_TYPE_INTEGER, partMethodTag, MB_TAG_SPARSE | MB_TAG_CREAT ); 00409 if( MB_SUCCESS != rval ) return 0; 00410 return partMethodTag; 00411 } 00412 00413 Tag ScdInterface::box_set_tag( bool create_if_missing ) 00414 { 00415 // Reset boxSetTag in case it has been deleted (e.g. by Core::clean_up_failed_read) 00416 if( boxSetTag ) 00417 { 00418 std::string tag_name; 00419 if( MB_TAG_NOT_FOUND == mbImpl->tag_get_name( boxSetTag, tag_name ) ) boxSetTag = NULL; 00420 } 00421 00422 if( boxSetTag || !create_if_missing ) return boxSetTag; 00423 00424 ErrorCode rval = mbImpl->tag_get_handle( "__BOX_SET", sizeof( ScdBox* ), MB_TYPE_OPAQUE, boxSetTag, 00425 MB_TAG_SPARSE | MB_TAG_CREAT ); 00426 if( MB_SUCCESS != rval ) return 0; 00427 return boxSetTag; 00428 } 00429 00430 //! Remove the box from the list on ScdInterface 00431 ErrorCode ScdInterface::remove_box( ScdBox* box ) 00432 { 00433 std::vector< ScdBox* >::iterator vit = std::find( scdBoxes.begin(), scdBoxes.end(), box ); 00434 if( vit != scdBoxes.end() ) 00435 { 00436 scdBoxes.erase( vit ); 00437 return MB_SUCCESS; 00438 } 00439 else 00440 return MB_FAILURE; 00441 } 00442 00443 //! Add the box to the list on ScdInterface 00444 ErrorCode ScdInterface::add_box( ScdBox* box ) 00445 { 00446 scdBoxes.push_back( box ); 00447 return MB_SUCCESS; 00448 } 00449 00450 ErrorCode ScdInterface::get_boxes( std::vector< ScdBox* >& boxes ) 00451 { 00452 std::copy( scdBoxes.begin(), scdBoxes.end(), std::back_inserter( boxes ) ); 00453 return MB_SUCCESS; 00454 } 00455 00456 ScdBox::ScdBox( ScdInterface* impl, EntityHandle bset, EntitySequence* seq1, EntitySequence* seq2 ) 00457 : scImpl( impl ), boxSet( bset ), vertDat( NULL ), elemSeq( NULL ), startVertex( 0 ), startElem( 0 ) 00458 { 00459 for( int i = 0; i < 6; i++ ) 00460 boxDims[i] = 0; 00461 for( int i = 0; i < 3; i++ ) 00462 locallyPeriodic[i] = false; 00463 VertexSequence* vseq = dynamic_cast< VertexSequence* >( seq1 ); 00464 if( vseq ) vertDat = dynamic_cast< ScdVertexData* >( vseq->data() ); 00465 if( vertDat ) 00466 { 00467 // retrieve the parametric space 00468 for( int i = 0; i < 3; i++ ) 00469 { 00470 boxDims[i] = vertDat->min_params()[i]; 00471 boxDims[3 + i] = vertDat->max_params()[i]; 00472 } 00473 startVertex = vertDat->start_handle(); 00474 } 00475 else if( impl->boxDimsTag ) 00476 { 00477 // look for parametric space info on set 00478 ErrorCode rval = impl->mbImpl->tag_get_data( impl->boxDimsTag, &bset, 1, boxDims ); 00479 if( MB_SUCCESS == rval ) 00480 { 00481 Range verts; 00482 impl->mbImpl->get_entities_by_dimension( bset, 0, verts ); 00483 if( !verts.empty() ) startVertex = *verts.begin(); 00484 } 00485 } 00486 00487 elemSeq = dynamic_cast< StructuredElementSeq* >( seq2 ); 00488 if( !elemSeq ) elemSeq = dynamic_cast< StructuredElementSeq* >( seq1 ); 00489 00490 if( elemSeq ) 00491 { 00492 if( vertDat ) 00493 { 00494 // check the parametric space to make sure it's consistent 00495 assert( elemSeq->sdata()->min_params() == HomCoord( boxDims, 3 ) && 00496 ( elemSeq->sdata()->max_params() + HomCoord( 1, 1, 1 ) ) == HomCoord( boxDims, 3 ) ); 00497 } 00498 else 00499 { 00500 // get the parametric space from the element sequence 00501 for( int i = 0; i < 3; i++ ) 00502 { 00503 boxDims[i] = elemSeq->sdata()->min_params()[i]; 00504 boxDims[3 + i] = elemSeq->sdata()->max_params()[i]; 00505 } 00506 } 00507 00508 startElem = elemSeq->start_handle(); 00509 } 00510 else 00511 { 00512 Range elems; 00513 impl->mbImpl->get_entities_by_dimension( 00514 bset, ( boxDims[2] == boxDims[5] ? ( boxDims[1] == boxDims[4] ? 1 : 2 ) : 3 ), elems ); 00515 if( !elems.empty() ) startElem = *elems.begin(); 00516 // call the following w/o looking at return value, since it doesn't really need to be there 00517 if( impl->boxPeriodicTag ) impl->mbImpl->tag_get_data( impl->boxPeriodicTag, &bset, 1, locallyPeriodic ); 00518 } 00519 00520 assert( vertDat || elemSeq || boxDims[0] != boxDims[3] || boxDims[1] != boxDims[4] || boxDims[2] != boxDims[5] ); 00521 00522 boxSize = HomCoord( boxDims + 3, 3 ) - HomCoord( boxDims, 3 ) + HomCoord( 1, 1, 1 ); 00523 boxSizeIJ = ( boxSize[1] ? boxSize[1] : 1 ) * boxSize[0]; 00524 boxSizeIM1 = boxSize[0] - ( locallyPeriodic[0] ? 0 : 1 ); 00525 boxSizeIJM1 = ( boxSize[1] ? ( boxSize[1] - ( locallyPeriodic[1] ? 0 : 1 ) ) : 1 ) * boxSizeIM1; 00526 00527 scImpl->add_box( this ); 00528 } 00529 00530 ScdBox::~ScdBox() 00531 { 00532 // Reset the tag on the set 00533 if( boxSet ) 00534 { 00535 // It is possible that the box set entity has been deleted (e.g. by 00536 // Core::clean_up_failed_read) 00537 Core* mbcore = dynamic_cast< Core* >( scImpl->mbImpl ); 00538 assert( mbcore != NULL ); 00539 if( mbcore->is_valid( boxSet ) ) 00540 { 00541 ScdBox* tmp_ptr = NULL; 00542 scImpl->mbImpl->tag_set_data( scImpl->box_set_tag(), &boxSet, 1, &tmp_ptr ); 00543 } 00544 else 00545 boxSet = 0; 00546 } 00547 00548 scImpl->remove_box( this ); 00549 } 00550 00551 EntityHandle ScdBox::get_vertex_from_seq( int i, int j, int k ) const 00552 { 00553 assert( elemSeq ); 00554 return elemSeq->get_vertex( i, j, k ); 00555 } 00556 00557 int ScdBox::box_dimension() const 00558 { 00559 return ( startElem ? scImpl->mbImpl->dimension_from_handle( startElem ) : -1 ); 00560 } 00561 00562 ErrorCode ScdBox::add_vbox( ScdBox* vbox, HomCoord from1, HomCoord to1, HomCoord from2, HomCoord to2, HomCoord from3, 00563 HomCoord to3, bool bb_input, const HomCoord& bb_min, const HomCoord& bb_max ) 00564 { 00565 if( !vbox->vertDat ) return MB_FAILURE; 00566 ScdVertexData* dum_data = dynamic_cast< ScdVertexData* >( vbox->vertDat ); 00567 ErrorCode rval = 00568 elemSeq->sdata()->add_vsequence( dum_data, from1, to1, from2, to2, from3, to3, bb_input, bb_min, bb_max ); 00569 return rval; 00570 } 00571 00572 bool ScdBox::boundary_complete() const 00573 { 00574 return elemSeq->boundary_complete(); 00575 } 00576 00577 ErrorCode ScdBox::get_coordinate_arrays( double*& xc, double*& yc, double*& zc ) 00578 { 00579 if( !vertDat ) return MB_FAILURE; 00580 00581 xc = reinterpret_cast< double* >( vertDat->get_sequence_data( 0 ) ); 00582 yc = reinterpret_cast< double* >( vertDat->get_sequence_data( 1 ) ); 00583 zc = reinterpret_cast< double* >( vertDat->get_sequence_data( 2 ) ); 00584 return MB_SUCCESS; 00585 } 00586 00587 ErrorCode ScdBox::get_coordinate_arrays( const double*& xc, const double*& yc, const double*& zc ) const 00588 { 00589 if( !vertDat ) return MB_FAILURE; 00590 xc = reinterpret_cast< const double* >( vertDat->get_sequence_data( 0 ) ); 00591 yc = reinterpret_cast< const double* >( vertDat->get_sequence_data( 1 ) ); 00592 zc = reinterpret_cast< const double* >( vertDat->get_sequence_data( 2 ) ); 00593 return MB_SUCCESS; 00594 } 00595 00596 ErrorCode ScdBox::vert_dat( ScdVertexData* vert_dt ) 00597 { 00598 vertDat = vert_dt; 00599 return MB_SUCCESS; 00600 } 00601 00602 ErrorCode ScdBox::elem_seq( EntitySequence* elem_sq ) 00603 { 00604 elemSeq = dynamic_cast< StructuredElementSeq* >( elem_sq ); 00605 if( elemSeq ) elemSeq->is_periodic( locallyPeriodic ); 00606 00607 if( locallyPeriodic[0] ) boxSizeIM1 = boxSize[0] - ( locallyPeriodic[0] ? 0 : 1 ); 00608 if( locallyPeriodic[0] || locallyPeriodic[1] ) 00609 boxSizeIJM1 = ( boxSize[1] ? ( boxSize[1] - ( locallyPeriodic[1] ? 0 : 1 ) ) : 1 ) * boxSizeIM1; 00610 00611 return ( elemSeq ? MB_SUCCESS : MB_FAILURE ); 00612 } 00613 00614 ErrorCode ScdBox::get_params( EntityHandle ent, HomCoord& ijkd ) const 00615 { 00616 // check first whether this is an intermediate entity, so we know what to do 00617 int dimension = box_dimension(); 00618 int this_dim = scImpl->impl()->dimension_from_handle( ent ); 00619 00620 if( ( 0 == this_dim && !vertDat ) || ( this_dim && this_dim == dimension ) ) 00621 { 00622 assert( elemSeq ); 00623 return elemSeq->get_params( ent, ijkd[0], ijkd[1], ijkd[2] ); 00624 } 00625 00626 else if( !this_dim && vertDat ) 00627 return vertDat->get_params( ent, ijkd[0], ijkd[1], ijkd[2] ); 00628 00629 else 00630 return MB_NOT_IMPLEMENTED; 00631 } 00632 00633 //! Get the entity of specified dimension adjacent to parametric element 00634 /** 00635 * \param dim Dimension of adjacent entity being requested 00636 * \param i Parametric coordinates of cell being evaluated 00637 * \param j Parametric coordinates of cell being evaluated 00638 * \param k Parametric coordinates of cell being evaluated 00639 * \param dir Direction (0, 1, or 2), for getting adjacent edges (2d, 3d) or faces (3d) 00640 * \param ent EntityHandle of adjacent entity 00641 * \param create_if_missing If true, creates the entity if it doesn't already exist 00642 */ 00643 ErrorCode ScdBox::get_adj_edge_or_face( int dim, int i, int j, int k, int dir, EntityHandle& ent, 00644 bool create_if_missing ) const 00645 { 00646 // describe connectivity of sub-element in static array 00647 // subconnect[dim-1][dir][numv][ijk] where dimensions are: 00648 // [dim-1]: dim=1 or 2, so this is 0 or 1 00649 // [dir]: one of 0..2, for ijk directions in a hex 00650 // [numv]: number of vertices describing sub entity = 2*dim <= 4 00651 // [ijk]: 3 values for i, j, k 00652 int subconnect[2][3][4][3] = { { { { 0, 0, 0 }, { 1, 0, 0 }, { -1, -1, -1 }, { -1, -1, -1 } }, // i edge 00653 { { 0, 0, 0 }, { 0, 1, 0 }, { -1, -1, -1 }, { -1, -1, -1 } }, // j edge 00654 { { 0, 0, 0 }, { 0, 0, 1 }, { -1, -1, -1 }, { -1, -1, -1 } } }, // k edge 00655 00656 { { { 0, 0, 0 }, { 0, 1, 0 }, { 0, 1, 1 }, { 0, 0, 1 } }, // i face 00657 { { 0, 0, 0 }, { 1, 0, 0 }, { 1, 0, 1 }, { 0, 0, 1 } }, // j face 00658 { { 0, 0, 0 }, { 1, 0, 0 }, { 1, 1, 0 }, { 0, 1, 0 } } } }; // k face 00659 00660 // check proper input dimensions and lower bound 00661 if( dim < 1 || dim > 2 || i < boxDims[0] || j < boxDims[1] || k < boxDims[2] ) return MB_FAILURE; 00662 00663 // now check upper bound; parameters must be <= upper corner, since edges/faces 00664 // follow element parameterization, not vertex parameterization 00665 else if( ( boxDims[3] != boxDims[0] && i > ( locallyPeriodic[0] ? boxDims[3] + 1 : boxDims[3] ) ) || 00666 ( boxDims[4] != boxDims[1] && j > ( locallyPeriodic[1] ? boxDims[4] + 1 : boxDims[4] ) ) || 00667 ( boxDims[5] != boxDims[2] && k > boxDims[5] ) ) 00668 return MB_FAILURE; 00669 00670 // get the vertices making up this entity 00671 EntityHandle verts[4]; 00672 for( int ind = 0; ind < 2 * dim; ind++ ) 00673 { 00674 int i1 = i + subconnect[dim - 1][dir][ind][0]; 00675 int j1 = j + subconnect[dim - 1][dir][ind][1]; 00676 // if periodic in i and i1 is boxDims[3]+1, wrap around 00677 if( locallyPeriodic[0] && i1 == boxDims[3] + 1 ) i1 = boxDims[0]; 00678 // if periodic in j and j1 is boxDims[4]+1, wrap around 00679 if( locallyPeriodic[1] && j1 == boxDims[4] + 1 ) j1 = boxDims[1]; 00680 verts[ind] = get_vertex( i1, j1, k + subconnect[dim - 1][dir][ind][2] ); 00681 if( !verts[ind] ) return MB_FAILURE; 00682 } 00683 00684 Range ents; 00685 ErrorCode rval = scImpl->impl()->get_adjacencies( verts, 2 * dim, dim, false, ents ); 00686 if( MB_SUCCESS != rval ) return rval; 00687 00688 if( ents.size() > 1 ) 00689 return MB_FAILURE; 00690 00691 else if( ents.size() == 1 ) 00692 { 00693 ent = *ents.begin(); 00694 } 00695 else if( create_if_missing ) 00696 rval = scImpl->impl()->create_element( ( 1 == dim ? MBEDGE : MBQUAD ), verts, 2 * dim, ent ); 00697 00698 return rval; 00699 } 00700 00701 #ifndef MOAB_HAVE_MPI 00702 ErrorCode ScdInterface::tag_shared_vertices( ParallelComm*, ScdBox* ) 00703 { 00704 return MB_FAILURE; 00705 #else 00706 ErrorCode ScdInterface::tag_shared_vertices( ParallelComm* pcomm, ScdBox* box ) 00707 { 00708 EntityHandle seth = box->box_set(); 00709 00710 // check the # ents in the box against the num in the set, to make sure it's only 1 box; 00711 // reuse tmp_range 00712 Range tmp_range; 00713 ErrorCode rval = mbImpl->get_entities_by_dimension( seth, box->box_dimension(), tmp_range ); 00714 if( MB_SUCCESS != rval ) return rval; 00715 if( box->num_elements() != (int)tmp_range.size() ) return MB_FAILURE; 00716 00717 const int* gdims = box->par_data().gDims; 00718 if( ( gdims[0] == gdims[3] && gdims[1] == gdims[4] && gdims[2] == gdims[5] ) || -1 == box->par_data().partMethod ) 00719 return MB_FAILURE; 00720 00721 // ok, we have a partitioned box; get the vertices shared with other processors 00722 std::vector< int > procs, offsets, shared_indices; 00723 rval = get_shared_vertices( pcomm, box, procs, offsets, shared_indices ); 00724 if( MB_SUCCESS != rval ) return rval; 00725 00726 // post receives for start handles once we know how many to look for 00727 std::vector< MPI_Request > recv_reqs( procs.size(), MPI_REQUEST_NULL ), send_reqs( procs.size(), MPI_REQUEST_NULL ); 00728 std::vector< EntityHandle > rhandles( 4 * procs.size() ), shandles( 4 ); 00729 for( unsigned int i = 0; i < procs.size(); i++ ) 00730 { 00731 int success = MPI_Irecv( (void*)&rhandles[4 * i], 4 * sizeof( EntityHandle ), MPI_UNSIGNED_CHAR, procs[i], 1, 00732 pcomm->proc_config().proc_comm(), &recv_reqs[i] ); 00733 if( success != MPI_SUCCESS ) return MB_FAILURE; 00734 } 00735 00736 // send our own start handles 00737 shandles[0] = box->start_vertex(); 00738 shandles[1] = 0; 00739 if( box->box_dimension() == 1 ) 00740 { 00741 shandles[1] = box->start_element(); 00742 shandles[2] = 0; 00743 shandles[3] = 0; 00744 } 00745 else if( box->box_dimension() == 2 ) 00746 { 00747 shandles[2] = box->start_element(); 00748 shandles[3] = 0; 00749 } 00750 else 00751 { 00752 shandles[2] = 0; 00753 shandles[3] = box->start_element(); 00754 } 00755 for( unsigned int i = 0; i < procs.size(); i++ ) 00756 { 00757 int success = MPI_Isend( (void*)&shandles[0], 4 * sizeof( EntityHandle ), MPI_UNSIGNED_CHAR, procs[i], 1, 00758 pcomm->proc_config().proc_comm(), &send_reqs[i] ); 00759 if( success != MPI_SUCCESS ) return MB_FAILURE; 00760 } 00761 00762 // receive start handles and save info to a tuple list 00763 int incoming = procs.size(); 00764 int p, j, k; 00765 MPI_Status status; 00766 TupleList shared_data; 00767 shared_data.initialize( 1, 0, 2, 0, shared_indices.size() / 2 ); 00768 shared_data.enableWriteAccess(); 00769 00770 j = 0; 00771 k = 0; 00772 while( incoming ) 00773 { 00774 int success = MPI_Waitany( procs.size(), &recv_reqs[0], &p, &status ); 00775 if( MPI_SUCCESS != success ) return MB_FAILURE; 00776 unsigned int num_indices = ( offsets[p + 1] - offsets[p] ) / 2; 00777 int *lh = &shared_indices[offsets[p]], *rh = lh + num_indices; 00778 for( unsigned int i = 0; i < num_indices; i++ ) 00779 { 00780 shared_data.vi_wr[j++] = procs[p]; 00781 shared_data.vul_wr[k++] = shandles[0] + lh[i]; 00782 shared_data.vul_wr[k++] = rhandles[4 * p] + rh[i]; 00783 shared_data.inc_n(); 00784 } 00785 incoming--; 00786 } 00787 00788 // still need to wait for the send requests 00789 std::vector< MPI_Status > mult_status( procs.size() ); 00790 int success = MPI_Waitall( procs.size(), &send_reqs[0], &mult_status[0] ); 00791 if( MPI_SUCCESS != success ) { MB_SET_ERR( MB_FAILURE, "Failed in waitall in ScdInterface::tag_shared_vertices" ); } 00792 // sort by local handle 00793 TupleList::buffer sort_buffer; 00794 sort_buffer.buffer_init( shared_indices.size() / 2 ); 00795 shared_data.sort( 1, &sort_buffer ); 00796 sort_buffer.reset(); 00797 00798 // process into sharing data 00799 std::map< std::vector< int >, std::vector< EntityHandle > > proc_nvecs; 00800 Range dum; 00801 rval = pcomm->tag_shared_verts( shared_data, proc_nvecs, dum, 0 ); 00802 if( MB_SUCCESS != rval ) return rval; 00803 00804 // create interface sets 00805 rval = pcomm->create_interface_sets( proc_nvecs ); 00806 if( MB_SUCCESS != rval ) return rval; 00807 00808 // add the box to the PComm's partitionSets 00809 pcomm->partition_sets().insert( box->box_set() ); 00810 00811 // make sure buffers are allocated for communicating procs 00812 for( std::vector< int >::iterator pit = procs.begin(); pit != procs.end(); ++pit ) 00813 pcomm->get_buffers( *pit ); 00814 00815 if( pcomm->get_debug_verbosity() > 1 ) pcomm->list_entities( NULL, 1 ); 00816 00817 #ifndef NDEBUG 00818 rval = pcomm->check_all_shared_handles(); 00819 if( MB_SUCCESS != rval ) return rval; 00820 #endif 00821 00822 return MB_SUCCESS; 00823 00824 #endif 00825 } 00826 00827 ErrorCode ScdInterface::get_neighbor_alljkbal( int np, int pfrom, const int* const gdims, const int* const gperiodic, 00828 const int* const dijk, int& pto, int* rdims, int* facedims, 00829 int* across_bdy ) 00830 { 00831 if( dijk[0] != 0 ) 00832 { 00833 pto = -1; 00834 return MB_SUCCESS; 00835 } 00836 00837 pto = -1; 00838 across_bdy[0] = across_bdy[1] = across_bdy[2] = 0; 00839 00840 int ldims[6], pijk[3], lperiodic[3]; 00841 ErrorCode rval = compute_partition_alljkbal( np, pfrom, gdims, gperiodic, ldims, lperiodic, pijk ); 00842 if( MB_SUCCESS != rval ) return rval; 00843 assert( pijk[1] * pijk[2] == np ); 00844 pto = -1; 00845 bool bot_j = pfrom< pijk[2], top_j = pfrom > np - pijk[2]; 00846 if( ( 1 == pijk[2] && dijk[2] ) || // 1d in j means no neighbors with dk != 0 00847 ( !( pfrom % pijk[2] ) && -1 == dijk[2] ) || // at -k bdy 00848 ( pfrom % pijk[2] == pijk[2] - 1 && 1 == dijk[2] ) || // at +k bdy 00849 ( pfrom < pijk[2] && -1 == dijk[1] && !gperiodic[1] ) || // down and not periodic 00850 ( pfrom >= np - pijk[2] && 1 == dijk[1] && !gperiodic[1] ) ) // up and not periodic 00851 return MB_SUCCESS; 00852 00853 pto = pfrom; 00854 std::copy( ldims, ldims + 6, rdims ); 00855 std::copy( ldims, ldims + 6, facedims ); 00856 00857 if( 0 != dijk[1] ) 00858 { 00859 pto = ( pto + dijk[1] * pijk[2] + np ) % np; 00860 assert( pto >= 0 && pto < np ); 00861 int dj = ( gdims[4] - gdims[1] ) / pijk[1], extra = ( gdims[4] - gdims[1] ) % pijk[1]; 00862 if( -1 == dijk[1] ) 00863 { 00864 facedims[4] = facedims[1]; 00865 if( bot_j ) 00866 { 00867 // going across periodic lower bdy in j 00868 rdims[4] = gdims[4]; 00869 across_bdy[1] = -1; 00870 } 00871 else 00872 { 00873 rdims[4] = ldims[1]; 00874 } 00875 rdims[1] = rdims[4] - dj; 00876 if( pto < extra ) rdims[1]--; 00877 } 00878 else 00879 { 00880 if( pfrom > np - pijk[2] ) facedims[4] = gdims[1]; 00881 facedims[1] = facedims[4]; 00882 if( top_j ) 00883 { 00884 // going across periodic upper bdy in j 00885 rdims[1] = gdims[1]; 00886 across_bdy[1] = 1; 00887 } 00888 else 00889 { 00890 rdims[1] = ldims[4]; 00891 } 00892 rdims[4] = rdims[1] + dj; 00893 if( pto < extra ) rdims[4]++; 00894 } 00895 } 00896 if( 0 != dijk[2] ) 00897 { 00898 pto = ( pto + dijk[2] ) % np; 00899 assert( pto >= 0 && pto < np ); 00900 facedims[2] = facedims[5] = ( -1 == dijk[2] ? facedims[2] : facedims[5] ); 00901 int dk = ( gdims[5] - gdims[2] ) / pijk[2]; 00902 if( -1 == dijk[2] ) 00903 { 00904 facedims[5] = facedims[2]; 00905 rdims[5] = ldims[2]; 00906 rdims[2] = rdims[5] - dk; // never any kextra for alljkbal 00907 } 00908 else 00909 { 00910 facedims[2] = facedims[5]; 00911 rdims[2] = ldims[5]; 00912 rdims[5] = rdims[2] + dk; // never any kextra for alljkbal 00913 } 00914 } 00915 00916 assert( -1 == pto || ( rdims[0] >= gdims[0] && rdims[3] <= gdims[3] ) ); 00917 assert( -1 == pto || ( rdims[1] >= gdims[1] && ( rdims[4] <= gdims[4] || ( across_bdy[1] && bot_j ) ) ) ); 00918 assert( -1 == pto || ( rdims[2] >= gdims[2] && rdims[5] <= gdims[5] ) ); 00919 assert( -1 == pto || ( ( facedims[0] >= rdims[0] || 00920 ( gperiodic[0] && rdims[3] == gdims[3] + 1 && facedims[0] == gdims[0] ) ) ) ); 00921 assert( -1 == pto || ( facedims[3] <= rdims[3] ) ); 00922 assert( -1 == pto || ( ( facedims[1] >= rdims[1] || 00923 ( gperiodic[1] && rdims[4] == gdims[4] + 1 && facedims[1] == gdims[1] ) ) ) ); 00924 assert( -1 == pto || ( facedims[4] <= rdims[4] ) ); 00925 assert( -1 == pto || ( facedims[2] >= rdims[2] ) ); 00926 assert( -1 == pto || ( facedims[5] <= rdims[5] ) ); 00927 assert( -1 == pto || ( facedims[0] >= ldims[0] ) ); 00928 assert( -1 == pto || ( facedims[3] <= ldims[3] ) ); 00929 assert( -1 == pto || ( facedims[1] >= ldims[1] ) ); 00930 assert( -1 == pto || ( facedims[4] <= ldims[4] ) ); 00931 assert( -1 == pto || ( facedims[2] >= ldims[2] ) ); 00932 assert( -1 == pto || ( facedims[5] <= ldims[5] ) ); 00933 00934 return MB_SUCCESS; 00935 } 00936 00937 ErrorCode ScdInterface::get_neighbor_sqij( int np, int pfrom, const int* const gdims, const int* const gperiodic, 00938 const int* const dijk, int& pto, int* rdims, int* facedims, int* across_bdy ) 00939 { 00940 if( dijk[2] != 0 ) 00941 { 00942 // for sqij, there is no k neighbor, ever 00943 pto = -1; 00944 return MB_SUCCESS; 00945 } 00946 00947 pto = -1; 00948 across_bdy[0] = across_bdy[1] = across_bdy[2] = 0; 00949 int lperiodic[3], pijk[3], ldims[6]; 00950 ErrorCode rval = compute_partition_sqij( np, pfrom, gdims, gperiodic, ldims, lperiodic, pijk ); 00951 if( MB_SUCCESS != rval ) return rval; 00952 assert( pijk[0] * pijk[1] == np ); 00953 pto = -1; 00954 bool top_i = 0, top_j = 0, bot_i = 0, bot_j = 0; 00955 int ni = pfrom % pijk[0], nj = pfrom / pijk[0]; // row / column number of me 00956 if( ni == pijk[0] - 1 ) top_i = 1; 00957 if( nj == pijk[1] - 1 ) top_j = 1; 00958 if( !ni ) bot_i = 1; 00959 if( !nj ) bot_j = 1; 00960 if( ( !gperiodic[0] && bot_i && -1 == dijk[0] ) || // left and not periodic 00961 ( !gperiodic[0] && top_i && 1 == dijk[0] ) || // right and not periodic 00962 ( !gperiodic[1] && bot_j && -1 == dijk[1] ) || // bottom and not periodic 00963 ( !gperiodic[1] && top_j && 1 == dijk[1] ) ) // top and not periodic 00964 return MB_SUCCESS; 00965 00966 std::copy( ldims, ldims + 6, facedims ); 00967 std::copy( ldims, ldims + 6, rdims ); 00968 pto = pfrom; 00969 int j = gdims[4] - gdims[1], dj = j / pijk[1], jextra = ( gdims[4] - gdims[1] ) % dj, i = gdims[3] - gdims[0], 00970 di = i / pijk[0], iextra = ( gdims[3] - gdims[0] ) % di; 00971 00972 if( 0 != dijk[0] ) 00973 { 00974 pto = ( ni + dijk[0] + pijk[0] ) % pijk[0]; // get pto's ni value 00975 pto = nj * pijk[0] + pto; // then convert to pto 00976 assert( pto >= 0 && pto < np ); 00977 if( -1 == dijk[0] ) 00978 { 00979 facedims[3] = facedims[0]; 00980 if( bot_i ) 00981 { 00982 // going across lower periodic bdy in i 00983 across_bdy[0] = -1; 00984 rdims[3] = gdims[3] + 1; // +1 because ldims[3] on remote proc is gdims[3]+1 00985 rdims[0] = rdims[3] - di - 1; // -1 to account for rdims[3] being one larger 00986 } 00987 else 00988 { 00989 rdims[3] = ldims[0]; 00990 rdims[0] = rdims[3] - di; 00991 } 00992 00993 if( pto % pijk[0] < iextra ) rdims[0]--; 00994 } 00995 else 00996 { 00997 if( top_i ) 00998 { 00999 // going across lower periodic bdy in i 01000 facedims[3] = gdims[0]; 01001 across_bdy[0] = 1; 01002 } 01003 facedims[0] = facedims[3]; 01004 rdims[0] = ( top_i ? gdims[0] : ldims[3] ); 01005 rdims[3] = rdims[0] + di; 01006 if( pto % pijk[0] < iextra ) rdims[3]++; 01007 if( gperiodic[0] && ni == pijk[0] - 2 ) rdims[3]++; // remote proc is top_i and periodic 01008 } 01009 } 01010 if( 0 != dijk[1] ) 01011 { 01012 pto = ( pto + dijk[1] * pijk[0] + np ) % np; 01013 assert( pto >= 0 && pto < np ); 01014 if( -1 == dijk[1] ) 01015 { 01016 facedims[4] = facedims[1]; 01017 if( bot_j ) 01018 { 01019 // going across lower periodic bdy in j 01020 rdims[4] = gdims[4] + 1; // +1 because ldims[4] on remote proc is gdims[4]+1 01021 rdims[1] = rdims[4] - dj - 1; // -1 to account for gdims[4] being one larger 01022 across_bdy[1] = -1; 01023 } 01024 else 01025 { 01026 rdims[4] = ldims[1]; 01027 rdims[1] = rdims[4] - dj; 01028 } 01029 if( pto / pijk[0] < jextra ) rdims[1]--; 01030 } 01031 else 01032 { 01033 if( top_j ) 01034 { 01035 // going across lower periodic bdy in j 01036 facedims[4] = gdims[1]; 01037 rdims[1] = gdims[1]; 01038 across_bdy[1] = 1; 01039 } 01040 else 01041 { 01042 rdims[1] = ldims[4]; 01043 } 01044 facedims[1] = facedims[4]; 01045 rdims[4] = rdims[1] + dj; 01046 if( nj + 1 < jextra ) rdims[4]++; 01047 if( gperiodic[1] && nj == pijk[1] - 2 ) rdims[4]++; // remote proc is top_j and periodic 01048 } 01049 } 01050 01051 // rdims within gdims 01052 assert( -1 == pto || ( rdims[0] >= gdims[0] && 01053 ( rdims[3] <= gdims[3] + ( gperiodic[0] && pto % pijk[0] == pijk[0] - 1 ? 1 : 0 ) ) ) ); 01054 assert( -1 == pto || ( rdims[1] >= gdims[1] && 01055 ( rdims[4] <= gdims[4] + ( gperiodic[1] && pto / pijk[0] == pijk[1] - 1 ? 1 : 0 ) ) ) ); 01056 assert( -1 == pto || ( rdims[2] >= gdims[2] && rdims[5] <= gdims[5] ) ); 01057 // facedims within rdims 01058 assert( -1 == pto || ( ( facedims[0] >= rdims[0] || 01059 ( gperiodic[0] && pto % pijk[0] == pijk[0] - 1 && facedims[0] == gdims[0] ) ) ) ); 01060 assert( -1 == pto || ( facedims[3] <= rdims[3] ) ); 01061 assert( -1 == pto || ( ( facedims[1] >= rdims[1] || 01062 ( gperiodic[1] && pto / pijk[0] == pijk[1] - 1 && facedims[1] == gdims[1] ) ) ) ); 01063 assert( -1 == pto || ( facedims[4] <= rdims[4] ) ); 01064 assert( -1 == pto || ( facedims[2] >= rdims[2] && facedims[5] <= rdims[5] ) ); 01065 // facedims within ldims 01066 assert( -1 == pto || ( ( facedims[0] >= ldims[0] || ( top_i && facedims[0] == gdims[0] ) ) ) ); 01067 assert( -1 == pto || ( facedims[3] <= ldims[3] ) ); 01068 assert( -1 == pto || ( ( facedims[1] >= ldims[1] || ( gperiodic[1] && top_j && facedims[1] == gdims[1] ) ) ) ); 01069 assert( -1 == pto || ( facedims[4] <= ldims[4] ) ); 01070 assert( -1 == pto || ( facedims[2] >= ldims[2] && facedims[5] <= ldims[5] ) ); 01071 01072 return MB_SUCCESS; 01073 } 01074 01075 ErrorCode ScdInterface::get_neighbor_sqjk( int np, int pfrom, const int* const gdims, const int* const gperiodic, 01076 const int* const dijk, int& pto, int* rdims, int* facedims, int* across_bdy ) 01077 { 01078 if( dijk[0] != 0 ) 01079 { 01080 pto = -1; 01081 return MB_SUCCESS; 01082 } 01083 01084 pto = -1; 01085 across_bdy[0] = across_bdy[1] = across_bdy[2] = 0; 01086 int pijk[3], lperiodic[3], ldims[6]; 01087 ErrorCode rval = compute_partition_sqjk( np, pfrom, gdims, gperiodic, ldims, lperiodic, pijk ); 01088 if( MB_SUCCESS != rval ) return rval; 01089 assert( pijk[1] * pijk[2] == np ); 01090 pto = -1; 01091 bool top_j = 0, top_k = 0, bot_j = 0, bot_k = 0; 01092 int nj = pfrom % pijk[1], nk = pfrom / pijk[1]; 01093 if( nj == pijk[1] - 1 ) top_j = 1; 01094 if( nk == pijk[2] - 1 ) top_k = 1; 01095 if( !nj ) bot_j = 1; 01096 if( !nk ) bot_k = 1; 01097 if( ( !gperiodic[1] && bot_j && -1 == dijk[1] ) || // down and not periodic 01098 ( !gperiodic[1] && top_j && 1 == dijk[1] ) || // up and not periodic 01099 ( bot_k && -1 == dijk[2] ) || // k- bdy 01100 ( top_k && 1 == dijk[2] ) ) // k+ bdy 01101 return MB_SUCCESS; 01102 01103 std::copy( ldims, ldims + 6, facedims ); 01104 std::copy( ldims, ldims + 6, rdims ); 01105 pto = pfrom; 01106 int dj = ( gdims[4] - gdims[1] ) / pijk[1], jextra = ( gdims[4] - gdims[1] ) % dj, 01107 dk = ( gdims[5] == gdims[2] ? 0 : ( gdims[5] - gdims[2] ) / pijk[2] ), 01108 kextra = ( gdims[5] - gdims[2] ) - dk * pijk[2]; 01109 assert( ( dj * pijk[1] + jextra == ( gdims[4] - gdims[1] ) ) && 01110 ( dk * pijk[2] + kextra == ( gdims[5] - gdims[2] ) ) ); 01111 if( 0 != dijk[1] ) 01112 { 01113 pto = ( nj + dijk[1] + pijk[1] ) % pijk[1]; // get pto's ni value 01114 pto = nk * pijk[1] + pto; // then convert to pto 01115 assert( pto >= 0 && pto < np ); 01116 if( -1 == dijk[1] ) 01117 { 01118 facedims[4] = facedims[1]; 01119 if( bot_j ) 01120 { 01121 // going across lower periodic bdy in j 01122 rdims[4] = gdims[4] + 1; // +1 because ldims[4] on remote proc is gdims[4]+1 01123 across_bdy[1] = -1; 01124 } 01125 else 01126 { 01127 rdims[4] = ldims[1]; 01128 } 01129 rdims[1] = rdims[4] - dj; 01130 if( nj < jextra ) rdims[1]--; 01131 } 01132 else 01133 { 01134 if( top_j ) 01135 { 01136 // going across upper periodic bdy in j 01137 rdims[1] = gdims[1]; 01138 facedims[4] = gdims[1]; 01139 across_bdy[1] = 1; 01140 } 01141 else 01142 { 01143 rdims[1] = ldims[4]; 01144 } 01145 facedims[1] = facedims[4]; 01146 rdims[4] = rdims[1] + dj; 01147 if( nj < jextra ) rdims[4]++; 01148 if( gperiodic[1] && nj == dijk[1] - 2 ) rdims[4]++; // +1 because next proc is on periodic bdy 01149 } 01150 } 01151 if( 0 != dijk[2] ) 01152 { 01153 pto = ( pto + dijk[2] * pijk[1] + np ) % np; 01154 assert( pto >= 0 && pto < np ); 01155 if( -1 == dijk[2] ) 01156 { 01157 facedims[5] = facedims[2]; 01158 rdims[5] = ldims[2]; 01159 rdims[2] -= dk; 01160 if( pto / pijk[1] < kextra ) rdims[2]--; 01161 } 01162 else 01163 { 01164 facedims[2] = facedims[5]; 01165 rdims[2] = ldims[5]; 01166 rdims[5] += dk; 01167 if( pto / pijk[1] < kextra ) rdims[5]++; 01168 } 01169 } 01170 01171 assert( -1 == pto || ( rdims[0] >= gdims[0] && rdims[3] <= gdims[3] ) ); 01172 assert( -1 == pto || ( rdims[1] >= gdims[1] && ( rdims[4] <= gdims[4] || ( across_bdy[1] && bot_j ) ) ) ); 01173 assert( -1 == pto || ( rdims[2] >= gdims[2] && rdims[5] <= gdims[5] ) ); 01174 assert( -1 == pto || ( facedims[0] >= rdims[0] && facedims[3] <= rdims[3] ) ); 01175 assert( -1 == pto || 01176 ( ( facedims[1] >= rdims[1] || ( gperiodic[1] && rdims[4] == gdims[4] && facedims[1] == gdims[1] ) ) ) ); 01177 assert( -1 == pto || ( facedims[4] <= rdims[4] ) ); 01178 assert( -1 == pto || ( facedims[2] >= rdims[2] && facedims[5] <= rdims[5] ) ); 01179 assert( -1 == pto || ( facedims[0] >= ldims[0] && facedims[3] <= ldims[3] ) ); 01180 assert( -1 == pto || ( facedims[1] >= ldims[1] && facedims[4] <= ldims[4] ) ); 01181 assert( -1 == pto || ( facedims[2] >= ldims[2] && facedims[5] <= ldims[5] ) ); 01182 01183 return MB_SUCCESS; 01184 } 01185 01186 ErrorCode ScdInterface::get_neighbor_sqijk( int np, int pfrom, const int* const gdims, const int* const gperiodic, 01187 const int* const dijk, int& pto, int* rdims, int* facedims, 01188 int* across_bdy ) 01189 { 01190 if( gperiodic[0] || gperiodic[1] || gperiodic[2] ) return MB_FAILURE; 01191 01192 pto = -1; 01193 across_bdy[0] = across_bdy[1] = across_bdy[2] = 0; 01194 int pijk[3], lperiodic[3], ldims[6]; 01195 ErrorCode rval = compute_partition_sqijk( np, pfrom, gdims, gperiodic, ldims, lperiodic, pijk ); 01196 if( MB_SUCCESS != rval ) return rval; 01197 assert( pijk[0] * pijk[1] * pijk[2] == np ); 01198 pto = -1; 01199 bool top[3] = { false, false, false }, bot[3] = { false, false, false }; 01200 // nijk: rank in i/j/k direction 01201 int nijk[3] = { pfrom % pijk[0], ( pfrom % ( pijk[0] * pijk[1] ) ) / pijk[0], pfrom / ( pijk[0] * pijk[1] ) }; 01202 01203 for( int i = 0; i < 3; i++ ) 01204 { 01205 if( nijk[i] == pijk[i] - 1 ) top[i] = true; 01206 if( !nijk[i] ) bot[i] = true; 01207 if( ( !gperiodic[i] && bot[i] && -1 == dijk[i] ) || // downward && not periodic 01208 ( !gperiodic[i] && top[i] && 1 == dijk[i] ) ) // upward && not periodic 01209 return MB_SUCCESS; 01210 } 01211 01212 std::copy( ldims, ldims + 6, facedims ); 01213 std::copy( ldims, ldims + 6, rdims ); 01214 pto = pfrom; 01215 int delijk[3], extra[3]; 01216 // nijk_to: rank of pto in i/j/k direction 01217 int nijk_to[3]; 01218 for( int i = 0; i < 3; i++ ) 01219 { 01220 delijk[i] = ( gdims[i + 3] == gdims[i] ? 0 : ( gdims[i + 3] - gdims[i] ) / pijk[i] ); 01221 extra[i] = ( gdims[i + 3] - gdims[i] ) % delijk[i]; 01222 nijk_to[i] = ( nijk[i] + dijk[i] + pijk[i] ) % pijk[i]; 01223 } 01224 pto = nijk_to[2] * pijk[0] * pijk[1] + nijk_to[1] * pijk[0] + nijk_to[0]; 01225 assert( pto >= 0 && pto < np ); 01226 for( int i = 0; i < 3; i++ ) 01227 { 01228 if( 0 != dijk[i] ) 01229 { 01230 if( -1 == dijk[i] ) 01231 { 01232 facedims[i + 3] = facedims[i]; 01233 if( bot[i] ) 01234 { 01235 // going across lower periodic bdy in i 01236 rdims[i + 3] = gdims[i + 3] + 1; // +1 because ldims[4] on remote proc is gdims[4]+1 01237 across_bdy[i] = -1; 01238 } 01239 else 01240 { 01241 rdims[i + 3] = ldims[i]; 01242 } 01243 rdims[i] = rdims[i + 3] - delijk[i]; 01244 if( nijk[i] < extra[i] ) rdims[i]--; 01245 } 01246 else 01247 { 01248 if( top[i] ) 01249 { 01250 // going across upper periodic bdy in i 01251 rdims[i] = gdims[i]; 01252 facedims[i + 3] = gdims[i]; 01253 across_bdy[i] = 1; 01254 } 01255 else 01256 { 01257 rdims[i] = ldims[i + 3]; 01258 } 01259 facedims[i] = facedims[i + 3]; 01260 rdims[i + 3] = rdims[i] + delijk[i]; 01261 if( nijk[i] < extra[i] ) rdims[i + 3]++; 01262 if( gperiodic[i] && nijk[i] == dijk[i] - 2 ) rdims[i + 3]++; // +1 because next proc is on periodic bdy 01263 } 01264 } 01265 } 01266 01267 assert( -1 != pto ); 01268 #ifndef NDEBUG 01269 for( int i = 0; i < 3; i++ ) 01270 { 01271 assert( ( rdims[i] >= gdims[i] && ( rdims[i + 3] <= gdims[i + 3] || ( across_bdy[i] && bot[i] ) ) ) ); 01272 assert( ( ( facedims[i] >= rdims[i] || 01273 ( gperiodic[i] && rdims[i + 3] == gdims[i + 3] && facedims[i] == gdims[i] ) ) ) ); 01274 assert( ( facedims[i] >= ldims[i] && facedims[i + 3] <= ldims[i + 3] ) ); 01275 } 01276 #endif 01277 01278 return MB_SUCCESS; 01279 } 01280 01281 ErrorCode ScdInterface::get_neighbor_alljorkori( int np, int pfrom, const int* const gdims, const int* const gperiodic, 01282 const int* const dijk, int& pto, int* rdims, int* facedims, 01283 int* across_bdy ) 01284 { 01285 ErrorCode rval = MB_SUCCESS; 01286 pto = -1; 01287 if( np == 1 ) return MB_SUCCESS; 01288 01289 int pijk[3], lperiodic[3], ldims[6]; 01290 rval = compute_partition_alljorkori( np, pfrom, gdims, gperiodic, ldims, lperiodic, pijk ); 01291 if( MB_SUCCESS != rval ) return rval; 01292 01293 int ind = -1; 01294 across_bdy[0] = across_bdy[1] = across_bdy[2] = 0; 01295 01296 for( int i = 0; i < 3; i++ ) 01297 { 01298 if( pijk[i] > 1 ) 01299 { 01300 ind = i; 01301 break; 01302 } 01303 } 01304 01305 assert( -1 < ind ); 01306 01307 if( !dijk[ind] ) 01308 // no neighbor, pto is already -1, return 01309 return MB_SUCCESS; 01310 01311 bool is_periodic = ( ( gperiodic[0] && ind == 0 ) || ( gperiodic[1] && ind == 1 ) ); 01312 if( dijk[( ind + 1 ) % 3] || dijk[( ind + 2 ) % 3] || // stepping in either other two directions 01313 ( !is_periodic && ldims[ind] == gdims[ind] && dijk[ind] == -1 ) || // lower side and going lower 01314 ( !is_periodic && ldims[3 + ind] >= gdims[3 + ind] && 01315 dijk[ind] == 1 ) ) // not >= because ldims is only > gdims when periodic; 01316 // higher side and going higher 01317 return MB_SUCCESS; 01318 01319 std::copy( ldims, ldims + 6, facedims ); 01320 std::copy( ldims, ldims + 6, rdims ); 01321 01322 int dind = ( gdims[ind + 3] - gdims[ind] ) / np; 01323 int extra = ( gdims[ind + 3] - gdims[ind] ) % np; 01324 if( -1 == dijk[ind] && pfrom ) 01325 { 01326 // actual left neighbor 01327 pto = pfrom - 1; // no need for %np, because pfrom > 0 01328 facedims[ind + 3] = facedims[ind]; 01329 rdims[ind + 3] = ldims[ind]; 01330 rdims[ind] = rdims[ind + 3] - dind - ( pto < extra ? 1 : 0 ); 01331 } 01332 else if( 1 == dijk[ind] && pfrom < np - 1 ) 01333 { 01334 // actual right neighbor 01335 pto = pfrom + 1; 01336 facedims[ind] = facedims[ind + 3]; 01337 rdims[ind] = ldims[ind + 3]; 01338 rdims[ind + 3] = rdims[ind] + dind + ( pto < extra ? 1 : 0 ); 01339 if( is_periodic && pfrom == np - 2 ) rdims[ind + 3]++; // neighbor is on periodic bdy 01340 } 01341 else if( -1 == dijk[ind] && !pfrom && gperiodic[ind] ) 01342 { 01343 // downward across periodic bdy 01344 pto = np - 1; 01345 facedims[ind + 3] = facedims[ind] = gdims[ind]; // by convention, facedims is within gdims, so lower value 01346 rdims[ind + 3] = 01347 gdims[ind + 3] + 1; // by convention, local dims one greater than gdims to indicate global lower value 01348 rdims[ind] = rdims[ind + 3] - dind - 1; 01349 across_bdy[ind] = -1; 01350 } 01351 else if( 1 == dijk[ind] && pfrom == np - 1 && is_periodic ) 01352 { 01353 // right across periodic bdy 01354 pto = 0; 01355 facedims[ind + 3] = facedims[ind] = gdims[ind]; // by convention, facedims is within gdims, so lowest value 01356 rdims[ind] = gdims[ind]; 01357 rdims[ind + 3] = rdims[ind] + dind + ( pto < extra ? 1 : 0 ); 01358 across_bdy[ind] = 1; 01359 } 01360 01361 assert( -1 == pto || ( rdims[0] >= gdims[0] && ( rdims[3] <= gdims[3] || ( across_bdy[0] && !pfrom ) ) ) ); 01362 assert( -1 == pto || ( rdims[1] >= gdims[1] && ( rdims[4] <= gdims[4] || ( across_bdy[1] && !pfrom ) ) ) ); 01363 assert( -1 == pto || ( rdims[2] >= gdims[2] && rdims[5] <= gdims[5] ) ); 01364 assert( -1 == pto || ( facedims[0] >= rdims[0] && facedims[3] <= rdims[3] ) ); 01365 assert( -1 == pto || ( facedims[1] >= rdims[1] && facedims[4] <= rdims[4] ) ); 01366 assert( -1 == pto || ( facedims[2] >= rdims[2] && facedims[5] <= rdims[5] ) ); 01367 assert( -1 == pto || ( facedims[0] >= ldims[0] && facedims[3] <= ldims[3] ) ); 01368 assert( -1 == pto || ( facedims[1] >= ldims[1] && facedims[4] <= ldims[4] ) ); 01369 assert( -1 == pto || ( facedims[2] >= ldims[2] && facedims[5] <= ldims[5] ) ); 01370 01371 return rval; 01372 } 01373 01374 //! get shared vertices for alljorkori partition scheme 01375 #ifndef MOAB_HAVE_MPI 01376 ErrorCode ScdInterface::get_shared_vertices( ParallelComm*, ScdBox*, std::vector< int >&, std::vector< int >&, 01377 std::vector< int >& ) 01378 { 01379 return MB_FAILURE; 01380 #else 01381 ErrorCode ScdInterface::get_shared_vertices( ParallelComm* pcomm, ScdBox* box, std::vector< int >& procs, 01382 std::vector< int >& offsets, std::vector< int >& shared_indices ) 01383 { 01384 // get index of partitioned dimension 01385 const int* ldims = box->box_dims(); 01386 ErrorCode rval; 01387 int ijkrem[6], ijkface[6], across_bdy[3]; 01388 01389 for( int k = -1; k <= 1; k++ ) 01390 { 01391 for( int j = -1; j <= 1; j++ ) 01392 { 01393 for( int i = -1; i <= 1; i++ ) 01394 { 01395 if( !i && !j && !k ) continue; 01396 int pto; 01397 int dijk[] = { i, j, k }; 01398 rval = get_neighbor( pcomm->proc_config().proc_size(), pcomm->proc_config().proc_rank(), 01399 box->par_data(), dijk, pto, ijkrem, ijkface, across_bdy ); 01400 if( MB_SUCCESS != rval ) return rval; 01401 if( -1 != pto ) 01402 { 01403 if( procs.empty() || pto != *procs.rbegin() ) 01404 { 01405 procs.push_back( pto ); 01406 offsets.push_back( shared_indices.size() ); 01407 } 01408 rval = get_indices( ldims, ijkrem, across_bdy, ijkface, shared_indices ); 01409 if( MB_SUCCESS != rval ) return rval; 01410 01411 // check indices against known #verts on local and remote 01412 // begin of this block is shared_indices[*offsets.rbegin()], end is 01413 // shared_indices.end(), halfway is 01414 // (shared_indices.size()-*offsets.rbegin())/2 01415 #ifndef NDEBUG 01416 int start_idx = *offsets.rbegin(), end_idx = shared_indices.size(), 01417 mid_idx = ( start_idx + end_idx ) / 2; 01418 01419 int num_local_verts = ( ldims[3] - ldims[0] + 1 ) * ( ldims[4] - ldims[1] + 1 ) * 01420 ( -1 == ldims[2] && -1 == ldims[5] ? 1 : ( ldims[5] - ldims[2] + 1 ) ), 01421 num_remote_verts = ( ijkrem[3] - ijkrem[0] + 1 ) * ( ijkrem[4] - ijkrem[1] + 1 ) * 01422 ( -1 == ijkrem[2] && -1 == ijkrem[5] ? 1 : ( ijkrem[5] - ijkrem[2] + 1 ) ); 01423 01424 assert( 01425 *std::min_element( &shared_indices[start_idx], &shared_indices[mid_idx] ) >= 0 && 01426 *std::max_element( &shared_indices[start_idx], &shared_indices[mid_idx] ) < num_local_verts && 01427 *std::min_element( &shared_indices[mid_idx], &shared_indices[end_idx] ) >= 0 && 01428 *std::max_element( &shared_indices[mid_idx], &shared_indices[end_idx] ) < num_remote_verts ); 01429 #endif 01430 } 01431 } 01432 } 01433 } 01434 01435 offsets.push_back( shared_indices.size() ); 01436 01437 return MB_SUCCESS; 01438 #endif 01439 } 01440 01441 std::ostream& operator<<( std::ostream& str, const ScdParData& pd ) 01442 { 01443 str << "Partition method = " << ScdParData::PartitionMethodNames[pd.partMethod] << ", gDims = (" << pd.gDims[0] 01444 << "," << pd.gDims[1] << "," << pd.gDims[2] << ")-(" << pd.gDims[3] << "," << pd.gDims[4] << "," << pd.gDims[5] 01445 << "), gPeriodic = (" << pd.gPeriodic[0] << "," << pd.gPeriodic[1] << "," << pd.gPeriodic[2] << "), pDims = (" 01446 << pd.pDims[0] << "," << pd.pDims[1] << "," << pd.pDims[2] << ")" << std::endl; 01447 return str; 01448 } 01449 01450 } // namespace moab