13828ba9db03d9f2b7dc0c392391eefe575ee1eb
[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 int nodeSetMap2d_tri[3][2] = {{0,1},{1,2},{2,0}};
12 int nodeSetMap2d_quad[4][2] = {{0,1},{1,2},{2,3},{3,0}};
13 int nodeSetMap3d_tet[4][3] = {{0,1,2},{0,3,1},{0,2,3},{1,3,2}};
14 int nodeSetMap3d_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 int nodeSetMap2d_cohquad[2][2] = {{0,1},{2,3}}; // Cohesive for 2D triangles
17 int nodeSetMap3d_cohprism[2][3] = {{0,1,2},{3,4,5}}; // Cohesive for 3D tets
18
19 inline void addSharedNodeData(int node,const IDXL_Rec *sharedChunkList,
20             adjNode *adaptAdjTable){
21   adaptAdjTable[node].numSharedPartitions = sharedChunkList->getShared();
22   adaptAdjTable[node].sharedWithPartition = 
23     new int [sharedChunkList->getShared()];
24   adaptAdjTable[node].sharedWithLocalIdx = 
25     new int [sharedChunkList->getShared()];
26   for(int j=0;j<sharedChunkList->getShared();j++){
27     int sharedChunk = sharedChunkList->getChk(j);
28     int sharedIdx = sharedChunkList->getIdx(j);
29     adaptAdjTable[node].sharedWithPartition[j] = sharedChunk;
30     adaptAdjTable[node].sharedWithLocalIdx[j] = sharedIdx;
31   }
32 }
33
34 inline void addElementNodeSetData(int elem,const int *conn,int numAdjElems,
35           const int nodesPerElem,int nodeSetSize,
36               int nodeSetMap[MAX_ADJELEMS][MAX_NODESET_SIZE],
37           adjNode *adaptAdjTable){
38   for (int j=0; j<numAdjElems; j++) { // one nodeSet per neighbor element
39     adjElem *e = new adjElem(nodeSetSize);
40     e->nodeSetID = j;
41     for (int k=0; k<nodeSetSize; k++) { // Build nodeSet for an element pairing
42       e->nodeSet[k] = conn[elem*nodesPerElem+nodeSetMap[j][k]];
43     }
44     // Add this element-nodeSet pair to table at min nodeID in the nodeSet
45     e->nodeSet.quickSort();
46     int minNode = e->nodeSet[0];
47     e->elemID = elem;
48     e->next = adaptAdjTable[minNode].adjElemList->next;
49     adaptAdjTable[minNode].adjElemList->next = e;
50     adaptAdjTable[minNode].adjElemCount++;
51   }
52 }
53
54
55 //Look for an adjElem (rover->next) whose nodeSet matches that
56 //specified in searchForNodeSet in the link list following adjStart
57 //and return rover such that rover->next matches with the
58 //searchForNodeSet. It also checks that the elemID of the element
59 //being searched does not match with that of searchForElemID.  *found
60 //is set to 1 if match is found .. else to 0
61 inline adjElem *searchAdjElemInList(adjElem *adjStart,int *searchForNodeSet,
62             int nodeSetSize,int searchForElemID,
63             int *found){
64   adjElem *rover = adjStart; // compare rover->next with adjStart
65   *found = 0;
66   
67   while (rover->next != NULL) {
68     if (rover->next->elemID != searchForElemID) {
69       *found = 1; // found an element that is not myself, possibly a match
70       for (int j=0; j<nodeSetSize; j++) {
71         if (rover->next->nodeSet[j] != searchForNodeSet[j]) {
72           *found = 0; // No, the nodeSets dont match
73           break;
74         }
75       }
76     }
77     if (*found) {
78       break; // We have found a nodeSet that matches adjStart
79     }else {
80       rover = rover->next; // Keep looking in adjElemList for matching nodeSet
81     }
82   }
83   return rover;
84 }
85
86 /** Create Adaptivity Adjacencies for elemType; dimension inferred. */
87 void CreateAdaptAdjacencies(int meshid, int elemType)
88 {
89   // Need to derive all of these from elemType;
90   int myRank,numChunks;
91   int numAdjElems;
92   int nodeSetSize; // number of nodes shared by two adjacent elems
93
94   MPI_Comm_rank(MPI_COMM_WORLD,&myRank);
95   MPI_Comm_size(MPI_COMM_WORLD,&numChunks);
96   FEM_Mesh *mesh = FEM_chunk::get("CreateAdaptAdjacencies")->lookup(meshid,"CreateAdaptAdjacencies");
97   FEM_Elem *elem = (FEM_Elem *)mesh->lookup(FEM_ELEM+elemType,"CreateAdaptAdjacencies");
98   FEM_Node *node = (FEM_Node *)mesh->lookup(FEM_NODE,"CreateAdaptAdjacencies");
99   const int numElems = elem->size();
100   const int numNodes = node->size();
101   const int nodesPerElem = (elem->getConn()).width();
102   assert(node->getCoord()!= NULL);
103   const int dim = (node->getCoord())->getWidth();
104   assert(dim == 2|| dim == 3);
105
106   // A nodeSet is a set of nodes that defines a pairing of two
107   // adjacent elements; For example, in 2D triangle meshes, the
108   // nodeSet is the nodes of an edge between two elements.  The
109   // nodeSetMap is an ordering of element-local node IDs that
110   // specifies all possible nodeSets for a particular element type
111   int nodeSetMap[MAX_ADJELEMS][MAX_NODESET_SIZE];
112   guessElementShape(dim,nodesPerElem,&numAdjElems,&nodeSetSize,nodeSetMap);
113   CkAssert(nodeSetSize <= MAX_NODESET_SIZE);
114   
115   // Add the FEM_ADAPT_ADJ attribute to the elements
116   // Set the correct width of the table and then get the pointer to the actual data
117   FEM_DataAttribute *adaptAdjAttr = (FEM_DataAttribute *) elem->lookup(FEM_ADAPT_ADJ,"CreateAdaptAdjacencies"); 
118   adaptAdjAttr->setWidth(sizeof(adaptAdj)*numAdjElems);
119   adaptAdj *adaptAdjacencies = (adaptAdj *)(adaptAdjAttr->getChar()).getData();
120
121   // Init adaptAdj array to at least have -1 as partID signifying no neighbor
122   for (int i=0; i<numElems*numAdjElems; i++) {
123     adaptAdjacencies[i].partID = -1;
124     adaptAdjacencies[i].localID = -1;
125   }
126
127   // Create an array of size equal to the number of local nodes, each
128   // entry has a pair: (shared partitions, list of elem-nodeSet pairs)
129   // Call this adaptAdjTable
130   adjNode *adaptAdjTable;
131   adaptAdjTable = new adjNode[numNodes];
132
133   // Loop through shared node list and add the partition ids to adaptAdjTable.
134   for(int i=0;i<numNodes;i++){
135     if(node->is_valid(i)){
136       const IDXL_Rec *sharedChunkList = node->shared.getRec(i);
137       if(sharedChunkList != NULL){
138         addSharedNodeData(i,sharedChunkList,adaptAdjTable);
139       }
140     }
141   }
142
143   // Pull out conn for elems of elemType
144   const int *conn = (elem->getConn()).getData();
145   
146   for (int i=0; i<numElems; i++) { // Add each element-nodeSet pair to table
147     if(elem->is_valid(i)){
148       addElementNodeSetData(i,conn,numAdjElems,nodesPerElem,nodeSetSize,
149           nodeSetMap,adaptAdjTable);
150     }
151   }
152
153   fillLocalAdaptAdjacencies(numNodes,node,adaptAdjTable,adaptAdjacencies,nodeSetSize,numAdjElems,myRank,elemType);
154
155
156   // Now all elements' local adjacencies are set; remainder in table are 
157   // nodeSets shared with other chunks or nodeSets on domain boundary
158
159   MPI_Barrier(MPI_COMM_WORLD);
160   MSA1DREQLIST *requestTable;
161   if(myRank == 0){
162     requestTable = new MSA1DREQLIST(numChunks,numChunks);
163   }else{
164     requestTable = new MSA1DREQLIST;
165   }
166   MPI_Bcast_pup(*requestTable,0,MPI_COMM_WORLD);
167   requestTable->enroll(numChunks);
168   requestTable->sync();
169
170   makeAdjacencyRequests(numNodes,node,adaptAdjTable,requestTable,nodeSetSize,myRank,elemType);
171
172   requestTable->sync();
173   printf("[%d] All requests made \n",myRank);
174
175   MSA1DREPLYLIST *replyTable;
176   if(myRank == 0){
177     replyTable = new MSA1DREPLYLIST(numChunks,numChunks);
178   }else{
179     replyTable = new MSA1DREPLYLIST;
180   }
181   MPI_Bcast_pup(*replyTable,0,MPI_COMM_WORLD);
182   replyTable->enroll(numChunks);
183   replyTable->sync();
184
185         
186   replyAdjacencyRequests(requestTable,replyTable,node,adaptAdjTable,adaptAdjacencies,nodeSetSize,numAdjElems,myRank,elemType);
187
188   requestTable->sync();
189   replyTable->sync();
190
191   //Once the replies are back, loop through each reply and update the
192   //adaptAdjacencies for each element in the reply
193   CkVec<adjReply> *receivedReplyVec = replyTable->get(myRank).vec;
194   for(int i=0;i< receivedReplyVec->size();i++){
195     adjReply *receivedReply = &(*receivedReplyVec)[i];
196     printf("[%d] Replies received for (%d,%d) (%d,%d,%d)\n",myRank,receivedReply->requestingElemID,receivedReply->requestingNodeSetID,receivedReply->replyingElem.partID,receivedReply->replyingElem.localID,receivedReply->replyingElem.elemType);
197     adaptAdjacencies[receivedReply->requestingElemID*numAdjElems + receivedReply->requestingNodeSetID] = receivedReply->replyingElem;
198   }
199
200   replyTable->sync();
201   dumpAdaptAdjacencies(adaptAdjacencies,numElems,numAdjElems,myRank);
202
203 /*  mesh->becomeSetting();  
204   //Register the adaptAdjacency with ParFUM
205   FEM_Register_array(meshid,FEM_ELEM+elemType,FEM_ADAPT_ADJ,(void *)adaptAdjacencies,FEM_BYTE,sizeof(adaptAdj)*numAdjElems);
206   //do not delete adaptAdjacencies. It will be used during the rest of adaptivity
207   mesh->becomeGetting();  
208         */
209 }
210
211 void fillLocalAdaptAdjacencies(int numNodes,FEM_Node *node,adjNode *adaptAdjTable,adaptAdj *adaptAdjacencies,int nodeSetSize,int numAdjElems,int myRank,int elemType){
212
213   for (int i=0; i<numNodes; i++) { 
214     // For each node, match up incident elements
215     // Each adjacency between two elements is represented by two adjElems
216     // We try to match those up
217     if(node->is_valid(i) && adaptAdjTable[i].adjElemList != NULL){  
218       // CkAssert(adaptAdjTable[i].adjElemList != NULL);
219       adjElem *preTarget = adaptAdjTable[i].adjElemList;
220       adjElem *target = adaptAdjTable[i].adjElemList->next;
221       while (target != NULL) { //loop over adjElemList of a node
222         int found = 0; 
223         //target represents an adjacency between two elements
224         //We search for the other adjElem corresponding to that adjancency:
225         //Look for an entry in adjElemList after target such that 
226         //the nodeset of that entry and that of target match but they 
227         //do not belong to the same element. 
228         adjElem *rover = searchAdjElemInList(preTarget,
229                                              target->nodeSet.getVec(),
230                                              nodeSetSize,target->elemID,
231                                              &found); 
232         
233         if (found) { // We found a local element adjacent to target->elemID
234           // Set adjacency of target->elemID corresponding to nodeSet to 
235           // rover->next->elemID, and vice versa
236           // Store adjacency info in adaptAdjacency of each one and
237           // use nodeSetID to index into adaptAdjacency
238           adaptAdjacencies[numAdjElems*target->elemID+target->nodeSetID] = 
239             adaptAdj(myRank,rover->next->elemID,elemType);
240           adaptAdjacencies[numAdjElems*rover->next->elemID+rover->next->nodeSetID] = adaptAdj(myRank,target->elemID,elemType);
241           // Remove both elem-nodeSet pairs from the list
242           adjElem *tmp = rover->next;
243           rover->next = rover->next->next;
244           delete tmp;
245           tmp = target;
246           preTarget->next = target->next;
247           target = target->next;
248           delete tmp;
249         }else { // No match for target was found in adjElemList
250           // This means that either target is on the domain boundary 
251           // or it is on a partition boundary and its neighbor is on another VP
252           // Move target to next entry in adjElemList
253           preTarget = target;
254           target = target->next;
255         }
256       }
257     }
258   }
259 }
260
261
262 void makeAdjacencyRequests(int numNodes,FEM_Node *node,adjNode *adaptAdjTable,MSA1DREQLIST *requestTable, int nodeSetSize,int myRank,int elemType){
263
264   for (int i=0; i<numNodes; i++) { // Examine each node's remaining elements
265     if(node->is_valid(i)){
266       if (adaptAdjTable[i].adjElemList->next!=NULL) {
267         adjElem *adjStart = adaptAdjTable[i].adjElemList->next;
268         while (adjStart !=NULL) {
269           // create and empty set, commonSharedChunks
270           std::set<int> commonSharedChunks;
271           for (int j=0; j<nodeSetSize; j++) {
272             // look up sharedWithPartitions for node: 
273             //    adaptAdjTable[i]->adjElemList->nodeset[j]
274             // intersect with commonSharedChunks
275             int sharedNode = adjStart->nodeSet[j];
276             adjNode *sharedNodeAdj = &adaptAdjTable[sharedNode];
277             std::set<int> sharedChunks(sharedNodeAdj->sharedWithPartition,
278                                        sharedNodeAdj->sharedWithPartition+sharedNodeAdj->numSharedPartitions);
279             if(j == 0){
280               commonSharedChunks = sharedChunks;
281             }else{
282               std::set<int> tmpIntersect;
283               set_intersection(commonSharedChunks.begin(), commonSharedChunks.end(), 
284              sharedChunks.begin(), sharedChunks.end(),inserter(tmpIntersect, tmpIntersect.begin()));
285               commonSharedChunks = tmpIntersect;              
286             }
287           }
288           // At this point commonSharedChunks contains the list of
289           // chunks with which the element pointed by adjStart might
290           // be shared
291           int numCommonSharedChunks = commonSharedChunks.size();
292           if(numCommonSharedChunks > 0){
293       //adjStart is possibly shared with these chunks. It is
294       //shared across the adjStart->nodeSet set of nodes.
295             adjRequest *adjRequestList = new adjRequest[numCommonSharedChunks];
296             //Translate the nodes in the nodeSet into the index in the
297             //idxl list for each chunk in the commonSharedChunks
298             for(int j=0;j<nodeSetSize;j++){
299               int sharedNode = adjStart->nodeSet[j];
300               const IDXL_Rec *recSharedNode = node->shared.getRec(sharedNode);
301               int countChunk=0;
302               for(std::set<int>::iterator chunkIterator = commonSharedChunks.begin();chunkIterator != commonSharedChunks.end();chunkIterator++){
303                 int chunk = *chunkIterator;
304                 // if(chunk > myRank){
305                 if(j == 0){
306                   // if this is the first node we need to initialize
307                   // the adjRequestobject
308                   adjRequestList[countChunk] =  adjRequest(adjStart->elemID,myRank,adjStart->nodeSetID,elemType);
309                 }
310                 int sharedNodeIdx=-1; // index of sharedNode in the idxl list of chunk
311                 //search for this chunk in the list of chunks shared
312                 //by this node
313                 for(int k=0;k<recSharedNode->getShared();k++){
314                   if(recSharedNode->getChk(k) == chunk){
315                     //found the correct chunk
316                     sharedNodeIdx = recSharedNode->getIdx(k);
317                     break;
318                   }
319                 }
320                 CkAssert(sharedNodeIdx != -1);
321                 //The index of sharedNode in the index list of chunk
322                 //has been found.  this needs to be saved in the
323                 //corresponding translatedNodeSet
324                 adjRequestList[countChunk].translatedNodeSet[j] = sharedNodeIdx;
325                 //                }
326                 countChunk++;
327               }
328             }
329             //Now the nodeNumbers for the nodeSets that might be along
330             //chunk boundaries have been translated into idxl indices
331             //between the two chunks We now need to write the requests
332             //into the msa array requestTable WARNING: This depends on
333             //sets getting enumerated in the same way always Might not
334             //be true
335             int countChunk=0;
336             for(std::set<int>::iterator chunkIterator = commonSharedChunks.begin();chunkIterator != commonSharedChunks.end();chunkIterator++){
337               int chunk = *chunkIterator;
338 //              if(chunk > myRank){
339                 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);
340                 (*requestTable).accumulate(chunk,adjRequestList[countChunk]);
341 //              requestVec.push_back(&adjRequestList[countChunk]);
342 //              }
343               countChunk++;
344             }
345             delete [] adjRequestList;
346           }
347           adjStart = adjStart->next;
348         }
349       }
350     }
351   }
352 }
353
354 void replyAdjacencyRequests(MSA1DREQLIST *requestTable,MSA1DREPLYLIST *replyTable,FEM_Node *node,adjNode *adaptAdjTable,adaptAdj *adaptAdjacencies,int nodeSetSize,int numAdjElems,int myRank,int elemType){
355   // look up request table, put answer in reply table 
356   // receive: for local elemID, the adjacency on this set of shared indices is (remote elemID, remote partition ID, remote elem type)
357   // add adjacencies to local table
358   // lots of error checking :)
359   
360   //Look at each request that in the requestTable for this chunk
361   //Put the data for the requests in our own table and then create replies
362   CkVec<adjRequest> *receivedRequestVec = requestTable->get(myRank).vec;
363   for(int i=0;i<receivedRequestVec->length();i++){    
364     adjRequest &receivedRequest = (*receivedRequestVec)[i];
365     const IDXL_List &sharedNodeList = node->shared.getList(receivedRequest.chunkID);
366     CkVec<int> sharedNodes(nodeSetSize);
367     //Translate all the nodes in the nodeSet of the request
368     for(int j=0;j<nodeSetSize;j++){
369       int sharedNode = sharedNodeList[receivedRequest.translatedNodeSet[j]];
370       sharedNodes[j] = sharedNode;
371     }
372     sharedNodes.quickSort();
373     //We need to find the matching nodeset for the nodeset in the request
374     //We look it up in the adaptAdjTable for the minimum node on this chunk
375     adjNode *minNode = &adaptAdjTable[sharedNodes[0]];
376     int found=0;
377     //search for the nodeSet in the list of adjacencies around minNode
378     adjElem *rover =  searchAdjElemInList(minNode->adjElemList,sharedNodes.getVec(),nodeSetSize,-1,&found);
379     printf("[%d] Received request (%d,%d,%d,%d) minNode %d found %d\n",myRank,receivedRequest.elemID,receivedRequest.chunkID,receivedRequest.elemType,receivedRequest.nodeSetID,sharedNodes[0],found);
380     if(found){ //we have found a matching adjElem for the requested nodeset
381        //we shall set the adjacency correctly in the adjacency Table
382        //for the elemID in the found adjElem. We need to send a reply
383        //to the requesting chunk
384        int matchingElemID;
385        matchingElemID = rover->next->elemID;
386        adaptAdj *matchingAdaptAdj = &adaptAdjacencies[matchingElemID*numAdjElems + rover->next->nodeSetID];
387        adjElem *tmp = rover->next;
388        rover->next = tmp->next;
389        delete tmp;
390        CkAssert(matchingAdaptAdj->localID == -1);
391        matchingAdaptAdj->partID = receivedRequest.chunkID;
392        matchingAdaptAdj->localID = receivedRequest.elemID;
393        matchingAdaptAdj->elemType = receivedRequest.elemType;
394
395        //Set requesting data in reply to that in receivedRequest
396        //Put in data from rover->next into the replyElem portion of data
397        adjReply reply;
398        reply.requestingElemID = receivedRequest.elemID;
399        reply.requestingNodeSetID = receivedRequest.nodeSetID;
400        reply.replyingElem.partID = myRank;
401        reply.replyingElem.localID = matchingElemID;
402        reply.replyingElem.elemType = elemType;
403        //Write into the replyTable
404        replyTable->accumulate(receivedRequest.chunkID,reply);
405     }else{
406       //we have no matching nodeset for this request.. hopefully some
407       //other chunk does; we can ignore this request
408     }
409   }
410 }
411
412
413
414 void dumpAdaptAdjacencies(adaptAdj *adaptAdjacencies,int numElems,int numAdjElems,int myRank){
415   for(int i=0;i<numElems;i++){
416     printf("[%d] %d  :",myRank,i);
417     for(int j=0;j<numAdjElems;j++){
418       adaptAdj *entry = &adaptAdjacencies[i*numAdjElems+j];
419       printf("(%d,%d,%d)",entry->partID,entry->localID,entry->elemType);
420     }
421     printf("\n");
422   }
423 }
424
425 inline void findNodeSet(int meshid,int elemType,int *numAdjElems,int *nodeSetSize,int  nodeSetMap[MAX_ADJELEMS][MAX_NODESET_SIZE]){
426         FEM_Mesh *mesh = FEM_chunk::get("GetAdaptAdj")->lookup(meshid,"GetAdaptAdj");
427   FEM_Elem *elem = (FEM_Elem *)mesh->lookup(FEM_ELEM+elemType,"GetAdaptAdj");
428   FEM_Node *node = (FEM_Node *)mesh->lookup(FEM_NODE,"GetAdaptAdj");
429   const int nodesPer = (elem->getConn()).width();
430   assert(node->getCoord()!= NULL);
431   const int dim = (node->getCoord())->getWidth();
432   assert(dim == 2|| dim == 3);
433   guessElementShape(dim,nodesPer,numAdjElems,nodeSetSize,nodeSetMap);
434 }
435
436 void getAndDumpAdaptAdjacencies(int meshid, int numElems, int elemType, int myRank){
437   int numAdjElems;
438   int nodeSetSize, nodeSetMap[MAX_ADJELEMS][MAX_NODESET_SIZE]; // not used
439         findNodeSet(meshid,elemType,&numAdjElems,&nodeSetSize,nodeSetMap);
440   
441   for(int i=0;i<numElems;i++){
442     printf("[%d] %d  :",myRank,i);
443     for(int j=0;j<numAdjElems;j++){
444       adaptAdj *entry = GetAdaptAdj(meshid, adaptAdj(myRank,i,elemType), j);
445       printf("(%d,%d,%d)",entry->partID,entry->localID,entry->elemType);
446     }
447     printf("\n");
448   }
449 }
450
451 // Access functions
452 inline adaptAdj *lookupAdaptAdjacencies(FEM_Mesh *mesh,int elemType,int *numAdjacencies){
453   FEM_Elem *elem = (FEM_Elem *)mesh->lookup(FEM_ELEM+elemType,"lookupAdaptAdjacencies");
454
455   FEM_DataAttribute *adaptAttr = (FEM_DataAttribute *)elem->lookup(FEM_ADAPT_ADJ,"lookupAdaptAdjacencies");
456         *numAdjacencies = adaptAttr->getWidth()/sizeof(adaptAdj);
457   AllocTable2d<unsigned char> &table = adaptAttr->getChar();
458
459   return (adaptAdj  *)table.getData();
460 }
461
462 inline adaptAdj *lookupAdaptAdjacencies(int meshid,int elemType,int *numAdjacencies){
463   FEM_Mesh *mesh = FEM_chunk::get("lookupAdaptAdjacencies")->lookup(meshid,"lookupAdaptAdjacencies");
464   return lookupAdaptAdjacencies(mesh, elemType, numAdjacencies);
465 }
466
467 /** Look up elemID in elemType array, access edgeFaceID-th adaptAdj. */
468 adaptAdj *GetAdaptAdj(int meshid, adaptAdj elem, int edgeFaceID)
469 {
470   int numAdjacencies;
471   adaptAdj *adaptAdjacencies = lookupAdaptAdjacencies(meshid, elem.elemType,
472                                                       &numAdjacencies);
473   return(&(adaptAdjacencies[elem.localID*numAdjacencies+edgeFaceID]));
474 }
475
476 adaptAdj *GetAdaptAdj(FEM_Mesh *meshPtr, adaptAdj elem, int edgeFaceID)
477 {
478   int numAdjacencies;
479   adaptAdj *adaptAdjacencies = lookupAdaptAdjacencies(meshPtr, elem.elemType,
480                                                       &numAdjacencies);
481   return(&(adaptAdjacencies[elem.localID*numAdjacencies+edgeFaceID]));
482 }
483
484 /** Look up elemID in elemType array, calculate edgeFaceID from
485     vertexList (with GetEdgeFace below), and access edgeFaceID-th
486     adaptAdj with GetAdaptAdj above. */
487 adaptAdj *GetAdaptAdj(int meshid, adaptAdj elem, int *vertexList)
488 {
489   int edgeFaceID = GetEdgeFace(meshid, elem, vertexList);
490   return GetAdaptAdj(meshid, elem, edgeFaceID);
491 }
492
493 adaptAdj *GetAdaptAdjOnEdge(int meshid, adaptAdj elem, int edgeID, int *size)
494 {
495   return NULL;
496 }
497
498 adaptAdj *GetAdaptAdjOnEdge(FEM_Mesh *meshPtr, adaptAdj elem, int edgeID, int *size)
499 {
500   return NULL;
501 }
502
503 adaptAdj *GetAdaptAdjOnFace(int meshid, adaptAdj elem, int faceID)
504 {
505   return NULL;
506 }
507
508 adaptAdj *GetAdaptAdjOnFace(FEM_Mesh *meshPtr, adaptAdj elem, int faceID)
509 {
510   return NULL;
511 }
512
513
514 /** Look up elemID in elemType array and determine the set of vertices
515     associated with the edge or face represented by edgeFaceID. */
516 void GetVertices(int meshid, adaptAdj elem, int edgeFaceID, int *vertexList)
517 {
518   int numAdjacencies;
519   int nodeSetSize, nodeSetMap[MAX_ADJELEMS][MAX_NODESET_SIZE];
520   findNodeSet(meshid,elem.elemType,&numAdjacencies,&nodeSetSize,nodeSetMap);
521   
522   for (int i=0; i<nodeSetSize; i++) {
523     vertexList[i] = nodeSetMap[edgeFaceID][i];
524   }
525 }
526
527 /** Look up elemID in elemType array and determine the edge or face ID
528     specified by the set of vertices in vertexList. */
529 int GetEdgeFace(int meshid, adaptAdj elem, int *vertexList)
530 {
531   int numAdjacencies;
532   int nodeSetSize, nodeSetMap[MAX_ADJELEMS][MAX_NODESET_SIZE];
533   findNodeSet(meshid,elem.elemType,&numAdjacencies,&nodeSetSize,nodeSetMap);
534   
535   std::set<int> vertexSet(vertexList, vertexList+nodeSetSize);
536   for (int i=0; i<numAdjacencies; i++) {
537     std::set<int> aNodeSet(nodeSetMap[i], nodeSetMap[i]+nodeSetSize);
538     if (vertexSet == aNodeSet) return i; // CHECK: does set equivalence exist?
539   }
540   CkAbort("ERROR: GetEdgeFace: vertexList is not a valid nodeSet./n");
541 }
542
543 // Update functions
544 /** Look up elemID in elemType array and set the adjacency on
545     edgeFaceID to nbr. */
546 void SetAdaptAdj(int meshID, adaptAdj elem, int edgeFaceID, adaptAdj nbr)
547 {
548   int numAdjacencies;
549   adaptAdj *adaptAdjTable = lookupAdaptAdjacencies(meshID, elem.elemType,
550                                                    &numAdjacencies);
551   adaptAdjTable[elem.localID*numAdjacencies + edgeFaceID] = nbr;
552 }
553
554 void AddToAdaptAdj(int meshid, adaptAdj elem, int edgeID, adaptAdj nbr)
555 {
556 }
557
558 void RemoveFromAdaptAdj(int meshid, adaptAdj elem, int edgeID, adaptAdj nbr)
559 {
560 }
561
562 /** Lookup elemID in elemType array and search for the edgeID which has originalNbr as
563  * a neighbor, then replace originalNbr with newNbr
564  */
565 void ReplaceAdaptAdj(FEM_Mesh *meshPtr, adaptAdj elem, adaptAdj originalNbr, 
566                      adaptAdj newNbr){
567   int numAdjacencies;
568   adaptAdj *adaptAdjTable = lookupAdaptAdjacencies(meshPtr, elem.elemType,
569                                                    &numAdjacencies);
570   for(int i=0;i<numAdjacencies;i++){
571     if(adaptAdjTable[elem.localID*numAdjacencies+i] == originalNbr){
572       adaptAdjTable[elem.localID*numAdjacencies+i] = newNbr;
573       return;
574     }
575   }
576   CkAbort("ReplaceAdaptAdj did not find the specified originalNbr");
577 }
578
579 void ReplaceAdaptAdj(int meshID, adaptAdj elem, adaptAdj originalNbr, 
580                      adaptAdj newNbr){
581   FEM_Mesh *meshPtr = FEM_chunk::get("ReplaceAdaptAdj")->lookup(meshID,"ReplaceAdaptAdj");
582   ReplaceAdaptAdj(meshPtr, elem, originalNbr, newNbr);
583 }
584
585 void ReplaceAdaptAdjOnEdge(int meshID, adaptAdj elem, adaptAdj originalNbr, 
586                      adaptAdj newNbr)
587 {
588 }
589
590 void ReplaceAdaptAdjOnEdge(FEM_Mesh *meshPtr, adaptAdj elem, adaptAdj originalNbr,
591                      adaptAdj newNbr)
592 {
593 }
594
595 //given the dimensions and nodes per element guess whether the element 
596 // is a triangle, quad, tet or hex. At the moment these are the 4 shapes
597 // that are handled
598   
599 #define copyNodeSetMap(numAdjElems,nodeSetSize,nodeSetMap,srcMap) {\
600   for (int i=0;i<numAdjElems;i++){ \
601     for(int j=0;j<nodeSetSize;j++){ \
602       nodeSetMap[i][j] = srcMap[i][j]; \
603     } \
604   } \
605 }
606
607 void guessElementShape(int dim,int nodesPerElem,int *numAdjElems,int *nodeSetSize,int nodeSetMap[MAX_ADJELEMS][MAX_NODESET_SIZE]){
608   switch(dim){
609   case 2:
610     {
611       //2 dimension
612       switch(nodesPerElem){
613       case 3:
614   //Triangles
615   *numAdjElems = 3;
616   *nodeSetSize = 2;
617   copyNodeSetMap(3,2,nodeSetMap,nodeSetMap2d_tri)
618     break;
619       case 4:
620   //quads
621   *numAdjElems = 4;
622   *nodeSetSize = 2;
623   copyNodeSetMap(4,2,nodeSetMap,nodeSetMap2d_quad)
624     break;
625       }
626     }
627     break;
628   case 3:
629     {
630       //3 dimension
631       switch(nodesPerElem){
632       case 4:
633   //Tetrahedra
634   *numAdjElems = 4;
635   *nodeSetSize = 3;
636   copyNodeSetMap(4,3,nodeSetMap,nodeSetMap3d_tet)
637     break;
638       case 6:
639   //Hexahedra
640   *numAdjElems = 6;
641   *nodeSetSize = 4;
642   copyNodeSetMap(6,4,nodeSetMap,nodeSetMap3d_hex)
643     break;
644       }
645     }
646     break;
647   }
648 }
649