 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 #include "moab/ScdInterface.hpp"
2 #include "moab/Core.hpp"
3 #include "SequenceManager.hpp"
4 #include "EntitySequence.hpp"
5 #include "StructuredElementSeq.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
33 ScdInterface::ScdInterface( Interface* imp, bool boxes )
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
41 ScdInterface::~ScdInterface()
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
52 Interface* ScdInterface::impl() const
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
75 ErrorCode ScdInterface::find_boxes( Range& scd_boxes )
76 {
77 ErrorCode rval = MB_SUCCESS;
78 box_dims_tag();
79 Range boxes;
80 if( !searchedBoxes )
81 {
82 rval = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &boxDimsTag, NULL, 1, boxes, Interface::UNION );
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
103 ScdBox* ScdInterface::get_scd_box( EntityHandle eh )
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
112 ErrorCode ScdInterface::construct_box( HomCoord low,
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
252 ErrorCode ScdInterface::assign_global_ids( ScdBox* box )
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
284 ErrorCode ScdInterface::create_scd_sequence( const HomCoord& low,
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
337 ErrorCode ScdInterface::create_box_set( const HomCoord& low,
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 =
376 mbImpl->tag_get_handle( "BOX_PERIODIC", 3, MB_TYPE_INTEGER, boxPeriodicTag, MB_TAG_SPARSE | MB_TAG_CREAT );
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
392 ErrorCode rval = mbImpl->tag_get_handle( "BOX_DIMS", 6, MB_TYPE_INTEGER, boxDimsTag, MB_TAG_SPARSE | MB_TAG_CREAT );
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;
403 if( MB_TAG_NOT_FOUND == mbImpl->tag_get_name( globalBoxDimsTag, tag_name ) ) globalBoxDimsTag = NULL;
404 }
405
406 if( globalBoxDimsTag || !create_if_missing ) return globalBoxDimsTag;
407
408 ErrorCode rval =
409 mbImpl->tag_get_handle( "GLOBAL_BOX_DIMS", 6, MB_TYPE_INTEGER, globalBoxDimsTag, MB_TAG_SPARSE | MB_TAG_CREAT );
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 =
426 mbImpl->tag_get_handle( "PARTITION_METHOD", 1, MB_TYPE_INTEGER, partMethodTag, MB_TAG_SPARSE | MB_TAG_CREAT );
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,
443 MB_TAG_SPARSE | MB_TAG_CREAT );
444 if( MB_SUCCESS != rval ) return 0;
445 return boxSetTag;
446 }
447
448 //! Remove the box from the list on ScdInterface
449 ErrorCode ScdInterface::remove_box( ScdBox* box )
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
462 ErrorCode ScdInterface::add_box( ScdBox* box )
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
474 ScdBox::ScdBox( ScdInterface* impl, EntityHandle bset, EntitySequence* seq1, EntitySequence* seq2 )
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 }
491 startVertex = vertDat->start_handle();
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
526 startElem = elemSeq->start_handle();
527 }
528 else
529 {
530 Range elems;
531 impl->mbImpl->get_entities_by_dimension(
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
548 ScdBox::~ScdBox()
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
575 int ScdBox::box_dimension() const
576 {
577 return ( startElem ? scImpl->mbImpl->dimension_from_handle( startElem ) : -1 );
578 }
579
580 ErrorCode ScdBox::add_vbox( ScdBox* vbox,
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
598 bool ScdBox::boundary_complete() const
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
622 ErrorCode ScdBox::vert_dat( ScdVertexData* vert_dt )
623 {
624 vertDat = vert_dt;
625 return MB_SUCCESS;
626 }
627
628 ErrorCode ScdBox::elem_seq( EntitySequence* elem_sq )
629 {
630 elemSeq = dynamic_cast< StructuredElementSeq* >( elem_sq );
631 if( elemSeq ) elemSeq->is_periodic( locallyPeriodic );
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
640 ErrorCode ScdBox::get_params( EntityHandle ent, HomCoord& ijkd ) const
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 */
669 ErrorCode ScdBox::get_adj_edge_or_face( int dim,
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
733 ErrorCode ScdInterface::tag_shared_vertices( ParallelComm*, ScdBox* )
734 {
735 return MB_FAILURE;
736 #else
737 ErrorCode ScdInterface::tag_shared_vertices( ParallelComm* pcomm, ScdBox* box )
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
861 ErrorCode ScdInterface::get_neighbor_alljkbal( int np,
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
977 ErrorCode ScdInterface::get_neighbor_sqij( int np,
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
1122 ErrorCode ScdInterface::get_neighbor_sqjk( int np,
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
1240 ErrorCode ScdInterface::get_neighbor_sqijk( int np,
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
1341 ErrorCode ScdInterface::get_neighbor_alljorkori( int np,
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
1442 ErrorCode ScdInterface::get_shared_vertices( ParallelComm*,
1443 ScdBox*,
1444 std::vector< int >&,
1445 std::vector< int >&,
1446 std::vector< int >& )
1447 {
1448 return MB_FAILURE;
1449 #else
1450 ErrorCode ScdInterface::get_shared_vertices( ParallelComm* pcomm,
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
1.9.1.