Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ScdInterface.cpp
Go to the documentation of this file.
1 #include "moab/ScdInterface.hpp"
2 #include "moab/Core.hpp"
3 #include "SequenceManager.hpp"
4 #include "EntitySequence.hpp"
6 #include "VertexSequence.hpp"
7 #include "ScdVertexData.hpp"
8 #include "MBTagConventions.hpp"
9 #ifdef MOAB_HAVE_MPI
10 #include "moab/ParallelComm.hpp"
11 #include "moab/TupleList.hpp"
12 #include "moab/gs.hpp"
13 #endif
14 #include <cassert>
15 #include <iostream>
16 #include <functional>
17 
18 #define ERRORR( rval, str ) \
19  { \
20  if( MB_SUCCESS != ( rval ) ) \
21  { \
22  std::cerr << ( str ); \
23  return rval; \
24  } \
25  }
26 
27 namespace moab
28 {
29 
30 const char* ScdParData::PartitionMethodNames[] = { "alljorkori", "alljkbal", "sqij", "sqjk",
31  "sqijk", "trivial", "rcbzoltan", "nopart" };
32 
34  : mbImpl( imp ), searchedBoxes( false ), boxPeriodicTag( 0 ), boxDimsTag( 0 ), globalBoxDimsTag( 0 ),
35  partMethodTag( 0 ), boxSetTag( 0 )
36 {
37  if( boxes ) find_boxes( scdBoxes );
38 }
39 
40 // Destructor
42 {
43  std::vector< ScdBox* > tmp_boxes;
44  tmp_boxes.swap( scdBoxes );
45 
46  for( std::vector< ScdBox* >::iterator rit = tmp_boxes.begin(); rit != tmp_boxes.end(); ++rit )
47  delete *rit;
48 
49  if( box_set_tag( false ) ) mbImpl->tag_delete( box_set_tag() );
50 }
51 
53 {
54  return mbImpl;
55 }
56 
57 ErrorCode ScdInterface::find_boxes( std::vector< ScdBox* >& scd_boxes )
58 {
59  Range tmp_boxes;
60  ErrorCode rval = find_boxes( tmp_boxes );
61  if( MB_SUCCESS != rval ) return rval;
62 
63  for( Range::iterator rit = tmp_boxes.begin(); rit != tmp_boxes.end(); ++rit )
64  {
65  ScdBox* tmp_box = get_scd_box( *rit );
66  if( tmp_box )
67  scd_boxes.push_back( tmp_box );
68  else
69  rval = MB_FAILURE;
70  }
71 
72  return rval;
73 }
74 
76 {
77  ErrorCode rval = MB_SUCCESS;
78  box_dims_tag();
79  Range boxes;
80  if( !searchedBoxes )
81  {
83  searchedBoxes = true;
84  if( !boxes.empty() )
85  {
86  scdBoxes.resize( boxes.size() );
87  rval = mbImpl->tag_get_data( boxSetTag, boxes, &scdBoxes[0] );
88  ScdBox* dum = nullptr;
89  // std::remove_if(scdBoxes.begin(), scdBoxes.end(),
90  // std::bind2nd(std::equal_to<ScdBox*>(), dum) ) ;
91  std::remove_if( scdBoxes.begin(), scdBoxes.end(),
92  std::bind( std::equal_to< ScdBox* >(), std::placeholders::_1, dum ) );
93  // https://stackoverflow.com/questions/32739018/a-replacement-for-stdbind2nd
94  }
95  }
96 
97  for( std::vector< ScdBox* >::iterator vit = scdBoxes.begin(); vit != scdBoxes.end(); ++vit )
98  scd_boxes.insert( ( *vit )->box_set() );
99 
100  return rval;
101 }
102 
104 {
105  ScdBox* scd_box = NULL;
106  if( !box_set_tag( false ) ) return scd_box;
107 
108  mbImpl->tag_get_data( box_set_tag(), &eh, 1, &scd_box );
109  return scd_box;
110 }
111 
113  HomCoord high,
114  const double* const coords,
115  unsigned int num_coords,
116  ScdBox*& new_box,
117  int* const lperiodic,
118  ScdParData* par_data,
119  bool assign_gids,
120  int tag_shared_ents )
121 {
122  // create a rectangular structured mesh block
123  ErrorCode rval;
124 
125  int tmp_lper[3] = { 0, 0, 0 };
126  if( lperiodic ) std::copy( lperiodic, lperiodic + 3, tmp_lper );
127 
128 #ifndef MOAB_HAVE_MPI
129  if( -1 != tag_shared_ents ) ERRORR( MB_FAILURE, "Parallel capability requested but MOAB not compiled parallel." );
130  if( -1 == tag_shared_ents && !assign_gids ) assign_gids = true; // need to assign gids in order to tag shared verts
131 #else
132  if( par_data && low == high && ScdParData::NOPART != par_data->partMethod )
133  {
134  // requesting creation of parallel mesh, so need to compute partition
135  if( !par_data->pComm )
136  {
137  // this is a really boneheaded way to have to create a PC
138  par_data->pComm = ParallelComm::get_pcomm( mbImpl, 0 );
139  if( NULL == par_data->pComm ) par_data->pComm = new ParallelComm( mbImpl, MPI_COMM_WORLD );
140  }
141  int ldims[6];
142  rval = compute_partition( par_data->pComm->size(), par_data->pComm->rank(), *par_data, ldims, tmp_lper,
143  par_data->pDims );ERRORR( rval, "Error returned from compute_partition." );
144  low.set( ldims[0], ldims[1], ldims[2] );
145  high.set( ldims[3], ldims[4], ldims[5] );
146  if( par_data->pComm->get_debug_verbosity() > 0 )
147  {
148  std::cout << "Proc " << par_data->pComm->rank() << ": " << *par_data;
149  std::cout << "Proc " << par_data->pComm->rank() << " local dims: " << low << "-" << high << std::endl;
150  }
151  }
152 #endif
153 
154  HomCoord tmp_size = high - low + HomCoord( 1, 1, 1, 0 );
155  if( ( tmp_size[1] && num_coords && (int)num_coords < tmp_size[0] ) ||
156  ( tmp_size[2] && num_coords && (int)num_coords < tmp_size[0] * tmp_size[1] ) )
157  return MB_FAILURE;
158 
159  rval = create_scd_sequence( low, high, MBVERTEX, 0, new_box );ERRORR( rval, "Trouble creating scd vertex sequence." );
160 
161  // set the vertex coordinates
162  double *xc, *yc, *zc;
163  rval = new_box->get_coordinate_arrays( xc, yc, zc );ERRORR( rval, "Couldn't get vertex coordinate arrays." );
164 
165  if( coords && num_coords )
166  {
167  unsigned int i = 0;
168  for( int kl = low[2]; kl <= high[2]; kl++ )
169  {
170  for( int jl = low[1]; jl <= high[1]; jl++ )
171  {
172  for( int il = low[0]; il <= high[0]; il++ )
173  {
174  xc[i] = coords[3 * i];
175  if( new_box->box_size()[1] ) yc[i] = coords[3 * i + 1];
176  if( new_box->box_size()[2] ) zc[i] = coords[3 * i + 2];
177  i++;
178  }
179  }
180  }
181  }
182  else
183  {
184  unsigned int i = 0;
185  for( int kl = low[2]; kl <= high[2]; kl++ )
186  {
187  for( int jl = low[1]; jl <= high[1]; jl++ )
188  {
189  for( int il = low[0]; il <= high[0]; il++ )
190  {
191  xc[i] = (double)il;
192  if( new_box->box_size()[1] )
193  yc[i] = (double)jl;
194  else
195  yc[i] = 0.0;
196  if( new_box->box_size()[2] )
197  zc[i] = (double)kl;
198  else
199  zc[i] = 0.0;
200  i++;
201  }
202  }
203  }
204  }
205 
206  // create element sequence
207  Core* mbcore = dynamic_cast< Core* >( mbImpl );
208  SequenceManager* seq_mgr = mbcore->sequence_manager();
209 
210  EntitySequence* tmp_seq;
211  EntityHandle start_ent;
212 
213  // construct the sequence
214  EntityType this_tp = MBHEX;
215  if( 1 >= tmp_size[2] ) this_tp = MBQUAD;
216  if( 1 >= tmp_size[2] && 1 >= tmp_size[1] ) this_tp = MBEDGE;
217  rval = seq_mgr->create_scd_sequence( low, high, this_tp, 0, start_ent, tmp_seq, tmp_lper );ERRORR( rval, "Trouble creating scd element sequence." );
218 
219  new_box->elem_seq( tmp_seq );
220  new_box->start_element( start_ent );
221 
222  // add vertex seq to element seq, forward orientation, unity transform
223  rval = new_box->add_vbox( new_box,
224  // p1: imin,jmin
225  low, low,
226  // p2: imax,jmin
227  low + HomCoord( 1, 0, 0 ), low + HomCoord( 1, 0, 0 ),
228  // p3: imin,jmax
229  low + HomCoord( 0, 1, 0 ), low + HomCoord( 0, 1, 0 ) );ERRORR( rval, "Error constructing structured element sequence." );
230 
231  // add the new hexes to the scd box set; vertices were added in call to create_scd_sequence
232  Range tmp_range( new_box->start_element(), new_box->start_element() + new_box->num_elements() - 1 );
233  rval = mbImpl->add_entities( new_box->box_set(), tmp_range );ERRORR( rval, "Couldn't add new hexes to box set." );
234 
235  if( par_data ) new_box->par_data( *par_data );
236 
237  if( assign_gids )
238  {
239  rval = assign_global_ids( new_box );ERRORR( rval, "Trouble assigning global ids" );
240  }
241 
242 #ifdef MOAB_HAVE_MPI
243  if( par_data && -1 != tag_shared_ents )
244  {
245  rval = tag_shared_vertices( par_data->pComm, new_box );
246  }
247 #endif
248 
249  return MB_SUCCESS;
250 }
251 
253 {
254  // Get a ptr to global id memory
255  void* data;
256  int count = 0;
257  Tag gid_tag = mbImpl->globalId_tag();
258  Range tmp_range( box->start_vertex(), box->start_vertex() + box->num_vertices() );
259  ErrorCode rval = mbImpl->tag_iterate( gid_tag, tmp_range.begin(), tmp_range.end(), count, data );ERRORR( rval, "Failed to get tag iterator." );
260  assert( count == box->num_vertices() );
261  int* gid_data = (int*)data;
262  int di = box->par_data().gDims[3] - box->par_data().gDims[0] + 1;
263  int dj = box->par_data().gDims[4] - box->par_data().gDims[1] + 1;
264 
265  for( int kl = box->box_dims()[2]; kl <= box->box_dims()[5]; kl++ )
266  {
267  for( int jl = box->box_dims()[1]; jl <= box->box_dims()[4]; jl++ )
268  {
269  for( int il = box->box_dims()[0]; il <= box->box_dims()[3]; il++ )
270  {
271  int itmp =
272  ( !box->locally_periodic()[0] && box->par_data().gPeriodic[0] && il == box->par_data().gDims[3]
273  ? box->par_data().gDims[0]
274  : il );
275  *gid_data = ( -1 != kl ? kl * di * dj : 0 ) + jl * di + itmp + 1;
276  gid_data++;
277  }
278  }
279  }
280 
281  return MB_SUCCESS;
282 }
283 
285  const HomCoord& high,
286  EntityType tp,
287  int starting_id,
288  ScdBox*& new_box,
289  int* is_periodic )
290 {
291  HomCoord tmp_size = high - low + HomCoord( 1, 1, 1, 0 );
292  if( ( tp == MBHEX && 1 >= tmp_size[2] ) || ( tp == MBQUAD && 1 >= tmp_size[1] ) ||
293  ( tp == MBEDGE && 1 >= tmp_size[0] ) )
294  return MB_TYPE_OUT_OF_RANGE;
295 
296  Core* mbcore = dynamic_cast< Core* >( mbImpl );
297  assert( mbcore != NULL );
298  SequenceManager* seq_mgr = mbcore->sequence_manager();
299 
300  EntitySequence* tmp_seq;
301  EntityHandle start_ent, scd_set;
302 
303  // construct the sequence
304  ErrorCode rval = seq_mgr->create_scd_sequence( low, high, tp, starting_id, start_ent, tmp_seq, is_periodic );
305  if( MB_SUCCESS != rval ) return rval;
306 
307  // create the set for this rectangle
308  rval = create_box_set( low, high, scd_set );
309  if( MB_SUCCESS != rval ) return rval;
310 
311  // make the ScdBox
312  new_box = new ScdBox( this, scd_set, tmp_seq );
313  if( !new_box ) return MB_FAILURE;
314 
315  // set the start vertex/element
316  Range new_range;
317  if( MBVERTEX == tp )
318  {
319  new_range.insert( start_ent, start_ent + new_box->num_vertices() - 1 );
320  }
321  else
322  {
323  new_range.insert( start_ent, start_ent + new_box->num_elements() - 1 );
324  }
325 
326  // put the entities in the box set
327  rval = mbImpl->add_entities( scd_set, new_range );
328  if( MB_SUCCESS != rval ) return rval;
329 
330  // tag the set with the box
331  rval = mbImpl->tag_set_data( box_set_tag(), &scd_set, 1, &new_box );
332  if( MB_SUCCESS != rval ) return rval;
333 
334  return MB_SUCCESS;
335 }
336 
338  const HomCoord& high,
339  EntityHandle& scd_set,
340  int* is_periodic )
341 {
342  // create the set and put the entities in it
343  ErrorCode rval = mbImpl->create_meshset( MESHSET_SET, scd_set );
344  if( MB_SUCCESS != rval ) return rval;
345 
346  // tag the set with parameter extents
347  int boxdims[6];
348  for( int i = 0; i < 3; i++ )
349  boxdims[i] = low[i];
350  for( int i = 0; i < 3; i++ )
351  boxdims[3 + i] = high[i];
352  rval = mbImpl->tag_set_data( box_dims_tag(), &scd_set, 1, boxdims );
353  if( MB_SUCCESS != rval ) return rval;
354 
355  if( is_periodic )
356  {
357  rval = mbImpl->tag_set_data( box_periodic_tag(), &scd_set, 1, is_periodic );
358  if( MB_SUCCESS != rval ) return rval;
359  }
360 
361  return rval;
362 }
363 
364 Tag ScdInterface::box_periodic_tag( bool create_if_missing )
365 {
366  // Reset boxPeriodicTag in case it has been deleted (e.g. by Core::clean_up_failed_read)
367  if( boxPeriodicTag )
368  {
369  std::string tag_name;
370  if( MB_TAG_NOT_FOUND == mbImpl->tag_get_name( boxPeriodicTag, tag_name ) ) boxPeriodicTag = NULL;
371  }
372 
373  if( boxPeriodicTag || !create_if_missing ) return boxPeriodicTag;
374 
375  ErrorCode rval =
377  if( MB_SUCCESS != rval ) return 0;
378  return boxPeriodicTag;
379 }
380 
381 Tag ScdInterface::box_dims_tag( bool create_if_missing )
382 {
383  // Reset boxDimsTag in case it has been deleted (e.g. by clean_up_failed_read)
384  if( boxDimsTag )
385  {
386  std::string tag_name;
387  if( MB_TAG_NOT_FOUND == mbImpl->tag_get_name( boxDimsTag, tag_name ) ) boxDimsTag = NULL;
388  }
389 
390  if( boxDimsTag || !create_if_missing ) return boxDimsTag;
391 
393  if( MB_SUCCESS != rval ) return 0;
394  return boxDimsTag;
395 }
396 
397 Tag ScdInterface::global_box_dims_tag( bool create_if_missing )
398 {
399  // Reset globalBoxDimsTag in case it has been deleted (e.g. by Core::clean_up_failed_read)
400  if( globalBoxDimsTag )
401  {
402  std::string tag_name;
404  }
405 
406  if( globalBoxDimsTag || !create_if_missing ) return globalBoxDimsTag;
407 
408  ErrorCode rval =
410  if( MB_SUCCESS != rval ) return 0;
411  return globalBoxDimsTag;
412 }
413 
414 Tag ScdInterface::part_method_tag( bool create_if_missing )
415 {
416  // Reset partMethodTag in case it has been deleted (e.g. by Core::clean_up_failed_read)
417  if( partMethodTag )
418  {
419  std::string tag_name;
420  if( MB_TAG_NOT_FOUND == mbImpl->tag_get_name( partMethodTag, tag_name ) ) partMethodTag = NULL;
421  }
422 
423  if( partMethodTag || !create_if_missing ) return partMethodTag;
424 
425  ErrorCode rval =
427  if( MB_SUCCESS != rval ) return 0;
428  return partMethodTag;
429 }
430 
431 Tag ScdInterface::box_set_tag( bool create_if_missing )
432 {
433  // Reset boxSetTag in case it has been deleted (e.g. by Core::clean_up_failed_read)
434  if( boxSetTag )
435  {
436  std::string tag_name;
437  if( MB_TAG_NOT_FOUND == mbImpl->tag_get_name( boxSetTag, tag_name ) ) boxSetTag = NULL;
438  }
439 
440  if( boxSetTag || !create_if_missing ) return boxSetTag;
441 
442  ErrorCode rval = mbImpl->tag_get_handle( "__BOX_SET", sizeof( ScdBox* ), MB_TYPE_OPAQUE, boxSetTag,
444  if( MB_SUCCESS != rval ) return 0;
445  return boxSetTag;
446 }
447 
448 //! Remove the box from the list on ScdInterface
450 {
451  std::vector< ScdBox* >::iterator vit = std::find( scdBoxes.begin(), scdBoxes.end(), box );
452  if( vit != scdBoxes.end() )
453  {
454  scdBoxes.erase( vit );
455  return MB_SUCCESS;
456  }
457  else
458  return MB_FAILURE;
459 }
460 
461 //! Add the box to the list on ScdInterface
463 {
464  scdBoxes.push_back( box );
465  return MB_SUCCESS;
466 }
467 
468 ErrorCode ScdInterface::get_boxes( std::vector< ScdBox* >& boxes )
469 {
470  std::copy( scdBoxes.begin(), scdBoxes.end(), std::back_inserter( boxes ) );
471  return MB_SUCCESS;
472 }
473 
475  : scImpl( impl ), boxSet( bset ), vertDat( NULL ), elemSeq( NULL ), startVertex( 0 ), startElem( 0 )
476 {
477  for( int i = 0; i < 6; i++ )
478  boxDims[i] = 0;
479  for( int i = 0; i < 3; i++ )
480  locallyPeriodic[i] = false;
481  VertexSequence* vseq = dynamic_cast< VertexSequence* >( seq1 );
482  if( vseq ) vertDat = dynamic_cast< ScdVertexData* >( vseq->data() );
483  if( vertDat )
484  {
485  // retrieve the parametric space
486  for( int i = 0; i < 3; i++ )
487  {
488  boxDims[i] = vertDat->min_params()[i];
489  boxDims[3 + i] = vertDat->max_params()[i];
490  }
492  }
493  else if( impl->boxDimsTag )
494  {
495  // look for parametric space info on set
496  ErrorCode rval = impl->mbImpl->tag_get_data( impl->boxDimsTag, &bset, 1, boxDims );
497  if( MB_SUCCESS == rval )
498  {
499  Range verts;
500  impl->mbImpl->get_entities_by_dimension( bset, 0, verts );
501  if( !verts.empty() ) startVertex = *verts.begin();
502  }
503  }
504 
505  elemSeq = dynamic_cast< StructuredElementSeq* >( seq2 );
506  if( !elemSeq ) elemSeq = dynamic_cast< StructuredElementSeq* >( seq1 );
507 
508  if( elemSeq )
509  {
510  if( vertDat )
511  {
512  // check the parametric space to make sure it's consistent
513  assert( elemSeq->sdata()->min_params() == HomCoord( boxDims, 3 ) &&
514  ( elemSeq->sdata()->max_params() + HomCoord( 1, 1, 1 ) ) == HomCoord( boxDims, 3 ) );
515  }
516  else
517  {
518  // get the parametric space from the element sequence
519  for( int i = 0; i < 3; i++ )
520  {
521  boxDims[i] = elemSeq->sdata()->min_params()[i];
522  boxDims[3 + i] = elemSeq->sdata()->max_params()[i];
523  }
524  }
525 
527  }
528  else
529  {
530  Range elems;
532  bset, ( boxDims[2] == boxDims[5] ? ( boxDims[1] == boxDims[4] ? 1 : 2 ) : 3 ), elems );
533  if( !elems.empty() ) startElem = *elems.begin();
534  // call the following w/o looking at return value, since it doesn't really need to be there
535  if( impl->boxPeriodicTag ) impl->mbImpl->tag_get_data( impl->boxPeriodicTag, &bset, 1, locallyPeriodic );
536  }
537 
538  assert( vertDat || elemSeq || boxDims[0] != boxDims[3] || boxDims[1] != boxDims[4] || boxDims[2] != boxDims[5] );
539 
540  boxSize = HomCoord( boxDims + 3, 3 ) - HomCoord( boxDims, 3 ) + HomCoord( 1, 1, 1 );
541  boxSizeIJ = ( boxSize[1] ? boxSize[1] : 1 ) * boxSize[0];
542  boxSizeIM1 = boxSize[0] - ( locallyPeriodic[0] ? 0 : 1 );
543  boxSizeIJM1 = ( boxSize[1] ? ( boxSize[1] - ( locallyPeriodic[1] ? 0 : 1 ) ) : 1 ) * boxSizeIM1;
544 
545  scImpl->add_box( this );
546 }
547 
549 {
550  // Reset the tag on the set
551  if( boxSet )
552  {
553  // It is possible that the box set entity has been deleted (e.g. by
554  // Core::clean_up_failed_read)
555  Core* mbcore = dynamic_cast< Core* >( scImpl->mbImpl );
556  assert( mbcore != NULL );
557  if( mbcore->is_valid( boxSet ) )
558  {
559  ScdBox* tmp_ptr = NULL;
560  scImpl->mbImpl->tag_set_data( scImpl->box_set_tag(), &boxSet, 1, &tmp_ptr );
561  }
562  else
563  boxSet = 0;
564  }
565 
566  scImpl->remove_box( this );
567 }
568 
569 EntityHandle ScdBox::get_vertex_from_seq( int i, int j, int k ) const
570 {
571  assert( elemSeq );
572  return elemSeq->get_vertex( i, j, k );
573 }
574 
576 {
577  return ( startElem ? scImpl->mbImpl->dimension_from_handle( startElem ) : -1 );
578 }
579 
581  HomCoord from1,
582  HomCoord to1,
583  HomCoord from2,
584  HomCoord to2,
585  HomCoord from3,
586  HomCoord to3,
587  bool bb_input,
588  const HomCoord& bb_min,
589  const HomCoord& bb_max )
590 {
591  if( !vbox->vertDat ) return MB_FAILURE;
592  ScdVertexData* dum_data = dynamic_cast< ScdVertexData* >( vbox->vertDat );
593  ErrorCode rval =
594  elemSeq->sdata()->add_vsequence( dum_data, from1, to1, from2, to2, from3, to3, bb_input, bb_min, bb_max );
595  return rval;
596 }
597 
599 {
600  return elemSeq->boundary_complete();
601 }
602 
603 ErrorCode ScdBox::get_coordinate_arrays( double*& xc, double*& yc, double*& zc )
604 {
605  if( !vertDat ) return MB_FAILURE;
606 
607  xc = reinterpret_cast< double* >( vertDat->get_sequence_data( 0 ) );
608  yc = reinterpret_cast< double* >( vertDat->get_sequence_data( 1 ) );
609  zc = reinterpret_cast< double* >( vertDat->get_sequence_data( 2 ) );
610  return MB_SUCCESS;
611 }
612 
613 ErrorCode ScdBox::get_coordinate_arrays( const double*& xc, const double*& yc, const double*& zc ) const
614 {
615  if( !vertDat ) return MB_FAILURE;
616  xc = reinterpret_cast< const double* >( vertDat->get_sequence_data( 0 ) );
617  yc = reinterpret_cast< const double* >( vertDat->get_sequence_data( 1 ) );
618  zc = reinterpret_cast< const double* >( vertDat->get_sequence_data( 2 ) );
619  return MB_SUCCESS;
620 }
621 
623 {
624  vertDat = vert_dt;
625  return MB_SUCCESS;
626 }
627 
629 {
630  elemSeq = dynamic_cast< StructuredElementSeq* >( elem_sq );
632 
633  if( locallyPeriodic[0] ) boxSizeIM1 = boxSize[0] - ( locallyPeriodic[0] ? 0 : 1 );
634  if( locallyPeriodic[0] || locallyPeriodic[1] )
635  boxSizeIJM1 = ( boxSize[1] ? ( boxSize[1] - ( locallyPeriodic[1] ? 0 : 1 ) ) : 1 ) * boxSizeIM1;
636 
637  return ( elemSeq ? MB_SUCCESS : MB_FAILURE );
638 }
639 
641 {
642  // check first whether this is an intermediate entity, so we know what to do
643  int dimension = box_dimension();
644  int this_dim = scImpl->impl()->dimension_from_handle( ent );
645 
646  if( ( 0 == this_dim && !vertDat ) || ( this_dim && this_dim == dimension ) )
647  {
648  assert( elemSeq );
649  return elemSeq->get_params( ent, ijkd[0], ijkd[1], ijkd[2] );
650  }
651 
652  else if( !this_dim && vertDat )
653  return vertDat->get_params( ent, ijkd[0], ijkd[1], ijkd[2] );
654 
655  else
656  return MB_NOT_IMPLEMENTED;
657 }
658 
659 //! Get the entity of specified dimension adjacent to parametric element
660 /**
661  * \param dim Dimension of adjacent entity being requested
662  * \param i Parametric coordinates of cell being evaluated
663  * \param j Parametric coordinates of cell being evaluated
664  * \param k Parametric coordinates of cell being evaluated
665  * \param dir Direction (0, 1, or 2), for getting adjacent edges (2d, 3d) or faces (3d)
666  * \param ent EntityHandle of adjacent entity
667  * \param create_if_missing If true, creates the entity if it doesn't already exist
668  */
670  int i,
671  int j,
672  int k,
673  int dir,
674  EntityHandle& ent,
675  bool create_if_missing ) const
676 {
677  // describe connectivity of sub-element in static array
678  // subconnect[dim-1][dir][numv][ijk] where dimensions are:
679  // [dim-1]: dim=1 or 2, so this is 0 or 1
680  // [dir]: one of 0..2, for ijk directions in a hex
681  // [numv]: number of vertices describing sub entity = 2*dim <= 4
682  // [ijk]: 3 values for i, j, k
683  int subconnect[2][3][4][3] = { { { { 0, 0, 0 }, { 1, 0, 0 }, { -1, -1, -1 }, { -1, -1, -1 } }, // i edge
684  { { 0, 0, 0 }, { 0, 1, 0 }, { -1, -1, -1 }, { -1, -1, -1 } }, // j edge
685  { { 0, 0, 0 }, { 0, 0, 1 }, { -1, -1, -1 }, { -1, -1, -1 } } }, // k edge
686 
687  { { { 0, 0, 0 }, { 0, 1, 0 }, { 0, 1, 1 }, { 0, 0, 1 } }, // i face
688  { { 0, 0, 0 }, { 1, 0, 0 }, { 1, 0, 1 }, { 0, 0, 1 } }, // j face
689  { { 0, 0, 0 }, { 1, 0, 0 }, { 1, 1, 0 }, { 0, 1, 0 } } } }; // k face
690 
691  // check proper input dimensions and lower bound
692  if( dim < 1 || dim > 2 || i < boxDims[0] || j < boxDims[1] || k < boxDims[2] ) return MB_FAILURE;
693 
694  // now check upper bound; parameters must be <= upper corner, since edges/faces
695  // follow element parameterization, not vertex parameterization
696  else if( ( boxDims[3] != boxDims[0] && i > ( locallyPeriodic[0] ? boxDims[3] + 1 : boxDims[3] ) ) ||
697  ( boxDims[4] != boxDims[1] && j > ( locallyPeriodic[1] ? boxDims[4] + 1 : boxDims[4] ) ) ||
698  ( boxDims[5] != boxDims[2] && k > boxDims[5] ) )
699  return MB_FAILURE;
700 
701  // get the vertices making up this entity
702  EntityHandle verts[4];
703  for( int ind = 0; ind < 2 * dim; ind++ )
704  {
705  int i1 = i + subconnect[dim - 1][dir][ind][0];
706  int j1 = j + subconnect[dim - 1][dir][ind][1];
707  // if periodic in i and i1 is boxDims[3]+1, wrap around
708  if( locallyPeriodic[0] && i1 == boxDims[3] + 1 ) i1 = boxDims[0];
709  // if periodic in j and j1 is boxDims[4]+1, wrap around
710  if( locallyPeriodic[1] && j1 == boxDims[4] + 1 ) j1 = boxDims[1];
711  verts[ind] = get_vertex( i1, j1, k + subconnect[dim - 1][dir][ind][2] );
712  if( !verts[ind] ) return MB_FAILURE;
713  }
714 
715  Range ents;
716  ErrorCode rval = scImpl->impl()->get_adjacencies( verts, 2 * dim, dim, false, ents );
717  if( MB_SUCCESS != rval ) return rval;
718 
719  if( ents.size() > 1 )
720  return MB_FAILURE;
721 
722  else if( ents.size() == 1 )
723  {
724  ent = *ents.begin();
725  }
726  else if( create_if_missing )
727  rval = scImpl->impl()->create_element( ( 1 == dim ? MBEDGE : MBQUAD ), verts, 2 * dim, ent );
728 
729  return rval;
730 }
731 
732 #ifndef MOAB_HAVE_MPI
734 {
735  return MB_FAILURE;
736 #else
738 {
739  EntityHandle seth = box->box_set();
740 
741  // check the # ents in the box against the num in the set, to make sure it's only 1 box;
742  // reuse tmp_range
743  Range tmp_range;
744  ErrorCode rval = mbImpl->get_entities_by_dimension( seth, box->box_dimension(), tmp_range );
745  if( MB_SUCCESS != rval ) return rval;
746  if( box->num_elements() != (int)tmp_range.size() ) return MB_FAILURE;
747 
748  const int* gdims = box->par_data().gDims;
749  if( ( gdims[0] == gdims[3] && gdims[1] == gdims[4] && gdims[2] == gdims[5] ) || -1 == box->par_data().partMethod )
750  return MB_FAILURE;
751 
752  // ok, we have a partitioned box; get the vertices shared with other processors
753  std::vector< int > procs, offsets, shared_indices;
754  rval = get_shared_vertices( pcomm, box, procs, offsets, shared_indices );
755  if( MB_SUCCESS != rval ) return rval;
756 
757  // post receives for start handles once we know how many to look for
758  std::vector< MPI_Request > recv_reqs( procs.size(), MPI_REQUEST_NULL ), send_reqs( procs.size(), MPI_REQUEST_NULL );
759  std::vector< EntityHandle > rhandles( 4 * procs.size() ), shandles( 4 );
760  for( unsigned int i = 0; i < procs.size(); i++ )
761  {
762  int success = MPI_Irecv( (void*)&rhandles[4 * i], 4 * sizeof( EntityHandle ), MPI_UNSIGNED_CHAR, procs[i], 1,
763  pcomm->proc_config().proc_comm(), &recv_reqs[i] );
764  if( success != MPI_SUCCESS ) return MB_FAILURE;
765  }
766 
767  // send our own start handles
768  shandles[0] = box->start_vertex();
769  shandles[1] = 0;
770  if( box->box_dimension() == 1 )
771  {
772  shandles[1] = box->start_element();
773  shandles[2] = 0;
774  shandles[3] = 0;
775  }
776  else if( box->box_dimension() == 2 )
777  {
778  shandles[2] = box->start_element();
779  shandles[3] = 0;
780  }
781  else
782  {
783  shandles[2] = 0;
784  shandles[3] = box->start_element();
785  }
786  for( unsigned int i = 0; i < procs.size(); i++ )
787  {
788  int success = MPI_Isend( (void*)&shandles[0], 4 * sizeof( EntityHandle ), MPI_UNSIGNED_CHAR, procs[i], 1,
789  pcomm->proc_config().proc_comm(), &send_reqs[i] );
790  if( success != MPI_SUCCESS ) return MB_FAILURE;
791  }
792 
793  // receive start handles and save info to a tuple list
794  int incoming = procs.size();
795  int p, j, k;
796  MPI_Status status;
797  TupleList shared_data;
798  shared_data.initialize( 1, 0, 2, 0, shared_indices.size() / 2 );
799  shared_data.enableWriteAccess();
800 
801  j = 0;
802  k = 0;
803  while( incoming )
804  {
805  int success = MPI_Waitany( procs.size(), &recv_reqs[0], &p, &status );
806  if( MPI_SUCCESS != success ) return MB_FAILURE;
807  unsigned int num_indices = ( offsets[p + 1] - offsets[p] ) / 2;
808  int *lh = &shared_indices[offsets[p]], *rh = lh + num_indices;
809  for( unsigned int i = 0; i < num_indices; i++ )
810  {
811  shared_data.vi_wr[j++] = procs[p];
812  shared_data.vul_wr[k++] = shandles[0] + lh[i];
813  shared_data.vul_wr[k++] = rhandles[4 * p] + rh[i];
814  shared_data.inc_n();
815  }
816  incoming--;
817  }
818 
819  // still need to wait for the send requests
820  std::vector< MPI_Status > mult_status( procs.size() );
821  int success = MPI_Waitall( procs.size(), &send_reqs[0], &mult_status[0] );
822  if( MPI_SUCCESS != success )
823  {
824  MB_SET_ERR( MB_FAILURE, "Failed in waitall in ScdInterface::tag_shared_vertices" );
825  }
826  // sort by local handle
827  TupleList::buffer sort_buffer;
828  sort_buffer.buffer_init( shared_indices.size() / 2 );
829  shared_data.sort( 1, &sort_buffer );
830  sort_buffer.reset();
831 
832  // process into sharing data
833  std::map< std::vector< int >, std::vector< EntityHandle > > proc_nvecs;
834  Range dum;
835  rval = pcomm->tag_shared_verts( shared_data, proc_nvecs, dum, 0 );
836  if( MB_SUCCESS != rval ) return rval;
837 
838  // create interface sets
839  rval = pcomm->create_interface_sets( proc_nvecs );
840  if( MB_SUCCESS != rval ) return rval;
841 
842  // add the box to the PComm's partitionSets
843  pcomm->partition_sets().insert( box->box_set() );
844 
845  // make sure buffers are allocated for communicating procs
846  for( std::vector< int >::iterator pit = procs.begin(); pit != procs.end(); ++pit )
847  pcomm->get_buffers( *pit );
848 
849  if( pcomm->get_debug_verbosity() > 1 ) pcomm->list_entities( NULL, 1 );
850 
851 #ifndef NDEBUG
852  rval = pcomm->check_all_shared_handles();
853  if( MB_SUCCESS != rval ) return rval;
854 #endif
855 
856  return MB_SUCCESS;
857 
858 #endif
859 }
860 
862  int pfrom,
863  const int* const gdims,
864  const int* const gperiodic,
865  const int* const dijk,
866  int& pto,
867  int* rdims,
868  int* facedims,
869  int* across_bdy )
870 {
871  if( dijk[0] != 0 )
872  {
873  pto = -1;
874  return MB_SUCCESS;
875  }
876 
877  pto = -1;
878  across_bdy[0] = across_bdy[1] = across_bdy[2] = 0;
879 
880  int ldims[6], pijk[3], lperiodic[3];
881  ErrorCode rval = compute_partition_alljkbal( np, pfrom, gdims, gperiodic, ldims, lperiodic, pijk );
882  if( MB_SUCCESS != rval ) return rval;
883  assert( pijk[1] * pijk[2] == np );
884  pto = -1;
885  bool bot_j = pfrom< pijk[2], top_j = pfrom > np - pijk[2];
886  if( ( 1 == pijk[2] && dijk[2] ) || // 1d in j means no neighbors with dk != 0
887  ( !( pfrom % pijk[2] ) && -1 == dijk[2] ) || // at -k bdy
888  ( pfrom % pijk[2] == pijk[2] - 1 && 1 == dijk[2] ) || // at +k bdy
889  ( pfrom < pijk[2] && -1 == dijk[1] && !gperiodic[1] ) || // down and not periodic
890  ( pfrom >= np - pijk[2] && 1 == dijk[1] && !gperiodic[1] ) ) // up and not periodic
891  return MB_SUCCESS;
892 
893  pto = pfrom;
894  std::copy( ldims, ldims + 6, rdims );
895  std::copy( ldims, ldims + 6, facedims );
896 
897  if( 0 != dijk[1] )
898  {
899  pto = ( pto + dijk[1] * pijk[2] + np ) % np;
900  assert( pto >= 0 && pto < np );
901  int dj = ( gdims[4] - gdims[1] ) / pijk[1], extra = ( gdims[4] - gdims[1] ) % pijk[1];
902  if( -1 == dijk[1] )
903  {
904  facedims[4] = facedims[1];
905  if( bot_j )
906  {
907  // going across periodic lower bdy in j
908  rdims[4] = gdims[4];
909  across_bdy[1] = -1;
910  }
911  else
912  {
913  rdims[4] = ldims[1];
914  }
915  rdims[1] = rdims[4] - dj;
916  if( pto < extra ) rdims[1]--;
917  }
918  else
919  {
920  if( pfrom > np - pijk[2] ) facedims[4] = gdims[1];
921  facedims[1] = facedims[4];
922  if( top_j )
923  {
924  // going across periodic upper bdy in j
925  rdims[1] = gdims[1];
926  across_bdy[1] = 1;
927  }
928  else
929  {
930  rdims[1] = ldims[4];
931  }
932  rdims[4] = rdims[1] + dj;
933  if( pto < extra ) rdims[4]++;
934  }
935  }
936  if( 0 != dijk[2] )
937  {
938  pto = ( pto + dijk[2] ) % np;
939  assert( pto >= 0 && pto < np );
940  facedims[2] = facedims[5] = ( -1 == dijk[2] ? facedims[2] : facedims[5] );
941  int dk = ( gdims[5] - gdims[2] ) / pijk[2];
942  if( -1 == dijk[2] )
943  {
944  facedims[5] = facedims[2];
945  rdims[5] = ldims[2];
946  rdims[2] = rdims[5] - dk; // never any kextra for alljkbal
947  }
948  else
949  {
950  facedims[2] = facedims[5];
951  rdims[2] = ldims[5];
952  rdims[5] = rdims[2] + dk; // never any kextra for alljkbal
953  }
954  }
955 
956  assert( -1 == pto || ( rdims[0] >= gdims[0] && rdims[3] <= gdims[3] ) );
957  assert( -1 == pto || ( rdims[1] >= gdims[1] && ( rdims[4] <= gdims[4] || ( across_bdy[1] && bot_j ) ) ) );
958  assert( -1 == pto || ( rdims[2] >= gdims[2] && rdims[5] <= gdims[5] ) );
959  assert( -1 == pto || ( ( facedims[0] >= rdims[0] ||
960  ( gperiodic[0] && rdims[3] == gdims[3] + 1 && facedims[0] == gdims[0] ) ) ) );
961  assert( -1 == pto || ( facedims[3] <= rdims[3] ) );
962  assert( -1 == pto || ( ( facedims[1] >= rdims[1] ||
963  ( gperiodic[1] && rdims[4] == gdims[4] + 1 && facedims[1] == gdims[1] ) ) ) );
964  assert( -1 == pto || ( facedims[4] <= rdims[4] ) );
965  assert( -1 == pto || ( facedims[2] >= rdims[2] ) );
966  assert( -1 == pto || ( facedims[5] <= rdims[5] ) );
967  assert( -1 == pto || ( facedims[0] >= ldims[0] ) );
968  assert( -1 == pto || ( facedims[3] <= ldims[3] ) );
969  assert( -1 == pto || ( facedims[1] >= ldims[1] ) );
970  assert( -1 == pto || ( facedims[4] <= ldims[4] ) );
971  assert( -1 == pto || ( facedims[2] >= ldims[2] ) );
972  assert( -1 == pto || ( facedims[5] <= ldims[5] ) );
973 
974  return MB_SUCCESS;
975 }
976 
978  int pfrom,
979  const int* const gdims,
980  const int* const gperiodic,
981  const int* const dijk,
982  int& pto,
983  int* rdims,
984  int* facedims,
985  int* across_bdy )
986 {
987  if( dijk[2] != 0 )
988  {
989  // for sqij, there is no k neighbor, ever
990  pto = -1;
991  return MB_SUCCESS;
992  }
993 
994  pto = -1;
995  across_bdy[0] = across_bdy[1] = across_bdy[2] = 0;
996  int lperiodic[3], pijk[3], ldims[6];
997  ErrorCode rval = compute_partition_sqij( np, pfrom, gdims, gperiodic, ldims, lperiodic, pijk );
998  if( MB_SUCCESS != rval ) return rval;
999  assert( pijk[0] * pijk[1] == np );
1000  pto = -1;
1001  bool top_i = 0, top_j = 0, bot_i = 0, bot_j = 0;
1002  int ni = pfrom % pijk[0], nj = pfrom / pijk[0]; // row / column number of me
1003  if( ni == pijk[0] - 1 ) top_i = 1;
1004  if( nj == pijk[1] - 1 ) top_j = 1;
1005  if( !ni ) bot_i = 1;
1006  if( !nj ) bot_j = 1;
1007  if( ( !gperiodic[0] && bot_i && -1 == dijk[0] ) || // left and not periodic
1008  ( !gperiodic[0] && top_i && 1 == dijk[0] ) || // right and not periodic
1009  ( !gperiodic[1] && bot_j && -1 == dijk[1] ) || // bottom and not periodic
1010  ( !gperiodic[1] && top_j && 1 == dijk[1] ) ) // top and not periodic
1011  return MB_SUCCESS;
1012 
1013  std::copy( ldims, ldims + 6, facedims );
1014  std::copy( ldims, ldims + 6, rdims );
1015  pto = pfrom;
1016  int j = gdims[4] - gdims[1], dj = j / pijk[1], jextra = ( gdims[4] - gdims[1] ) % dj, i = gdims[3] - gdims[0],
1017  di = i / pijk[0], iextra = ( gdims[3] - gdims[0] ) % di;
1018 
1019  if( 0 != dijk[0] )
1020  {
1021  pto = ( ni + dijk[0] + pijk[0] ) % pijk[0]; // get pto's ni value
1022  pto = nj * pijk[0] + pto; // then convert to pto
1023  assert( pto >= 0 && pto < np );
1024  if( -1 == dijk[0] )
1025  {
1026  facedims[3] = facedims[0];
1027  if( bot_i )
1028  {
1029  // going across lower periodic bdy in i
1030  across_bdy[0] = -1;
1031  rdims[3] = gdims[3] + 1; // +1 because ldims[3] on remote proc is gdims[3]+1
1032  rdims[0] = rdims[3] - di - 1; // -1 to account for rdims[3] being one larger
1033  }
1034  else
1035  {
1036  rdims[3] = ldims[0];
1037  rdims[0] = rdims[3] - di;
1038  }
1039 
1040  if( pto % pijk[0] < iextra ) rdims[0]--;
1041  }
1042  else
1043  {
1044  if( top_i )
1045  {
1046  // going across lower periodic bdy in i
1047  facedims[3] = gdims[0];
1048  across_bdy[0] = 1;
1049  }
1050  facedims[0] = facedims[3];
1051  rdims[0] = ( top_i ? gdims[0] : ldims[3] );
1052  rdims[3] = rdims[0] + di;
1053  if( pto % pijk[0] < iextra ) rdims[3]++;
1054  if( gperiodic[0] && ni == pijk[0] - 2 ) rdims[3]++; // remote proc is top_i and periodic
1055  }
1056  }
1057  if( 0 != dijk[1] )
1058  {
1059  pto = ( pto + dijk[1] * pijk[0] + np ) % np;
1060  assert( pto >= 0 && pto < np );
1061  if( -1 == dijk[1] )
1062  {
1063  facedims[4] = facedims[1];
1064  if( bot_j )
1065  {
1066  // going across lower periodic bdy in j
1067  rdims[4] = gdims[4] + 1; // +1 because ldims[4] on remote proc is gdims[4]+1
1068  rdims[1] = rdims[4] - dj - 1; // -1 to account for gdims[4] being one larger
1069  across_bdy[1] = -1;
1070  }
1071  else
1072  {
1073  rdims[4] = ldims[1];
1074  rdims[1] = rdims[4] - dj;
1075  }
1076  if( pto / pijk[0] < jextra ) rdims[1]--;
1077  }
1078  else
1079  {
1080  if( top_j )
1081  {
1082  // going across lower periodic bdy in j
1083  facedims[4] = gdims[1];
1084  rdims[1] = gdims[1];
1085  across_bdy[1] = 1;
1086  }
1087  else
1088  {
1089  rdims[1] = ldims[4];
1090  }
1091  facedims[1] = facedims[4];
1092  rdims[4] = rdims[1] + dj;
1093  if( nj + 1 < jextra ) rdims[4]++;
1094  if( gperiodic[1] && nj == pijk[1] - 2 ) rdims[4]++; // remote proc is top_j and periodic
1095  }
1096  }
1097 
1098  // rdims within gdims
1099  assert( -1 == pto || ( rdims[0] >= gdims[0] &&
1100  ( rdims[3] <= gdims[3] + ( gperiodic[0] && pto % pijk[0] == pijk[0] - 1 ? 1 : 0 ) ) ) );
1101  assert( -1 == pto || ( rdims[1] >= gdims[1] &&
1102  ( rdims[4] <= gdims[4] + ( gperiodic[1] && pto / pijk[0] == pijk[1] - 1 ? 1 : 0 ) ) ) );
1103  assert( -1 == pto || ( rdims[2] >= gdims[2] && rdims[5] <= gdims[5] ) );
1104  // facedims within rdims
1105  assert( -1 == pto || ( ( facedims[0] >= rdims[0] ||
1106  ( gperiodic[0] && pto % pijk[0] == pijk[0] - 1 && facedims[0] == gdims[0] ) ) ) );
1107  assert( -1 == pto || ( facedims[3] <= rdims[3] ) );
1108  assert( -1 == pto || ( ( facedims[1] >= rdims[1] ||
1109  ( gperiodic[1] && pto / pijk[0] == pijk[1] - 1 && facedims[1] == gdims[1] ) ) ) );
1110  assert( -1 == pto || ( facedims[4] <= rdims[4] ) );
1111  assert( -1 == pto || ( facedims[2] >= rdims[2] && facedims[5] <= rdims[5] ) );
1112  // facedims within ldims
1113  assert( -1 == pto || ( ( facedims[0] >= ldims[0] || ( top_i && facedims[0] == gdims[0] ) ) ) );
1114  assert( -1 == pto || ( facedims[3] <= ldims[3] ) );
1115  assert( -1 == pto || ( ( facedims[1] >= ldims[1] || ( gperiodic[1] && top_j && facedims[1] == gdims[1] ) ) ) );
1116  assert( -1 == pto || ( facedims[4] <= ldims[4] ) );
1117  assert( -1 == pto || ( facedims[2] >= ldims[2] && facedims[5] <= ldims[5] ) );
1118 
1119  return MB_SUCCESS;
1120 }
1121 
1123  int pfrom,
1124  const int* const gdims,
1125  const int* const gperiodic,
1126  const int* const dijk,
1127  int& pto,
1128  int* rdims,
1129  int* facedims,
1130  int* across_bdy )
1131 {
1132  if( dijk[0] != 0 )
1133  {
1134  pto = -1;
1135  return MB_SUCCESS;
1136  }
1137 
1138  pto = -1;
1139  across_bdy[0] = across_bdy[1] = across_bdy[2] = 0;
1140  int pijk[3], lperiodic[3], ldims[6];
1141  ErrorCode rval = compute_partition_sqjk( np, pfrom, gdims, gperiodic, ldims, lperiodic, pijk );
1142  if( MB_SUCCESS != rval ) return rval;
1143  assert( pijk[1] * pijk[2] == np );
1144  pto = -1;
1145  bool top_j = 0, top_k = 0, bot_j = 0, bot_k = 0;
1146  int nj = pfrom % pijk[1], nk = pfrom / pijk[1];
1147  if( nj == pijk[1] - 1 ) top_j = 1;
1148  if( nk == pijk[2] - 1 ) top_k = 1;
1149  if( !nj ) bot_j = 1;
1150  if( !nk ) bot_k = 1;
1151  if( ( !gperiodic[1] && bot_j && -1 == dijk[1] ) || // down and not periodic
1152  ( !gperiodic[1] && top_j && 1 == dijk[1] ) || // up and not periodic
1153  ( bot_k && -1 == dijk[2] ) || // k- bdy
1154  ( top_k && 1 == dijk[2] ) ) // k+ bdy
1155  return MB_SUCCESS;
1156 
1157  std::copy( ldims, ldims + 6, facedims );
1158  std::copy( ldims, ldims + 6, rdims );
1159  pto = pfrom;
1160  int dj = ( gdims[4] - gdims[1] ) / pijk[1], jextra = ( gdims[4] - gdims[1] ) % dj,
1161  dk = ( gdims[5] == gdims[2] ? 0 : ( gdims[5] - gdims[2] ) / pijk[2] ),
1162  kextra = ( gdims[5] - gdims[2] ) - dk * pijk[2];
1163  assert( ( dj * pijk[1] + jextra == ( gdims[4] - gdims[1] ) ) &&
1164  ( dk * pijk[2] + kextra == ( gdims[5] - gdims[2] ) ) );
1165  if( 0 != dijk[1] )
1166  {
1167  pto = ( nj + dijk[1] + pijk[1] ) % pijk[1]; // get pto's ni value
1168  pto = nk * pijk[1] + pto; // then convert to pto
1169  assert( pto >= 0 && pto < np );
1170  if( -1 == dijk[1] )
1171  {
1172  facedims[4] = facedims[1];
1173  if( bot_j )
1174  {
1175  // going across lower periodic bdy in j
1176  rdims[4] = gdims[4] + 1; // +1 because ldims[4] on remote proc is gdims[4]+1
1177  across_bdy[1] = -1;
1178  }
1179  else
1180  {
1181  rdims[4] = ldims[1];
1182  }
1183  rdims[1] = rdims[4] - dj;
1184  if( nj < jextra ) rdims[1]--;
1185  }
1186  else
1187  {
1188  if( top_j )
1189  {
1190  // going across upper periodic bdy in j
1191  rdims[1] = gdims[1];
1192  facedims[4] = gdims[1];
1193  across_bdy[1] = 1;
1194  }
1195  else
1196  {
1197  rdims[1] = ldims[4];
1198  }
1199  facedims[1] = facedims[4];
1200  rdims[4] = rdims[1] + dj;
1201  if( nj < jextra ) rdims[4]++;
1202  if( gperiodic[1] && nj == dijk[1] - 2 ) rdims[4]++; // +1 because next proc is on periodic bdy
1203  }
1204  }
1205  if( 0 != dijk[2] )
1206  {
1207  pto = ( pto + dijk[2] * pijk[1] + np ) % np;
1208  assert( pto >= 0 && pto < np );
1209  if( -1 == dijk[2] )
1210  {
1211  facedims[5] = facedims[2];
1212  rdims[5] = ldims[2];
1213  rdims[2] -= dk;
1214  if( pto / pijk[1] < kextra ) rdims[2]--;
1215  }
1216  else
1217  {
1218  facedims[2] = facedims[5];
1219  rdims[2] = ldims[5];
1220  rdims[5] += dk;
1221  if( pto / pijk[1] < kextra ) rdims[5]++;
1222  }
1223  }
1224 
1225  assert( -1 == pto || ( rdims[0] >= gdims[0] && rdims[3] <= gdims[3] ) );
1226  assert( -1 == pto || ( rdims[1] >= gdims[1] && ( rdims[4] <= gdims[4] || ( across_bdy[1] && bot_j ) ) ) );
1227  assert( -1 == pto || ( rdims[2] >= gdims[2] && rdims[5] <= gdims[5] ) );
1228  assert( -1 == pto || ( facedims[0] >= rdims[0] && facedims[3] <= rdims[3] ) );
1229  assert( -1 == pto ||
1230  ( ( facedims[1] >= rdims[1] || ( gperiodic[1] && rdims[4] == gdims[4] && facedims[1] == gdims[1] ) ) ) );
1231  assert( -1 == pto || ( facedims[4] <= rdims[4] ) );
1232  assert( -1 == pto || ( facedims[2] >= rdims[2] && facedims[5] <= rdims[5] ) );
1233  assert( -1 == pto || ( facedims[0] >= ldims[0] && facedims[3] <= ldims[3] ) );
1234  assert( -1 == pto || ( facedims[1] >= ldims[1] && facedims[4] <= ldims[4] ) );
1235  assert( -1 == pto || ( facedims[2] >= ldims[2] && facedims[5] <= ldims[5] ) );
1236 
1237  return MB_SUCCESS;
1238 }
1239 
1241  int pfrom,
1242  const int* const gdims,
1243  const int* const gperiodic,
1244  const int* const dijk,
1245  int& pto,
1246  int* rdims,
1247  int* facedims,
1248  int* across_bdy )
1249 {
1250  if( gperiodic[0] || gperiodic[1] || gperiodic[2] ) return MB_FAILURE;
1251 
1252  pto = -1;
1253  across_bdy[0] = across_bdy[1] = across_bdy[2] = 0;
1254  int pijk[3], lperiodic[3], ldims[6];
1255  ErrorCode rval = compute_partition_sqijk( np, pfrom, gdims, gperiodic, ldims, lperiodic, pijk );
1256  if( MB_SUCCESS != rval ) return rval;
1257  assert( pijk[0] * pijk[1] * pijk[2] == np );
1258  pto = -1;
1259  bool top[3] = { false, false, false }, bot[3] = { false, false, false };
1260  // nijk: rank in i/j/k direction
1261  int nijk[3] = { pfrom % pijk[0], ( pfrom % ( pijk[0] * pijk[1] ) ) / pijk[0], pfrom / ( pijk[0] * pijk[1] ) };
1262 
1263  for( int i = 0; i < 3; i++ )
1264  {
1265  if( nijk[i] == pijk[i] - 1 ) top[i] = true;
1266  if( !nijk[i] ) bot[i] = true;
1267  if( ( !gperiodic[i] && bot[i] && -1 == dijk[i] ) || // downward && not periodic
1268  ( !gperiodic[i] && top[i] && 1 == dijk[i] ) ) // upward && not periodic
1269  return MB_SUCCESS;
1270  }
1271 
1272  std::copy( ldims, ldims + 6, facedims );
1273  std::copy( ldims, ldims + 6, rdims );
1274  pto = pfrom;
1275  int delijk[3], extra[3];
1276  // nijk_to: rank of pto in i/j/k direction
1277  int nijk_to[3];
1278  for( int i = 0; i < 3; i++ )
1279  {
1280  delijk[i] = ( gdims[i + 3] == gdims[i] ? 0 : ( gdims[i + 3] - gdims[i] ) / pijk[i] );
1281  extra[i] = ( gdims[i + 3] - gdims[i] ) % delijk[i];
1282  nijk_to[i] = ( nijk[i] + dijk[i] + pijk[i] ) % pijk[i];
1283  }
1284  pto = nijk_to[2] * pijk[0] * pijk[1] + nijk_to[1] * pijk[0] + nijk_to[0];
1285  assert( pto >= 0 && pto < np );
1286  for( int i = 0; i < 3; i++ )
1287  {
1288  if( 0 != dijk[i] )
1289  {
1290  if( -1 == dijk[i] )
1291  {
1292  facedims[i + 3] = facedims[i];
1293  if( bot[i] )
1294  {
1295  // going across lower periodic bdy in i
1296  rdims[i + 3] = gdims[i + 3] + 1; // +1 because ldims[4] on remote proc is gdims[4]+1
1297  across_bdy[i] = -1;
1298  }
1299  else
1300  {
1301  rdims[i + 3] = ldims[i];
1302  }
1303  rdims[i] = rdims[i + 3] - delijk[i];
1304  if( nijk[i] < extra[i] ) rdims[i]--;
1305  }
1306  else
1307  {
1308  if( top[i] )
1309  {
1310  // going across upper periodic bdy in i
1311  rdims[i] = gdims[i];
1312  facedims[i + 3] = gdims[i];
1313  across_bdy[i] = 1;
1314  }
1315  else
1316  {
1317  rdims[i] = ldims[i + 3];
1318  }
1319  facedims[i] = facedims[i + 3];
1320  rdims[i + 3] = rdims[i] + delijk[i];
1321  if( nijk[i] < extra[i] ) rdims[i + 3]++;
1322  if( gperiodic[i] && nijk[i] == dijk[i] - 2 ) rdims[i + 3]++; // +1 because next proc is on periodic bdy
1323  }
1324  }
1325  }
1326 
1327  assert( -1 != pto );
1328 #ifndef NDEBUG
1329  for( int i = 0; i < 3; i++ )
1330  {
1331  assert( ( rdims[i] >= gdims[i] && ( rdims[i + 3] <= gdims[i + 3] || ( across_bdy[i] && bot[i] ) ) ) );
1332  assert( ( ( facedims[i] >= rdims[i] ||
1333  ( gperiodic[i] && rdims[i + 3] == gdims[i + 3] && facedims[i] == gdims[i] ) ) ) );
1334  assert( ( facedims[i] >= ldims[i] && facedims[i + 3] <= ldims[i + 3] ) );
1335  }
1336 #endif
1337 
1338  return MB_SUCCESS;
1339 }
1340 
1342  int pfrom,
1343  const int* const gdims,
1344  const int* const gperiodic,
1345  const int* const dijk,
1346  int& pto,
1347  int* rdims,
1348  int* facedims,
1349  int* across_bdy )
1350 {
1351  ErrorCode rval = MB_SUCCESS;
1352  pto = -1;
1353  if( np == 1 ) return MB_SUCCESS;
1354 
1355  int pijk[3], lperiodic[3], ldims[6];
1356  rval = compute_partition_alljorkori( np, pfrom, gdims, gperiodic, ldims, lperiodic, pijk );
1357  if( MB_SUCCESS != rval ) return rval;
1358 
1359  int ind = -1;
1360  across_bdy[0] = across_bdy[1] = across_bdy[2] = 0;
1361 
1362  for( int i = 0; i < 3; i++ )
1363  {
1364  if( pijk[i] > 1 )
1365  {
1366  ind = i;
1367  break;
1368  }
1369  }
1370 
1371  assert( -1 < ind );
1372 
1373  if( !dijk[ind] )
1374  // no neighbor, pto is already -1, return
1375  return MB_SUCCESS;
1376 
1377  bool is_periodic = ( ( gperiodic[0] && ind == 0 ) || ( gperiodic[1] && ind == 1 ) );
1378  if( dijk[( ind + 1 ) % 3] || dijk[( ind + 2 ) % 3] || // stepping in either other two directions
1379  ( !is_periodic && ldims[ind] == gdims[ind] && dijk[ind] == -1 ) || // lower side and going lower
1380  ( !is_periodic && ldims[3 + ind] >= gdims[3 + ind] &&
1381  dijk[ind] == 1 ) ) // not >= because ldims is only > gdims when periodic;
1382  // higher side and going higher
1383  return MB_SUCCESS;
1384 
1385  std::copy( ldims, ldims + 6, facedims );
1386  std::copy( ldims, ldims + 6, rdims );
1387 
1388  int dind = ( gdims[ind + 3] - gdims[ind] ) / np;
1389  int extra = ( gdims[ind + 3] - gdims[ind] ) % np;
1390  if( -1 == dijk[ind] && pfrom )
1391  {
1392  // actual left neighbor
1393  pto = pfrom - 1; // no need for %np, because pfrom > 0
1394  facedims[ind + 3] = facedims[ind];
1395  rdims[ind + 3] = ldims[ind];
1396  rdims[ind] = rdims[ind + 3] - dind - ( pto < extra ? 1 : 0 );
1397  }
1398  else if( 1 == dijk[ind] && pfrom < np - 1 )
1399  {
1400  // actual right neighbor
1401  pto = pfrom + 1;
1402  facedims[ind] = facedims[ind + 3];
1403  rdims[ind] = ldims[ind + 3];
1404  rdims[ind + 3] = rdims[ind] + dind + ( pto < extra ? 1 : 0 );
1405  if( is_periodic && pfrom == np - 2 ) rdims[ind + 3]++; // neighbor is on periodic bdy
1406  }
1407  else if( -1 == dijk[ind] && !pfrom && gperiodic[ind] )
1408  {
1409  // downward across periodic bdy
1410  pto = np - 1;
1411  facedims[ind + 3] = facedims[ind] = gdims[ind]; // by convention, facedims is within gdims, so lower value
1412  rdims[ind + 3] =
1413  gdims[ind + 3] + 1; // by convention, local dims one greater than gdims to indicate global lower value
1414  rdims[ind] = rdims[ind + 3] - dind - 1;
1415  across_bdy[ind] = -1;
1416  }
1417  else if( 1 == dijk[ind] && pfrom == np - 1 && is_periodic )
1418  {
1419  // right across periodic bdy
1420  pto = 0;
1421  facedims[ind + 3] = facedims[ind] = gdims[ind]; // by convention, facedims is within gdims, so lowest value
1422  rdims[ind] = gdims[ind];
1423  rdims[ind + 3] = rdims[ind] + dind + ( pto < extra ? 1 : 0 );
1424  across_bdy[ind] = 1;
1425  }
1426 
1427  assert( -1 == pto || ( rdims[0] >= gdims[0] && ( rdims[3] <= gdims[3] || ( across_bdy[0] && !pfrom ) ) ) );
1428  assert( -1 == pto || ( rdims[1] >= gdims[1] && ( rdims[4] <= gdims[4] || ( across_bdy[1] && !pfrom ) ) ) );
1429  assert( -1 == pto || ( rdims[2] >= gdims[2] && rdims[5] <= gdims[5] ) );
1430  assert( -1 == pto || ( facedims[0] >= rdims[0] && facedims[3] <= rdims[3] ) );
1431  assert( -1 == pto || ( facedims[1] >= rdims[1] && facedims[4] <= rdims[4] ) );
1432  assert( -1 == pto || ( facedims[2] >= rdims[2] && facedims[5] <= rdims[5] ) );
1433  assert( -1 == pto || ( facedims[0] >= ldims[0] && facedims[3] <= ldims[3] ) );
1434  assert( -1 == pto || ( facedims[1] >= ldims[1] && facedims[4] <= ldims[4] ) );
1435  assert( -1 == pto || ( facedims[2] >= ldims[2] && facedims[5] <= ldims[5] ) );
1436 
1437  return rval;
1438 }
1439 
1440 //! get shared vertices for alljorkori partition scheme
1441 #ifndef MOAB_HAVE_MPI
1443  ScdBox*,
1444  std::vector< int >&,
1445  std::vector< int >&,
1446  std::vector< int >& )
1447 {
1448  return MB_FAILURE;
1449 #else
1451  ScdBox* box,
1452  std::vector< int >& procs,
1453  std::vector< int >& offsets,
1454  std::vector< int >& shared_indices )
1455 {
1456  // get index of partitioned dimension
1457  const int* ldims = box->box_dims();
1458  ErrorCode rval;
1459  int ijkrem[6], ijkface[6], across_bdy[3];
1460 
1461  for( int k = -1; k <= 1; k++ )
1462  {
1463  for( int j = -1; j <= 1; j++ )
1464  {
1465  for( int i = -1; i <= 1; i++ )
1466  {
1467  if( !i && !j && !k ) continue;
1468  int pto;
1469  int dijk[] = { i, j, k };
1470  rval = get_neighbor( pcomm->proc_config().proc_size(), pcomm->proc_config().proc_rank(),
1471  box->par_data(), dijk, pto, ijkrem, ijkface, across_bdy );
1472  if( MB_SUCCESS != rval ) return rval;
1473  if( -1 != pto )
1474  {
1475  if( procs.empty() || pto != *procs.rbegin() )
1476  {
1477  procs.push_back( pto );
1478  offsets.push_back( shared_indices.size() );
1479  }
1480  rval = get_indices( ldims, ijkrem, across_bdy, ijkface, shared_indices );
1481  if( MB_SUCCESS != rval ) return rval;
1482 
1483  // check indices against known #verts on local and remote
1484  // begin of this block is shared_indices[*offsets.rbegin()], end is
1485  // shared_indices.end(), halfway is
1486  // (shared_indices.size()-*offsets.rbegin())/2
1487 #ifndef NDEBUG
1488  int start_idx = *offsets.rbegin(), end_idx = shared_indices.size(),
1489  mid_idx = ( start_idx + end_idx ) / 2;
1490 
1491  int num_local_verts = ( ldims[3] - ldims[0] + 1 ) * ( ldims[4] - ldims[1] + 1 ) *
1492  ( -1 == ldims[2] && -1 == ldims[5] ? 1 : ( ldims[5] - ldims[2] + 1 ) ),
1493  num_remote_verts = ( ijkrem[3] - ijkrem[0] + 1 ) * ( ijkrem[4] - ijkrem[1] + 1 ) *
1494  ( -1 == ijkrem[2] && -1 == ijkrem[5] ? 1 : ( ijkrem[5] - ijkrem[2] + 1 ) );
1495 
1496  assert(
1497  *std::min_element( &shared_indices[start_idx], &shared_indices[mid_idx] ) >= 0 &&
1498  *std::max_element( &shared_indices[start_idx], &shared_indices[mid_idx] ) < num_local_verts &&
1499  *std::min_element( &shared_indices[mid_idx], &shared_indices[end_idx] ) >= 0 &&
1500  *std::max_element( &shared_indices[mid_idx], &shared_indices[end_idx] ) < num_remote_verts );
1501 #endif
1502  }
1503  }
1504  }
1505  }
1506 
1507  offsets.push_back( shared_indices.size() );
1508 
1509  return MB_SUCCESS;
1510 #endif
1511 }
1512 
1513 std::ostream& operator<<( std::ostream& str, const ScdParData& pd )
1514 {
1515  str << "Partition method = " << ScdParData::PartitionMethodNames[pd.partMethod] << ", gDims = (" << pd.gDims[0]
1516  << "," << pd.gDims[1] << "," << pd.gDims[2] << ")-(" << pd.gDims[3] << "," << pd.gDims[4] << "," << pd.gDims[5]
1517  << "), gPeriodic = (" << pd.gPeriodic[0] << "," << pd.gPeriodic[1] << "," << pd.gPeriodic[2] << "), pDims = ("
1518  << pd.pDims[0] << "," << pd.pDims[1] << "," << pd.pDims[2] << ")" << std::endl;
1519  return str;
1520 }
1521 
1522 } // namespace moab