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