moab
|
00001 /* 00002 * Intx2Mesh.cpp 00003 * 00004 * Created on: Oct 2, 2012 00005 */ 00006 00007 #include "Intx2Mesh.hpp" 00008 #ifdef USE_MPI 00009 #include "moab/ParallelComm.hpp" 00010 #endif /* USE_MPI */ 00011 #include "moab/AdaptiveKDTree.hpp" 00012 #include "MBParallelConventions.h" 00013 #include "MBTagConventions.hpp" 00014 // this is for DBL_MAX 00015 #include <float.h> 00016 #include <queue> 00017 #include <sstream> 00018 #include "moab/GeomUtil.hpp" 00019 00020 namespace moab { 00021 00022 Intx2Mesh::Intx2Mesh(Interface * mbimpl): mb(mbimpl) 00023 #ifdef USE_MPI 00024 , parcomm(NULL), remote_cells(NULL) 00025 #endif 00026 { 00027 dbg_1=0; 00028 box_error=0; 00029 my_rank=0; 00030 BlueFlagTag=0; 00031 RedFlagTag=0; 00032 redParentTag =0; 00033 blueParentTag = 0; 00034 countTag = 0; 00035 } 00036 00037 Intx2Mesh::~Intx2Mesh() 00038 { 00039 // TODO Auto-generated destructor stub 00040 #ifdef USE_MPI 00041 if (remote_cells) 00042 { 00043 delete remote_cells; 00044 remote_cells=NULL; 00045 } 00046 #endif 00047 } 00048 void Intx2Mesh::createTags() 00049 { 00050 if (redParentTag) 00051 mb->tag_delete(redParentTag); 00052 if(blueParentTag) 00053 mb->tag_delete(blueParentTag); 00054 if (countTag) 00055 mb->tag_delete(countTag); 00056 /*RedEdges.clear(); 00057 localEnts.clear()*/ 00058 00059 unsigned char def_data_bit = 0; // unused by default 00060 ErrorCode rval = mb->tag_get_handle("blueFlag", 1, MB_TYPE_BIT, BlueFlagTag, 00061 MB_TAG_CREAT, &def_data_bit); 00062 if (MB_SUCCESS != rval) 00063 return; 00064 // maybe the red tag is better to be deleted every time, and recreated; 00065 // or is it easy to set all values to something again? like 0? 00066 rval = mb->tag_get_handle("redFlag", 1, MB_TYPE_BIT, RedFlagTag, MB_TAG_CREAT, 00067 &def_data_bit); 00068 ERRORV(rval, "can't get red flag tag"); 00069 00070 // assume that the edges are on the red triangles 00071 Range redElements; 00072 //Range redEdges; 00073 rval = mb->get_entities_by_dimension(mbs2, 2, redElements, false); // so all tri, quad, poly 00074 ERRORV(rval, "can't get ents by dimension"); 00075 00076 // create red edges if they do not exist yet; so when they are looked upon, they are found 00077 // this is the only call that is potentially NlogN, in the whole method 00078 rval = mb->get_adjacencies(redElements, 1, true, RedEdges, Interface::UNION); 00079 ERRORV(rval, "can't get adjacent red edges"); 00080 00081 // now, create a map from each edge to a list of potential new nodes on a red edge 00082 // this memory has to be cleaned up 00083 // change it to a vector, and use the index in range of red edges 00084 int indx = 0; 00085 extraNodesVec.reserve(RedEdges.size()); 00086 for (Range::iterator eit = RedEdges.begin(); eit != RedEdges.end(); 00087 eit++, indx++) 00088 { 00089 //EntityHandle edge = *eit; 00090 //extraNodesMap[edge] = new std::vector<EntityHandle>; 00091 std::vector<EntityHandle> * nv = new std::vector<EntityHandle>; 00092 extraNodesVec.push_back(nv); 00093 } 00094 00095 int defaultInt = 0; 00096 00097 rval = mb->tag_get_handle("RedParent", 1, MB_TYPE_INTEGER, redParentTag, 00098 MB_TAG_SPARSE | MB_TAG_CREAT, &defaultInt); 00099 ERRORV(rval, "can't create positive tag"); 00100 00101 rval = mb->tag_get_handle("BlueParent", 1, MB_TYPE_INTEGER, blueParentTag, 00102 MB_TAG_SPARSE | MB_TAG_CREAT, &defaultInt); 00103 ERRORV(rval, "can't create negative tag"); 00104 00105 rval = mb->tag_get_handle("Counting", 1, MB_TYPE_INTEGER, countTag, 00106 MB_TAG_SPARSE | MB_TAG_CREAT, &defaultInt); 00107 ERRORV(rval, "can't create Counting tag"); 00108 00109 return; 00110 } 00111 00112 00113 // specify also desired set; we are interested only in neighbors in the set! 00114 // we should always get manifold mesh, each edge is adjacent to 2 cell 00115 // maybe we should check that first, just in case 00116 ErrorCode Intx2Mesh::GetOrderedNeighbors(EntityHandle set, EntityHandle cell, 00117 EntityHandle neighbors[MAXEDGES]) 00118 { 00119 int nnodes = 3; 00120 // will get the nnodes ordered neighbors; 00121 // first cell is for nodes 0, 1, second to 1, 2, third to 2, 3, last to nnodes-1, 00122 const EntityHandle * conn4; 00123 ErrorCode rval = mb->get_connectivity(cell, conn4, nnodes); 00124 int nsides = nnodes; 00125 // account for possible padded polygons 00126 while (conn4[nsides-2]==conn4[nsides-1] && nsides>3) 00127 nsides--; 00128 ERRORR(rval, "can't get connectivity on an element"); 00129 for (int i = 0; i < nsides; i++) 00130 { 00131 EntityHandle v[2]; 00132 v[0] = conn4[i]; 00133 v[1] = conn4[(i + 1) % nsides]; 00134 // get quads adjacent to vertices 00135 std::vector<EntityHandle> cells; 00136 std::vector<EntityHandle> cellsInSet; 00137 rval = mb->get_adjacencies(v, 2, 2, false, cells, Interface::INTERSECT); 00138 ERRORR(rval, "can't get adjacencies on 2 nodes"); 00139 size_t siz = cells.size(); 00140 for (size_t j = 0; j < siz; j++) 00141 if (mb->contains_entities(set, &(cells[j]), 1)) 00142 cellsInSet.push_back(cells[j]); 00143 siz = cellsInSet.size(); 00144 00145 if (siz > 2) 00146 { 00147 std::cout << "non manifold mesh, error" 00148 << mb->list_entities(&(cellsInSet[0]), cellsInSet.size()) << "\n"; 00149 return MB_FAILURE; // non-manifold 00150 } 00151 if (siz == 1) 00152 { 00153 // it must be the border, 00154 neighbors[i] = 0; // we are guaranteed that ids are !=0; this is marking a border 00155 // borders do not appear for a sphere in serial, but they do appear for 00156 // parallel processing anyway 00157 continue; 00158 } 00159 // here siz ==2, it is either the first or second 00160 if (cell == cellsInSet[0]) 00161 neighbors[i] = cellsInSet[1]; 00162 else 00163 neighbors[i] = cellsInSet[0]; 00164 } 00165 return MB_SUCCESS; 00166 } 00167 // main interface; this will do the advancing front trick 00168 // some are triangles, some are quads, some are polygons ... 00169 ErrorCode Intx2Mesh::intersect_meshes(EntityHandle mbset1, EntityHandle mbset2, 00170 EntityHandle & outputSet) 00171 { 00172 00173 ErrorCode rval; 00174 mbs1 = mbset1; // set 1 is departure, and it is completely covering the euler set on proc 00175 mbs2 = mbset2; 00176 outSet = outputSet; 00177 00178 // really, should be something from t1 and t2; blue is 1 (lagrange), red is 2 (euler) 00179 createTags(); // 00180 EntityHandle startBlue, startRed; 00181 00182 mb->get_entities_by_dimension(mbs1, 2, rs1); 00183 mb->get_entities_by_dimension(mbs2, 2, rs2); 00184 Range rs22=rs2; // a copy of the initial range; we will remove from it elements as we 00185 // advance ; rs2 is needed for marking the polygon to the red parent 00186 while (!rs22.empty()) 00187 { 00188 for (Range::iterator it = rs1.begin(); it != rs1.end(); it++) 00189 { 00190 startBlue = *it; 00191 int found = 0; 00192 for (Range::iterator it2 = rs22.begin(); it2 != rs22.end() && !found; it2++) 00193 { 00194 startRed = *it2; 00195 double area = 0; 00196 // if area is > 0 , we have intersections 00197 double P[10*MAXEDGES]; // max 8 intx points + 8 more in the polygon 00198 // 00199 int nP = 0; 00200 int nb[MAXEDGES], nr[MAXEDGES]; // sides 3 or 4? also, check boxes first 00201 int nsRed, nsBlue; 00202 computeIntersectionBetweenRedAndBlue(startRed, startBlue, P, nP, area, nb, nr, 00203 nsBlue, nsRed, true); 00204 if (area > 0) 00205 { 00206 found = 1; 00207 break; // found 2 elements that intersect; these will be the seeds 00208 } 00209 } 00210 if (found) 00211 break; 00212 } 00213 00214 std::queue<EntityHandle> blueQueue; // these are corresponding to Ta, 00215 blueQueue.push(startBlue); 00216 std::queue<EntityHandle> redQueue; 00217 redQueue.push(startRed); 00218 00219 Range toResetBlues; // will be used to reset blue flags for every red quad 00220 // processed 00221 00222 /*if (my_rank==0) 00223 dbg_1 = 1;*/ 00224 unsigned char used = 1; 00225 unsigned char unused = 0; // for red flags 00226 // mark the start blue quad as used, so it will not come back again 00227 mb->tag_set_data(RedFlagTag, &startRed, 1, &used); 00228 while (!redQueue.empty()) 00229 { 00230 // flags for the side : 0 means a blue quad not found on side 00231 // a paired blue not found yet for the neighbors of red 00232 EntityHandle n[MAXEDGES] = { EntityHandle(0) }; 00233 00234 int nsidesRed; // will be initialized later, when we compute intersection 00235 EntityHandle currentRed = redQueue.front(); 00236 00237 redQueue.pop(); 00238 // for (k=0; k<m_numPos; k++) 00239 // redFlag[k] = 0; 00240 // redFlag[m_numPos] = 1; // to guard for the boundary 00241 // all reds that were tagged, are now cleared 00242 for (Range::iterator itr = toResetBlues.begin(); itr != toResetBlues.end(); 00243 itr++) 00244 { 00245 EntityHandle ttt = *itr; 00246 rval = mb->tag_set_data(BlueFlagTag, &ttt, 1, &unused); 00247 ERRORR(rval, "can't set blue unused tag"); 00248 } 00249 //rval = mb2->tag_set_data(RedFlagTag, toResetReds, &unused); 00250 /*if (dbg_1) 00251 { 00252 std::cout << "reset blues: "; 00253 mb->list_entities(toResetBlues); 00254 }*/ 00255 //rval = mb2->tag_set_data(RedFlagTag, toResetReds, &unused); 00256 if (dbg_1) 00257 { 00258 std::cout << "reset blues: "; 00259 for (Range::iterator itr = toResetBlues.begin(); itr != toResetBlues.end(); 00260 itr++) 00261 std::cout << mb->id_from_handle(*itr) << " "; 00262 std::cout << std::endl; 00263 } 00264 EntityHandle currentBlue = blueQueue.front(); // where do we check for redQueue???? 00265 // red and blue queues are parallel 00266 blueQueue.pop(); // mark the current red 00267 //redFlag[currentRed] = 1; // 00268 toResetBlues.clear(); // empty the range of used blues, will have to be set unused again, 00269 // at the end of red element processing 00270 toResetBlues.insert(currentBlue); 00271 rval = mb->tag_set_data(BlueFlagTag, ¤tBlue, 1, &used); 00272 ERRORR(rval, "can't set blue tag"); 00273 //mb2->set_tag_data 00274 std::queue<EntityHandle> localBlue; 00275 localBlue.push(currentBlue); 00276 while (!localBlue.empty()) 00277 { 00278 // 00279 EntityHandle blueT = localBlue.front(); 00280 localBlue.pop(); 00281 double P[10*MAXEDGES], area; // 00282 int nP = 0; 00283 int nb[MAXEDGES] = {0}; 00284 int nr[MAXEDGES] = {0}; 00285 00286 int nsidesBlue; 00287 // area is in 2d, points are in 3d (on a sphere), back-projected, or in a plane 00288 // intersection points could include the vertices of initial elements 00289 // nb [j] = 0 means no intersection on the side k for element blue (markers) 00290 // nb [j] = 1 means that the side j (from j to j+1) of blue poly intersects the 00291 // red poly. A potential next poly is the red poly that is adjacent to this side 00292 computeIntersectionBetweenRedAndBlue(/* red */currentRed, blueT, P, nP, 00293 area, nb, nr, nsidesBlue, nsidesRed); 00294 if (nP > 0) 00295 { 00296 if (dbg_1) 00297 { 00298 for (int k=0; k<3; k++) 00299 { 00300 std::cout << " nb, nr: " << k << " " << nb[k] << " " << nr[k] << "\n"; 00301 } 00302 } 00303 // intersection found: output P and original triangles if nP > 2 00304 00305 EntityHandle neighbors[MAXEDGES]; 00306 rval = GetOrderedNeighbors(mbs1, blueT, neighbors); 00307 if (rval != MB_SUCCESS) 00308 { 00309 std::cout << " can't get the neighbors for blue element " 00310 << mb->id_from_handle(blueT); 00311 return MB_FAILURE; 00312 } 00313 00314 // add neighbors to the localBlue queue, if they are not marked 00315 for (int nn = 0; nn < nsidesBlue; nn++) 00316 { 00317 EntityHandle neighbor = neighbors[nn]; 00318 if (neighbor > 0 && nb[nn]>0) // advance across blue boundary n 00319 { 00320 //n[nn] = redT; // start from 0!! 00321 unsigned char status = 0; 00322 mb->tag_get_data(BlueFlagTag, &neighbor, 1, &status); 00323 if (status == 0) 00324 { 00325 localBlue.push(neighbor); 00326 /*if (dbg_1) 00327 { 00328 std::cout << " local blue elem " 00329 << mb->list_entities(&neighbor, 1) << "\n for red:" 00330 << mb->list_entities(¤tRed, 1) << "\n"; 00331 }*/ 00332 if (dbg_1) 00333 { 00334 std::cout << " local blue elem " << mb->id_from_handle(neighbor) 00335 << " for red:" << mb->id_from_handle(currentRed) << "\n"; 00336 mb->list_entities(&neighbor, 1); 00337 } 00338 rval = mb->tag_set_data(BlueFlagTag, &neighbor, 1, &used); 00339 //redFlag[neighbor] = 1; // flag it to not be added anymore 00340 toResetBlues.insert(neighbor); // this is used to reset the red flag 00341 } 00342 } 00343 // n(find(nc>0))=ac; % ac is starting candidate for neighbor 00344 if (nr[nn] > 0) 00345 n[nn] = blueT; 00346 00347 } 00348 if (nP > 1) // this will also construct triangles/polygons in the new mesh, if needed 00349 findNodes(currentRed, nsidesRed, blueT, nsidesBlue, P, nP); 00350 } 00351 else if (dbg_1) 00352 /*{ 00353 std::cout << " red, blue, do not intersect: " 00354 << mb->list_entities(¤tRed, 1) << " " 00355 << mb->list_entities(&blueT, 1) << "\n"; 00356 }*/ 00357 { 00358 std::cout << " red, blue, do not intersect: " 00359 << mb->id_from_handle(currentRed) << " " 00360 << mb->id_from_handle(blueT) << "\n"; 00361 } 00362 00363 } // end while (!localBlue.empty()) 00364 // here, we are finished with redCurrent, take it out of the rs22 range (red, arrival mesh) 00365 rs22.erase(currentRed); 00366 // also, look at its neighbors, and add to the seeds a next one 00367 00368 EntityHandle redNeighbors[MAXEDGES]; 00369 rval = GetOrderedNeighbors(mbs2, currentRed, redNeighbors); 00370 ERRORR(rval, "can't get neighbors"); 00371 if (dbg_1) 00372 { 00373 std::cout << "Next: neighbors for current red "; 00374 for (int kk = 0; kk < nsidesRed; kk++) 00375 { 00376 if (redNeighbors[kk] > 0) 00377 std::cout << mb->id_from_handle(redNeighbors[kk]) << " "; 00378 else 00379 std::cout << 0 << " "; 00380 } 00381 std::cout << std::endl; 00382 } 00383 for (int j = 0; j < nsidesRed; j++) 00384 { 00385 EntityHandle redNeigh = redNeighbors[j]; 00386 unsigned char status = 1; 00387 if (redNeigh == 0) 00388 continue; 00389 mb->tag_get_data(RedFlagTag, &redNeigh, 1, &status); // status 0 is unused 00390 if (status == 0 && n[j] > 0) // not treated yet and marked as a neighbor 00391 { 00392 // we identified red quad n[j] as intersecting with neighbor j of the blue quad 00393 redQueue.push(redNeigh); 00394 blueQueue.push(n[j]); 00395 if (dbg_1) 00396 std::cout << "new polys pushed: blue, red:" 00397 << mb->id_from_handle(redNeigh) << " " 00398 << mb->id_from_handle(n[j]) << std::endl; 00399 mb->tag_set_data(RedFlagTag, &redNeigh, 1, &used); 00400 } 00401 } 00402 00403 } // end while (!redQueue.empty()) 00404 } 00405 if (dbg_1) 00406 { 00407 for (int k = 0; k < 6; k++) 00408 mout_1[k].close(); 00409 } 00410 // before cleaning up , we need to settle the position of the intersection points 00411 // on the boundary edges 00412 // this needs to be collective, so we should maybe wait something 00413 #ifdef USE_MPI 00414 rval = correct_intersection_points_positions(); 00415 if (rval!=MB_SUCCESS) 00416 { 00417 std::cout << "can't correct position, Intx2Mesh.cpp \n"; 00418 } 00419 #endif 00420 clean(); 00421 return MB_SUCCESS; 00422 } 00423 00424 // clean some memory allocated 00425 void Intx2Mesh::clean() 00426 { 00427 // 00428 int indx = 0; 00429 for (Range::iterator eit = RedEdges.begin(); eit != RedEdges.end(); 00430 eit++, indx++) 00431 { 00432 //EntityHandle edge = *eit; 00433 //delete extraNodesMap[edge]; 00434 delete extraNodesVec[indx]; 00435 } 00436 //extraNodesMap.clear(); 00437 extraNodesVec.clear(); 00438 // also, delete some bit tags, used to mark processed reds and blues 00439 mb->tag_delete(RedFlagTag);// to mark blue quads already considered 00440 mb->tag_delete(BlueFlagTag); 00441 00442 /*mb->tag_delete(redParentTag); 00443 mb->tag_delete(blueParentTag); 00444 mb->tag_delete(countTag); 00445 RedEdges.clear(); 00446 localEnts.clear();*/ 00447 00448 } 00449 // this method will reduce number of nodes, collapse edges that are of length 0 00450 // so a polygon like 428 431 431 will become a line 428 431 00451 // or something like 428 431 431 531 -> 428 431 531 00452 void Intx2Mesh::correct_polygon(EntityHandle * nodes, int & nP) 00453 { 00454 int i = 0; 00455 while(i<nP) 00456 { 00457 int nextIndex = (i+1)%nP; 00458 if (nodes[i]==nodes[nextIndex]) 00459 { 00460 // we need to reduce nP, and collapse nodes 00461 if (dbg_1) 00462 { 00463 std::cout<<" nodes duplicated in list: " ; 00464 for (int j=0; j<nP; j++) 00465 std::cout<<nodes[j] << " " ; 00466 std::cout<<"\n"; 00467 std::cout<<" node " << nodes[i] << " at index " << i << " is duplicated" << "\n"; 00468 } 00469 // this will work even if we start from 1 2 3 1; when i is 3, we find nextIndex is 0, then next thing does nothing 00470 // (nP-1 is 3, so k is already >= nP-1); it will result in nodes -> 1, 2, 3 00471 for (int k=i; k<nP-1; k++) 00472 nodes[k] = nodes[k+1]; 00473 nP--; // decrease the number of nodes; also, decrease i, just if we may need to check again 00474 i--; 00475 } 00476 i++; 00477 } 00478 return; 00479 } 00480 #if USE_MPI 00481 ErrorCode Intx2Mesh::build_processor_euler_boxes(EntityHandle euler_set, Range & local_verts) 00482 { 00483 localEnts.clear(); 00484 ErrorCode rval = mb->get_entities_by_dimension(euler_set, 2, localEnts); 00485 ERRORR(rval, "can't get ents by dimension"); 00486 00487 rval = mb->get_connectivity(localEnts, local_verts); 00488 int num_local_verts = (int) local_verts.size(); 00489 ERRORR(rval, "can't get local vertices"); 00490 00491 parcomm = ParallelComm::get_pcomm(mb, 0); 00492 if (NULL==parcomm) 00493 return MB_FAILURE; 00494 00495 // get the position of local vertices, and decide local boxes (allBoxes...) 00496 double bmin[3]={DBL_MAX, DBL_MAX, DBL_MAX}; 00497 double bmax[3] ={-DBL_MAX, -DBL_MAX, -DBL_MAX}; 00498 00499 std::vector<double> coords(3*num_local_verts); 00500 rval = mb->get_coords(local_verts, &coords[0]); 00501 ERRORR(rval, "can't get coords of vertices "); 00502 00503 for (int i=0; i< num_local_verts; i++) 00504 { 00505 for (int k=0; k<3; k++) 00506 { 00507 double val=coords[3*i+k]; 00508 if (val < bmin[k]) 00509 bmin[k]=val; 00510 if (val > bmax[k]) 00511 bmax[k] = val; 00512 } 00513 } 00514 int numprocs=parcomm->proc_config().proc_size(); 00515 allBoxes.resize(6*numprocs); 00516 00517 my_rank = parcomm->proc_config().proc_rank() ; 00518 for (int k=0; k<3; k++) 00519 { 00520 allBoxes[6*my_rank+k]=bmin[k]; 00521 allBoxes[6*my_rank+3+k] = bmax[k]; 00522 } 00523 00524 // now communicate to get all boxes 00525 int mpi_err; 00526 #if (MPI_VERSION >= 2) 00527 // use "in place" option 00528 mpi_err = MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, 00529 &allBoxes[0], 6, MPI_DOUBLE, 00530 parcomm->proc_config().proc_comm()); 00531 #else 00532 { 00533 std::vector<double> allBoxes_tmp(6*parcomm->proc_config().proc_size()); 00534 mpi_err = MPI_Allgather( &allBoxes[6*my_rank], 6, MPI_DOUBLE, 00535 &allBoxes_tmp[0], 6, MPI_DOUBLE, 00536 parcomm->proc_config().proc_comm()); 00537 allBoxes = allBoxes_tmp; 00538 } 00539 #endif 00540 if (MPI_SUCCESS != mpi_err) return MB_FAILURE; 00541 00542 // also process the max number of vertices per cell (4 for quads, but could be more for polygons) 00543 int local_max_edges = 3; 00544 for (Range::iterator it = localEnts.begin(); it!=localEnts.end(); it++) 00545 { 00546 const EntityHandle * conn; 00547 int num_nodes; 00548 rval = mb->get_connectivity(*it, conn, num_nodes); 00549 ERRORR(rval, "can't get connectivity"); 00550 if (num_nodes>local_max_edges) 00551 local_max_edges = num_nodes; 00552 } 00553 00554 // now reduce max_edges over all processors 00555 mpi_err = MPI_Allreduce(&local_max_edges, &max_edges, 1, MPI_INTEGER, MPI_MAX, parcomm->proc_config().proc_comm()); 00556 if (MPI_SUCCESS != mpi_err) return MB_FAILURE; 00557 00558 if (my_rank==0) 00559 { 00560 std::cout << " maximum number of vertices per cell is " << max_edges << "\n"; 00561 for (int i=0; i<numprocs; i++) 00562 { 00563 std::cout<<"proc: " << i << " box min: " << allBoxes[6*i ] << " " <<allBoxes[6*i+1] << " " << allBoxes[6*i+2] << " \n"; 00564 std::cout<< " box max: " << allBoxes[6*i+3] << " " <<allBoxes[6*i+4] << " " << allBoxes[6*i+5] << " \n"; 00565 } 00566 } 00567 00568 return MB_SUCCESS; 00569 } 00570 ErrorCode Intx2Mesh::create_departure_mesh_2nd_alg(EntityHandle & euler_set, EntityHandle & covering_lagr_set) 00571 { 00572 // compute the bounding box on each proc 00573 parcomm = ParallelComm::get_pcomm(mb, 0); 00574 if (NULL==parcomm) 00575 return MB_FAILURE; 00576 00577 localEnts.clear(); 00578 ErrorCode rval = mb->get_entities_by_dimension(euler_set, 2, localEnts); 00579 ERRORR(rval, "can't get ents by dimension"); 00580 00581 Tag dpTag = 0; 00582 std::string tag_name("DP"); 00583 rval = mb->tag_get_handle(tag_name.c_str(), 3, MB_TYPE_DOUBLE, dpTag, MB_TAG_DENSE); 00584 ERRORR(rval, "can't get DP tag"); 00585 00586 EntityHandle dum=0; 00587 Tag corrTag; 00588 rval = mb->tag_get_handle(CORRTAGNAME, 00589 1, MB_TYPE_HANDLE, corrTag, 00590 MB_TAG_DENSE|MB_TAG_CREAT, &dum); 00591 ERRORR(rval, "can't get CORR tag"); 00592 // get all local verts 00593 Range local_verts; 00594 rval = mb->get_connectivity(localEnts, local_verts); 00595 int num_local_verts = (int) local_verts.size(); 00596 ERRORR(rval, "can't get local vertices"); 00597 00598 rval = build_processor_euler_boxes(euler_set, local_verts); 00599 ERRORR(rval, "can't build processor boxes"); 00600 Tag gid; 00601 rval = mb->tag_get_handle(GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, gid, MB_TAG_DENSE); 00602 ERRORR(rval,"can't get global ID tag" ); 00603 std::vector<int> gids(num_local_verts); 00604 rval = mb->tag_get_data(gid, local_verts, &gids[0]); 00605 ERRORR(rval, "can't get local vertices gids"); 00606 00607 // now see the departure points; to what boxes should we send them? 00608 std::vector<double> dep_points(3*num_local_verts); 00609 rval = mb->tag_get_data(dpTag, local_verts, (void*)&dep_points[0]); 00610 ERRORR(rval, "can't get DP tag values"); 00611 // ranges to send to each processor; will hold vertices and elements (quads?) 00612 // will look if the box of the dep quad covers box of euler mesh on proc (with tolerances) 00613 std::map<int, Range> Rto; 00614 int numprocs=parcomm->proc_config().proc_size(); 00615 00616 for (Range::iterator eit = localEnts.begin(); eit!=localEnts.end(); eit++) 00617 { 00618 EntityHandle q=*eit; 00619 const EntityHandle * conn4; 00620 int num_nodes; 00621 rval= mb->get_connectivity(q, conn4, num_nodes); 00622 ERRORR(rval, "can't get DP tag values"); 00623 CartVect qbmin(DBL_MAX); 00624 CartVect qbmax(-DBL_MAX); 00625 for (int i=0; i<num_nodes; i++) 00626 { 00627 EntityHandle v=conn4[i]; 00628 size_t index=local_verts.find(v)-local_verts.begin(); 00629 CartVect dp( &dep_points[3*index] ); // will use constructor 00630 for (int j=0; j<3; j++) 00631 { 00632 if (qbmin[j]>dp[j]) 00633 qbmin[j]=dp[j]; 00634 if (qbmax[j]<dp[j]) 00635 qbmax[j]=dp[j]; 00636 } 00637 } 00638 for (int p=0; p<numprocs; p++) 00639 { 00640 CartVect bbmin(&allBoxes[6*p]); 00641 CartVect bbmax(&allBoxes[6*p+3]); 00642 if ( GeomUtil::boxes_overlap( bbmin, bbmax, qbmin, qbmax, box_error) ) 00643 { 00644 Rto[p].insert(q); 00645 } 00646 } 00647 } 00648 00649 // now, build TLv and TLq, for each p 00650 size_t numq=0; 00651 size_t numv=0; 00652 for (int p=0; p<numprocs; p++) 00653 { 00654 if (p==(int)my_rank) 00655 continue; // do not "send" it, because it is already here 00656 Range & range_to_P = Rto[p]; 00657 // add the vertices to it 00658 if (range_to_P.empty()) 00659 continue;// nothing to send to proc p 00660 Range vertsToP; 00661 rval = mb->get_connectivity(range_to_P, vertsToP); 00662 ERRORR(rval, "can't get connectivity"); 00663 numq=numq+range_to_P.size(); 00664 numv=numv+vertsToP.size(); 00665 range_to_P.merge(vertsToP); 00666 } 00667 TupleList TLv; 00668 TupleList TLq; 00669 TLv.initialize(2, 0, 0, 3, numv); // to proc, GLOBAL ID, DP points 00670 TLv.enableWriteAccess(); 00671 00672 int sizeTuple = 2+max_edges; // determined earlier 00673 TLq.initialize(2+max_edges, 0, 1, 0, numq); // to proc, elem GLOBAL ID, connectivity[10] (global ID v), local eh 00674 TLq.enableWriteAccess(); 00675 std::cout << "from proc " << my_rank << " send " << numv << " vertices and " << numq << " elements\n"; 00676 00677 for (int to_proc=0; to_proc<numprocs; to_proc++) 00678 { 00679 if (to_proc==(int)my_rank) 00680 continue; 00681 Range & range_to_P = Rto[to_proc]; 00682 Range V = range_to_P.subset_by_type(MBVERTEX); 00683 00684 for (Range::iterator it=V.begin(); it!=V.end(); it++) 00685 { 00686 EntityHandle v=*it; 00687 unsigned int index = local_verts.find(v)-local_verts.begin(); 00688 int n=TLv.get_n(); 00689 TLv.vi_wr[2*n] = to_proc; // send to processor 00690 TLv.vi_wr[2*n+1] = gids[index]; // global id needs index in the local_verts range 00691 TLv.vr_wr[3*n] = dep_points[3*index]; // departure position, of the node local_verts[i] 00692 TLv.vr_wr[3*n+1] = dep_points[3*index+1]; 00693 TLv.vr_wr[3*n+2] = dep_points[3*index+2]; 00694 TLv.inc_n(); 00695 } 00696 // also, prep the quad for sending ... 00697 Range Q = range_to_P.subset_by_dimension(2); 00698 for (Range::iterator it=Q.begin(); it!=Q.end(); it++) 00699 { 00700 EntityHandle q=*it; 00701 int global_id; 00702 rval = mb->tag_get_data(gid, &q, 1, &global_id); 00703 ERRORR(rval, "can't get gid for polygon"); 00704 int n=TLq.get_n(); 00705 TLq.vi_wr[sizeTuple*n] = to_proc; // 00706 TLq.vi_wr[sizeTuple*n+1] = global_id; // global id of element, used to identify it ... 00707 const EntityHandle * conn4; 00708 int num_nodes; 00709 rval = mb->get_connectivity(q, conn4, num_nodes);// could be up to MAXEDGES, but it is limited by max_edges 00710 ERRORR(rval, "can't get connectivity for cell"); 00711 if (num_nodes > MAXEDGES) 00712 ERRORR(MB_FAILURE, "too many nodes in a polygon"); 00713 for (int i=0; i<num_nodes; i++) 00714 { 00715 EntityHandle v = conn4[i]; 00716 unsigned int index = local_verts.find(v)-local_verts.begin(); 00717 TLq.vi_wr[sizeTuple*n+2+i] = gids[index]; 00718 } 00719 for (int k=num_nodes; k<max_edges; k++) 00720 { 00721 TLq.vi_wr[sizeTuple*n+2+k] = 0; // fill the rest of node ids with 0; we know that the node ids start from 1! 00722 } 00723 TLq.vul_wr[n]=q; // save here the entity handle, it will be communicated back 00724 // mabe we should forget about global ID 00725 TLq.inc_n(); 00726 00727 } 00728 00729 } 00730 00731 // now we are done populating the tuples; route them to the appropriate processors 00732 (parcomm->proc_config().crystal_router())->gs_transfer(1, TLv, 0); 00733 (parcomm->proc_config().crystal_router())->gs_transfer(1, TLq, 0); 00734 // the elements are already in localEnts; 00735 00736 // maps from global ids to new vertex and quad handles, that are added 00737 std::map<int, EntityHandle> globalID_to_handle; 00738 /*std::map<int, EntityHandle> globalID_to_eh;*/ 00739 globalID_to_eh.clear();// need for next iteration 00740 // now, look at every TLv, and see if we have to create a vertex there or not 00741 int n=TLv.get_n();// the size of the points received 00742 for (int i=0; i<n; i++) 00743 { 00744 int globalId = TLv.vi_rd[2*i+1]; 00745 if (globalID_to_handle.find(globalId)==globalID_to_handle.end()) 00746 { 00747 EntityHandle new_vert; 00748 double dp_pos[3]= {TLv.vr_wr[3*i], TLv.vr_wr[3*i+1], TLv.vr_wr[3*i+2]}; 00749 rval = mb->create_vertex(dp_pos, new_vert); 00750 ERRORR(rval, "can't create new vertex "); 00751 globalID_to_handle[globalId]= new_vert; 00752 } 00753 } 00754 00755 // now, all dep points should be at their place 00756 // look in the local list of q for this proc, and create all those quads and vertices if needed 00757 // it may be an overkill, but because it does not involve communication, we do it anyway 00758 Range & local=Rto[my_rank]; 00759 Range local_q = local.subset_by_dimension(2); 00760 // the local should have all the vertices in local_verts 00761 for (Range::iterator it=local_q.begin(); it!=local_q.end(); it++) 00762 { 00763 EntityHandle q=*it; 00764 int nnodes; 00765 const EntityHandle * conn4; 00766 rval = mb->get_connectivity(q, conn4, nnodes); 00767 ERRORR(rval, "can't get connectivity of local q "); 00768 EntityHandle new_conn[MAXEDGES]; 00769 for (int i=0; i<nnodes; i++) 00770 { 00771 EntityHandle v1=conn4[i]; 00772 unsigned int index = local_verts.find(v1)-local_verts.begin(); 00773 int globalId=gids[index]; 00774 if(globalID_to_handle.find(globalId)==globalID_to_handle.end()) 00775 { 00776 // we need to create that vertex, at this position dep_points 00777 double dp_pos[3]={dep_points[3*index], dep_points[3*index+1], dep_points[3*index+2]}; 00778 EntityHandle new_vert; 00779 rval = mb->create_vertex(dp_pos, new_vert); 00780 ERRORR(rval, "can't create new vertex "); 00781 globalID_to_handle[globalId]= new_vert; 00782 } 00783 new_conn[i] = globalID_to_handle[gids[index]]; 00784 } 00785 EntityHandle new_element; 00786 // 00787 EntityType entType = MBQUAD; 00788 if (nnodes >4) 00789 entType = MBPOLYGON; 00790 if (nnodes <4) 00791 entType = MBTRI; 00792 00793 rval = mb->create_element(entType, new_conn, nnodes, new_element); 00794 ERRORR(rval, "can't create new quad "); 00795 rval = mb->add_entities(covering_lagr_set, &new_element, 1); 00796 ERRORR(rval, "can't add new element to dep set"); 00797 int gid_el; 00798 // get the global ID of the initial quad 00799 rval=mb->tag_get_data(gid, &q, 1, &gid_el); 00800 ERRORR(rval, "can't get element global ID "); 00801 globalID_to_eh[gid_el]=new_element; 00802 // is this redundant or not? 00803 rval = mb->tag_set_data(corrTag, &new_element, 1, &q); 00804 ERRORR(rval, "can't set corr tag on new el"); 00805 // set the global id on new elem 00806 rval = mb->tag_set_data(gid, &new_element, 1, &gid_el); 00807 ERRORR(rval, "can't set global id tag on new el"); 00808 } 00809 // now look at all elements received through; we do not want to duplicate them 00810 n=TLq.get_n();// number of elements received by this processor 00811 // form the remote cells, that will be used to send the tracer info back to the originating proc 00812 remote_cells = new TupleList(); 00813 remote_cells->initialize(2, 0, 1, 1, n); 00814 remote_cells->enableWriteAccess(); 00815 for (int i=0; i<n; i++) 00816 { 00817 int globalIdEl = TLq.vi_rd[sizeTuple*i+1]; 00818 int from_proc = TLq.vi_wr[sizeTuple*i]; 00819 // do we already have a quad with this global ID, represented? 00820 if (globalID_to_eh.find(globalIdEl)==globalID_to_eh.end()) 00821 { 00822 // construct the conn quad 00823 EntityHandle new_conn[MAXEDGES]; 00824 int nnodes = -1; 00825 for (int j=0; j<max_edges; j++) 00826 { 00827 int vgid = TLq.vi_rd[sizeTuple*i+2+j];// vertex global ID 00828 if (vgid==0) 00829 new_conn[j] = 0; 00830 else 00831 { 00832 assert(globalID_to_handle.find(vgid)!=globalID_to_handle.end()); 00833 new_conn[j]=globalID_to_handle[vgid]; 00834 nnodes = j+1;// nodes are at the beginning, and are variable number 00835 } 00836 } 00837 EntityHandle new_element; 00838 // 00839 EntityType entType = MBQUAD; 00840 if (nnodes >4) 00841 entType = MBPOLYGON; 00842 if (nnodes <4) 00843 entType = MBTRI; 00844 rval = mb->create_element(entType, new_conn, nnodes, new_element); 00845 ERRORR(rval, "can't create new element "); 00846 globalID_to_eh[globalIdEl]=new_element; 00847 rval = mb->add_entities(covering_lagr_set, &new_element, 1); 00848 ERRORR(rval, "can't add new element to dep set"); 00849 /* rval = mb->tag_set_data(corrTag, &new_element, 1, &q); 00850 ERRORR(rval, "can't set corr tag on new el");*/ 00851 remote_cells->vi_wr[2*i]=from_proc; 00852 remote_cells->vi_wr[2*i+1]=globalIdEl; 00853 remote_cells->vr_wr[i] = 0.; // no contribution yet sent back 00854 remote_cells->vul_wr[i]= TLq.vul_rd[i];// this is the corresponding red cell (arrival) 00855 remote_cells->inc_n(); 00856 // set the global id on new elem 00857 rval = mb->tag_set_data(gid, &new_element, 1, &globalIdEl); 00858 ERRORR(rval, "can't set global id tag on new el"); 00859 } 00860 } 00861 // order the remote cells tuple list, with the global id, because we will search in it 00862 //remote_cells->print("remote_cells before sorting"); 00863 moab::TupleList::buffer sort_buffer; 00864 sort_buffer.buffer_init(n); 00865 remote_cells->sort(1, &sort_buffer); 00866 sort_buffer.reset(); 00867 return MB_SUCCESS; 00868 } 00869 00870 // this algorithm assumes lagr set is already created, and some elements will be coming from 00871 // other procs, and populate the covering_set 00872 // we need to keep in a tuple list the remote cells from other procs, because we need to send back 00873 // the intersection info (like area of the intx polygon, and the current concentration) maybe total 00874 // mass in that intx 00875 ErrorCode Intx2Mesh::create_departure_mesh_3rd_alg(EntityHandle & lagr_set, 00876 EntityHandle & covering_set) 00877 { 00878 EntityHandle dum = 0; 00879 00880 Tag corrTag; 00881 ErrorCode rval = mb->tag_get_handle(CORRTAGNAME, 00882 1, MB_TYPE_HANDLE, corrTag, 00883 MB_TAG_DENSE, &dum); 00884 //start copy from 2nd alg 00885 // compute the bounding box on each proc 00886 parcomm = ParallelComm::get_pcomm(mb, 0); 00887 if (NULL == parcomm || ( 1==parcomm->proc_config().proc_size())) 00888 { 00889 covering_set = lagr_set; // nothing to communicate, it must be serial 00890 return MB_SUCCESS; 00891 } 00892 00893 // get all local verts 00894 Range local_verts; 00895 rval = mb->get_connectivity(localEnts, local_verts); 00896 int num_local_verts = (int) local_verts.size(); 00897 ERRORR(rval, "can't get local vertices"); 00898 00899 Tag gid; 00900 rval = mb->tag_get_handle(GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, gid, 00901 MB_TAG_DENSE); 00902 ERRORR(rval, "can't get global ID tag"); 00903 std::vector<int> gids(num_local_verts); 00904 rval = mb->tag_get_data(gid, local_verts, &gids[0]); 00905 ERRORR(rval, "can't get local vertices gids"); 00906 00907 Range localDepCells; 00908 rval = mb->get_entities_by_dimension(lagr_set, 2, localDepCells); 00909 ERRORR(rval, "can't get ents by dimension from lagr set"); 00910 00911 // get all lagr verts (departure vertices) 00912 Range lagr_verts; 00913 rval = mb->get_connectivity(localDepCells, lagr_verts);// they should be created in 00914 // the same order as the euler vertices 00915 int num_lagr_verts = (int) lagr_verts.size(); 00916 ERRORR(rval, "can't get local lagr vertices"); 00917 00918 // now see the departure points position; to what boxes should we send them? 00919 std::vector<double> dep_points(3 * num_lagr_verts); 00920 rval = mb->get_coords(lagr_verts, &dep_points[0]); 00921 ERRORR(rval, "can't get departure points position"); 00922 // ranges to send to each processor; will hold vertices and elements (quads?) 00923 // will look if the box of the dep quad covers box of euler mesh on proc (with tolerances) 00924 std::map<int, Range> Rto; 00925 int numprocs = parcomm->proc_config().proc_size(); 00926 00927 for (Range::iterator eit = localDepCells.begin(); eit != localDepCells.end(); eit++) 00928 { 00929 EntityHandle q = *eit; 00930 const EntityHandle * conn4; 00931 int num_nodes; 00932 rval = mb->get_connectivity(q, conn4, num_nodes); 00933 ERRORR(rval, "can't get DP tag values"); 00934 CartVect qbmin(DBL_MAX); 00935 CartVect qbmax(-DBL_MAX); 00936 for (int i = 0; i < num_nodes; i++) 00937 { 00938 EntityHandle v = conn4[i]; 00939 int index = lagr_verts.index(v); 00940 assert(-1!=index); 00941 CartVect dp(&dep_points[3 * index]); // will use constructor 00942 for (int j = 0; j < 3; j++) 00943 { 00944 if (qbmin[j] > dp[j]) 00945 qbmin[j] = dp[j]; 00946 if (qbmax[j] < dp[j]) 00947 qbmax[j] = dp[j]; 00948 } 00949 } 00950 for (int p = 0; p < numprocs; p++) 00951 { 00952 CartVect bbmin(&allBoxes[6 * p]); 00953 CartVect bbmax(&allBoxes[6 * p + 3]); 00954 if (GeomUtil::boxes_overlap(bbmin, bbmax, qbmin, qbmax, box_error)) 00955 { 00956 Rto[p].insert(q); 00957 } 00958 } 00959 } 00960 00961 // now, build TLv and TLq, for each p 00962 size_t numq = 0; 00963 size_t numv = 0; 00964 for (int p = 0; p < numprocs; p++) 00965 { 00966 if (p == (int) my_rank) 00967 continue; // do not "send" it, because it is already here 00968 Range & range_to_P = Rto[p]; 00969 // add the vertices to it 00970 if (range_to_P.empty()) 00971 continue; // nothing to send to proc p 00972 Range vertsToP; 00973 rval = mb->get_connectivity(range_to_P, vertsToP); 00974 ERRORR(rval, "can't get connectivity"); 00975 numq = numq + range_to_P.size(); 00976 numv = numv + vertsToP.size(); 00977 range_to_P.merge(vertsToP); 00978 } 00979 TupleList TLv; 00980 TupleList TLq; 00981 TLv.initialize(2, 0, 0, 3, numv); // to proc, GLOBAL ID, DP points 00982 TLv.enableWriteAccess(); 00983 00984 int sizeTuple = 2 + max_edges; // max edges could be up to MAXEDGES :) for polygons 00985 TLq.initialize(2+max_edges, 0, 1, 0, numq); // to proc, elem GLOBAL ID, connectivity[max_edges] (global ID v) 00986 // send also the corresponding red cell it will come to 00987 TLq.enableWriteAccess(); 00988 std::cout << "from proc " << my_rank << " send " << numv << " vertices and " 00989 << numq << " elements\n"; 00990 00991 for (int to_proc = 0; to_proc < numprocs; to_proc++) 00992 { 00993 if (to_proc == (int) my_rank) 00994 continue; 00995 Range & range_to_P = Rto[to_proc]; 00996 Range V = range_to_P.subset_by_type(MBVERTEX); 00997 00998 for (Range::iterator it = V.begin(); it != V.end(); it++) 00999 { 01000 EntityHandle v = *it; 01001 int index = lagr_verts.index(v);// will be the same index as the corresponding vertex in euler verts 01002 assert(-1!=index); 01003 int n = TLv.get_n(); 01004 TLv.vi_wr[2 * n] = to_proc; // send to processor 01005 TLv.vi_wr[2 * n + 1] = gids[index]; // global id needs index in the local_verts range 01006 TLv.vr_wr[3 * n] = dep_points[3 * index]; // departure position, of the node local_verts[i] 01007 TLv.vr_wr[3 * n + 1] = dep_points[3 * index + 1]; 01008 TLv.vr_wr[3 * n + 2] = dep_points[3 * index + 2]; 01009 TLv.inc_n(); 01010 } 01011 // also, prep the 2d cells for sending ... 01012 Range Q = range_to_P.subset_by_dimension(2); 01013 for (Range::iterator it = Q.begin(); it != Q.end(); it++) 01014 { 01015 EntityHandle q = *it; // this is a blue cell 01016 int global_id; 01017 rval = mb->tag_get_data(gid, &q, 1, &global_id); 01018 ERRORR(rval, "can't get gid for polygon"); 01019 int n = TLq.get_n(); 01020 TLq.vi_wr[sizeTuple * n] = to_proc; // 01021 TLq.vi_wr[sizeTuple * n + 1] = global_id; // global id of element, used to identify it ... 01022 const EntityHandle * conn4; 01023 int num_nodes; 01024 rval = mb->get_connectivity(q, conn4, num_nodes); // could be up to 10; 01025 ERRORR(rval, "can't get connectivity for quad"); 01026 if (num_nodes > MAXEDGES) 01027 ERRORR(MB_FAILURE, "too many nodes in a polygon"); 01028 for (int i = 0; i < num_nodes; i++) 01029 { 01030 EntityHandle v = conn4[i]; 01031 int index = lagr_verts.index(v); 01032 assert(-1!=index); 01033 TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index]; 01034 } 01035 for (int k = num_nodes; k < max_edges; k++) 01036 { 01037 TLq.vi_wr[sizeTuple * n + 2 + k] = 0; // fill the rest of node ids with 0; we know that the node ids start from 1! 01038 } 01039 EntityHandle redCell; 01040 rval = mb->tag_get_data(corrTag, &q, 1, &redCell); 01041 ERRORR(rval, "can't get corresponding red cell for dep cell"); 01042 TLq.vul_wr[n]=redCell; // this will be sent to remote_cells, to be able to come back 01043 TLq.inc_n(); 01044 01045 } 01046 01047 } 01048 // now we can route them to each processor 01049 // now we are done populating the tuples; route them to the appropriate processors 01050 (parcomm->proc_config().crystal_router())->gs_transfer(1, TLv, 0); 01051 (parcomm->proc_config().crystal_router())->gs_transfer(1, TLq, 0); 01052 // the elements are already in localEnts; 01053 01054 // maps from global ids to new vertex and quad handles, that are added 01055 std::map<int, EntityHandle> globalID_to_handle; 01056 // we already have vertices from lagr set; they are already in the processor, even before receiving other 01057 // verts from neighbors 01058 int k=0; 01059 for (Range::iterator vit=lagr_verts.begin(); vit!=lagr_verts.end(); vit++, k++) 01060 { 01061 globalID_to_handle[gids[k]] = *vit; // a little bit of overkill 01062 // we do know that the global ids between euler and lagr verts are parallel 01063 } 01064 /*std::map<int, EntityHandle> globalID_to_eh;*/ // do we need this one? 01065 globalID_to_eh.clear(); 01066 // now, look at every TLv, and see if we have to create a vertex there or not 01067 int n = TLv.get_n(); // the size of the points received 01068 for (int i = 0; i < n; i++) 01069 { 01070 int globalId = TLv.vi_rd[2 * i + 1]; 01071 if (globalID_to_handle.find(globalId) == globalID_to_handle.end()) 01072 { 01073 EntityHandle new_vert; 01074 double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 01075 * i + 2] }; 01076 rval = mb->create_vertex(dp_pos, new_vert); 01077 ERRORR(rval, "can't create new vertex "); 01078 globalID_to_handle[globalId] = new_vert; 01079 } 01080 } 01081 01082 // now, all dep points should be at their place 01083 // look in the local list of 2d cells for this proc, and create all those cells if needed 01084 // it may be an overkill, but because it does not involve communication, we do it anyway 01085 Range & local = Rto[my_rank]; 01086 Range local_q = local.subset_by_dimension(2); 01087 // the local should have all the vertices in lagr_verts 01088 for (Range::iterator it = local_q.begin(); it != local_q.end(); it++) 01089 { 01090 EntityHandle q = *it;// these are from lagr cells, local 01091 int gid_el; 01092 rval = mb->tag_get_data(gid, &q, 1, &gid_el); 01093 ERRORR(rval, "can't get element global ID "); 01094 globalID_to_eh[gid_el] = q; // do we need this? maybe to just mark the ones on this processor 01095 // maybe a range of global cell ids is fine? 01096 } 01097 // now look at all elements received through; we do not want to duplicate them 01098 n = TLq.get_n(); // number of elements received by this processor 01099 // a cell should be received from one proc only; so why are we so worried about duplicated elements? 01100 // a vertex can be received from multiple sources, that is fine 01101 01102 remote_cells = new TupleList(); 01103 remote_cells->initialize(2, 0, 1, 1, n); 01104 remote_cells->enableWriteAccess(); 01105 for (int i = 0; i < n; i++) 01106 { 01107 int globalIdEl = TLq.vi_rd[sizeTuple * i + 1]; 01108 int from_proc=TLq.vi_rd[sizeTuple * i ]; 01109 // do we already have a quad with this global ID, represented? 01110 if (globalID_to_eh.find(globalIdEl) == globalID_to_eh.end()) 01111 { 01112 // construct the conn quad 01113 EntityHandle new_conn[MAXEDGES]; 01114 int nnodes = -1; 01115 for (int j = 0; j < max_edges; j++) 01116 { 01117 int vgid = TLq.vi_rd[sizeTuple * i + 2 + j]; // vertex global ID 01118 if (vgid == 0) 01119 new_conn[j] = 0; 01120 else 01121 { 01122 assert(globalID_to_handle.find(vgid)!=globalID_to_handle.end()); 01123 new_conn[j] = globalID_to_handle[vgid]; 01124 nnodes = j + 1; // nodes are at the beginning, and are variable number 01125 } 01126 } 01127 EntityHandle new_element; 01128 // 01129 EntityType entType = MBQUAD; 01130 if (nnodes > 4) 01131 entType = MBPOLYGON; 01132 if (nnodes < 4) 01133 entType = MBTRI; 01134 rval = mb->create_element(entType, new_conn, nnodes, new_element); 01135 ERRORR(rval, "can't create new element "); 01136 globalID_to_eh[globalIdEl] = new_element; 01137 local_q.insert(new_element); 01138 rval = mb->tag_set_data(gid, &new_element, 1, &globalIdEl); 01139 ERRORR(rval, "can't set gid on new element "); 01140 } 01141 remote_cells->vi_wr[2*i]=from_proc; 01142 remote_cells->vi_wr[2*i+1]=globalIdEl; 01143 remote_cells->vr_wr[i] = 0.; // no contribution yet sent back 01144 remote_cells->vul_wr[i]= TLq.vul_rd[i];// this is the corresponding red cell (arrival) 01145 remote_cells->inc_n(); 01146 } 01147 // now, create a new set, covering_set 01148 rval = mb->create_meshset(MESHSET_SET, covering_set); 01149 ERRORR(rval, "can't create new mesh set "); 01150 rval = mb->add_entities(covering_set, local_q); 01151 ERRORR(rval, "can't add entities to new mesh set "); 01152 // order the remote cells tuple list, with the global id, because we will search in it 01153 //remote_cells->print("remote_cells before sorting"); 01154 moab::TupleList::buffer sort_buffer; 01155 sort_buffer.buffer_init(n); 01156 remote_cells->sort(1, &sort_buffer); 01157 sort_buffer.reset(); 01158 return MB_SUCCESS; 01159 //end copy 01160 } 01161 01162 ErrorCode Intx2Mesh::correct_intersection_points_positions() 01163 { 01164 if (parcomm) 01165 { 01166 // first, find out the edges that are shared between processors, and owned by the current processor 01167 Range shared_edges_owned; 01168 ErrorCode rval = parcomm->get_shared_entities(-1, // all other proc 01169 shared_edges_owned, 01170 1, 01171 true, // only on the interface 01172 true); // only the edges owned by the current processor 01173 ERRORR(rval, "can't get shared edges owned"); 01174 01175 rval = parcomm->settle_intersection_points(RedEdges, shared_edges_owned, extraNodesVec, epsilon_1); 01176 ERRORR(rval, "can't settle intx points"); 01177 } 01178 return MB_SUCCESS; 01179 } 01180 #endif /* USE_MPI */ 01181 } /* namespace moab */