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