Charm++: Charm Source Code Documentation

00001 /***************************************************
00002  * fem_mesh_adjacency.C
00003  *
00004  * All the adjacency creation and update functions.
00005  * The first half of the file is for initial creation
00006  * of all needed tables. The second half provides 
00007  * easy accessor and modifier functions.
00008  *
00009  * Authors include: Isaac, Sayantan, Terry, Nilesh
00010  */
00011 
00012 #include "ParFUM.h"
00013 #include "ParFUM_internals.h"
00014 
00015 
00016 //#define DEBUG
00017 
00018 #ifdef DEBUG
00019 FORTRAN_AS_C(CMIMEMORYCHECK,
00020              CmiMemoryCheck,
00021              cmimemorycheck, 
00022              (void),  () )
00023 #endif
00024 
00025 
00026 CLINKAGE void
00027 FEM_Mesh_create_node_elem_adjacency(int fem_mesh){
00028     const char *caller="FEM_Mesh_create_node_elem_adjacency"; FEMAPI(caller);
00029     FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,caller);
00030     m->createNodeElemAdj();
00031 }
00032 FORTRAN_AS_C(FEM_MESH_CREATE_NODE_ELEM_ADJACENCY,
00033              FEM_Mesh_create_node_elem_adjacency,
00034              fem_mesh_create_node_elem_adjacency, 
00035              (int *fem_mesh),  (*fem_mesh) )
00036 
00037 CLINKAGE void
00038 FEM_Mesh_create_node_node_adjacency(int fem_mesh){
00039     const char *caller="FEM_Mesh_create_node_node_adjacency"; FEMAPI(caller);
00040     FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,caller);
00041     m->createNodeNodeAdj();
00042 }
00043 FORTRAN_AS_C(FEM_MESH_CREATE_NODE_NODE_ADJACENCY,
00044              FEM_Mesh_create_node_node_adjacency,
00045              fem_mesh_create_node_node_adjacency, 
00046              (int *fem_mesh),  (*fem_mesh) )
00047 
00048 CLINKAGE void
00049 FEM_Mesh_create_elem_elem_adjacency(int fem_mesh){
00050     const char *caller="FEM_Mesh_create_elem_elem_adjacency"; FEMAPI(caller);
00051     FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,caller);
00052     m->createElemElemAdj();
00053 }
00054 FORTRAN_AS_C(FEM_MESH_CREATE_ELEM_ELEM_ADJACENCY,
00055              FEM_Mesh_create_elem_elem_adjacency,
00056              fem_mesh_create_elem_elem_adjacency, 
00057              (int *fem_mesh),  (*fem_mesh) )
00058 
00059 CLINKAGE void
00060 FEM_Mesh_create_elem_node_adjacency(int fem_mesh){
00061   CkPrintf("WARNING: Do Not Call FEM_Mesh_create_elem_node_adjacency(), as the connectivity table should already exist\n");
00062 }
00063 
00064 
00065 CLINKAGE void
00066 FEM_Mesh_get2ElementsOnEdge(int fem_mesh, int n1, int n2, int *e1, int *e2){
00067     const char *caller="FEM_Mesh_get2ElementsOnEdge"; FEMAPI(caller);
00068     FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,caller);
00069     m->get2ElementsOnEdge(n1, n2, e1, e2);
00070 }
00071 FORTRAN_AS_C(FEM_MESH_GET2ELEMENTSONEDGE,
00072              FEM_Mesh_get2ElementsOnEdge,
00073              fem_mesh_get2elementsonedge, 
00074              (int *fem_mesh, int *n1, int *n2, int *e1, int *e2),  
00075              (*fem_mesh,*n1,*n2,e1,e2) )
00076 
00077 
00078 CLINKAGE void
00079 FEM_Mesh_get2ElementsOnEdgeSorted(int fem_mesh, int n1, int n2, int *e1, int *e2){
00080   const char *caller="FEM_Mesh_get2ElementsOnEdge"; FEMAPI(caller);
00081   FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,caller);
00082   m->get2ElementsOnEdgeSorted(n1, n2, e1, e2);
00083 }
00084 FORTRAN_AS_C(FEM_MESH_GET2ELEMENTSONEDGESORTED,
00085              FEM_Mesh_get2ElementsOnEdgeSorted,
00086              fem_mesh_get2elementsonedgesorted,
00087              (int *fem_mesh, int *n1, int *n2, int *e1, int *e2),
00088              (*fem_mesh,*n1,*n2,e1,e2) )
00089 
00090 
00091 void FEM_Node::allocateElemAdjacency(){
00092     if(elemAdjacency){
00093         delete elemAdjacency;
00094     }
00095     elemAdjacency = new FEM_VarIndexAttribute(this,FEM_NODE_ELEM_ADJACENCY);
00096     add(elemAdjacency);
00097 
00098 }
00099 
00100 
00101 void FEM_Node::allocateNodeAdjacency(){
00102     if(nodeAdjacency){
00103         delete nodeAdjacency;
00104     }
00105     nodeAdjacency = new FEM_VarIndexAttribute(this,FEM_NODE_NODE_ADJACENCY);
00106     add(nodeAdjacency);
00107 }
00108 
00109 
00113 void FEM_Node::setElemAdjacency(int type, const FEM_Elem &elem){
00114     int nodesPerElem = elem.getNodesPer();
00115     FEM_VarIndexAttribute *adjacencyAttr = elemAdjacency;
00116     CkVec<CkVec<ElemID> > &adjacencyTable = elemAdjacency->get();
00117     FEM_VarIndexAttribute *ghostAdjacencyAttr = ((FEM_Node *)getGhost())->elemAdjacency;
00118     CkVec<CkVec<ElemID> > &ghostAdjacencyTable = ghostAdjacencyAttr->get();
00119 
00120     // Scan through elements
00121     for(int i=0;i<elem.size();i++){
00122       const int *conn = elem.connFor(i);
00123       for(int j=0;j<nodesPerElem;j++){
00124         int node = conn[j];
00125         if (node!=-1){
00126           if(FEM_Is_ghost_index(node)){
00127             int idx = ghostAdjacencyAttr->findInRow(FEM_To_ghost_index(node),var_id(type,i));
00128             if(idx == -1) {// If not currently in the adjacency list, push onto list
00129               ghostAdjacencyTable[FEM_To_ghost_index(node)].push_back(var_id(type,i));
00130             }
00131           }
00132           else{
00133             int idx = adjacencyAttr->findInRow(node,var_id(type,i));
00134             if(idx == -1) {// If not currently in the adjacency list, push onto list
00135               adjacencyTable[node].push_back(var_id(type,i));
00136             }
00137           }
00138         }
00139       }
00140     }
00141     
00142     // Scan through ghost elements
00143     if(elem.getGhost()){
00144       for(int i=0;i<((FEM_Elem*)elem.getGhost())->size();i++){
00145         const int *conn = ((FEM_Elem*)elem.getGhost())->connFor(i);
00146         for(int j=0;j<nodesPerElem;j++){
00147           int node = conn[j];
00148           if (node!=-1){
00149             if(FEM_Is_ghost_index(node)){
00150               int idx = ghostAdjacencyAttr->findInRow(FEM_To_ghost_index(node),var_id(type,FEM_From_ghost_index(i)));
00151               if(idx == -1){ // If not currently in the adjacency list, push onto list
00152                 ghostAdjacencyTable[FEM_To_ghost_index(node)].push_back(var_id(type,FEM_From_ghost_index(i)));
00153               }
00154             }
00155             
00156             else{
00157               int idx = adjacencyAttr->findInRow(node,var_id(type,FEM_From_ghost_index(i)));
00158               if(idx == -1){// If not currently in the adjacency list, push onto list
00159                 adjacencyTable[node].push_back(var_id(type,FEM_From_ghost_index(i)));
00160               }
00161             }
00162           }
00163         }
00164       }
00165     }
00166 }
00167 
00168 
00179 void FEM_Node::setNodeAdjacency(const FEM_Elem &elem){
00180   //CkPrintf("In FEM_Node::setNodeAdjacency()\n");
00181   int nodesPerElem = elem.getNodesPer();
00182   CkVec<CkVec<var_id> > &adjacencyTable = nodeAdjacency->get();
00183   FEM_VarIndexAttribute *ghostAdjacencyAttr = ((FEM_Node *)getGhost())->nodeAdjacency;
00184   CkVec<CkVec<var_id> > &ghostAdjacencyTable = ghostAdjacencyAttr->get();
00185   
00186   // Add the adjacencies defined by the non-ghost elements
00187   for(int i=0;i<elem.size();i++) {        // for each element of the given type
00188     const int *conn = elem.connFor(i);
00189     for(int j=0;j<nodesPerElem;j++){   // for each node adjacent to the element
00190       const int nodej = conn[j];
00191       if (nodej!=-1){
00192         if(FEM_Is_ghost_index(nodej)) { // A ghost node
00193           for(int k=0;k<nodesPerElem;k++){
00194             const int nodek=conn[k];
00195             if(nodek != nodej){
00196               var_id nodeID = var_id::createNodeID(1,nodek);
00197               int idx = ghostAdjacencyAttr->findInRow(FEM_To_ghost_index(nodej),nodeID);
00198               if(idx == -1){
00199                 //if(nodej==-5|| nodek==-5) CkPrintf("G %d->%d not found adding\n", nodej, nodek);
00200                 ghostAdjacencyTable[FEM_To_ghost_index(nodej)].push_back(nodeID);
00201               }
00202               {
00203                 //if(nodej==-5|| nodek==-5) CkPrintf("G %d->%d found already\n", nodej, nodek);
00204               }
00205             }
00206           }
00207         }
00208         else { // A non-ghost node, almost same as for ghost nodes
00209           for(int k=0;k<nodesPerElem;k++){
00210             const int nodek=conn[k];
00211             if(nodek != nodej){
00212               var_id nodeID = var_id::createNodeID(1,nodek);
00213               int idx = nodeAdjacency->findInRow(nodej,nodeID);
00214               if(idx == -1){
00215                 //if(nodej==-5|| nodek==-5) CkPrintf("NG %d->%d not found--adding\n", nodej, nodek);
00216                 adjacencyTable[nodej].push_back(nodeID);
00217               }
00218               {
00219                 //if(nodej==-5 || nodek==-5) CkPrintf("NG %d->%d found already\n", nodej, nodek);
00220               }
00221             }
00222           }
00223         }
00224       }
00225     }
00226   }
00227 
00228 
00229   for(int i=0;i<((FEM_Elem*)elem.getGhost())->size();i++) {        // for each element of the given type
00230     const int *conn = ((FEM_Elem*)elem.getGhost())->connFor(i);
00231     for(int j=0;j<nodesPerElem;j++){   // for each node adjacent to the element
00232       const int nodej = conn[j];
00233       if (nodej!=-1){
00234         if(FEM_Is_ghost_index(nodej)) { // A ghost node
00235           for(int k=0;k<nodesPerElem;k++){
00236             const int nodek=conn[k];
00237             if(nodek != nodej){
00238               var_id nodeID = var_id::createNodeID(1,nodek);
00239               int idx = ghostAdjacencyAttr->findInRow(FEM_To_ghost_index(nodej),nodeID);
00240               if(idx == -1){
00241                 //if(nodej==-5|| nodek==-5) CkPrintf("G-G %d->%d not found adding\n", nodej, nodek);
00242                 ghostAdjacencyTable[FEM_To_ghost_index(nodej)].push_back(nodeID);
00243               }
00244               {
00245                 //if(nodej==-5|| nodek==-5) CkPrintf("G-G %d->%d found already\n", nodej, nodek);
00246               }
00247             }
00248           }
00249         }
00250         else { // A non-ghost node, almost same as for ghost nodes
00251           for(int k=0;k<nodesPerElem;k++){
00252             const int nodek=conn[k];
00253             if(nodek != nodej){
00254               var_id nodeID = var_id::createNodeID(1,nodek);
00255               int idx = nodeAdjacency->findInRow(nodej,nodeID);
00256               if(idx == -1){
00257                 //if(nodej==-5|| nodek==-5) CkPrintf("G-NG %d->%d not found--adding\n", nodej, nodek);
00258                 adjacencyTable[nodej].push_back(nodeID);
00259               }
00260               {
00261                 //if (nodej==-5 || nodek==-5) CkPrintf("G-NG %d->%d found already\n", nodej, nodek);
00262               }
00263             }
00264           }
00265         }
00266       }
00267     }
00268   }
00269 
00270 }
00271 
00272 
00273 
00277 void FEM_Elem::allocateElemAdjacency(){
00278     
00279   if(elemAdjacency){
00280     CkPrintf("FEM> WARNING: Deleting previously allocated(?) elemAdjacency array. allocateElemAdjacency() should probably only be called once.\n");
00281     delete elemAdjacency;
00282   }
00283   if(elemAdjacencyTypes){
00284     CkPrintf("FEM> WARNING: Deleting previously allocated(?) elemAdjacencyTypes array. allocateElemAdjacency() should probably only be called once.\n");
00285     delete elemAdjacencyTypes;
00286   }
00287 
00288   elemAdjacency = new FEM_IndexAttribute(this,FEM_ELEM_ELEM_ADJACENCY);
00289   elemAdjacencyTypes = new FEM_IndexAttribute(this,FEM_ELEM_ELEM_ADJ_TYPES);
00290 
00291   elemAdjacency->setLength(size());
00292   elemAdjacency->setWidth(conn->getWidth());
00293   elemAdjacencyTypes->setLength(size());
00294   elemAdjacencyTypes->setWidth(conn->getWidth());
00295     
00296   add(elemAdjacency);
00297   add(elemAdjacencyTypes);
00298 
00299 }
00300 
00301 
00302 void FEM_Mesh::createNodeElemAdj(){
00303     node.lookup(FEM_NODE_ELEM_ADJACENCY,"FEM_Mesh::createNodeElemAdj");
00304     for(int i=0;i<elem.size();i++){
00305         node.setElemAdjacency(i,elem[i]);
00306     }
00307 }
00308 
00309 void FEM_Mesh::createNodeNodeAdj(){
00310     node.lookup(FEM_NODE_NODE_ADJACENCY,"FEM_Mesh::createNodeNodeAdj");
00311     for(int i=0;i<elem.size();i++){
00312         node.setNodeAdjacency(elem[i]);
00313     }
00314 }
00315 
00316 
00323 FEM_ElemAdj_Layer* FEM_Mesh::getElemAdjLayer(void) {
00324   if (! lastElemAdjLayer) { 
00325     lastElemAdjLayer=new FEM_ElemAdj_Layer();
00326     lastLayerSet = true;
00327   }
00328   return lastElemAdjLayer;
00329 }
00330 
00331 
00364 void FEM_Mesh::createElemElemAdj()
00365 {
00366     //CkPrintf("createElemElemAdj()\n");
00367     
00368   FEM_ElemAdj_Layer *g = getElemAdjLayer();
00369   if(! g->initialized)
00370     CkAbort("FEM> Cannot call FEM_Mesh_create_elem_elem_adjacency() before registering tuples with FEM_Add_elem2face_tuples()\n");
00371 
00372   
00373   int nodesPerTuple = g->nodesPerTuple;
00374   tupleTable table(nodesPerTuple);
00375 
00376   // Put tuples into table
00377   for (int t=0;t<elem.size();t++){ // for each element type
00378       if(elem.hasNonEmpty(t)) {
00379           const int tuplesPerElem = g->elem[t].tuplesPerElem;
00380           const int numElements = elem[t].size();
00381           // for every element of  type t:
00382           for (int elemIndex=0;elemIndex<numElements;elemIndex++)   {
00383               // insert element into the tuple table
00384               const int *conn=elem[t].connFor(elemIndex);
00385               int tuple[tupleTable::MAX_TUPLE];
00386               FEM_Symmetries_t allSym;
00387               // copy node numbers into tuple
00388               for (int u=0;u<tuplesPerElem;u++) {
00389                   for (int i=0;i<nodesPerTuple;i++) {
00390                       int eidx=g->elem[t].elem2tuple[i+u*g->nodesPerTuple];
00391                       if (eidx==-1)    // "not-there" node
00392                           tuple[i]=-1;   // Don't map via connectivity
00393                       else             // Ordinary node
00394                           tuple[i]=conn[eidx]; 
00395                   }
00396                   // add tuple to table
00397                   table.addTuple(tuple,new elemList(0,elemIndex,t,allSym,u)); 
00398               }
00399           }
00400 
00401           // Put corresponding ghost elements into tuple table
00402           if(elem[t].getGhost() != NULL){
00403               FEM_Elem *ghostElem = (FEM_Elem *)elem[t].getGhost();
00404               const int numElements = ghostElem->size();
00405               // for every element of  type t:
00406               for (int elemIndex=0;elemIndex<numElements;elemIndex++)   {
00407                   // insert element into the tuple table
00408                   const int *conn=ghostElem->connFor(elemIndex);
00409                   int tuple[tupleTable::MAX_TUPLE];
00410                   FEM_Symmetries_t allSym;
00411                   // copy node numbers into tuple
00412                   for (int u=0;u<tuplesPerElem;u++) {
00413                       for (int i=0;i<nodesPerTuple;i++) {
00414 
00415                           int eidx=g->elem[t].elem2tuple[i+u*g->nodesPerTuple];
00416                           if (eidx==-1) { //"not-there" node--
00417                               tuple[i]=-1; //Don't map via connectivity
00418                           } else { //Ordinary node
00419                               int n=conn[eidx];
00420                               tuple[i]=n; 
00421                           }
00422                       }
00423                       // add tuple to table
00424                       table.addTuple(tuple,new elemList(0,FEM_From_ghost_index(elemIndex),t,allSym,u)); 
00425                   }
00426               }
00427           }
00428       }
00429   }
00430 
00431   /* Extract adjacencies from table and store into 
00432    * FEM_ELEM_ELEM_ADJACENCY and FEM_ELEM_ELEM_ADJ_TYPES 
00433    * attribute fields
00434    */
00435   for (int t=0;t<elem.size();t++) { // for each element type t
00436       if (elem.hasNonEmpty(t)) {
00437           elemList *l;
00438           const int tuplesPerElem = g->elem[t].tuplesPerElem;
00439           const int numElements = elem[t].size();
00440           const int numGhostElements = ((FEM_Elem*)(elem[t].getGhost()))->size();
00441 
00442           // directly modify the element adjacency table for element type t
00443           FEM_IndexAttribute *elemAdjTypesAttr = (FEM_IndexAttribute *)elem[t].lookup(FEM_ELEM_ELEM_ADJ_TYPES,"createElemElemAdj");
00444           FEM_IndexAttribute *elemAdjAttr = (FEM_IndexAttribute *)elem[t].lookup(FEM_ELEM_ELEM_ADJACENCY,"createElemElemAdj");
00445           FEM_IndexAttribute *elemAdjTypesAttrGhost = (FEM_IndexAttribute *)elem[t].getGhost()->lookup(FEM_ELEM_ELEM_ADJ_TYPES,"createElemElemAdj");
00446           FEM_IndexAttribute *elemAdjAttrGhost = (FEM_IndexAttribute *)elem[t].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"createElemElemAdj");
00447           CkAssert(elemAdjTypesAttr && elemAdjAttr);
00448 
00449           AllocTable2d<int> &adjTable = elemAdjAttr->get();
00450           int *adjs = adjTable.getData();
00451           AllocTable2d<int> &adjTypesTable = elemAdjTypesAttr->get();
00452           int *adjTypes = adjTypesTable.getData();
00453 
00454           AllocTable2d<int> &adjTableGhost = elemAdjAttrGhost->get();
00455           int *adjsGhost = adjTableGhost.getData();
00456           AllocTable2d<int> &adjTypesTableGhost = elemAdjTypesAttrGhost->get();
00457           int *adjTypesGhost = adjTypesTableGhost.getData();
00458 
00459           // initialize tables
00460           for(int i=0;i<numElements*tuplesPerElem;i++){
00461               adjs[i]=-1;
00462               adjTypes[i]=-1;
00463           }
00464           for(int i=0;i<numGhostElements*tuplesPerElem;i++){
00465               adjsGhost[i]=-1;
00466               adjTypesGhost[i]=-1;
00467           }
00468 
00469           // look through each elemList that is returned by the tuple table
00470           table.beginLookup();
00471           while (NULL!=(l=table.lookupNext())) {
00472               if (l->next==NULL) { 
00473                   // One-entry list: must be a symmetry
00474                   // UNHANDLED CASE: not sure exactly what this means
00475               }
00476               else { /* Several elements in list: normal case */
00477                   //CkPrintf("Found in table list: ");
00478                   //for(const elemList *c=l;c!=NULL;c=c->next){
00479                                   //      CkPrintf("     %d,%d", c->type, c->localNo);
00480                   //}
00481                   //CkPrintf("\n");
00482                   
00483                   // for each a,b from the list
00484                   for (const elemList *a=l;a!=NULL;a=a->next){
00485                       for (const elemList *b=l;b!=NULL;b=b->next){
00486                           // if a and b are different elements
00487                           if((a->localNo != b->localNo) || (a->type != b->type)){
00488                               int j;
00489                               
00490                               if(a->type == t){ // only update the entries for element type t
00491                                   if(FEM_Is_ghost_index(a->localNo)){
00492                                       j = FEM_To_ghost_index(a->localNo)*tuplesPerElem + a->tupleNo;
00493                                       CkAssert(j<numGhostElements*tuplesPerElem);
00494                                       adjsGhost[j] = b->localNo;
00495                                       adjTypesGhost[j] = b->type;
00496                                   }
00497                                   else {
00498                                       j = a->localNo*tuplesPerElem + a->tupleNo;
00499                                       CkAssert(j<numElements*tuplesPerElem);
00500                                       //CkPrintf("Recording that %d,%d at position %d has neighbor %d,%d \n", a->type,a->localNo, a->tupleNo, b->type, b->localNo);
00501                                       adjs[j] = b->localNo;
00502                                       adjTypes[j] = b->type;
00503                                   }
00504                               }
00505                           }
00506                       }
00507                   }
00508               }
00509           }
00510       }
00511   }
00512   delete g;
00513 }
00514 
00515 
00516 
00521 void FEM_Mesh::e2e_getAll(int e, int *neighbors, int etype) 
00522 {
00523   if (e == -1) return; // non existent element
00524   FEM_IndexAttribute *eAdj;
00525   if(FEM_Is_ghost_index(e)){
00526     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getAll");
00527     AllocTable2d<int> &eAdjs = eAdj->get();
00528     for (int i=0; i<eAdjs.width(); i++) {
00529       neighbors[i] = eAdjs[FEM_To_ghost_index(e)][i];
00530     }
00531   }
00532   else {
00533     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getAll");
00534     AllocTable2d<int> &eAdjs = eAdj->get();
00535     for (int i=0; i<eAdjs.width(); i++) {
00536       neighbors[i] = eAdjs[e][i];
00537     }
00538   }
00539 }
00540 
00543 int FEM_Mesh::e2e_getNbr(int e, short idx, int etype) 
00544 {     
00545   if (e == -1) return -1;
00546   FEM_IndexAttribute *eAdj;
00547   if(FEM_Is_ghost_index(e)){
00548     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getNbr");
00549     AllocTable2d<int> &eAdjs = eAdj->get();
00550     return eAdjs[FEM_To_ghost_index(e)][idx];
00551   }
00552   else {
00553     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getNbr");
00554     AllocTable2d<int> &eAdjs = eAdj->get();
00555     return eAdjs[e][idx];
00556   }
00557 }
00558 
00559 
00560 ElemID FEM_Mesh::e2e_getElem(ElemID elem, int idx){
00561     return e2e_getElem(elem.id, idx, elem.type);
00562 }
00563 
00564 ElemID FEM_Mesh::e2e_getElem(ElemID *elem, int idx){
00565     return e2e_getElem(elem->id, idx, elem->type);
00566 }
00567 
00568 
00569 
00571 ElemID FEM_Mesh::e2e_getElem(int e, int idx, int etype){
00572   if (e == -1){
00573     ElemID ele(-1,-1);    
00574     return ele;
00575   }
00576   
00577   if(FEM_Is_ghost_index(e)){
00578 
00579       FEM_IndexAttribute *eAdj;
00580       FEM_IndexAttribute *eAdjType;
00581 
00582       eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getElem");
00583       AllocTable2d<int> &eAdjs = eAdj->get();
00584       eAdjType = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJ_TYPES,"e2e_getElem");
00585       AllocTable2d<int> &eAdjTypes = eAdjType->get();
00586 
00587       int t =  eAdjTypes[FEM_To_ghost_index(e)][idx];
00588       int id =  eAdjs[FEM_To_ghost_index(e)][idx];
00589       ElemID ele(t,id);
00590 
00591       return ele;
00592   }
00593   else {
00594 
00595       CkAssert(elem.has(etype) && e < elem[etype].size());
00596 
00597       FEM_IndexAttribute *eAdj;
00598       FEM_IndexAttribute *eAdjType;
00599 
00600       eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getElem");
00601       AllocTable2d<int> &eAdjs = eAdj->get();
00602       eAdjType = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJ_TYPES,"e2e_getElem");
00603       AllocTable2d<int> &eAdjTypes = eAdjType->get();
00604 
00605       int t =  eAdjTypes[e][idx];
00606       int id =  eAdjs[e][idx];
00607       ElemID ele(t,id);
00608 
00609       return ele;
00610   }
00611   
00612 }
00613 
00614 
00618 int FEM_Mesh::e2e_getIndex(int e, int nbr, int etype) 
00619 { 
00620   if (e == -1) return -1;
00621   FEM_IndexAttribute *eAdj;
00622   if(FEM_Is_ghost_index(e)){
00623     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getIndex");
00624     AllocTable2d<int> &eAdjs = eAdj->get();
00625     for (int i=0; i<eAdjs.width(); i++) {
00626       if (eAdjs[FEM_To_ghost_index(e)][i] == nbr) {
00627         return i;
00628       }
00629     }
00630   }
00631   else{
00632     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getIndex");
00633     AllocTable2d<int> &eAdjs = eAdj->get();
00634     for (int i=0; i<eAdjs.width(); i++) {
00635       if (eAdjs[e][i] == nbr) {
00636         return i;
00637       }
00638     }
00639   }
00640   return -1;
00641 }
00642 
00646 void FEM_Mesh::e2e_setAll(int e, int *neighbors, int etype) 
00647 {
00648   if (e == -1) return;
00649   FEM_IndexAttribute *eAdj;
00650   if(FEM_Is_ghost_index(e)){
00651     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_setAll");
00652     AllocTable2d<int> &eAdjs = eAdj->get();
00653     //CkPrintf("e2e_setAll: Setting element %d's neighbors to %d,%d,%d\n", e, 
00654     //neighbors[0], neighbors[1], neighbors[2]);
00655     for (int i=0; i<eAdjs.width(); i++) {
00656       eAdjs[FEM_To_ghost_index(e)][i] = neighbors[i];
00657     }
00658   }
00659   else{
00660     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_setAll");
00661     AllocTable2d<int> &eAdjs = eAdj->get();
00662     //CkPrintf("e2e_setAll: Setting element %d's neighbors to %d,%d,%d\n", e, 
00663     //neighbors[0], neighbors[1], neighbors[2]);
00664     for (int i=0; i<eAdjs.width(); i++) {
00665       eAdjs[e][i] = neighbors[i];
00666     }
00667   }
00668 }
00669 
00672 void FEM_Mesh::e2e_setIndex(int e, short idx, int newElem, int etype) 
00673 {
00674   if (e == -1) return;
00675   FEM_IndexAttribute *eAdj;
00676   if(FEM_Is_ghost_index(e)){
00677     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_setIndex");
00678     AllocTable2d<int> &eAdjs = eAdj->get();
00679     eAdjs[FEM_To_ghost_index(e)][idx] = newElem;
00680   }
00681   else {
00682     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_setIndex");    
00683     AllocTable2d<int> &eAdjs = eAdj->get();
00684     eAdjs[e][idx] = newElem;
00685   }
00686 }
00687 
00690 void FEM_Mesh::e2e_replace(int e, int oldNbr, int newNbr, int etype) 
00691 {
00692   if (e == -1) return;
00693   FEM_IndexAttribute *eAdj;
00694   if(FEM_Is_ghost_index(e)){
00695     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_replace");
00696     AllocTable2d<int> &eAdjs = eAdj->get();
00697     for (int i=0; i<eAdjs.width(); i++) {
00698       if (eAdjs[FEM_To_ghost_index(e)][i] == oldNbr) {
00699         eAdjs[FEM_To_ghost_index(e)][i] = newNbr;
00700         break;
00701       }
00702     }
00703   }
00704   else{
00705     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_replace");
00706     AllocTable2d<int> &eAdjs = eAdj->get();
00707     for (int i=0; i<eAdjs.width(); i++) {
00708       if (eAdjs[e][i] == oldNbr) {
00709         eAdjs[e][i] = newNbr;
00710         break;
00711       }
00712     }
00713   }
00714 }
00715 
00716 
00717 
00720 void FEM_Mesh::e2e_replace(ElemID e, ElemID oldNbr, ElemID newNbr) 
00721 {
00722   if (e.id == -1) return;
00723 
00724 //  CkPrintf("replacing %d,%d with %d,%d\n", oldNbr.type, oldNbr.id, newNbr.type, newNbr.id);
00725 
00726   
00727   FEM_IndexAttribute *eAdj;
00728   if(FEM_Is_ghost_index(e.id)){
00729       FEM_IndexAttribute *elemAdjTypesAttrGhost = (FEM_IndexAttribute *)elem[e.type].getGhost()->lookup(FEM_ELEM_ELEM_ADJ_TYPES,"e2e_replace");
00730       FEM_IndexAttribute *elemAdjAttrGhost = (FEM_IndexAttribute *)elem[e.type].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_replace");
00731    
00732       AllocTable2d<int> &eAdjs = elemAdjAttrGhost->get();
00733       AllocTable2d<int> &eAdjTypes = elemAdjTypesAttrGhost->get();
00734       
00735       for (int i=0; i<eAdjs.width(); i++) {
00736           if (eAdjs[FEM_To_ghost_index(e.id)][i] == oldNbr.id && eAdjTypes[FEM_To_ghost_index(e.id)][i] == oldNbr.type ) {
00737               eAdjs[FEM_To_ghost_index(e.id)][i] = newNbr.id;
00738               eAdjTypes[FEM_To_ghost_index(e.id)][i] = newNbr.type;
00739               break;
00740           }
00741       }
00742   }
00743   else{
00744     eAdj = (FEM_IndexAttribute *)elem[e.type].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_replace");
00745   
00746     FEM_IndexAttribute *elemAdjTypesAttr = (FEM_IndexAttribute *)elem[e.type].lookup(FEM_ELEM_ELEM_ADJ_TYPES,"e2e_replace");
00747     FEM_IndexAttribute *elemAdjAttr = (FEM_IndexAttribute *)elem[e.type].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_replace");
00748     
00749     AllocTable2d<int> &eAdjs = elemAdjAttr->get();
00750     AllocTable2d<int> &eAdjTypes = elemAdjTypesAttr->get();
00751 
00752     for (int i=0; i<eAdjs.width(); i++) {
00753         if (eAdjs[e.id][i] == oldNbr.id && eAdjTypes[e.id][i] == oldNbr.type ) {
00754 //          CkPrintf("eAdjTypes[e.id][i] was %d eAdjs[e.id][i] was %d\n", eAdjTypes[e.id][i], eAdjs[e.id][i]);
00755             eAdjs[e.id][i] = newNbr.id;
00756             eAdjTypes[e.id][i] = newNbr.type;
00757 //          CkPrintf("newNbr.type = %d\n", newNbr.type);
00758 //          CkPrintf("eAdjTypes[e.id][i] is %d eAdjs[e.id][i] is %d\n", eAdjTypes[e.id][i], eAdjs[e.id][i]);
00759             break;
00760         }
00761     }
00762     
00763   }
00764 }
00765 
00766 
00769 void FEM_Mesh::e2e_removeAll(int e, int etype)
00770 {
00771   if (e == -1) return;
00772   FEM_IndexAttribute *eAdj;
00773   if(FEM_Is_ghost_index(e)) {
00774     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_removeAll");
00775     AllocTable2d<int> &eAdjs = eAdj->get();
00776     for (int i=0; i<eAdjs.width(); i++) {
00777       eAdjs[FEM_To_ghost_index(e)][i] = -1;
00778     }
00779   }
00780   else {
00781     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_removeAll");
00782     AllocTable2d<int> &eAdjs = eAdj->get();
00783     for (int i=0; i<eAdjs.width(); i++) {
00784       eAdjs[e][i] = -1;
00785     }
00786   }
00787 }
00788 
00789 
00792 void FEM_Mesh::e2e_printAll(ElemID e) 
00793 {
00794     if (e.id == -1) return; // non existent element
00795     
00796     FEM_IndexAttribute *eAdj;
00797     FEM_IndexAttribute *eAdjType;
00798     
00799     if(FEM_Is_ghost_index(e.id)){
00800         eAdj = (FEM_IndexAttribute *)elem[e.type].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_printAll");
00801         eAdjType = (FEM_IndexAttribute *)elem[e.type].getGhost()->lookup(FEM_ELEM_ELEM_ADJ_TYPES,"e2e_printAll");
00802     }
00803     else {
00804         eAdj = (FEM_IndexAttribute *)elem[e.type].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_printAll");
00805         eAdjType = (FEM_IndexAttribute *)elem[e.type].lookup(FEM_ELEM_ELEM_ADJ_TYPES,"e2e_printAll");
00806     }
00807     
00808     
00809     AllocTable2d<int> &eAdjs = eAdj->get();
00810     AllocTable2d<int> &eAdjTypes = eAdjType->get();
00811     CkAssert(eAdjs.width() == eAdjTypes.width());
00812 
00813     CkAssert(e.getSignedId()>=0);
00814 
00815     for (int i=0; i<eAdjs.width(); i++) {
00816         CkPrintf("Element %d,%d is adjacent to %d,%d\n", e.type, e.id, eAdjTypes[e.getSignedId()][i], eAdjs[e.getSignedId()][i]);
00817 
00818     }
00819 
00820 }
00821 
00822 
00823 
00828 void FEM_Mesh::e2n_getAll(int e, int *adjnodes, int etype) 
00829 {
00830   if (e == -1) return;
00831   if(FEM_Is_ghost_index(e)){
00832     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00833     AllocTable2d<int> &conn = eConn->get();
00834     for (int i=0; i<conn.width(); i++) {
00835       adjnodes[i] = conn[FEM_To_ghost_index(e)][i];
00836     }
00837   }
00838   else{
00839     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00840     AllocTable2d<int> &conn = eConn->get();
00841     for (int i=0; i<conn.width(); i++) {
00842       adjnodes[i] = conn[e][i];
00843     }
00844   }
00845 }
00846 
00849 int FEM_Mesh::e2n_getNode(int e, short idx, int etype) 
00850 { 
00851   if (e == -1) return -1;
00852   if(FEM_Is_ghost_index(e)){
00853     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00854     AllocTable2d<int> &conn = eConn->get();
00855     return conn[FEM_To_ghost_index(e)][idx];
00856   }
00857   else{
00858     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00859     AllocTable2d<int> &conn = eConn->get();
00860     return conn[e][idx];
00861   }
00862 }
00863 
00867 short FEM_Mesh::e2n_getIndex(int e, int n, int etype) 
00868 {
00869   if (e == -1) return -1;
00870   if(FEM_Is_ghost_index(e)){
00871     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00872     AllocTable2d<int> &conn = eConn->get();
00873     for (int i=0; i<conn.width(); i++)
00874       if (conn[FEM_To_ghost_index(e)][i] == n) 
00875         return i;
00876   }
00877   else{
00878     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00879     AllocTable2d<int> &conn = eConn->get();
00880     for (int i=0; i<conn.width(); i++)
00881       if (conn[e][i] == n)
00882         return i;
00883   }
00884   return -1;
00885 }
00886 
00890 void FEM_Mesh::e2n_setAll(int e, int *adjnodes, int etype) 
00891 {
00892   if (e == -1) return;
00893   if(FEM_Is_ghost_index(e)){
00894     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00895     AllocTable2d<int> &conn = eConn->get();
00896     for (int i=0; i<conn.width(); i++) {
00897       conn[FEM_To_ghost_index(e)][i] = adjnodes[i];
00898     }
00899   }
00900   else{
00901     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00902     AllocTable2d<int> &conn = eConn->get();
00903     for (int i=0; i<conn.width(); i++) {
00904       conn[e][i] = adjnodes[i];
00905     }
00906   }  
00907 }
00908 
00911 void FEM_Mesh::e2n_setIndex(int e, short idx, int newNode, int etype) 
00912 {
00913   if (e == -1) return;
00914   if(FEM_Is_ghost_index(e)){
00915     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00916     AllocTable2d<int> &conn = eConn->get();
00917     conn[FEM_To_ghost_index(e)][idx] = newNode;
00918   }
00919   else{
00920     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00921     AllocTable2d<int> &conn = eConn->get();
00922     conn[e][idx] = newNode;
00923   }
00924 }
00925 
00928 void FEM_Mesh::e2n_replace(int e, int oldNode, int newNode, int etype) 
00929 {
00930   if (e == -1) return;
00931   if(FEM_Is_ghost_index(e)){
00932     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00933     AllocTable2d<int> &conn = eConn->get();
00934     for (int i=0; i<conn.width(); i++) {
00935       if (conn[FEM_To_ghost_index(e)][i] == oldNode) {
00936         conn[FEM_To_ghost_index(e)][i] = newNode;
00937         break;
00938       }
00939     }
00940   }
00941   else{
00942     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00943     AllocTable2d<int> &conn = eConn->get();
00944     for (int i=0; i<conn.width(); i++) {
00945       if (conn[e][i] == oldNode) {
00946         conn[e][i] = newNode;
00947         break;
00948       }
00949     }
00950   }  
00951 }
00952 
00953 void FEM_Mesh::e2n_removeAll(int e, int etype)
00954 {
00955   if (e == -1) return;
00956   if(FEM_Is_ghost_index(e)){
00957     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00958     AllocTable2d<int> &conn = eConn->get();
00959     for (int i=0; i<conn.width(); i++)
00960       conn[FEM_To_ghost_index(e)][i] = -1;
00961   }
00962   else{
00963     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00964     AllocTable2d<int> &conn = eConn->get();
00965     for (int i=0; i<conn.width(); i++)
00966       conn[e][i] = -1;
00967   }  
00968 }
00969 
00970 
00971 
00972 
00974 int FEM_Mesh::n2n_getLength(int n) {
00975   if (n == -1) {
00976     return 0;
00977   }
00978   FEM_VarIndexAttribute *eAdj;
00979   if(FEM_Is_ghost_index(n)){
00980     eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_getLength");
00981     n = FEM_To_ghost_index(n);
00982   }
00983   else {
00984     eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_getLength");
00985   }
00986   CkVec<CkVec<ElemID> > &eVec = eAdj->get();
00987   CkVec<ElemID> &nsVec = eVec[n];
00988   return nsVec.length();
00989 }
00990 
00994 void FEM_Mesh::n2n_getAll(int n, int *&adjnodes, int &sz) 
00995 {
00996   if (n == -1) {
00997     sz = 0;
00998     return;
00999   }
01000   if(FEM_Is_ghost_index(n)){
01001     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_getAll");
01002     CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01003     CkVec<ElemID> &nsVec = nVec[FEM_To_ghost_index(n)];
01004     sz = nsVec.length();
01005     if(sz > 0) adjnodes = new int[sz];
01006     for (int i=0; i<sz; i++) {
01007       adjnodes[i] = nsVec[i].getSignedId();
01008     }
01009   }
01010   else{
01011     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_getAll");
01012     CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01013     CkVec<ElemID> &nsVec = nVec[n];
01014     sz = nsVec.length();
01015     if(sz > 0) adjnodes = new int[sz];
01016     for (int i=0; i<sz; i++) {
01017       adjnodes[i] = nsVec[i].getSignedId();
01018     }
01019   }
01020 }
01021  
01024 void FEM_Mesh::n2n_add(int n, int newNode) 
01025 {
01026   if (n == -1) return;
01027   if(FEM_Is_ghost_index(n)){
01028     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_add");
01029     CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01030     ElemID nn(0, newNode);
01031     CkVec<ElemID> &nsVec = nVec[FEM_To_ghost_index(n)];
01032     nsVec.push_back(nn);
01033   }
01034   else{
01035     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_add");
01036     CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01037     ElemID nn(0, newNode);
01038     CkVec<ElemID> &nsVec = nVec[n];
01039     nsVec.push_back(nn);
01040   }
01041 }
01042 
01043 
01044 
01045 
01048 void FEM_Mesh::n2n_remove(int n, int oldNode) 
01049 {
01050   if (n == -1) return;
01051   if(FEM_Is_ghost_index(n)){
01052     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_remove");
01053     CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01054     CkVec<ElemID> &nsVec = nVec[FEM_To_ghost_index(n)];
01055     for (int i=0; i<nsVec.length(); i++) {
01056       if (nsVec[i].getSignedId() == oldNode) {
01057         nsVec.remove(i);
01058         break;
01059       }
01060     }
01061   }
01062   else {
01063     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_remove");
01064     CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01065     CkVec<ElemID> &nsVec = nVec[n];
01066     for (int i=0; i<nsVec.length(); i++) {
01067       if (nsVec[i].getSignedId() == oldNode) {
01068         nsVec.remove(i);
01069         break;
01070       }
01071     }
01072   }
01073 }
01074 
01077 int FEM_Mesh::n2n_exists(int n, int queryNode) 
01078 {
01079   if (n == -1) return 0;
01080   if(FEM_Is_ghost_index(n)){
01081     CkAssert(node.getGhost());
01082     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_exists");
01083     CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01084     CkVec<ElemID> &nsVec = nVec[FEM_To_ghost_index(n)];
01085     for (int i=0; i<nsVec.length(); i++)
01086       if (nsVec[i].getSignedId() == queryNode) 
01087         return 1;
01088   }
01089   else {
01090     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_exists");
01091     CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01092     CkVec<ElemID> &nsVec = nVec[n];
01093     for (int i=0; i<nsVec.length(); i++)
01094       if (nsVec[i].getSignedId() == queryNode) 
01095     return 1;
01096   }
01097   return 0;
01098 }
01099 
01102 void FEM_Mesh::n2n_replace(int n, int oldNode, int newNode) 
01103 {
01104   if (n == -1) return;
01105   if(FEM_Is_ghost_index(n)){
01106     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_replace");
01107     CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01108     CkVec<ElemID> &nsVec = nVec[FEM_To_ghost_index(n)];
01109     for (int i=0; i<nsVec.length(); i++) {
01110       if (nsVec[i].getSignedId() == oldNode) {
01111         nsVec[i] = ElemID(0,newNode);
01112         break;
01113       }
01114     }
01115   }
01116   else {
01117     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_replace");
01118     CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01119     CkVec<ElemID> &nsVec = nVec[n];
01120     for (int i=0; i<nsVec.length(); i++) {
01121       if (nsVec[i].getSignedId() == oldNode) {
01122         nsVec[i] = ElemID(0,newNode);
01123         break;
01124       }
01125     }
01126 
01127   }
01128 }
01129 
01132 void FEM_Mesh::n2n_removeAll(int n)
01133 {
01134   if (n == -1) return;
01135   if(FEM_Is_ghost_index(n)){
01136     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_removeAll");  
01137     CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01138     CkVec<ElemID> &nsVec = nVec[FEM_To_ghost_index(n)];
01139     nsVec.free();
01140   }
01141   else{
01142     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_removeAll");  
01143     CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01144     CkVec<ElemID> &nsVec = nVec[n];
01145     nsVec.free();
01146   }
01147 }
01148 
01150 int FEM_Mesh::n2e_getLength(int n) {
01151   if (n == -1) {
01152     return 0;
01153   }
01154   FEM_VarIndexAttribute *eAdj;
01155   if(FEM_Is_ghost_index(n)){
01156     eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_getLength");
01157     n = FEM_To_ghost_index(n);
01158   }
01159   else {
01160     eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_getLength");
01161   }
01162   CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01163   CkVec<ElemID> &nsVec = eVec[n];
01164   return nsVec.length();
01165 }
01166 
01170 ElemID FEM_Mesh::n2e_getElem(int n, int whichIdx){
01171 
01172   if(FEM_Is_ghost_index(n)){
01173     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_getElem");
01174     CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01175     CkVec<ElemID> &nsVec = eVec[FEM_To_ghost_index(n)];
01176     assert(whichIdx < nsVec.length());
01177     return  nsVec[whichIdx];
01178     
01179   }
01180   else {
01181     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_getAll");
01182     CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01183     CkVec<ElemID> &nsVec = eVec[n];
01184     assert(whichIdx < nsVec.length());
01185     return  nsVec[whichIdx];
01186   }
01187   
01188 }
01189 
01190 
01191 
01196 void FEM_Mesh::n2e_getAll(int n, int *&adjelements, int &sz) 
01197 {
01198   if (n == -1) {
01199     sz = 0;
01200     return;
01201   }
01202   if(FEM_Is_ghost_index(n)){
01203     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_getAll");  
01204     CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01205     CkVec<ElemID> &nsVec = eVec[FEM_To_ghost_index(n)];
01206     sz = nsVec.length();
01207     if(sz > 0) adjelements = new int[sz];
01208     for (int i=0; i<sz; i++) {
01209       adjelements[i] = nsVec[i].getSignedId();
01210     }
01211   }
01212   else {
01213     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_getAll");  
01214     CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01215     CkVec<ElemID> &nsVec = eVec[n];
01216     int len = nsVec.length();
01217     sz = len;
01218     if(sz > 0) adjelements = new int[sz];
01219     for (int i=0; i<sz; i++) {
01220       adjelements[i] = nsVec[i].getSignedId();
01221     }
01222   }
01223 }
01224  
01225 
01232 const CkVec<ElemID> &  FEM_Mesh::n2e_getAll(int n) 
01233 {
01234     assert(n!=-1); 
01235 
01236     if(FEM_Is_ghost_index(n)){
01237         FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_getAll");  
01238         CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01239         return eVec[FEM_To_ghost_index(n)];
01240     }
01241     else {
01242         FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_getAll");  
01243         CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01244         return eVec[n];
01245     }
01246 }
01247  
01248 
01249 
01252 void FEM_Mesh::n2e_add(int n, int newElem) 
01253 {
01254   if (n == -1) return;
01255   if(FEM_Is_ghost_index(n)){
01256     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_add");     
01257     CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01258     CkVec<ElemID> &nsVec = eVec[FEM_To_ghost_index(n)];
01259     ElemID ne(0, newElem);
01260     nsVec.push_back(ne);
01261     int *testn2e, testn2ec;
01262     n2e_getAll(n,testn2e,testn2ec);
01263     for(int i=0; i<testn2ec; i++) {
01264       if(FEM_Is_ghost_index(testn2e[i]))
01265         CkAssert(elem[0].ghost->is_valid(FEM_From_ghost_index(testn2e[i])));
01266       else 
01267         CkAssert(elem[0].is_valid(testn2e[i]));
01268     }
01269     if(testn2ec!=0) delete[] testn2e;
01270   }
01271   else {
01272     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_add");     
01273     CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01274     CkVec<ElemID> &nsVec = eVec[n];
01275     ElemID ne(0, newElem);
01276     nsVec.push_back(ne);
01277   }
01278 }
01279 
01282 void FEM_Mesh::n2e_remove(int n, int oldElem) 
01283 {
01284   if (n == -1) return;
01285   if(FEM_Is_ghost_index(n)){
01286     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_remove");
01287     CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01288     CkVec<ElemID> &nsVec = eVec[FEM_To_ghost_index(n)];
01289     for (int i=0; i<nsVec.length(); i++) {
01290       if (nsVec[i].getSignedId() == oldElem) {
01291         nsVec.remove(i);
01292         break;
01293       }
01294     }
01295   }
01296   else {
01297     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_remove");
01298     CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01299     CkVec<ElemID> &nsVec = eVec[n];
01300     for (int i=0; i<nsVec.length(); i++) {
01301       if (nsVec[i].getSignedId() == oldElem) {
01302         nsVec.remove(i);
01303         break;
01304       }
01305     }
01306   }
01307 }
01308  
01311 void FEM_Mesh::n2e_replace(int n, int oldElem, int newElem) 
01312 {
01313   if (n == -1) return;
01314   if(FEM_Is_ghost_index(n)){
01315     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_replace");
01316     CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01317     CkVec<ElemID> &nsVec = eVec[FEM_To_ghost_index(n)];
01318     for (int i=0; i<nsVec.length(); i++) {
01319       if (nsVec[i].getSignedId() == oldElem) {
01320         nsVec[i] = ElemID(0,newElem);
01321         break;
01322       }
01323     }
01324     int *testn2e, testn2ec;
01325     n2e_getAll(n,testn2e,testn2ec);
01326     for(int i=0; i<testn2ec; i++) {
01327       if(FEM_Is_ghost_index(testn2e[i]))
01328         CkAssert(elem[0].ghost->is_valid(FEM_From_ghost_index(testn2e[i])));
01329       else 
01330         CkAssert(elem[0].is_valid(testn2e[i]));
01331     }
01332     if(testn2ec!=0) delete[] testn2e;
01333   }
01334   else{
01335     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_replace");
01336     CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01337     CkVec<ElemID> &nsVec = eVec[n];
01338     for (int i=0; i<nsVec.length(); i++) {
01339       if (nsVec[i].getSignedId() == oldElem) {
01340         nsVec[i] = ElemID(0,newElem);
01341         break;
01342       }
01343     }
01344   }
01345 }
01346 
01349 void FEM_Mesh::n2e_removeAll(int n)
01350 {
01351   if (n == -1) return;
01352   if(FEM_Is_ghost_index(n)){
01353     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_removeAll");
01354     CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01355     CkVec<ElemID> &nsVec = eVec[FEM_To_ghost_index(n)];
01356     nsVec.free();
01357   }
01358   else{
01359     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_removeAll");
01360     CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01361     CkVec<ElemID> &nsVec = eVec[n];
01362     nsVec.free();
01363   }
01364 }
01365 
01374 int FEM_Mesh::getElementOnEdge(int n1, int n2) 
01375 {
01376   int *n1AdjElems, *n2AdjElems;
01377   int n1NumElems, n2NumElems;
01378   n2e_getAll(n1, n1AdjElems, n1NumElems);
01379   n2e_getAll(n2, n2AdjElems, n2NumElems);
01380   int ret = -1;
01381   //CkPrintf("%d has %d neighboring elements, %d has %d\n", n1, n1NumElems, n2, n2NumElems);
01382 
01383   bool flag = false;
01384   for (int i=0; i<n1NumElems; i++) {
01385     for (int j=0; j<n2NumElems; j++) {
01386       if (n1AdjElems[i] == n2AdjElems[j]) {
01387         if(n1AdjElems[i] >= 0) {
01388           ret = n1AdjElems[i];
01389       flag = true;
01390           break;
01391         }
01392         else {
01393           ret = n1AdjElems[i];
01394         }
01395       }
01396     }
01397     if(flag) break;
01398   }
01399   delete[] n1AdjElems;
01400   delete[] n2AdjElems;
01401   return ret; //preferably return a local element, otherwise return a ghost 
01402 }
01403 
01404 
01409 void FEM_Mesh::get2ElementsOnEdge(int n1, int n2, int *result_e1, int *result_e2) 
01410 {
01411   int *n1AdjElems=0, *n2AdjElems=0;
01412   int n1NumElems, n2NumElems;
01413 
01414   if(n1==n2){
01415     CkPrintf("ERROR: You called get2ElementsOnEdge() with two identical nodes %d, and %d \n", n1, n2);
01416     CkExit();
01417   }
01418 
01419   n2e_getAll(n1, n1AdjElems, n1NumElems);
01420   n2e_getAll(n2, n2AdjElems, n2NumElems);
01421   CkAssert(n1AdjElems!=0);
01422   CkAssert(n2AdjElems!=0);
01423   int found=0;
01424 
01425   *result_e1=-1;
01426   *result_e2=-1;
01427 
01428   for (int i=0; i<n1NumElems; i++) {
01429     for (int j=0; j<n2NumElems; j++) {
01430       if (n1AdjElems[i] == n2AdjElems[j]) {
01431         if(found==0){
01432           //          CkPrintf("found element1 %d\n", n1AdjElems[i]);
01433           *result_e1 = n1AdjElems[i];
01434           found++;
01435         }
01436         else if(found==1){
01437           //   CkPrintf("found element2 %d\n", n1AdjElems[i]);
01438           *result_e2 = n1AdjElems[i];
01439           found++;
01440         }
01441         else {
01442           CkPrintf("ERROR: Found a third element(%d) on edge %d,%d \n", n1AdjElems[i], n1, n2);
01443            CkExit();
01444         }
01445       }
01446     }
01447   }
01448   delete[] n1AdjElems;
01449   delete[] n2AdjElems;
01450 }
01451 
01455 int getID(FEM_Mesh *mesh, int elem){
01456     int id = 0;
01457 
01458     if(elem<0){
01459       AllocTable2d<int> &idTable = ((FEM_DataAttribute *)mesh->elem[0].getGhost()->lookup(FEM_DATA+3,"getID"))->getInt();
01460       id = idTable(0,FEM_To_ghost_index(elem));
01461     } else {
01462       AllocTable2d<int> &idTable = ((FEM_DataAttribute *)mesh->elem[0].lookup(FEM_DATA+3,"getID"))->getInt();
01463       id = idTable(0,elem);
01464     }
01465     
01466   return id;
01467 }
01468 
01469 
01477 void FEM_Mesh::get2ElementsOnEdgeSorted(int n1, int n2, int *result_e1, int *result_e2)
01478 {
01479   int *n1AdjElems=0, *n2AdjElems=0;
01480   int n1NumElems, n2NumElems;
01481   
01482   get2ElementsOnEdge(n1, n2, result_e1, result_e2);
01483 
01484   //   Swap the two elements if the first has greater id than the second 
01485     if(getID(this,*result_e1) > getID(this,*result_e2)){
01486     int temp = *result_e1;
01487     *result_e1 = *result_e2;
01488     *result_e2 = temp;
01489   }
01490 }
01491 
01492 
01493 
01494 
01496 int FEM_Mesh::countElementsOnEdge(int n1, int n2) {
01497   if (n1==n2) {
01498     CkPrintf("ERROR: You called countElementsOnEdge() with two identical nodes %d, and %d \n", n1, n2);
01499     CkExit();
01500   }
01501   
01502   int *n1AdjElems=0, *n2AdjElems=0;
01503   int n1NumElems, n2NumElems;
01504 
01505   CkAssert(node.is_valid_any_idx(n1));
01506   CkAssert(node.is_valid_any_idx(n2));
01507 
01508   n2e_getAll(n1, n1AdjElems, n1NumElems);
01509   n2e_getAll(n2, n2AdjElems, n2NumElems);
01510   CkAssert(n1AdjElems!=0);
01511   CkAssert(n2AdjElems!=0);
01512   int count=0;
01513 
01514 
01515   for (int i=0; i<n1NumElems; i++) {
01516     for (int j=0; j<n2NumElems; j++) {
01517       if (n1AdjElems[i] == n2AdjElems[j]) {
01518         count++;
01519       }
01520     }
01521   }
01522   delete[] n1AdjElems;
01523   delete[] n2AdjElems;
01524 
01525   return count;
01526 }
01527 
01528
libs/ck-libs/ParFUM/mesh_adjacency.C
