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