Got rid of double free, though I left a memory leak in its place
[charm.git] / src / libs / ck-libs / ParFUM / adapt_adj.C
1 /* Adaptivity Adjacencies: element-to-element adjacencies for use by
2    adaptivity codes only.  Adaptivity codes should keep them
3    up-to-date for each mesh modification primitive.
4    
5    Created 11 Sept 2006 - Terry L. Wilmarth
6 */
7 #include "ParFUM.h"
8 #include "ParFUM_internals.h"
9 using namespace std;
10
11 const int faceMap2d_tri[3][2] = {{0,1},{1,2},{2,0}};
12 const int faceMap2d_quad[4][2] = {{0,1},{1,2},{2,3},{3,0}};
13 const int faceMap3d_tet[4][3] = {{0,1,2},{0,3,1},{0,2,3},{1,3,2}};
14 const int faceMap3d_hex[6][4] = {{0,1,2,3},{1,5,6,2},{2,6,7,3},{3,7,4,0},{0,4,5,1},{5,4,6,7}};
15
16 const int edgeMap3d_tet[6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};
17 const int edgeMap3d_hex[12][2] = {{0,1},{0,3},{0,4},{1,2},{1,5},{2,3},{2,6},{3,7},{4,5},{4,7},{5,6},{6,7}};
18
19 const int faceMap2d_cohquad[2][2] = {{0,1},{2,3}}; // Cohesive for 2D triangles
20 const int faceMap3d_cohprism[2][3] = {{0,1,2},{3,4,5}}; // Cohesive for 3D tets
21
22
23 inline void addSharedNodeData(int node, const IDXL_Rec *sharedChunkList,
24             adjNode *adaptAdjTable)
25 {
26     adaptAdjTable[node].numSharedPartitions = sharedChunkList->getShared();
27     adaptAdjTable[node].sharedWithPartition = 
28         new int [sharedChunkList->getShared()];
29     adaptAdjTable[node].sharedWithLocalIdx = 
30         new int [sharedChunkList->getShared()];
31     for(int j=0; j<sharedChunkList->getShared(); j++){
32         int sharedChunk = sharedChunkList->getChk(j);
33         int sharedIdx = sharedChunkList->getIdx(j);
34         adaptAdjTable[node].sharedWithPartition[j] = sharedChunk;
35         adaptAdjTable[node].sharedWithLocalIdx[j] = sharedIdx;
36     }
37 }
38
39
40 inline void addElementNodeSetData(
41         const int elem, 
42         const int *conn, 
43         const int nodesPerElem,
44         const int numFaces, 
45         const int numEdges, 
46         const int faceSize,
47         const int faceMap[MAX_ADJELEMS][MAX_FACE_SIZE],
48         const int edgeMap[MAX_EDGES][2],
49         adjNode *faceTable, 
50         adjNode* edgeTable)
51 {
52     // first add face adjacencies
53     for (int j=0; j<numFaces; j++) { 
54         adjElem *e = new adjElem(faceSize);
55         e->nodeSetID = j;
56         for (int k=0; k<faceSize; k++) { 
57             e->nodeSet[k] = conn[elem * nodesPerElem + faceMap[j][k]];
58         }
59         // Add this element-nodeSet pair to table at min nodeID in the nodeSet
60         e->nodeSet.quickSort();
61         int minNode = e->nodeSet[0];
62         e->elemID = elem;
63         e->next = faceTable[minNode].adjElemList->next;
64         faceTable[minNode].adjElemList->next = e;
65         faceTable[minNode].adjElemCount++;
66     }
67
68     // then add edge adjacencies
69     if (edgeTable) {
70         for (int j=0; j<numEdges; j++) { 
71             adjElem *e = new adjElem(2);
72             for (int k=0; k<2; k++) { 
73                 e->nodeSet[k] = conn[elem * nodesPerElem + edgeMap[j][k]];
74             }
75             // Add this element-nodeSet pair to table at min nodeID in the 
76             // nodeSet
77             e->nodeSet.quickSort();
78             int minNode = e->nodeSet[0];
79             e->elemID = elem;
80             e->next = edgeTable[minNode].adjElemList->next;
81             edgeTable[minNode].adjElemList->next = e;
82             edgeTable[minNode].adjElemCount++;
83         }
84     }
85 }
86
87
88 //Look for an adjElem (rover->next) whose nodeSet matches that
89 //specified in searchForNodeSet in the link list following adjStart
90 //and return rover such that rover->next matches with the
91 //searchForNodeSet. It also checks that the elemID of the element
92 //being searched does not match with that of searchForElemID.  *found
93 //is set to 1 if match is found .. else to 0
94 inline adjElem *searchAdjElemInList(adjElem *adjStart, int *searchForNodeSet,
95             int nodeSetSize, int searchForElemID, int *found)
96 {
97     adjElem *rover = adjStart; // compare rover->next with adjStart
98     *found = 0;
99
100     while (rover->next != NULL) {
101         if (rover->next->elemID != searchForElemID) {
102             *found = 1; // found an element that is not myself, 
103                         // possibly a match
104
105             if (nodeSetSize != rover->next->nodeSet.size()) {
106                 *found = 0; // not a match
107                 continue; 
108             }
109             for (int j=0; j<nodeSetSize; j++) {
110                 if (rover->next->nodeSet[j] != searchForNodeSet[j]) {
111                     *found = 0; // not a match
112                     break;
113                 }
114             }
115         }
116         if (*found) {
117             break; // We have found a nodeSet that matches adjStart
118         } else {
119             rover = rover->next; // Keep looking for matching nodeSet
120         }
121     }
122     return rover;
123 }
124
125
126 /** Create Adaptivity Adjacencies for elemType; dimension inferred. */
127 void CreateAdaptAdjacencies(int meshid, int elemType)
128 {
129     // Need to derive all of these from elemType;
130     int myRank,numChunks;
131     int faceSize;    // number of nodes shared by two adjacent elems
132     int faceMapSize; // number of faces per element
133     int edgeMapSize; // number of edges per element
134
135     MPI_Comm_rank(MPI_COMM_WORLD,&myRank);
136     MPI_Comm_size(MPI_COMM_WORLD,&numChunks);
137     FEM_Mesh* mesh = FEM_chunk::get("CreateAdaptAdjacencies")->lookup(
138             meshid,"CreateAdaptAdjacencies");
139     FEM_Elem* elem = (FEM_Elem *)mesh->lookup(
140             FEM_ELEM+elemType,"CreateAdaptAdjacencies");
141     FEM_Node* node = (FEM_Node *)mesh->lookup(
142             FEM_NODE,"CreateAdaptAdjacencies");
143
144     const int numElems = elem->size();
145     const int numNodes = node->size();
146     const int nodesPerElem = (elem->getConn()).width();
147     CkAssert(node->getCoord() != NULL);
148     const int dim = (node->getCoord())->getWidth();
149     CkAssert(dim == 2|| dim == 3);
150
151     // A nodeSet is a set of nodes that defines a pairing of two
152     // adjacent elements; For example, in 2D triangle meshes, the
153     // nodeSet is the nodes of an edge between two elements.
154     // These adjacency relationships are defined by face and edge maps that
155     // define which set of nodes constitute faces and edges for each element
156     // type.
157     int faceMap[MAX_ADJELEMS][MAX_FACE_SIZE];
158     int edgeMap[MAX_EDGES][2];
159
160     guessElementShape(dim, nodesPerElem, 
161             &faceSize, &faceMapSize, &edgeMapSize, 
162             faceMap, edgeMap);
163     CkAssert(faceMapSize <= MAX_ADJELEMS);
164     CkAssert(edgeMapSize <= MAX_EDGES);
165     CkAssert(faceSize <= MAX_FACE_SIZE);
166
167     // Add the FEM_ADAPT_FACE_ADJ and FEM_ADAPT_EDGE_ADJ attributes
168     FEM_DataAttribute* adaptAdjAttr = (FEM_DataAttribute*) elem->lookup(
169             FEM_ADAPT_FACE_ADJ, "CreateAdaptAdjacencies");      
170     adaptAdjAttr->setWidth(faceMapSize*sizeof(adaptAdj));
171     adaptAdj* adaptFaceAdjacencies = 
172         reinterpret_cast<adaptAdj*>((adaptAdjAttr->getChar()).getData());
173
174     FEM_DataAttribute* adaptAdjEdgeAttr = (FEM_DataAttribute*) elem->lookup(
175             FEM_ADAPT_EDGE_ADJ, "CreateAdaptAdjacencies");      
176     adaptAdjEdgeAttr->setWidth(edgeMapSize*sizeof(CkVec<adaptAdj>*));
177     CkVec<adaptAdj>** adaptEdgeAdjacencies = edgeMapSize == 0 ?
178         NULL :
179         reinterpret_cast<CkVec<adaptAdj>**>(
180                 (adaptAdjEdgeAttr->getChar()).getData());
181
182     // Initialize adaptAdj arrays
183     for (int i=0; i<numElems; i++) {
184         for (int j=0; j<faceMapSize; ++j) {
185             adaptFaceAdjacencies[i*faceMapSize + j].partID = -1;
186             adaptFaceAdjacencies[i*faceMapSize + j].localID = -1;
187         }
188         if (adaptEdgeAdjacencies) {
189             for (int j=0; j<edgeMapSize; ++j) {
190                 adaptEdgeAdjacencies[i*edgeMapSize + j] = new CkVec<adaptAdj>;
191                 assert(adaptEdgeAdjacencies[i*edgeMapSize + j] != NULL);
192             }
193         }
194     }
195
196     // Create an array of size equal to the number of local nodes, each
197     // entry has a pair: (shared partitions, list of elem-nodeSet pairs)
198     // make one array for face adjacencies and one for edge adjacencies
199
200     adjNode *faceTable;
201     adjNode *edgeTable;
202     faceTable = new adjNode[numNodes];
203     edgeTable = adaptEdgeAdjacencies ? new adjNode[numNodes] : NULL;
204
205     // Loop through the shared node list and add shared partition ids
206     for (int i=0; i<numNodes; i++) {
207         if (node->is_valid(i)) {
208             const IDXL_Rec *sharedChunkList = node->shared.getRec(i);
209             if (sharedChunkList != NULL) {
210                 addSharedNodeData(i, sharedChunkList, faceTable);
211                 if (edgeTable)
212                     addSharedNodeData(i, sharedChunkList, edgeTable);
213             } else {
214                 faceTable[i].sharedWithPartition = NULL;
215                 faceTable[i].sharedWithLocalIdx= NULL;
216                 if (edgeTable) {
217                     edgeTable[i].sharedWithPartition = NULL;
218                     edgeTable[i].sharedWithLocalIdx= NULL;
219                 }
220             }
221         } else {
222             faceTable[i].sharedWithPartition = NULL;
223             faceTable[i].sharedWithLocalIdx= NULL;
224             if (edgeTable) {
225                 edgeTable[i].sharedWithPartition = NULL;
226                 edgeTable[i].sharedWithLocalIdx= NULL;
227             }
228         }
229     }
230
231     // add local node adjacency info via connectivity
232     const int *conn = (elem->getConn()).getData();
233     for (int i=0; i<numElems; i++) {
234         if(elem->is_valid(i)){
235             addElementNodeSetData(i, conn, nodesPerElem, faceMapSize, 
236                     edgeMapSize, faceSize, faceMap, edgeMap, faceTable, 
237                     edgeTable);
238         }
239     }
240
241     fillLocalAdaptAdjacencies(
242             numNodes,
243             node,
244             faceTable,
245             edgeTable,
246             faceMapSize,
247             edgeMapSize,
248             adaptFaceAdjacencies,
249             adaptEdgeAdjacencies,
250             myRank,
251             elemType);
252
253     // Now all elements' local adjacencies are set; remainder in table are 
254     // nodeSets shared with other chunks or nodeSets on domain boundary
255     // We handle face nodeSets first, then do the edges.
256
257     MPI_Barrier(MPI_COMM_WORLD);
258     MSA1DREQLIST *requestTable;
259     if (myRank == 0) {
260         requestTable = new MSA1DREQLIST(numChunks,numChunks);
261     } else {
262         requestTable = new MSA1DREQLIST;
263     }
264     MPI_Bcast_pup(*requestTable,0,MPI_COMM_WORLD);
265     requestTable->enroll(numChunks);
266     requestTable->sync();
267
268     makeAdjacencyRequests(
269             numNodes,
270             node,
271             faceTable,
272             requestTable,
273             faceSize,
274             myRank,
275             elemType);
276
277     requestTable->sync();
278     printf("[%d] All face requests made \n",myRank);
279
280     MSA1DREPLYLIST *replyTable;
281     if (myRank == 0) {
282         replyTable = new MSA1DREPLYLIST(numChunks,numChunks);
283     } else {
284         replyTable = new MSA1DREPLYLIST;
285     }
286     MPI_Bcast_pup(*replyTable,0,MPI_COMM_WORLD);
287     replyTable->enroll(numChunks);
288     replyTable->sync();
289
290     replyAdjacencyRequests(
291             requestTable,
292             replyTable,
293             node,
294             faceTable,
295             adaptFaceAdjacencies,
296             adaptEdgeAdjacencies,
297             faceSize,
298             faceMapSize,
299             myRank,
300             elemType,
301             false);
302
303     requestTable->sync();
304     replyTable->sync();
305
306     // Once the replies are back, loop through each reply and update the
307     // adjacencies for each element in the reply
308     CkVec<adjReply> *receivedReplyVec = replyTable->get(myRank).vec;
309     for(int i=0;i< receivedReplyVec->size();i++){
310         adjReply *receivedReply = &(*receivedReplyVec)[i];
311         printf("[%d] Replies received for (%d,%d) (%d,%d,%d)\n",
312                 myRank,receivedReply->requestingElemID,
313                 receivedReply->requestingNodeSetID,
314                 receivedReply->replyingElem.partID,
315                 receivedReply->replyingElem.localID,
316                 receivedReply->replyingElem.elemType);
317         adaptFaceAdjacencies[receivedReply->requestingElemID*faceMapSize + 
318             receivedReply->requestingNodeSetID] = receivedReply->replyingElem;
319     }
320     replyTable->sync();
321
322     if (adaptEdgeAdjacencies != NULL) {
323         //delete requestTable;
324         //delete replyTable;
325
326         // do the same thing for the edges
327         if (myRank == 0) {
328             requestTable = new MSA1DREQLIST(numChunks,numChunks);
329         } else {
330             requestTable = new MSA1DREQLIST;
331         }
332         MPI_Bcast_pup(*requestTable,0,MPI_COMM_WORLD);
333         requestTable->enroll(numChunks);
334         requestTable->sync();
335
336         makeAdjacencyRequests(
337                 numNodes,
338                 node,
339                 edgeTable,
340                 requestTable,
341                 2,
342                 myRank,
343                 elemType);
344
345         requestTable->sync();
346         printf("[%d] All edge requests made \n",myRank);
347
348         MSA1DREPLYLIST *replyTable;
349         if (myRank == 0) {
350             replyTable = new MSA1DREPLYLIST(numChunks,numChunks);
351         } else {
352             replyTable = new MSA1DREPLYLIST;
353         }
354         MPI_Bcast_pup(*replyTable,0,MPI_COMM_WORLD);
355         replyTable->enroll(numChunks);
356         replyTable->sync();
357
358         replyAdjacencyRequests(
359                 requestTable,
360                 replyTable,
361                 node,
362                 edgeTable,
363                 adaptFaceAdjacencies,
364                 adaptEdgeAdjacencies,
365                 2,
366                 edgeMapSize,
367                 myRank,
368                 elemType,
369                 true);
370
371         requestTable->sync();
372         replyTable->sync();
373
374         // Once the replies are back, loop through each reply and update the
375         // adjacencies for each element in the reply
376         CkVec<adjReply> *receivedReplyVec = replyTable->get(myRank).vec;
377         for(int i=0;i< receivedReplyVec->size();i++){
378             adjReply *receivedReply = &(*receivedReplyVec)[i];
379             printf("[%d] Replies received for (%d,%d) (%d,%d,%d)\n",
380                     myRank,receivedReply->requestingElemID,
381                     receivedReply->requestingNodeSetID,
382                     receivedReply->replyingElem.partID,
383                     receivedReply->replyingElem.localID,
384                     receivedReply->replyingElem.elemType);
385             adaptEdgeAdjacencies[receivedReply->requestingElemID*edgeMapSize + 
386                 receivedReply->requestingNodeSetID]->push_back(
387                         receivedReply->replyingElem);
388         }
389         replyTable->sync();
390     }
391
392     // cleanup
393     //delete requestTable;
394     //delete replyTable;
395
396     //for (int i=0; i<numNodes; ++i) {
397     //    delete[] faceTable[i].sharedWithPartition;
398     //    delete[] faceTable[i].sharedWithLocalIdx;
399     //    adjElem* e = faceTable[i].adjElemList->next;
400     //    adjElem* dummy;
401     //    //while (e != NULL) {
402     //    //    dummy = e;
403     //    //    e = e->next;
404     //    //    delete dummy;
405     //    //}
406
407     //    if (edgeTable) {
408     //        delete[] edgeTable[i].sharedWithPartition;
409     //        delete[] edgeTable[i].sharedWithLocalIdx;
410     //        e = edgeTable[i].adjElemList->next;
411     //        //while (e != NULL) {
412     //        //    dummy = e;
413     //        //    e = e->next;
414     //        //    delete dummy;
415     //        //}
416     //    }
417     //}
418     //delete[] faceTable;
419     //delete[] edgeTable;
420 }
421
422
423 void fillLocalAdaptAdjacencies(
424         const int numNodes, 
425         FEM_Node* node, 
426         adjNode* faceTable, 
427         adjNode* edgeTable,
428         const int faceMapSize,
429         const int edgeMapSize,
430         adaptAdj* adaptFaceAdjacencies,
431         CkVec<adaptAdj>** adaptEdgeAdjacencies,
432         const int myRank, 
433         const int elemType)
434 {
435     adjNode* adjTables[2] = { faceTable, edgeTable };
436     adjNode* adaptAdjTable;
437
438     for (int table=0; table<2; ++table) {
439         adaptAdjTable = adjTables[table];
440         for (int i=0; i<numNodes; i++) { 
441             if (!adaptAdjTable) continue;
442             // For each node, match up incident elements
443             // Each adjacency between two elements is represented by two 
444             // adjElems. We try to match those up
445             if(node->is_valid(i) && adaptAdjTable[i].adjElemList != NULL){  
446                 // CkAssert(adaptAdjTable[i].adjElemList != NULL);
447                 adjElem *preTarget = adaptAdjTable[i].adjElemList;
448                 adjElem *target = adaptAdjTable[i].adjElemList->next;
449                 while (target != NULL) { //loop over adjElemList of a node
450                     int found = 0; 
451                     // target represents an adjacency between two elements
452                     // We search for the other adjElem corresponding to that 
453                     // adjancency:
454                     // Look for an entry in adjElemList after target such that 
455                     // the nodeset of that entry and that of target match but 
456                     // they do not belong to the same element. 
457                     adjElem *rover = searchAdjElemInList(
458                             preTarget,
459                             target->nodeSet.getVec(),
460                             target->nodeSet.size(),
461                             target->elemID,
462                             &found); 
463
464                     if (found) { 
465                         // We found a local element adjacent to target->elemID
466                         // Set adjacency of target->elemID corresponding to 
467                         // nodeSet to rover->next->elemID, and vice versa
468                         // Store adjacency info in adaptAdjacency of each one 
469                         // and use nodeSetID to index into adaptAdjacency
470                         
471                         if (table == 0)  { // working on face nodesets
472                             adaptFaceAdjacencies[faceMapSize*target->elemID + 
473                                 target->nodeSetID] = adaptAdj(
474                                         myRank, rover->next->elemID, elemType);
475                             adaptFaceAdjacencies[
476                                 faceMapSize*rover->next->elemID + 
477                                 rover->next->nodeSetID] = adaptAdj(
478                                         myRank, target->elemID, elemType);
479                         } else if (table == 1) { // working on edge nodesets
480                             adaptEdgeAdjacencies[edgeMapSize*target->elemID + 
481                                 target->nodeSetID]->push_back(adaptAdj(
482                                             myRank, 
483                                             rover->next->elemID, 
484                                             elemType));
485                             adaptEdgeAdjacencies[
486                                 edgeMapSize*rover->next->elemID + 
487                                 rover->next->nodeSetID]->push_back(adaptAdj(
488                                             myRank,
489                                             target->elemID,
490                                             elemType));
491                         }
492                         
493                         // Remove both elem-nodeSet pairs from the list
494                         adjElem *tmp = rover->next;
495                         rover->next = rover->next->next;
496                         //delete tmp;
497                         tmp = target;
498                         preTarget->next = target->next;
499                         target = target->next;
500                         //delete tmp;
501                     } else { 
502                         // No match for target was found in adjElemList
503                         // This means that either target is on the domain 
504                         // boundary or it is on a partition boundary and its 
505                         // neighbor is on another VP. Move target to next 
506                         // entry in adjElemList
507                         preTarget = target;
508                         target = target->next;
509                     }
510                 }
511             }
512         }
513     }
514 }
515
516
517 void makeAdjacencyRequests(
518         const int numNodes,
519         FEM_Node *node,
520         adjNode *adaptAdjTable,
521         MSA1DREQLIST *requestTable, 
522         const int nodeSetSize,
523         const int myRank,
524         const int elemType)
525 {
526     for (int i=0; i<numNodes; i++) { 
527         // Examine each node's remaining elements
528         if(node->is_valid(i)){
529             if (adaptAdjTable[i].adjElemList->next!=NULL) {
530                 adjElem *adjStart = adaptAdjTable[i].adjElemList->next;
531                 while (adjStart !=NULL) {
532                     // create and empty set, commonSharedChunks
533                     std::set<int> commonSharedChunks;
534                     for (int j=0; j<nodeSetSize; j++) {
535                         // look up sharedWithPartitions for node: 
536                         //    adaptAdjTable[i]->adjElemList->nodeset[j]
537                         // intersect with commonSharedChunks
538                         int sharedNode = adjStart->nodeSet[j];
539                         adjNode *sharedNodeAdj = &adaptAdjTable[sharedNode];
540                         std::set<int> sharedChunks(
541                                 sharedNodeAdj->sharedWithPartition,
542                                 sharedNodeAdj->sharedWithPartition + 
543                                 sharedNodeAdj->numSharedPartitions);
544                         if(j == 0){
545                             commonSharedChunks = sharedChunks;
546                         }else{
547                             std::set<int> tmpIntersect;
548                             set_intersection(
549                                     commonSharedChunks.begin(), 
550                                     commonSharedChunks.end(), 
551                                     sharedChunks.begin(), 
552                                     sharedChunks.end(),
553                                     inserter(
554                                         tmpIntersect, 
555                                         tmpIntersect.begin()));
556                             commonSharedChunks = tmpIntersect;              
557                         }
558                     }
559                     // At this point commonSharedChunks contains the list of
560                     // chunks with which the element pointed by adjStart might
561                     // be shared
562                     int numCommonSharedChunks = commonSharedChunks.size();
563                     if(numCommonSharedChunks > 0){
564                         //adjStart is possibly shared with these chunks. It is
565                         //shared across the adjStart->nodeSet set of nodes.
566                         adjRequest *adjRequestList = 
567                             new adjRequest[numCommonSharedChunks];
568                         // Translate the nodes in the nodeSet into the index in 
569                         // the idxl list for each chunk in the 
570                         // commonSharedChunks
571                         for(int j=0;j<nodeSetSize;j++){
572                             int sharedNode = adjStart->nodeSet[j];
573                             const IDXL_Rec *recSharedNode = 
574                                 node->shared.getRec(sharedNode);
575                             int countChunk=0;
576                             for(std::set<int>::iterator chunkIterator = 
577                                     commonSharedChunks.begin();
578                                     chunkIterator != commonSharedChunks.end();
579                                     chunkIterator++){
580                                 int chunk = *chunkIterator;
581                                 if(j == 0){
582                                     // if this is the first node we need to 
583                                     // initialize the adjRequestobject
584                                     adjRequestList[countChunk] = 
585                                         adjRequest(adjStart->elemID,
586                                                 myRank,
587                                                 adjStart->nodeSetID,
588                                                 elemType);
589                                 }
590
591                                 // index of sharedNode in the idxl list of 
592                                 // chunk search for this chunk in the list of 
593                                 // chunks shared by this node
594                                 int sharedNodeIdx=-1; 
595                                 for(int k=0;k<recSharedNode->getShared();k++){
596                                     if(recSharedNode->getChk(k) == chunk){
597                                         //found the correct chunk
598                                         sharedNodeIdx = 
599                                             recSharedNode->getIdx(k);
600                                         break;
601                                     }
602                                 }
603                                 CkAssert(sharedNodeIdx != -1);
604                                 //The index of sharedNode in the index list of
605                                 //chunks has been found.  This needs to be 
606                                 //saved in the corresponding translatedNodeSet
607                                 adjRequestList[countChunk].
608                                     translatedNodeSet[j] = sharedNodeIdx;
609                                 countChunk++;
610                             }
611                         }
612
613                         // Now the nodeNumbers for the nodeSets that might be 
614                         // along chunk boundaries have been translated into 
615                         // idxl indices between the two chunks We now need to 
616                         // write the requests into the msa array requestTable 
617                         int countChunk=0;
618                         for(std::set<int>::iterator chunkIterator = 
619                                 commonSharedChunks.begin();
620                                 chunkIterator != commonSharedChunks.end();
621                                 chunkIterator++){
622                             int chunk = *chunkIterator;
623                             printf("[%d] Sending to chunk %d request (%d,%d,%d,%d) \n",myRank,chunk,adjRequestList[countChunk].elemID,adjRequestList[countChunk].chunkID,adjRequestList[countChunk].elemType,adjRequestList[countChunk].nodeSetID);
624                             (*requestTable).accumulate(
625                                     chunk,adjRequestList[countChunk]);
626                             countChunk++;
627                         }
628                         //delete [] adjRequestList;
629                     }
630                     adjStart = adjStart->next;
631                 }
632             }
633         }
634     }
635 }
636
637
638 void replyAdjacencyRequests(
639         MSA1DREQLIST* requestTable,
640         MSA1DREPLYLIST* replyTable,
641         FEM_Node* node,
642         adjNode* adaptAdjTable,
643         adaptAdj* adaptFaceAdjacencies,
644         CkVec<adaptAdj>** adaptEdgeAdjacencies,
645         const int nodeSetSize,
646         const int numAdjElems,
647         const int myRank,
648         const int elemType,
649         bool isEdgeRequest)
650 {
651     // look up request table, put answer in reply table 
652     // receive: for local elemID, the adjacency on this set of shared indices 
653     // is (remote elemID, remote partition ID, remote elem type)
654     // add adjacencies to local table
655     // lots of error checking :)
656
657     //Look at each request that in the requestTable for this chunk
658     //Put the data for the requests in our own table and then create replies
659     CkVec<adjRequest> *receivedRequestVec = requestTable->get(myRank).vec;
660     for (int i=0;i<receivedRequestVec->length();i++) {    
661         adjRequest &receivedRequest = (*receivedRequestVec)[i];
662         const IDXL_List &sharedNodeList = 
663             node->shared.getList(receivedRequest.chunkID);
664         CkVec<int> sharedNodes(nodeSetSize);
665         //Translate all the nodes in the nodeSet of the request
666         for (int j=0;j<nodeSetSize;j++) {
667             int sharedNode = sharedNodeList[
668                 receivedRequest.translatedNodeSet[j]];
669             sharedNodes[j] = sharedNode;
670         }
671         sharedNodes.quickSort();
672         //We need to find the matching nodeset for the nodeset in the request
673         //We look it up in the adaptAdjTable for the minimum node on this chunk
674         adjNode *minNode = &adaptAdjTable[sharedNodes[0]];
675         int found=0;
676         //search for the nodeSet in the list of adjacencies around minNode
677         adjElem *rover =  searchAdjElemInList(
678                 minNode->adjElemList, 
679                 sharedNodes.getVec(), 
680                 nodeSetSize, 
681                 -1, 
682                 &found);
683         printf("[%d] Received request (%d,%d,%d,%d) minNode %d found %d\n",
684                 myRank,receivedRequest.elemID,receivedRequest.chunkID,
685                 receivedRequest.elemType,receivedRequest.nodeSetID,
686                 sharedNodes[0],found);
687         if (found) { 
688             //we have found a matching adjElem for the requested nodeset
689             //we shall set the adjacency correctly in the adjacency Table
690             //for the elemID in the found adjElem. We need to send a reply
691             //to the requesting chunk
692             int matchingElemID;
693             matchingElemID = rover->next->elemID;
694
695             if (isEdgeRequest) {
696                 CkVec<adaptAdj>* matchingAdaptAdj = adaptEdgeAdjacencies[
697                     matchingElemID*numAdjElems + rover->next->nodeSetID];
698                 matchingAdaptAdj->push_back(adaptAdj(
699                             receivedRequest.chunkID,
700                             receivedRequest.elemID, 
701                             receivedRequest.elemType));
702             } else {
703                 adaptAdj *matchingAdaptAdj = &adaptFaceAdjacencies[
704                     matchingElemID*numAdjElems + rover->next->nodeSetID];
705                 CkAssert(matchingAdaptAdj->localID == -1);
706                 matchingAdaptAdj->partID = receivedRequest.chunkID;
707                 matchingAdaptAdj->localID = receivedRequest.elemID;
708                 matchingAdaptAdj->elemType = receivedRequest.elemType;
709             }
710             
711             adjElem *tmp = rover->next;
712             rover->next = tmp->next;
713             //delete tmp;
714
715             //Set requesting data in reply to that in receivedRequest
716             //Put in data from rover->next into the replyElem portion of data
717             adjReply reply;
718             reply.requestingElemID = receivedRequest.elemID;
719             reply.requestingNodeSetID = receivedRequest.nodeSetID;
720             reply.replyingElem.partID = myRank;
721             reply.replyingElem.localID = matchingElemID;
722             reply.replyingElem.elemType = elemType;
723             //Write into the replyTable
724             replyTable->accumulate(receivedRequest.chunkID,reply);
725         } else {
726             //we have no matching nodeset for this request.. hopefully some
727             //other chunk does; we can ignore this request
728         }
729     }
730 }
731
732
733 inline void findNodeSet(
734         int meshid,
735         int elemType,
736         int* faceSize, 
737         int* faceMapSize, 
738         int* edgeMapSize,
739         int nodeSetMap[MAX_ADJELEMS][MAX_NODESET_SIZE],
740         int edgeMap[MAX_EDGES][2])
741 {
742     FEM_Mesh *mesh = FEM_chunk::get("GetAdaptAdj")->
743         lookup(meshid,"GetAdaptAdj");
744     FEM_Elem *elem = (FEM_Elem *)mesh->
745         lookup(FEM_ELEM+elemType,"GetAdaptAdj");
746     FEM_Node *node = (FEM_Node *)mesh->
747         lookup(FEM_NODE,"GetAdaptAdj");
748     const int nodesPer = (elem->getConn()).width();
749     assert(node->getCoord()!= NULL);
750     assert(nodesPer > 0 && nodesPer < MAX_FACE_SIZE);
751     const int dim = (node->getCoord())->getWidth();
752     assert(dim == 2|| dim == 3);
753
754     guessElementShape(dim, nodesPer, 
755             faceSize, faceMapSize, edgeMapSize, 
756             nodeSetMap, edgeMap);
757 }
758
759 void getAndDumpAdaptAdjacencies(
760         const int meshid, 
761         const int numElems, 
762         const int elemType, 
763         const int myRank)
764 {
765   int faceSize;
766   int faceMapSize, edgeMapSize;
767   int nodeSetMap[MAX_ADJELEMS][MAX_NODESET_SIZE]; // not used
768   int edgeMap[MAX_EDGES][2]; // not used
769   
770   findNodeSet(
771           meshid, elemType,
772           &faceSize, &faceMapSize, &edgeMapSize,
773           nodeSetMap, edgeMap);
774   
775   for (int i=0; i<numElems; i++) {
776       printf("[%d] %d  :",myRank,i);
777       for (int j=0; j<faceMapSize; j++) {
778           adaptAdj *entry = getAdaptAdj(meshid, i, elemType, j);
779           printf("(%d,%d,%d)", entry->partID, entry->localID, entry->elemType);
780       }
781       printf("\n");
782   }
783
784   if (edgeMapSize == 0) return;
785   for (int i=0; i<numElems; i++) {
786       printf("[%d] %d  :", myRank, i);
787       for (int j=0; j<edgeMapSize; j++) {
788           CkVec<adaptAdj>* entry = getEdgeAdaptAdj(meshid, i, elemType, j);
789           for (int k=0; k<entry->size(); ++k) {
790               printf("(%d,%d,%d)", (*entry)[k].partID, (*entry)[k].localID, 
791                       (*entry)[k].elemType);
792           }
793       }
794       printf("\n");
795   }
796
797 }
798
799 // Access functions
800 inline adaptAdj* lookupAdaptAdjacencies(
801         const FEM_Mesh* const mesh,
802         const int elemType,
803         int* numAdjacencies)
804 {
805     FEM_Elem *elem = (FEM_Elem*)mesh->lookup(
806             FEM_ELEM+elemType, "lookupAdaptAdjacencies");
807
808     FEM_DataAttribute* adaptAttr = (FEM_DataAttribute*)elem->lookup(
809             FEM_ADAPT_FACE_ADJ, "lookupAdaptAdjacencies");
810     *numAdjacencies = adaptAttr->getWidth()/sizeof(adaptAdj);
811     AllocTable2d<unsigned char> &table = adaptAttr->getChar();
812
813     return (adaptAdj*)(adaptAttr->getChar().getData());
814 }
815
816 inline adaptAdj* lookupAdaptAdjacencies(
817         const int meshid,
818         const int elemType,
819         int* numAdjacencies)
820 {
821     FEM_Mesh* mesh = FEM_chunk::get("lookupAdaptAdjacencies")->lookup(
822             meshid, "lookupAdaptAdjacencies");
823     return lookupAdaptAdjacencies(mesh, elemType, numAdjacencies);
824 }
825
826 inline CkVec<adaptAdj>** lookupEdgeAdaptAdjacencies(
827         const FEM_Mesh* const mesh,
828         const int elemType,
829         int* numAdjacencies)
830 {
831     FEM_Elem *elem = (FEM_Elem*)mesh->lookup(
832             FEM_ELEM+elemType,"lookupAdaptAdjacencies");
833     FEM_DataAttribute* adaptAdjAttr = (FEM_DataAttribute*)elem->lookup(
834             FEM_ADAPT_EDGE_ADJ, "CreateAdaptAdjacencies");      
835     *numAdjacencies = adaptAdjAttr->getWidth()/sizeof(CkVec<adaptAdj>**);
836     CkVec<adaptAdj>** adaptAdjacencies = 
837         reinterpret_cast<CkVec<adaptAdj>**>((adaptAdjAttr->getInt()).getData());
838     return adaptAdjacencies;
839 }
840
841 inline CkVec<adaptAdj>** lookupEdgeAdaptAdjacencies(
842         const int meshID,
843         const int elemType,
844         int* numAdjacencies)
845 {
846     FEM_Mesh *mesh = FEM_chunk::get("lookupAdaptAdjacencies")->
847         lookup(meshID,"lookupAdaptAdjacencies");
848     return lookupEdgeAdaptAdjacencies(mesh, elemType, numAdjacencies);
849 }
850
851 /** Look up elemID in elemType array, access edgeFaceID-th adaptAdj. */
852 adaptAdj* getAdaptAdj(
853         const int meshID, 
854         const int localID,
855         const int elemType, 
856         const int faceID)
857 {
858     int numAdjacencies;
859     adaptAdj* adaptAdjacencies = 
860         lookupAdaptAdjacencies(meshID, elemType, &numAdjacencies);
861     return &(adaptAdjacencies[localID*numAdjacencies+faceID]);
862 }
863
864 adaptAdj* getAdaptAdj(
865         const FEM_Mesh* const meshPtr, 
866         const int localID,
867         const int elemType, 
868         const int faceID)
869 {
870     int numAdjacencies;
871     adaptAdj* adaptAdjacencies = 
872         lookupAdaptAdjacencies(meshPtr, elemType, &numAdjacencies);
873     return &(adaptAdjacencies[localID*numAdjacencies+faceID]);
874 }
875
876 CkVec<adaptAdj>* getEdgeAdaptAdj(
877         const int meshID, 
878         const int localID,
879         const int elemType, 
880         const int edgeID)
881 {
882     int numAdjacencies;
883     CkVec<adaptAdj>** adaptAdjacencies = 
884         lookupEdgeAdaptAdjacencies(meshID, elemType, &numAdjacencies);
885     return adaptAdjacencies[localID*numAdjacencies+edgeID];
886 }
887
888 CkVec<adaptAdj>* getEdgeAdaptAdj(
889         const FEM_Mesh* const meshPtr, 
890         const int localID,
891         const int elemType, 
892         const int edgeID)
893 {
894     int numAdjacencies;
895     CkVec<adaptAdj>** adaptAdjacencies = 
896         lookupEdgeAdaptAdjacencies(meshPtr, elemType, &numAdjacencies);
897     return adaptAdjacencies[localID*numAdjacencies+edgeID];
898 }
899
900 /** Look up elemID in elemType array, calculate edgeFaceID from
901     vertexList (with GetEdgeFace below), and access edgeFaceID-th
902     adaptAdj with GetAdaptAdj above. */
903 adaptAdj* getAdaptAdj(
904         const int meshID, 
905         const int localID,
906         const int elemType, 
907         const int* const vertexList)
908 {
909   int faceID = getElemFace(meshID, elemType, vertexList);
910   return getAdaptAdj(meshID, localID, elemType, faceID);
911 }
912
913 CkVec<adaptAdj>* getEdgeAdaptAdj(
914         const int meshID, 
915         const int localID,
916         const int elemType, 
917         const int* const vertexList)
918 {
919   int edgeID = getElemEdge(meshID, elemType, vertexList);
920   return getEdgeAdaptAdj(meshID, localID, elemType, edgeID);
921 }
922
923 adaptAdj* getFaceAdaptAdj(
924         const int meshID, 
925         const int localID,
926         const int elemType, 
927         const int* const vertexList)
928 {
929   return getAdaptAdj(meshID, localID, elemType, vertexList);
930 }
931
932 adaptAdj* getFaceAdaptAdj(
933         const int meshID, 
934         const int localID,
935         const int elemType, 
936         const int faceID)
937 {
938     return getAdaptAdj(meshID, localID, elemType, faceID);
939 }
940
941 adaptAdj* getFaceAdaptAdj(
942         const FEM_Mesh* const meshPtr, 
943         const int localID,
944         const int elemType, 
945         const int faceID)
946 {
947     return getAdaptAdj(meshPtr, localID, elemType, faceID);
948 }
949
950
951
952 /** Look up elemID in elemType array and determine the edge or face ID
953     specified by the set of vertices in vertexList. */
954 int getElemFace(
955         const int meshID, 
956         const int type, 
957         const int* vertexList)
958 {
959   int faceSize;
960   int faceMapSize;
961   int faceMap[MAX_ADJELEMS][MAX_NODESET_SIZE];
962   int edgeMapSize; // not used
963   int edgeMap[MAX_EDGES][2]; // not used
964   findNodeSet(meshID, type, &faceSize, &faceMapSize, &edgeMapSize,
965          faceMap, edgeMap);
966   
967   // look for vertexList in the face map
968   std::set<int> vertexSet(vertexList, vertexList+faceSize);
969   for (int i=0; i<faceMapSize; i++) {
970     std::set<int> aNodeSet(faceMap[i], faceMap[i]+faceSize);
971     if (vertexSet == aNodeSet) return i;
972   }
973
974   // uh-oh, didn't find it
975   CkAbort("ERROR: GetEdgeFace: vertexList is not a valid face./n");
976 }
977
978 int getElemEdge(
979         const int meshID, 
980         const int type, 
981         const int* vertexList)
982 {
983   int faceSize; // not used
984   int faceMapSize; // not used
985   int faceMap[MAX_ADJELEMS][MAX_NODESET_SIZE]; // not used
986   int edgeMapSize;
987   int edgeMap[MAX_EDGES][2];
988   findNodeSet(meshID, type, &faceSize, &faceMapSize, &edgeMapSize,
989          faceMap, edgeMap);
990   
991   // look for vertexList in the edge map
992   std::set<int> vertexSet(vertexList, vertexList+2);
993   for (int i=0; i<edgeMapSize; i++) {
994     std::set<int> aNodeSet(edgeMap[i], edgeMap[i]+2);
995     if (vertexSet == aNodeSet) return i;
996   }
997
998   // uh-oh, didn't find it
999   CkAbort("ERROR: GetEdgeFace: vertexList is not a valid edge./n");
1000 }
1001
1002
1003 // Update functions
1004 /** Look up elemID in elemType array and set the adjacency on
1005     edgeFaceID to nbr. */
1006 void setAdaptAdj(
1007         const int meshID, 
1008         const adaptAdj elem, 
1009         const int faceID, 
1010         const adaptAdj nbr)
1011 {
1012   int numAdjacencies;
1013   adaptAdj *adaptAdjTable = lookupAdaptAdjacencies(meshID, elem.elemType,
1014                                                    &numAdjacencies);
1015   adaptAdjTable[elem.localID*numAdjacencies + faceID] = nbr;
1016 }
1017
1018 void addToAdaptAdj(
1019         const int meshid, 
1020         const adaptAdj elem, 
1021         const int edgeID, 
1022         const adaptAdj nbr)
1023 {
1024     CkVec<adaptAdj>* adjVec = getEdgeAdaptAdj(meshid, elem.localID, 
1025             elem.elemType, edgeID);
1026     adjVec->push_back(nbr);
1027 }
1028
1029 void removeFromAdaptAdj(
1030         const int meshid, 
1031         const adaptAdj elem, 
1032         const int edgeID, 
1033         const adaptAdj nbr)
1034 {
1035     CkVec<adaptAdj>* adjVec = getEdgeAdaptAdj(meshid, elem.localID, 
1036             elem.elemType, edgeID);
1037     for (int i=0; i<adjVec->size(); ++i) {
1038         if ((*adjVec)[i] == nbr) {
1039             adjVec->remove(i);
1040             return;
1041         }
1042     }
1043     CkAbort("removeFromAdaptAdj did not find the specified nbr");
1044 }
1045
1046 /** Lookup elemID in elemType array and search for the face which has 
1047  *  originalNbr as a neighbor, then replace originalNbr with newNbr
1048  */
1049 void replaceAdaptAdj(
1050         const FEM_Mesh* const meshPtr, 
1051         const adaptAdj elem, 
1052         const adaptAdj originalNbr, 
1053         const adaptAdj newNbr)
1054 {
1055     int numAdjacencies;
1056     adaptAdj *adaptAdjTable = lookupAdaptAdjacencies(meshPtr, elem.elemType,
1057             &numAdjacencies);
1058     for(int i=0;i<numAdjacencies;i++){
1059         if(adaptAdjTable[elem.localID*numAdjacencies+i] == originalNbr){
1060             adaptAdjTable[elem.localID*numAdjacencies+i] = newNbr;
1061             return;
1062         }
1063     }
1064     CkAbort("replaceAdaptAdj did not find the specified originalNbr");
1065 }
1066
1067
1068 void replaceAdaptAdj(
1069         const int meshID, 
1070         const adaptAdj elem, 
1071         const adaptAdj originalNbr, 
1072         const adaptAdj newNbr)
1073 {
1074     FEM_Mesh* meshPtr = FEM_chunk::get("ReplaceAdaptAdj")->lookup(
1075             meshID,"ReplaceAdaptAdj");
1076     replaceAdaptAdj(meshPtr, elem, originalNbr, newNbr);
1077 }
1078
1079
1080 void replaceAdaptAdjOnEdge(
1081         const FEM_Mesh* const meshPtr, 
1082         const adaptAdj elem, 
1083         const adaptAdj originalNbr,
1084         const adaptAdj newNbr,
1085         const int edgeID)
1086 {
1087     int numAdjacencies;
1088     CkVec<adaptAdj>** adaptAdjTable = lookupEdgeAdaptAdjacencies(
1089             meshPtr, elem.elemType, &numAdjacencies);
1090     CkVec<adaptAdj>* innerTable = 
1091         adaptAdjTable[elem.localID*numAdjacencies + edgeID];
1092     for (int n=0; n<innerTable->size(); ++n) {
1093         if((*innerTable)[n] == originalNbr){
1094             (*innerTable)[n] = newNbr;
1095             return;
1096         }
1097     }
1098     CkAbort("replaceAdaptAdjOnEdge did not find the specified originalNbr");
1099 }
1100
1101
1102 void replaceAdaptAdjOnEdge(
1103         const int meshID, 
1104         const adaptAdj elem, 
1105         const adaptAdj originalNbr, 
1106         const adaptAdj newNbr,
1107         const int edgeID)
1108 {
1109     FEM_Mesh* meshPtr = FEM_chunk::get("ReplaceAdaptAdj")->lookup(
1110             meshID,"ReplaceAdaptAdj");
1111     replaceAdaptAdjOnEdge(meshPtr, elem, originalNbr, newNbr, edgeID);
1112 }
1113
1114
1115 #define copyNodeMap(numEntries,entrySize,map,srcMap) \
1116 do { \
1117   for (int i=0;i<numEntries;i++){ \
1118     for(int j=0;j<entrySize;j++){ \
1119       map[i][j] = srcMap[i][j]; \
1120     } \
1121   } \
1122 } while (0) 
1123
1124
1125 // Given the dimensions and nodes per element guess whether the element 
1126 // is a triangle, quad, tet or hex. At the moment these are the 4 shapes
1127 // that are handled  
1128 void guessElementShape(
1129         const int dim, 
1130         const int nodesPerElem, 
1131         int* faceSize, 
1132         int *faceMapSize, 
1133         int* edgeMapSize,
1134         int faceMap[MAX_ADJELEMS][MAX_FACE_SIZE], 
1135         int edgeMap[MAX_EDGES][2])
1136 {
1137     switch(dim){
1138         case 2:
1139             {
1140                 //2 dimension
1141                 switch (nodesPerElem) {
1142                     case 3:
1143                         //Triangles
1144                         *faceSize = 2;
1145                         *faceMapSize = 3;
1146                         *edgeMapSize = 0;
1147                         copyNodeMap(3,2,faceMap,faceMap2d_tri);
1148                         break;
1149                     case 4:
1150                         //quads
1151                         *faceSize = 2;
1152                         *faceMapSize = 4;
1153                         *edgeMapSize = 0;
1154                         copyNodeMap(4,2,faceMap,faceMap2d_quad);
1155                         break;
1156                     default:
1157                         CkPrintf("Unknown element type\n");
1158                         CkAssert(false);
1159                 }
1160             }
1161             break;
1162         case 3:
1163             {
1164                 //3 dimension
1165                 switch(nodesPerElem){
1166                     case 4:
1167                         //Tetrahedra
1168                         *faceSize = 3;
1169                         *faceMapSize = 4;
1170                         *edgeMapSize = 6;
1171                         copyNodeMap(4,3,faceMap,faceMap3d_tet);
1172                         copyNodeMap(6,2,edgeMap,edgeMap3d_tet);
1173                         break;
1174                     case 8:
1175                         //Hexahedra
1176                         *faceSize = 4;
1177                         *faceMapSize = 6;
1178                         *edgeMapSize = 12;
1179                         copyNodeMap(6,4,faceMap,faceMap3d_hex);
1180                         copyNodeMap(12,2,edgeMap,edgeMap3d_hex);
1181                         break;
1182                     default:
1183                         CkPrintf("Unknown element type\n");
1184                         CkAssert(false);
1185                 }
1186             }
1187             break;
1188         default:
1189             CkPrintf("Unknown element type\n");
1190             CkAssert(false);
1191     }
1192 }
1193