Charm++: Charm Source Code Documentation

00001 
00005 #include "bulk_adapt_ops.h"
00006 #include "import.h"
00007 
00008 #define BULK_DEBUG(x) x
00009 
00011 BulkAdapt::BulkAdapt(int meshid,FEM_Mesh *mPtr, int partID, 
00012              CProxy_ParFUMShadowArray sa_proxy)
00013 {
00014   meshID = meshid;
00015   meshPtr = mPtr;
00016   partitionID = partID;
00017   shadowProxy = sa_proxy;
00018   localShadow = meshPtr->parfumSA;
00019   for (int i=0; i<10; i++) { freeTable[i] = 1; }
00020   firstFree = 0;
00021 }
00022 
00024 BulkAdapt::~BulkAdapt()
00025 {
00026 }
00027 
00029 void BulkAdapt::pup(PUP::er &p)
00030 {
00031 // MIGRATION NOTES:
00032 // * will need to fix meshPtr when this partition migrates
00033 // * will need to fix localShadow too.  shadowProxy? should be fine?
00034 }
00035 
00036 /* BULK MESH OPERATIONS: These are all called locally, but may invoke
00037    remote operations. */
00038 
00040 
00045 bool BulkAdapt::edge_bisect(int elemID, int elemType, int edgeID, int dim, RegionID lockRegionID)
00046 {
00047   if (dim == 2) {
00048     return edge_bisect_2D(elemID, elemType, edgeID);
00049   }
00050   else if (dim == 3) {
00051     return edge_bisect_3D(elemID, elemType, edgeID, lockRegionID);
00052   }
00053   else return false;
00054 }
00055 
00057 /*
00058                o                                     o   
00059               /|\                           /|\  
00060  startElem   / | \   nbrElem           startElem   / | \   nbrElem   
00061             /  |  \                       /  |  \    
00062            /   |   \                         /   |   \   
00063           /    |    \                       /    |    \  
00064          o     |     o                     o-----o-----o     
00065           \    |    /                       \    |    /  
00066            \   |   /                         \   |   /   
00067             \  |  /            splitElem  \  |  /  nbrSplitElem
00068              \ | /                         \ | /     
00069               \|/                           \|/  
00070                o                                 o         
00071  */
00072 bool BulkAdapt::edge_bisect_2D(int elemID, int elemType, int edgeID)
00073 {
00074   //BULK_DEBUG(CkPrintf("[%d] BulkAdapt::edge_bisect_2D starts elemID %d \n",partitionID,elemID));
00075   // lock partitions for the two involved elements
00076   adaptAdj elemsToLock[2];
00077   adaptAdj startElem(partitionID, elemID, elemType);
00078   adaptAdj nbrElem = *getAdaptAdj(meshID, elemID, elemType, edgeID);
00079   //BULK_DEBUG(printf("[%d] neighbor of elem %d is elem (%d,%d) \n",partitionID,elemID,nbrElem.partID,nbrElem.localID);)
00080   elemsToLock[0] = startElem;
00081   elemsToLock[1] = nbrElem;
00082   RegionID lockRegionID;
00083   if (localShadow->lockRegion(2, elemsToLock, &lockRegionID,0.0)) {
00084     //BULK_DEBUG(CkPrintf("[%d] Lock obtained.\n",partitionID););
00085   }
00086   else {
00087     //BULK_DEBUG(CkPrintf("[%d] Lock not obtained.\n",partitionID););
00088     return false;
00089   }
00090 
00091   // ******** LOCAL OPS *********
00092   int node1idx, node2idx, newNodeID;
00093   adaptAdj splitElem;
00094   // split the local element, i.e. the first "side"
00095   one_side_split_2D(startElem, splitElem, edgeID, &node1idx, &node2idx, &newNodeID, true);
00096 
00097   if ((nbrElem.partID > -1) && (nbrElem.partID == partitionID)) { // if nbrElem exists and is local...
00098     // PRE: neighbor-side operations
00099     FEM_Elem &elem = meshPtr->elem[elemType]; // elem is local elements
00100     int *nbrConn = elem.connFor(nbrElem.localID);
00101     int relNode1 = getRelNode(node1idx, nbrConn, 3);
00102     int relNode2 = getRelNode(node2idx, nbrConn, 3);
00103     int nbrEdgeID = getEdgeID(relNode1, relNode2, 3, 2);
00104 
00105     int nbrNode1, nbrNode2;
00106     adaptAdj nbrSplitElem;
00107     // split the local neighbor element, i.e. the second "side"
00108     one_side_split_2D(nbrElem, nbrSplitElem, nbrEdgeID, &nbrNode1, &nbrNode2, &newNodeID, false);
00109 
00110     // now fix the adjacencies across the new edge to the two new elements
00111     adaptAdj *splitElemAdaptAdj = getAdaptAdj(meshID, splitElem.localID, splitElem.elemType, 0);
00112     adaptAdj *nbrSplitElemAdaptAdj = getAdaptAdj(meshID, nbrSplitElem.localID, nbrSplitElem.elemType, 0);
00113     nbrSplitElemAdaptAdj[nbrEdgeID] = splitElem;
00114     //BULK_DEBUG(printf("[%d] For nbrSplitElem %d set adjacency to %d across splitEdge\n",partitionID,nbrSplitElem.localID,splitElem.localID);)
00115     // POST: start-side operations
00116     splitElemAdaptAdj[edgeID] = nbrSplitElem;
00117     //BULK_DEBUG(printf("[%d] For splitElem %d set adjacency to %d across splitEdge\n",partitionID,splitElem.localID,nbrSplitElem.localID);)
00118   }
00119   else if (nbrElem.partID == -1) { // startElem's edgeID is on domain boundary
00120     adaptAdj *splitElemAdaptAdj = getAdaptAdj(meshID, splitElem.localID, splitElem.elemType, 0);
00121     splitElemAdaptAdj[edgeID] = adaptAdj(-1, -1, 0);
00122     //BULK_DEBUG(printf("[%d] For splitElem %d splitEdge is on the domain boundary.\n",partitionID,splitElem.localID);)
00123   }
00124   else { // nbrElem exists and is remote
00125     int chunks[1] = {nbrElem.partID};
00126     make_node_shared(newNodeID, 1, &chunks[0]);
00127     int new_idxl, n1_idxl, n2_idxl;
00128     new_idxl = get_idxl_for_node(newNodeID, nbrElem.partID);
00129     n1_idxl = get_idxl_for_node(node1idx, nbrElem.partID);
00130     n2_idxl = get_idxl_for_node(node2idx, nbrElem.partID);
00131 
00132     // make sync call on partition nbrElem.partID
00133     // SEND: nbrElem, splitElem, new_idxl, n1_idxl, n2_idxl, partitionID
00134     // RETURNS: nbrSplitElem
00135     adaptAdjMsg *am = shadowProxy[nbrElem.partID].remote_bulk_edge_bisect_2D(nbrElem, splitElem, new_idxl, n1_idxl, n2_idxl, partitionID);
00136     
00137     adaptAdj nbrSplitElem = am->elem;
00138     // now fix the adjacencies across the new edge to remote new element
00139     adaptAdj *splitElemAdaptAdj = getAdaptAdj(meshID, splitElem.localID, splitElem.elemType, 0);
00140     splitElemAdaptAdj[edgeID] = nbrSplitElem;
00141     //BULK_DEBUG(printf("[%d] For splitElem %d set adjacency to %d across splitEdge\n",partitionID,splitElem.localID,nbrSplitElem.localID);)
00142   }    
00143 
00144   // unlock the two partitions
00145   localShadow->unlockRegion(lockRegionID);
00146   getAndDumpAdaptAdjacencies(meshID, meshPtr->nElems(), elemType, partitionID);
00147   return true;
00148 }
00149 
00150 void BulkAdapt::dumpConn()
00151 {
00152   FEM_Elem &elem = meshPtr->elem[0]; // elem is local elements
00153   int nelems=FEM_Mesh_get_length(meshID, FEM_ELEM+0); // Get number of elements
00154   for (int i=0; i<nelems; i++) {
00155     int *conn = elem.connFor(i);
00156     CkPrintf("[%d] %d: %d %d %d %d\n", partitionID, i, conn[0], conn[1], conn[2], conn[3]);
00157   }
00158 }
00159 
00160 adaptAdj BulkAdapt::remote_edge_bisect_2D(adaptAdj nbrElem, adaptAdj splitElem, int new_idxl, int n1_idxl, int n2_idxl, int remotePartID)
00161 {
00162   int node1idx, node2idx, newNodeID;
00163   node1idx = get_node_from_idxl(n1_idxl, remotePartID);
00164   node2idx = get_node_from_idxl(n2_idxl, remotePartID);
00165   
00166   FEM_Elem &elem = meshPtr->elem[nbrElem.elemType]; // elem is local elements
00167   int *nbrConn = elem.connFor(nbrElem.localID);
00168   int relNode1 = getRelNode(node1idx, nbrConn, 3);
00169   int relNode2 = getRelNode(node2idx, nbrConn, 3);
00170   int nbrEdgeID = getEdgeID(relNode1, relNode2, 3, 2);
00171 
00172   FEM_DataAttribute *coord = meshPtr->node.getCoord(); // entire local coords
00173   double *node1coords = (coord->getDouble()).getRow(node1idx); // ptrs to ACTUAL coords!
00174   double *node2coords = (coord->getDouble()).getRow(node2idx); // ptrs to ACTUAL coords!
00175   double bisectCoords[2];
00176   midpoint(node1coords, node2coords, 2, &bisectCoords[0]);
00177   newNodeID = add_node(2, &bisectCoords[0]);
00178   int chunks[2] = {nbrElem.partID, remotePartID};
00179   make_node_shared(newNodeID, 2, chunks);
00180   int local_new_idxl = get_idxl_for_node(newNodeID, remotePartID);
00181   if (new_idxl != local_new_idxl)
00182     printf("ERROR: Partition %d added shared node at different idxl index %d than other copy at %d on partition %d!", 
00183        nbrElem.partID, local_new_idxl, new_idxl, remotePartID);
00184 
00185   int nbrNode1, nbrNode2;
00186   adaptAdj nbrSplitElem;
00187   // split the local neighbor element, i.e. the second "side"
00188   one_side_split_2D(nbrElem, nbrSplitElem, nbrEdgeID, &nbrNode1, &nbrNode2, &newNodeID, false);
00189 
00190   adaptAdj *nbrSplitElemAdaptAdj = getAdaptAdj(meshPtr, nbrSplitElem.localID, nbrSplitElem.elemType, 0);
00191   nbrSplitElemAdaptAdj[nbrEdgeID] = splitElem;
00192   //BULK_DEBUG(printf("[%d] For nbrSplitElem %d set adjacency to %d across splitEdge\n",partitionID,nbrSplitElem.localID,splitElem.localID);)
00193   
00194   return nbrSplitElem;
00195 }
00196 
00197 int BulkAdapt::lock_3D_region(int elemID, int elemType, int edgeID, double prio, RegionID *lockRegionID)
00198 { // Lock region defined by all elements surrounding edgeID, i.e. the "star"
00199   int numElemsToLock=0, success=-1;
00200   adaptAdj *elemsToLock;
00201   adaptAdj startElem(partitionID, elemID, elemType);
00202 
00203   BULK_DEBUG(CkPrintf("[%d] BulkAdapt::lock_3D_region acquiring list of elements to build locked region.\n",
00204               partitionID));
00205   get_elemsToLock(startElem, &elemsToLock, edgeID, &numElemsToLock);
00206 
00207   success = localShadow->lockRegion(numElemsToLock, elemsToLock, lockRegionID, prio);
00208   free(elemsToLock);
00209   if (success==2) { // got the lock free and clear; refine the element
00210     BULK_DEBUG(CkPrintf("[%d] BulkAdapt::lock_3D_region: Lock (%d,%d,%6.4f) obtained.\n", partitionID,
00211             lockRegionID->chunkID, lockRegionID->localID, lockRegionID->prio));
00212   }
00213   else if (success==1) { // lock is pending
00214     BULK_DEBUG(CkPrintf("[%d] BulkAdapt::lock_3D_region: Lock (%d,%d,%6.4f) pending.\n", partitionID,
00215             lockRegionID->chunkID, lockRegionID->localID, lockRegionID->prio));
00216   }
00217   else if (success==0) { // lock failed
00218     BULK_DEBUG(CkPrintf("[%d] BulkAdapt::edge_bisect_3D: Lock (%d,%d,%6.4f) not obtained.\n", partitionID,
00219             lockRegionID->chunkID, lockRegionID->localID, lockRegionID->prio));
00220   }
00221   else {
00222     CkPrintf("Lock status=%d\n", success);
00223     CkAbort("Invalid lock return status!\n");
00224   }
00225   return success;
00226 }
00227 
00228 void BulkAdapt::unlock_3D_region(RegionID lockRegionID)
00229 { // unlock the region
00230   BULK_DEBUG(CkPrintf("[%d] BulkAdapt::unlock_3D_region: Unlocking lock (%d,%d,%6.4f).\n", partitionID,
00231               lockRegionID.chunkID, lockRegionID.localID, lockRegionID.prio));
00232   localShadow->unlockRegion(lockRegionID);
00233 }
00234 
00235 void BulkAdapt::unpend_3D_region(RegionID lockRegionID)
00236 { // unlock the region
00237   BULK_DEBUG(CkPrintf("[%d] BulkAdapt::unpend_3D_region: Unlocking lock (%d,%d,%6.4f).\n", partitionID,
00238               lockRegionID.chunkID, lockRegionID.localID, lockRegionID.prio));
00239   localShadow->unpendRegion(lockRegionID);
00240 }
00241 
00242 
00244 bool BulkAdapt::edge_bisect_3D(int elemID, int elemType, int edgeID, RegionID lockRegionID)
00245 { // ASSERT: An edge can only be on one surface. Region around edge is locked, along with appropriate 
00246   // partitions and IDXL lists.
00247   BULK_DEBUG(CkPrintf("[%d] BulkAdapt::edge_bisect_3D starts at elemID %d \n",partitionID,elemID));
00248   //dumpConn();
00249   //getAndDumpAdaptAdjacencies(meshID, meshPtr->nElems(), elemType, partitionID);
00250 
00251   int numElemsToLock=0;
00252   adaptAdj *elemsToLock;
00253   adaptAdj startElem(partitionID, elemID, elemType);
00254 
00255   CkAssert(lockRegionID == localShadow->holdingLock);
00256 
00257   BULK_DEBUG(CkPrintf("[%d] BulkAdapt::edge_bisect_3D acquiring list of elements surrounding edge.\n",
00258               partitionID));
00259   get_elemsToLock(startElem, &elemsToLock, edgeID, &numElemsToLock);
00260 
00261   // Find the nodes on the edge to be bisected
00262   int localNodes[3], localRelNodes[2]; 
00263   double nodeCoords[9];
00264   FEM_Elem &elems = meshPtr->elem[elemType]; // elems is all local elements
00265   FEM_DataAttribute *coord = meshPtr->node.getCoord(); // all local coords
00266   int *conn = elems.connFor(elemID); // conn points at elemID's ACTUAL data!
00267 
00268   // edgeID has element-relative node indices localRelNodes[0], localRelNodes[1]
00269   getRelNodes(edgeID, 4, &(localRelNodes[0]), &(localRelNodes[1]));
00270   localNodes[0] = conn[localRelNodes[0]];
00271   localNodes[1] = conn[localRelNodes[1]];
00272   double *n0co = (coord->getDouble()).getRow(localNodes[0]);
00273   double *n1co = (coord->getDouble()).getRow(localNodes[1]);
00274   nodeCoords[0] = n0co[0];
00275   nodeCoords[1] = n0co[1];
00276   nodeCoords[2] = n0co[2];
00277   nodeCoords[3] = n1co[0];
00278   nodeCoords[4] = n1co[1];
00279   nodeCoords[5] = n1co[2];
00280 
00281   // Create the new shared node that will be added
00282   midpoint(&(nodeCoords[0]), &(nodeCoords[3]), 3, &(nodeCoords[6]));
00283   localNodes[2] = add_node(3, &(nodeCoords[6]));
00284   int *chunks = (int *)malloc(numElemsToLock*sizeof(int));
00285   int numParts=0;
00286 
00287   for (int i=0; i<numElemsToLock; i++) {
00288     chunks[i] = -1;
00289     int j;
00290     for (j=0; j<numParts; j++) {
00291       if (chunks[j] == elemsToLock[i].partID) {
00292     break;
00293       }
00294     }
00295     chunks[j] = elemsToLock[i].partID;
00296     if (j==numParts) numParts++;
00297   }
00298 
00299   if (numParts > 1)
00300     make_node_shared(localNodes[2], numParts, &chunks[0]);
00301 
00302   free(chunks);
00303 
00304   // Perform the splits on all affected elements
00305   int numRemote=0;
00306   int tableID = getTableID();
00307   elemPairs[tableID] = (adaptAdj *)malloc(2*numElemsToLock*sizeof(adaptAdj));
00308 
00309   for (int i=0; i<numElemsToLock; i++) {
00310     if (elemsToLock[i].partID != partitionID) {
00311       elemPairs[tableID][2*i] = elemsToLock[i];
00312       numRemote++;
00313       shadowProxy[elemsToLock[i].partID].
00314     handle_split_3D(partitionID, i, tableID, elemsToLock[i], lockRegionID,
00315             get_idxl_for_node(localNodes[0], elemsToLock[i].partID), 
00316             get_idxl_for_node(localNodes[1], elemsToLock[i].partID), 
00317             get_idxl_for_node(localNodes[2], elemsToLock[i].partID));
00318     }
00319   }
00320 
00321   for (int i=0; i<numElemsToLock; i++) {
00322     if (elemsToLock[i].partID == partitionID) {
00323       elemPairs[tableID][2*i] = elemsToLock[i];
00324       adaptAdj *splitElem = local_split_3D(elemsToLock[i], localNodes[0], 
00325                        localNodes[1], localNodes[2]);
00326       elemPairs[tableID][2*i+1] = (*splitElem);
00327       delete splitElem;
00328     }
00329   }
00330 
00331   if (numRemote > 0) {
00332     localShadow->recv_splits(tableID, numRemote);
00333   }
00334 
00335   // Now we need to send out the elemPairs and make sure all
00336   // adjacencies are updated
00337   for (int i=0; i<numElemsToLock; i++) {
00338     if (elemsToLock[i].partID != partitionID) {
00339       shadowProxy[elemsToLock[i].partID].
00340     update_asterisk_3D(partitionID, i, elemsToLock[i], numElemsToLock, 
00341                elemPairs[tableID], lockRegionID,
00342                get_idxl_for_node(localNodes[0], elemsToLock[i].partID), 
00343                get_idxl_for_node(localNodes[1], elemsToLock[i].partID),
00344                get_idxl_for_node(localNodes[2], elemsToLock[i].partID));
00345     }
00346   }
00347   for (int i=0; i<numElemsToLock; i++) {
00348     if (elemsToLock[i].partID == partitionID) {
00349       local_update_asterisk_3D(i, elemsToLock[i], numElemsToLock, elemPairs[tableID],
00350                    localNodes[0], localNodes[1], localNodes[2]);
00351     }
00352   }
00353 
00354   // clean up elemPairs
00355   freeTableID(tableID);
00356   free(elemsToLock);
00357 
00358   //dumpConn();
00359   //getAndDumpAdaptAdjacencies(meshID, meshPtr->nElems(), elemType, partitionID);
00360 
00361   BULK_DEBUG(CkPrintf("[%d] BulkAdapt::edge_bisect_3D of elem %d successful.\n",partitionID, elemID));  
00362   return true;
00363 }
00364 
00366 
00371 int BulkAdapt::edge_flip(int elemID, int edgeID)
00372 {
00373   CkPrintf("BulkAdapt::edge_flip not yet implemented.\n");
00374   return 0;
00375 }
00376 
00378 
00383 int BulkAdapt::flip_23(int elemID, int faceID)
00384 {
00385   CkPrintf("BulkAdapt::flip_23 not yet implemented.\n");
00386   return 0;
00387 }
00388 
00390 
00395 int BulkAdapt::flip_32(int elemID, int edgeID)
00396 {
00397   CkPrintf("BulkAdapt::flip_32 not yet implemented.\n");
00398   return 0;
00399 }
00400 
00402 
00407 int BulkAdapt::edge_collapse(int elemID, int edgeID)
00408 {
00409   CkPrintf("BulkAdapt::edge_collapse not yet implemented.\n");
00410   return 0;
00411 }
00412 
00413 void BulkAdapt::one_side_split_2D(adaptAdj &startElem, adaptAdj &splitElem, 
00414                   int edgeID, int *node1idx, int *node2idx, 
00415                   int *newNodeID, bool startSide)
00416 {
00417   // get startElem.localID's conn.
00418   FEM_Elem &elem = meshPtr->elem[startElem.elemType]; // elem is local elements
00419   int *startConn = elem.connFor(startElem.localID); // startConn points at startElem.localID's ACTUAL data!
00420   // let edgeID have element-relative node indices relNode1, relNode2.
00421   int relNode1 = edgeID, relNode2 = (edgeID+1)%3;
00422   if (!startSide) {
00423     relNode1 = (edgeID+1)%3;
00424     relNode2 = edgeID;
00425   }
00426   *node1idx = startConn[relNode1];
00427   *node2idx = startConn[relNode2];
00428 
00429 
00430   //BULK_DEBUG(printf("[%d] one_side_split_2D called for elem %d edge %d nodes %d %d \n",partitionID,startElem.localID,edgeID,*node1idx,*node2idx);)
00431   // get node1coords and node2Coords. find midpoint: bisectCoords
00432   FEM_DataAttribute *coord = meshPtr->node.getCoord(); // entire local coords
00433   double *node1coords = (coord->getDouble()).getRow(*node1idx); // ptrs to ACTUAL coords!
00434   double *node2coords = (coord->getDouble()).getRow(*node2idx); // ptrs to ACTUAL coords!
00435   double bisectCoords[2];
00436   midpoint(node1coords, node2coords, 2, &bisectCoords[0]);
00437   // add node at bisectCoords to get bisectNode
00438   int bisectNode;
00439   if (startSide) { // For local ops, startSide creates the node, nbr uses it
00440     bisectNode = add_node(2, &bisectCoords[0]);
00441     *newNodeID = bisectNode;
00442   }
00443   else { // local nbr just uses node created on startSide
00444     // if neighbor is remote, the newNodeID is set in remote_edge_bisect_2D before calling this function
00445     bisectNode = *newNodeID;
00446   }
00447   // duplicate conn in conn2.  
00448   int splitConn[3];
00449   memcpy(&splitConn[0], startConn, 3*sizeof(int));
00450   // in startConn, replace node2idx with bisectNode
00451   startConn[relNode2] = bisectNode; // ACTUAL data was changed!
00452   //BULK_DEBUG(printf("[%d] conn[%d] of elem %d changed to node %d\n",partitionID,relNode2,startElem.localID,bisectNode);)
00453   // in splitConn, replace node1idx with bisectNode
00454   splitConn[relNode1] = bisectNode;
00455   // add element split with splitConn
00456   int splitElemID = add_element(startElem.elemType, 3, &splitConn[0], elem.getMeshSizing(startElem.localID));
00457   //BULK_DEBUG(printf("[%d] new element %d with conn %d %d %d added \n", partitionID, splitElemID, splitConn[0], splitConn[1], splitConn[2]);)
00458   // copy startElem.localID's adapt adj for all edges.
00459   splitElem = adaptAdj(partitionID, splitElemID, startElem.elemType);
00460   adaptAdj *startElemAdaptAdj = getAdaptAdj(meshPtr, startElem.localID, startElem.elemType, 0);
00461   adaptAdj *splitElemAdaptAdj = getAdaptAdj(meshPtr, splitElem.localID, startElem.elemType, 0);
00462   memcpy(splitElemAdaptAdj, startElemAdaptAdj, 3*sizeof(adaptAdj));
00463   adaptAdj startElemNbr;  // startElem's original nbr on the edge that will now border with splitElem
00464   if (startSide) {
00465     // update startElemAdaptAdj for edge (edgeID+1)%3 to local splitElem
00466     startElemNbr = startElemAdaptAdj[(edgeID+1)%3];
00467     startElemAdaptAdj[(edgeID+1)%3] = splitElem;
00468     // update splitElemAdaptAdj for edge (edgeID+2)%3 to local startElem
00469     splitElemAdaptAdj[(edgeID+2)%3] = startElem;
00470     //BULK_DEBUG(printf("[%d] For startElem %d edge %d is now set to %d\n",partitionID,startElem.localID,(edgeID+1)%3,splitElem.localID));
00471     //BULK_DEBUG(printf("[%d] For splitElem %d edge %d is now set to %d\n",partitionID,splitElem.localID,(edgeID+2)%3,startElem.localID));
00472   }
00473   else {
00474     // update startElemAdaptAdj for edge (edgeID+1)%3 to local splitElem
00475     startElemNbr = startElemAdaptAdj[(edgeID+2)%3];
00476     startElemAdaptAdj[(edgeID+2)%3] = splitElem;
00477     // update splitElemAdaptAdj for edge (edgeID+2)%3 to local startElem
00478     splitElemAdaptAdj[(edgeID+1)%3] = startElem;
00479     //BULK_DEBUG(printf("[%d] For startElem %d edge %d is now set to %d\n",partitionID,startElem.localID,(edgeID+2)%3,splitElem.localID));
00480     //BULK_DEBUG(printf("[%d] For splitElem %d edge %d is now set to %d\n",partitionID,splitElem.localID,(edgeID+1)%3,startElem.localID));
00481   }
00482   if (startElemNbr.partID == startElem.partID) {
00483     replaceAdaptAdj(meshPtr, startElemNbr, startElem, splitElem);
00484     //BULK_DEBUG(printf("[%d] For startElemNbr %d replaced startElem %d with splitElem %d\n",partitionID,startElemNbr.localID,startElem.localID,splitElem.localID);)
00485   }
00486   else if (startElemNbr.partID != -1) { // startElemNbr exists and is remote
00487     // need to call replaceAdaptAdj on startElemNbr.partID
00488     shadowProxy[startElemNbr.partID].remote_adaptAdj_replace(startElemNbr, startElem, splitElem); 
00489   }
00490   // interpolate nodal data, copy startElem data to splitElem
00491   CkPrintf("WARNING: Data transfer not yet implemented.\n");
00492 }
00493 
00494 
00495 
00496 /* COMMUNICATION HELPERS FOR BULK ADAPTIVITY OPERATIONS 
00497    ARE LOCATED IN ParFUM_SA. */
00498 
00499 void BulkAdapt::remote_adaptAdj_replace(adaptAdj elem, adaptAdj oldElem, adaptAdj newElem)
00500 {
00501   replaceAdaptAdj(meshPtr, elem, oldElem, newElem);
00502 }
00503 
00504 void BulkAdapt::remote_edgeAdj_replace(int remotePartID, adaptAdj elem, adaptAdj oldElem, adaptAdj newElem, int n1_idxl, int n2_idxl)
00505 {
00506   int n1 = get_node_from_idxl(n1_idxl, remotePartID), 
00507     n2 = get_node_from_idxl(n2_idxl, remotePartID);
00508   // Find the relative nodes on the edge to be bisected
00509   int relNodes[4]; 
00510   FEM_Elem &elems = meshPtr->elem[elem.elemType]; // elems is all local elements
00511   int *conn = elems.connFor(elem.localID); // conn points at elemID's ACTUAL data!
00512   fillNodes(relNodes, n1, n2, conn);
00513   // find which edgeID is bisected
00514   int edgeID = getEdgeID(relNodes[0], relNodes[1], 4, 3);
00515 
00516   replaceAdaptAdjOnEdge(meshPtr, elem, oldElem, newElem, edgeID);
00517 }
00518 
00519 void BulkAdapt::remote_edgeAdj_add(int remotePartID, adaptAdj adj, adaptAdj splitElem, int n1_idxl, int n2_idxl)
00520 {
00521   int n1 = get_node_from_idxl(n1_idxl, remotePartID), 
00522     n2 = get_node_from_idxl(n2_idxl, remotePartID);
00523   // Find the relative nodes on the edge to be bisected
00524   int relNodes[4]; 
00525   FEM_Elem &elems = meshPtr->elem[adj.elemType]; // elems is all local elements
00526   int *conn = elems.connFor(adj.localID); // conn points at elemID's ACTUAL data!
00527   fillNodes(relNodes, n1, n2, conn);
00528   // find which edgeID is bisected
00529   int edgeID = getEdgeID(relNodes[0], relNodes[1], 4, 3);
00530 
00531   addToAdaptAdj(meshPtr, adj, edgeID, splitElem);
00532 }
00533 
00534 void BulkAdapt::handle_split_3D(int remotePartID, int pos, int tableID, adaptAdj elem, 
00535                 RegionID lockRegionID, int n1_idxl, int n2_idxl, int n5_idxl)
00536 {
00537   int n1 = get_node_from_idxl(n1_idxl, remotePartID), 
00538     n2 = get_node_from_idxl(n2_idxl, remotePartID), 
00539     n5;
00540 
00541   if (is_node_in_idxl(n5_idxl, remotePartID)) {
00542     n5 = get_node_from_idxl(n5_idxl, remotePartID);
00543   }
00544   else {
00545     FEM_DataAttribute *coord = meshPtr->node.getCoord(); // all local coords
00546     double *n0co = (coord->getDouble()).getRow(n1);
00547     double *n1co = (coord->getDouble()).getRow(n2);
00548     double nodeCoords[9];
00549     nodeCoords[0] = n0co[0];
00550     nodeCoords[1] = n0co[1];
00551     nodeCoords[2] = n0co[2];
00552     nodeCoords[3] = n1co[0];
00553     nodeCoords[4] = n1co[1];
00554     nodeCoords[5] = n1co[2];
00555     midpoint(&(nodeCoords[0]), &(nodeCoords[3]), 3, &(nodeCoords[6]));
00556     n5 = add_node(3, &(nodeCoords[6]));  
00557     
00558     // Find the relative nodes on the edge to be bisected
00559     int relNodes[4]; 
00560     FEM_Elem &elems = meshPtr->elem[elem.elemType]; // elems is all local elements
00561     int *conn = elems.connFor(elem.localID); // conn points at elemID's ACTUAL data!
00562     
00563     fillNodes(relNodes, n1, n2, conn);
00564     // find which edgeID is bisected
00565     int edgeID = getEdgeID(relNodes[0], relNodes[1], 4, 3);
00566     
00567     // find elements on edge to be bisected
00568     int numElemsToLock = 0;
00569     adaptAdj *elemsToLock;
00570     get_elemsToLock(elem, &elemsToLock, edgeID, &numElemsToLock);
00571     
00572     // find chunks that share the edge to be bisected
00573     int *chunks = (int *)malloc(numElemsToLock*sizeof(int));
00574     int numParts=0;
00575     for (int i=0; i<numElemsToLock; i++) {
00576       chunks[i] = -1;
00577       int j;
00578       for (j=0; j<numParts; j++) {
00579     if (chunks[j] == elemsToLock[i].partID) {
00580       break;
00581     }
00582       }
00583       chunks[j] = elemsToLock[i].partID;
00584       if (j==numParts) numParts++;
00585     }
00586     
00587     if (numParts > 1)
00588       make_node_shared(n5, numParts, &chunks[0]);
00589     
00590     free(chunks);
00591     free(elemsToLock);
00592   }
00593 
00594   adaptAdj *splitElem = local_split_3D(elem, n1, n2, n5);
00595   shadowProxy[remotePartID].recv_split_3D(pos, tableID, elem, *splitElem);
00596   delete splitElem;
00597 }
00598 
00599 void BulkAdapt::recv_split_3D(int pos, int tableID, adaptAdj elem,
00600                   adaptAdj splitElem) 
00601 {
00602   assert(elemPairs[tableID][2*pos] == elem);
00603   elemPairs[tableID][2*pos+1] = splitElem;
00604   numGathered[tableID]++;
00605 }
00606 
00607 bool BulkAdapt::all_splits_received(int tableID, int expectedSplits)
00608 {
00609   return (expectedSplits == numGathered[tableID]);
00610 }
00611 
00612 void BulkAdapt::update_asterisk_3D(int remotePartID, int i, adaptAdj elem, 
00613                    int numElemPairs, adaptAdj *elemPairs, RegionID lockRegionID,
00614                    int n1_idxl, int n2_idxl, int n5_idxl)
00615 {
00616   CkAssert(remotePartID != elem.partID);
00617   int n1 = get_node_from_idxl(n1_idxl, remotePartID), 
00618     n2 = get_node_from_idxl(n2_idxl, remotePartID), 
00619     n5 = get_node_from_idxl(n5_idxl, remotePartID);
00620   //printf("[%d] update_asterisk_3D(from %d): elem %d n1=%d n2=%d n5=%d\n",
00621   //elem.partID, remotePartID, elem.localID, n1, n2, n5);
00622   local_update_asterisk_3D(i, elem, numElemPairs, elemPairs, n1, n2, n5);
00623 }
00624 
00625 
00626 /* LOCAL HELPERS FOR BULK ADAPTIVITY OPERATIONS */
00627 
00632 int BulkAdapt::add_element(int elemType,int nodesPerElem,int *conn, double sizing){ 
00633   //since the element array might be resized we need to set the correct thread
00634   //context before we call get_next_invalid
00635   localShadow->setRunningTCharm();
00636   FEM_Elem &elem = meshPtr->elem[elemType];
00637   int newElem = elem.get_next_invalid();
00638   elem.connIs(newElem,conn);
00639   elem.setMeshSizing(newElem, sizing);
00640 
00641   int nAdj;
00642   adaptAdj* adaptAdjacencies = lookupAdaptAdjacencies(
00643         meshPtr, elemType, &nAdj);
00644   for (int a = 0; a<nAdj; ++a) {
00645     //adaptAdjacencies[newElem*nAdj + a].partID = TCHARM_Element();
00646     adaptAdjacencies[newElem*nAdj + a].partID = -1;
00647     adaptAdjacencies[newElem*nAdj + a].localID = -1;
00648     adaptAdjacencies[newElem*nAdj + a].elemType = elemType;
00649   }
00650   CkVec<adaptAdj>** adaptEdgeAdjacencies = 
00651     lookupEdgeAdaptAdjacencies(meshPtr,elemType,&nAdj);
00652   for (int a=0; a<nAdj; ++a) {
00653     adaptEdgeAdjacencies[newElem*nAdj + a] = new CkVec<adaptAdj>;
00654   }
00655   return newElem;
00656 }
00657 
00659 void BulkAdapt::update_element_conn(int elemType,int elemID,int nodesPerElem,int *conn){
00660   FEM_Elem &elem = meshPtr->elem[elemType];
00661   elem.connIs(elemID,conn);
00662 }
00663 
00664 bool BulkAdapt::isLongest(int elem, int elemType, double len) {
00665   FEM_Elem &elems = meshPtr->elem[elemType];
00666   int *conn = elems.connFor(elem);
00667   if ((len < length(conn[0], conn[1], 3)) || (len < length(conn[0], conn[2], 3)) || (len < length(conn[0], conn[3], 3)) || (len < length(conn[1], conn[2], 3)) || (len < length(conn[1], conn[3], 3)) || (len < length(conn[2], conn[3], 3))) {
00668     return false;
00669   }
00670   return true;
00671 }
00672 
00673 
00674 
00680 int BulkAdapt::add_node(int dim,double *coords){ 
00681   //since the node array might be resized we need to set the correct thread
00682   //context before we call get_next_invalid
00683   localShadow->setRunningTCharm();
00684   int newNode = meshPtr->node.get_next_invalid();
00685   FEM_DataAttribute *coord = meshPtr->node.getCoord();
00686   (coord->getDouble()).setRow(newNode,coords);
00687   AllocTable2d<int> &lockTable = ((FEM_DataAttribute *)(meshPtr->node.lookup(FEM_ADAPT_LOCK,"lockRegion")))->getInt();
00688   
00689   lockTable[newNode][0] = -1;
00690   lockTable[newNode][1] = -1;
00691   return newNode;
00692 }
00693 
00694 
00695 double BulkAdapt::length(int n1, int n2, int dim) {
00696   FEM_DataAttribute *coord = meshPtr->node.getCoord();
00697   double *coordsn1 = (coord->getDouble()).getRow(n1);
00698   double *coordsn2 = (coord->getDouble()).getRow(n2);
00699   return length(coordsn1, coordsn2, dim);
00700 }
00701 
00702 double BulkAdapt::length(double *n1, double *n2, int dim) {
00703   double d, ds_sum=0.0;
00704   for (int i=0; i<dim; i++) {
00705     if(n1[i]<-1.0 || n2[i]<-1.0) return -2.0;
00706     d = n1[i] - n2[i];
00707     ds_sum += d*d;
00708   }
00709   return (sqrt(ds_sum));
00710 }
00711 
00712 
00713 
00716 void BulkAdapt::update_node_coord(int nodeID,int dim,double *coords){
00717   FEM_DataAttribute *coord = meshPtr->node.getCoord();
00718   (coord->getDouble()).setRow(nodeID,coords);
00719 }
00720 
00721 void BulkAdapt::make_node_shared(int nodeID,int numSharedChunks,int *sharedChunks){
00722   for(int i=0;i<numSharedChunks;i++){
00723     IDXL_List &sharedList = meshPtr->node.shared.addList(sharedChunks[i]);
00724     sharedList.push_back(nodeID);
00725   }
00726   meshPtr->node.shared.flushMap();
00727 }
00728 
00729 int BulkAdapt::get_idxl_for_node(int nodeID, int partID) 
00730 {
00731   IDXL_List *list = meshPtr->node.shared.getIdxlListN(partID);
00732   CkAssert(list!=NULL);
00733   for (int i=0; i<list->size(); i++) {
00734     if ((*list)[i] == nodeID) {
00735       return i;
00736     }
00737   }
00738   CkAssert(0);
00739   return -1;
00740 }
00741 
00742 int BulkAdapt::get_node_from_idxl(int node_idxl, int partID)
00743 {
00744   IDXL_List *list = meshPtr->node.shared.getIdxlListN(partID);
00745   CkAssert(list!=NULL);
00746   CkAssert(list->size()>node_idxl);
00747   return (*list)[node_idxl];
00748 }
00749 
00750 bool BulkAdapt::is_node_in_idxl(int node_idxl, int partID)
00751 {
00752   IDXL_List *list = meshPtr->node.shared.getIdxlListN(partID);
00753   return(list->size()>node_idxl);
00754 }
00755 
00756 void BulkAdapt::get_elemsToLock(adaptAdj startElem, adaptAdj **elemsToLock, int edgeID, int *count)
00757 {
00758   CkVec<adaptAdj>* nbrElems;
00759   // find the elements adjacent to startElem along the edge edgeID
00760   //BULK_DEBUG(CkPrintf("[%d] BulkAdapt::get_elemsToLock: calling getEdgeAdaptAdj on elem %d\n",partitionID,startElem.localID));
00761 
00762   nbrElems = getEdgeAdaptAdj(meshID, startElem.localID, startElem.elemType, 
00763                  edgeID);
00764   // extract adjacencies from CkVec into array needed by the locking code 
00765   (*count) = nbrElems->size();
00766   (*elemsToLock) = (adaptAdj *)malloc((*count + 1) * sizeof(adaptAdj));
00767 
00768   for (int i=0; i<*count; i++) {
00769     (*elemsToLock)[i] = (*nbrElems)[i];
00770   }
00771   // add the start element
00772   (*elemsToLock)[*count] = startElem;
00773   (*count)++;
00774 
00775   /*
00776   printf("Elems to lock: ");
00777   for (int i=0; i<*count; i++) {
00778     printf("(%d, %d, %d) ", (*elemsToLock)[i].partID, (*elemsToLock)[i].localID,
00779        (*elemsToLock)[i].elemType);
00780   }
00781   printf("\n");
00782   */
00783 }
00784 
00785 
00786 
00787 void midpoint(double *n1, double *n2, int dim, double *result) {
00788   for(int i=0;i<dim;i++){
00789     result[i] = (n1[i]+n2[i])/2.0;
00790   }
00791 }
00792 
00793 int getRelNode(int nodeIdx, int *conn, int nodesPerElem) {
00794   for(int i=0;i<nodesPerElem;i++){
00795     if(conn[i] == nodeIdx){
00796       return i;
00797     }
00798   }
00799   CkAbort("Could not find node in given connectivity");
00800   return -1;
00801 }
00802 
00803 
00804 void getRelNodes(int edgeID, int nodesPerElem, int *r1, int *r2)
00805 {
00806   if (nodesPerElem == 3) {
00807     (*r1) = edgeID;
00808     (*r2) = (edgeID + 1) % 3;
00809   }
00810   else if (nodesPerElem == 4) {
00811     if (edgeID < 3) {
00812       (*r1) = 0;
00813       (*r2) = edgeID+1;
00814     }
00815     else if (edgeID < 5) {
00816       (*r1) = 1;
00817       (*r2) = edgeID-1;
00818     }
00819     else {
00820       (*r1) = 2; 
00821       (*r2) = 3;
00822     }
00823   }
00824 }
00825 
00826 int getEdgeID(int node1, int node2, int nodesPerElem, int dim) {
00827   switch(dim){
00828   case 2: {
00829     switch(nodesPerElem){
00830     case 3:
00831       {
00832     static int edgeFromNode[3][3]={{-1,0,2},{0,-1,1},{2,1,-1}};
00833     return edgeFromNode[node1][node2];
00834       }
00835       break;
00836     default:
00837       CkAbort("This shape is not yet implemented");
00838     };
00839     break;
00840   }
00841   case 3: {
00842     switch(nodesPerElem){
00843     case 4:
00844       {
00845     static int edgeFromNode[4][4]={{-1,0,1,2},{0,-1,3,4},{1,3,-1,5},{2,4,5,-1}};
00846     return edgeFromNode[node1][node2];
00847       }
00848       break;
00849     default:
00850       CkAbort("This shape is not yet implemented");
00851     };
00852     break;
00853   }
00854   default:
00855     CkAbort("This dimension not yet implemented");
00856   };
00857   return -1;
00858 }
00859 
00860 
00861 int getFaceID(int node1, int node2, int node3, int nodesPerElem)
00862 {
00863   switch(nodesPerElem){
00864   case 4: {
00865     static int faceFromNode[4][4][4]=
00866       {{{-1,-1,-1,-1},  {-1,-1, 0, 1},  {-1, 0,-1, 2},  {-1, 1, 2,-1}}, 
00867        {{-1,-1, 0, 1},  {-1,-1,-1,-1},   {0,-1,-1, 3},   {1,-1, 3,-1}}, 
00868        {{-1, 0,-1, 2},   {0,-1,-1, 3},  {-1,-1,-1,-1},   {2, 3,-1,-1}},
00869        {{-1, 1, 2,-1},   {1,-1, 3,-1},   {2, 3,-1,-1},  {-1,-1,-1,-1}}};
00870     return faceFromNode[node1][node2][node3];
00871     break;
00872   }
00873   default:
00874     CkAbort("This shape is not yet implemented");
00875   };
00876   return -1;
00877 }
00878 
00879 void fillNodes(int *relNode, int *nodeIDs, int *conn)
00880 { // ASSERT: this is only for tets.  Assumes positions 0 and 1 are filled.
00881   if ((relNode[0] != 0) && (relNode[1] != 0))
00882     relNode[2] = 0;
00883   else if ((relNode[0] != 1) && (relNode[1] != 1))
00884     relNode[2] = 1;
00885   else relNode[2] = 2;
00886   relNode[3] = 6 - relNode[0] - relNode[1] - relNode[2];
00887   nodeIDs[2] = conn[relNode[2]];
00888   nodeIDs[3] = conn[relNode[3]];
00889 }
00890 
00891 void fillNodes(int *relNode, int n1, int n2, int *conn)
00892 { // Given real node IDs n1 and n2 representing, respectively, the
00893   // unchanged and the changed ends of the edge to bisect, and given
00894   // the conn of a tet with edge (n1, n2), fill out a set of relative
00895   // node numberings that put n1 and n2 in the first two relNode
00896   // positions, respectively, and the remaining nodes in the remaining
00897   // positions. i.e.:
00898   // conn[relNode[0]] = n1, conn[relNode[1] = n2, and the rest map to 
00899   // either of the remaining nodes.
00900   for (int i=0; i<4; i++) {
00901     if (conn[i] == n1)
00902       relNode[0] = i;
00903     else if (conn[i] == n2)
00904       relNode[1] = i;
00905   }
00906   if ((relNode[0] != 0) && (relNode[1] != 0))
00907     relNode[2] = 0;
00908   else if ((relNode[0] != 1) && (relNode[1] != 1))
00909     relNode[2] = 1;
00910   else relNode[2] = 2;
00911   relNode[3] = 6 - relNode[0] - relNode[1] - relNode[2];
00912 }
00913 
00915 adaptAdj *BulkAdapt::local_split_3D(adaptAdj elem, int n1, int n2, int n5)
00916 {
00917   CkPrintf("[%d] BEGIN local_split_3D\n", partitionID);
00918   int nElems = meshPtr->lookup(FEM_ELEM+elem.elemType,"BulkAdapt::local_split_3D")->size();  // Get number of elements
00919   int nNodes=meshPtr->lookup(FEM_NODE,"BulkAdapt::local_split_3D")->size();    // Get number of nodes
00920   CkAssert(elem.localID < nElems);
00921   CkAssert((n1 < nNodes) && (n1 >= 0));
00922   CkAssert((n2 < nNodes) && (n2 >= 0));
00923   CkAssert((n5 < nNodes) && (n5 >= 0));
00924 
00925   // get elem's conn
00926   FEM_Elem &elems = meshPtr->elem[elem.elemType];
00927   int *conn = elems.connFor(elem.localID);
00928   // make splitElem's conn and init with elem's conn
00929   int splitConn[4];
00930   memcpy(&splitConn[0], conn, 4*sizeof(int));
00931   // find relative node numbers for n1 and n2
00932   int rel_n1 = -1, rel_n2 = -1;
00933   for (int i=0; i<4; i++) {
00934     if (conn[i] == n1) rel_n1 = i;
00935     if (conn[i] == n2) rel_n2 = i;
00936   }
00937   CkAssert((rel_n1 != -1) && (rel_n2 != -1));
00938   // adjust elem's conn to reflect the split
00939   CkPrintf("[%d] elem %d conn was %d,%d,%d,%d\n", partitionID, elem.localID, conn[0], conn[1], conn[2], conn[3]);
00940   CkPrintf("[%d] old edge is %d,%d, inserting node %d in between\n", partitionID, n1, n2, n5);
00941   conn[rel_n2] = n5;
00942   CkPrintf("[%d] modifying elem %d to be %d,%d,%d,%d\n", partitionID, elem.localID, conn[0], conn[1], conn[2], conn[3]);
00943   // adjust splitElem's conn to reflect the split
00944   splitConn[rel_n1] = n5;
00945   // add splitElem
00946   int splitElemID = add_element(elem.elemType, 4, &splitConn[0], elems.getMeshSizing(elem.localID));
00947   CkPrintf("[%d] Adding elem %d with conn %d,%d,%d,%d\n", partitionID, splitElemID, splitConn[0], splitConn[1], splitConn[2], splitConn[3]); 
00948   adaptAdj *splitElem = new adaptAdj(partitionID, splitElemID, elem.elemType);
00949   // call updates here
00950   update_local_face_adj(elem, *splitElem, n1, n2, n5);
00951   update_local_edge_adj(elem, *splitElem, n1, n2, n5);
00952   CkPrintf("[%d] END local_split_3D\n", partitionID);
00953   return splitElem;
00954 }
00955 
00956 void BulkAdapt::local_update_asterisk_3D(int i, adaptAdj elem, int numElemPairs,
00957                      adaptAdj *elemPairs, 
00958                      int n1, int n2, int n5)
00959 {
00960   FEM_Elem &elems = meshPtr->elem[elem.elemType];
00961   int *conn = elems.connFor(elem.localID); 
00962   int n3=-1, n4=-1;
00963   for (int i=0; i<4; i++) {
00964     if ((conn[i] != n1) && (conn[i] != n5) && (n3 == -1))
00965       n3 = conn[i]; 
00966     else if ((conn[i] != n1) && (conn[i] != n5))
00967       n4 = conn[i];
00968   }
00969   // derive relative nodes
00970   int relNode[4];
00971   relNode[0] = getRelNode(n1, conn, 4);
00972   relNode[1] = getRelNode(n5, conn, 4);
00973   relNode[2] = getRelNode(n3, conn, 4);
00974   relNode[3] = getRelNode(n4, conn, 4);
00975   int face[4]; // face[0] and [1] are split; others are not
00976   face[0] = (relNode[0] + relNode[1] + relNode[2]) - 3;
00977   face[1] = (relNode[0] + relNode[1] + relNode[3]) - 3;
00978   face[2] = (relNode[0] + relNode[3] + relNode[2]) - 3;
00979   face[3] = (relNode[1] + relNode[3] + relNode[2]) - 3;
00980   adaptAdj neighbors[4]; // elem's neighbors
00981   neighbors[0] = *getFaceAdaptAdj(meshPtr, elem.localID, elem.elemType, face[0]);
00982   neighbors[1] = *getFaceAdaptAdj(meshPtr, elem.localID, elem.elemType, face[1]);
00983   neighbors[2] = *getFaceAdaptAdj(meshPtr, elem.localID, elem.elemType, face[2]);
00984   neighbors[3] = *getFaceAdaptAdj(meshPtr, elem.localID, elem.elemType, face[3]);
00985 
00986   adaptAdj splitElem = neighbors[3], nbr1 = neighbors[0], nbr2 = neighbors[1];
00987   adaptAdj nbr1split(-1,-1,0), nbr2split(-1,-1,0);
00988   int found=0;
00989   if (nbr1.localID != -1) found++;
00990   if (nbr2.localID != -1) found++;
00991   for (int i=0; i<numElemPairs; i++) {
00992     if (elemPairs[2*i] == nbr1) {
00993       nbr1split = elemPairs[2*i+1];
00994       found--;
00995     }
00996     else if (elemPairs[2*i] == nbr2) {
00997       nbr2split = elemPairs[2*i+1];
00998       found--;
00999     }
01000     if (found == 0) break;
01001   }
01002   int *splitConn = elems.connFor(splitElem.localID); 
01003   int splitRelNode[4];
01004   splitRelNode[0] = getRelNode(n5, splitConn, 4);
01005   splitRelNode[1] = getRelNode(n2, splitConn, 4);
01006   splitRelNode[2] = getRelNode(n3, splitConn, 4);
01007   splitRelNode[3] = getRelNode(n4, splitConn, 4);
01008   int splitFace[4]; // splitFace[0] and [1] are split; others are not
01009   splitFace[0] = (splitRelNode[0] + splitRelNode[1] + splitRelNode[2]) - 3;
01010   splitFace[1] = (splitRelNode[0] + splitRelNode[1] + splitRelNode[3]) - 3;
01011   splitFace[2] = (splitRelNode[0] + splitRelNode[3] + splitRelNode[2]) - 3;
01012   splitFace[3] = (splitRelNode[1] + splitRelNode[3] + splitRelNode[2]) - 3;
01013 
01014   // Update faces first
01015   replaceAdaptAdj(meshPtr, splitElem, nbr1, nbr1split);
01016   replaceAdaptAdj(meshPtr, splitElem, nbr2, nbr2split);
01017   // Now, update edges
01018   // Here are the edges to update for elem:
01019   int n3_n5, n4_n5;
01020   // Here are the edges to update for splitElem:
01021   int n5_n2, n5_n3, n5_n4, n2_n3, n2_n4;
01022   // set the edge IDs
01023   n3_n5 = getEdgeID(relNode[2], relNode[1], 4, 3);
01024   n4_n5 = getEdgeID(relNode[3], relNode[1], 4, 3);
01025   n5_n2 = getEdgeID(splitRelNode[0], splitRelNode[1], 4, 3);
01026   n5_n3 = getEdgeID(splitRelNode[0], splitRelNode[2], 4, 3);
01027   n5_n4 = getEdgeID(splitRelNode[0], splitRelNode[3], 4, 3);
01028   n2_n3 = getEdgeID(splitRelNode[1], splitRelNode[2], 4, 3);
01029   n2_n4 = getEdgeID(splitRelNode[1], splitRelNode[3], 4, 3);
01030 
01031   clearEdgeAdjacency(meshPtr, elem.localID, elem.elemType, n3_n5);  
01032   clearEdgeAdjacency(meshPtr, elem.localID, elem.elemType, n4_n5);  
01033   clearEdgeAdjacency(meshPtr, splitElem.localID, splitElem.elemType, n5_n2);  
01034   clearEdgeAdjacency(meshPtr, splitElem.localID, splitElem.elemType, n5_n3);  
01035   clearEdgeAdjacency(meshPtr, splitElem.localID, splitElem.elemType, n5_n4);  
01036 
01037   addEdgeAdjacency(meshPtr, elem.localID, elem.elemType, n3_n5, splitElem);
01038   if (nbr1.localID != -1) {
01039     addEdgeAdjacency(meshPtr, elem.localID, elem.elemType, n3_n5, nbr1);
01040     addEdgeAdjacency(meshPtr, elem.localID, elem.elemType, n3_n5, nbr1split);
01041   }
01042 
01043   addEdgeAdjacency(meshPtr, elem.localID, elem.elemType, n4_n5, splitElem);
01044   if (nbr2.localID != -1) {
01045     addEdgeAdjacency(meshPtr, elem.localID, elem.elemType, n4_n5, nbr2);
01046     addEdgeAdjacency(meshPtr, elem.localID, elem.elemType, n4_n5, nbr2split);
01047   }
01048 
01049   addEdgeAdjacency(meshPtr, splitElem.localID, splitElem.elemType, n5_n3, elem);
01050   if (nbr1.localID != -1) {
01051     addEdgeAdjacency(meshPtr, splitElem.localID, splitElem.elemType, n5_n3, nbr1);
01052     addEdgeAdjacency(meshPtr, splitElem.localID, splitElem.elemType, n5_n3, nbr1split);
01053   }
01054 
01055   addEdgeAdjacency(meshPtr, splitElem.localID, splitElem.elemType, n5_n4, elem);
01056   if (nbr2.localID != -1) {
01057     addEdgeAdjacency(meshPtr, splitElem.localID, splitElem.elemType, n5_n4, nbr2);
01058     addEdgeAdjacency(meshPtr, splitElem.localID, splitElem.elemType, n5_n4, nbr2split);
01059   }
01060 
01061   for (int i=0; i<numElemPairs; i++) {
01062     if (splitElem != elemPairs[2*i+1]) {
01063       addEdgeAdjacency(meshPtr, splitElem.localID, splitElem.elemType, n5_n2,
01064                elemPairs[2*i+1]);
01065     }
01066   }
01067 
01068   if(nbr1.localID != -1)
01069     replaceAdaptAdjOnEdge(meshPtr, splitElem, nbr1, nbr1split, n2_n3);
01070   if(nbr2.localID != -1)
01071     replaceAdaptAdjOnEdge(meshPtr, splitElem, nbr2, nbr2split, n2_n4);
01072 }
01073 
01074 
01076 void BulkAdapt::update_local_face_adj(adaptAdj elem, adaptAdj splitElem, 
01077                       int n1, int n2, int n5)
01078 {
01079   // init splitElem's face and edge adjacencies to elem's
01080   CmiMemoryCheck();
01081   copyAdaptAdj(meshPtr, &elem, &splitElem);
01082   // get n3 and n4 from conn
01083   FEM_Elem &elems = meshPtr->elem[elem.elemType];
01084   int *conn = elems.connFor(elem.localID); 
01085   int n3=-1, n4=-1;
01086   for (int i=0; i<4; i++) {
01087     if ((conn[i] != n1) && (conn[i] != n5) && (n3 == -1))
01088       n3 = conn[i]; 
01089     else if ((conn[i] != n1) && (conn[i] != n5))
01090       n4 = conn[i];
01091   }
01092   // derive relative nodes
01093   int relNode[4];
01094   relNode[0] = getRelNode(n1, conn, 4);
01095   relNode[1] = getRelNode(n5, conn, 4);
01096   relNode[2] = getRelNode(n3, conn, 4);
01097   relNode[3] = getRelNode(n4, conn, 4);
01098   int face[4]; // face[0] and [1] are split; others are not
01099   face[0] = (relNode[0] + relNode[1] + relNode[2]) - 3;
01100   face[1] = (relNode[0] + relNode[1] + relNode[3]) - 3;
01101   face[2] = (relNode[0] + relNode[3] + relNode[2]) - 3;
01102   face[3] = (relNode[1] + relNode[3] + relNode[2]) - 3;
01103   adaptAdj neighbors[4]; // elem's neighbors
01104   neighbors[0] = *getFaceAdaptAdj(meshPtr,elem.localID, elem.elemType, face[0]);
01105   neighbors[1] = *getFaceAdaptAdj(meshPtr,elem.localID, elem.elemType, face[1]);
01106   neighbors[2] = *getFaceAdaptAdj(meshPtr,elem.localID, elem.elemType, face[2]);
01107   neighbors[3] = *getFaceAdaptAdj(meshPtr,elem.localID, elem.elemType, face[3]);
01108   // elem's non-split faces have neighbors[2] and [3]
01109   // update elem's new nbr from neighbors[3] to splitElem
01110   setAdaptAdj(meshPtr, elem, face[3], splitElem);
01111   //replaceAdaptAdj(meshPtr, elem, neighbors[3], splitElem);
01112   // update splitElem's neighbor from neighbors[2] to elem
01113   setAdaptAdj(meshPtr, splitElem, face[2], elem);
01114   //replaceAdaptAdj(meshPtr, splitElem, neighbors[2], elem);
01115   // update elem's nbr's back-adjacency from elem to splitElem
01116   if (neighbors[3].partID == elem.partID) {
01117     replaceAdaptAdj(meshPtr, neighbors[3], elem, splitElem);
01118   }
01119   else if (neighbors[3].partID != -1) { // elem's nbr exists and is remote
01120     shadowProxy[neighbors[3].partID].remote_adaptAdj_replace(neighbors[3], elem, splitElem); 
01121   }
01122 }
01123 
01125 void BulkAdapt::update_local_edge_adj(adaptAdj elem, adaptAdj splitElem, 
01126                       int n1, int n2, int n5)
01127 {
01128   copyEdgeAdaptAdj(meshPtr, &elem, &splitElem);
01129   // first, extract the modifiable data structures in question
01130   CkVec<adaptAdj> *elemEdgeAdj[6];
01131   CkVec<adaptAdj> *splitElemEdgeAdj[6];
01132   for (int i=0; i<6; i++) {
01133     elemEdgeAdj[i] = getEdgeAdaptAdj(meshPtr, elem.localID, elem.elemType, i);
01134     splitElemEdgeAdj[i] = getEdgeAdaptAdj(meshPtr, splitElem.localID, splitElem.elemType, i);
01135   }
01136   // the edges modified in this operation are the old edges of elem,
01137   // now on splitElem: (n2,n3) and (n2,n4), the old edge (n3,n4) which
01138   // was on elem, but is now on both elem and splitElem, and the new
01139   // edges that get added incident on n5: (n3,n5) and (n4,n5).  The
01140   // new edge (n2,n5) must be updated in a different context.
01141   
01142   // Here's how each of these needs to be modified:
01143   // (n2,n3): originally on elem, for splitElem, it will be identical, but we
01144   // need to replace elem with splitElem on all the other elements along that 
01145   // edge.
01146   // (n2,n4): identical to (n2,n3)
01147   // (n3,n4): elem needs to take it's original list and add splitElem to it; 
01148   // splitElem takes elem's original list and adds element to it; then
01149   // all the elements in that original list need to add splitElem.
01150   // (n3,n5): start with an empty list on both elem and splitElem.
01151   // elem adds splitElem, splitElem adds elem.  Other elements
01152   // surrounding the split edge (n1,n2) must be added outside of this
01153   // operation
01154   // (n4,n5): identical to (n3,n5)
01155 
01156   // Start with splitElem: get relNodes and edgeIDs
01157   FEM_Elem &elems = meshPtr->elem[splitElem.elemType];
01158   int *conn = elems.connFor(splitElem.localID); 
01159   int n3=-1, n4=-1;
01160   for (int i=0; i<4; i++) {
01161     if ((conn[i] != n5) && (conn[i] != n2) && (n3 == -1))
01162       n3 = conn[i]; 
01163     else if ((conn[i] != n5) && (conn[i] != n2))
01164       n4 = conn[i];
01165   }
01166   // derive relative nodes
01167   int relNode[4];
01168   relNode[0] = getRelNode(n5, conn, 4);
01169   relNode[1] = getRelNode(n2, conn, 4);
01170   relNode[2] = getRelNode(n3, conn, 4);
01171   relNode[3] = getRelNode(n4, conn, 4);
01172   
01173   // edgeIDs on splitElem
01174   int n2_n3 = getEdgeID(relNode[1], relNode[2], 4, 3);
01175   int n2_n4 = getEdgeID(relNode[1], relNode[3], 4, 3);
01176   int n3_n4 = getEdgeID(relNode[2], relNode[3], 4, 3);
01177   int n3_n5 = getEdgeID(relNode[2], relNode[0], 4, 3);
01178   int n4_n5 = getEdgeID(relNode[3], relNode[0], 4, 3);
01179 
01180   // Get similar info for elem
01181   int *elemConn = elems.connFor(elem.localID); 
01182   // derive relative nodes
01183   int elem_relNode[4];
01184   elem_relNode[0] = getRelNode(n1, elemConn, 4);
01185   elem_relNode[1] = getRelNode(n5, elemConn, 4);
01186   elem_relNode[2] = getRelNode(n3, elemConn, 4);
01187   elem_relNode[3] = getRelNode(n4, elemConn, 4);
01188 
01189   // edgeIDs on elem
01190   int elem_n3_n4 = getEdgeID(relNode[2], relNode[3], 4, 3);
01191   int elem_n3_n5 = getEdgeID(relNode[2], relNode[1], 4, 3);
01192   int elem_n4_n5 = getEdgeID(relNode[3], relNode[1], 4, 3);
01193 
01194   // split face IDs on elem
01195   int face[4]; // face[0] and [1] are split; others are not
01196   face[0] = (elem_relNode[0] + elem_relNode[1] + elem_relNode[2]) - 3;
01197   face[1] = (elem_relNode[0] + elem_relNode[1] + elem_relNode[3]) - 3;
01198   // split neighbors on elem
01199   adaptAdj elem_neighbors[4]; // elem's neighbors
01200   elem_neighbors[0] = *getFaceAdaptAdj(meshPtr,elem.localID, elem.elemType, face[0]);
01201   elem_neighbors[1] = *getFaceAdaptAdj(meshPtr,elem.localID, elem.elemType, face[1]);
01202 
01203   // Now we're ready to go to town
01204   // n2_n3
01205   for (int i=0; i<elemEdgeAdj[elem_n3_n5]->size(); i++) {
01206     adaptAdj adj = (*elemEdgeAdj[elem_n3_n5])[i];
01207     if ((adj != elem_neighbors[0]) && (adj != elem_neighbors[1])) {
01208       if (adj.partID == splitElem.partID) { // do the local replace on adj
01209     int *adjConn = elems.connFor(adj.localID); 
01210     // derive relative nodes
01211     int r2, r3;
01212     r2 = getRelNode(n2, adjConn, 4);
01213     r3 = getRelNode(n3, adjConn, 4);
01214     // edgeID on adj
01215     int edgeID = getEdgeID(r2, r3, 4, 3);
01216     replaceAdaptAdjOnEdge(meshPtr, adj, elem, splitElem, edgeID);
01217       }
01218       else if (adj.partID != -1) { // call remote replacement
01219     int n2_idxl = get_idxl_for_node(n2, adj.partID);
01220     int n3_idxl = get_idxl_for_node(n3, adj.partID);
01221     //printf("[%d] A:calling remote edgeAdj replace on %don%d, to replace %don%d with %don%d, with edge at %d,%d thinking elem's forbidden face neighbors are %don%d and %don%d\n", partitionID, adj.localID, adj.partID, elem.localID, elem.partID, splitElem.localID, splitElem.partID, n2, n3, elem_neighbors[0].localID, elem_neighbors[0].partID, elem_neighbors[1].localID, elem_neighbors[1].partID);
01222     shadowProxy[adj.partID].remote_edgeAdj_replace(partitionID, adj, elem, 
01223                                splitElem, n2_idxl, n3_idxl);
01224       }
01225     }
01226   }
01227   // wow!  all that was just for the (n2,n3) edge!  This sucks!
01228   // n2_n4 -- easy, just like the previous one
01229   for (int i=0; i<elemEdgeAdj[elem_n4_n5]->size(); i++) {
01230     adaptAdj adj = (*elemEdgeAdj[elem_n4_n5])[i];
01231     if ((adj != elem_neighbors[0]) && (adj != elem_neighbors[1])) {
01232       if (adj.partID == splitElem.partID) { // do the local replace on adj
01233     int *adjConn = elems.connFor(adj.localID); 
01234     // derive relative nodes
01235     int r2, r4;
01236     r2 = getRelNode(n2, adjConn, 4);
01237     r4 = getRelNode(n4, adjConn, 4);
01238     // edgeID on adj
01239     int edgeID = getEdgeID(r2, r4, 4, 3);
01240     replaceAdaptAdjOnEdge(meshPtr, adj, elem, splitElem, edgeID);
01241       }
01242       else if (adj.partID != -1) { // call remote replacement
01243     int n2_idxl = get_idxl_for_node(n2, adj.partID);
01244     int n4_idxl = get_idxl_for_node(n4, adj.partID);
01245     //printf("[%d] B:calling remote edgeAdj replace on %don%d, to replace %don%d with %don%d, with edge at %d,%d, thinking elem's forbidden face neighbors are %don%d and %don%d\n", partitionID, adj.localID, adj.partID, elem.localID, elem.partID, splitElem.localID, splitElem.partID, n2, n4, elem_neighbors[0].localID, elem_neighbors[0].partID, elem_neighbors[1].localID, elem_neighbors[1].partID);
01246     shadowProxy[adj.partID].remote_edgeAdj_replace(partitionID, adj, elem, 
01247                                splitElem, n2_idxl, n4_idxl);
01248       }
01249     }
01250   }
01251   // n3_n4: elem needs to take it's original list and add splitElem to it; 
01252   // splitElem takes elem's original list and adds element to it; then
01253   // all the elements in that original list need to add splitElem.
01254   // This is *similar* to the previous case, with the exceptions that
01255   // we are *adding* instead of replacing on the other elems on the
01256   // edge, and the addition of elem and splitElem to each other's edgeAdj.
01257   for (int i=0; i<elemEdgeAdj[elem_n3_n4]->size(); i++) {
01258     adaptAdj adj = (*elemEdgeAdj[elem_n3_n4])[i];
01259     if (adj.partID == splitElem.partID) { // do the local replace on adj
01260       int *adjConn = elems.connFor(adj.localID); 
01261       // derive relative nodes
01262       int r3, r4;
01263       r3 = getRelNode(n3, adjConn, 4);
01264       r4 = getRelNode(n4, adjConn, 4);
01265       // edgeID on adj
01266       int edgeID = getEdgeID(r3, r4, 4, 3);
01267       addToAdaptAdj(meshPtr, adj, edgeID, splitElem);
01268     }
01269     else if (adj.partID != -1) { // call remote replacement
01270       int n3_idxl = get_idxl_for_node(n3, adj.partID);
01271       int n4_idxl = get_idxl_for_node(n4, adj.partID);
01272       shadowProxy[adj.partID].remote_edgeAdj_add(partitionID, adj, splitElem, 
01273                          n3_idxl, n4_idxl);
01274     }
01275   }
01276   addToAdaptAdj(meshPtr, splitElem, n3_n4, elem);
01277   addToAdaptAdj(meshPtr, elem, elem_n3_n4, splitElem);
01278   // The last two cases are performed elsewhere
01279 }
01280
libs/ck-libs/ParFUM/bulk_adapt_ops.C
