MOAB: Mesh Oriented datABase  (version 5.3.1)
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, 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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines