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 "fem.h"
00013 #include "fem_impl.h"
00014 #include "charm-api.h" /*for CLINKAGE, FTN_NAME*/
00015 
00016 
00017 
00018 //#define DEBUG
00019 
00020 
00021 #ifdef DEBUG
00022 FORTRAN_AS_C(CMIMEMORYCHECK,
00023              CmiMemoryCheck,
00024              cmimemorycheck, 
00025              (void),  () )
00026 #endif
00027 
00028 
00029 
00030 
00031 CLINKAGE void
00032 FEM_Mesh_create_node_elem_adjacency(int fem_mesh){
00033     const char *caller="FEM_Mesh_create_node_elem_adjacency"; FEMAPI(caller);
00034     FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,caller);
00035     m->createNodeElemAdj();
00036 }
00037 FORTRAN_AS_C(FEM_MESH_CREATE_NODE_ELEM_ADJACENCY,
00038              FEM_Mesh_create_node_elem_adjacency,
00039              fem_mesh_create_node_elem_adjacency, 
00040              (int *fem_mesh),  (*fem_mesh) )
00041 
00042 
00043 CLINKAGE void
00044 FEM_Mesh_create_node_node_adjacency(int fem_mesh){
00045     const char *caller="FEM_Mesh_create_node_node_adjacency"; FEMAPI(caller);
00046     FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,caller);
00047     m->createNodeNodeAdj();
00048 }
00049 FORTRAN_AS_C(FEM_MESH_CREATE_NODE_NODE_ADJACENCY,
00050              FEM_Mesh_create_node_node_adjacency,
00051              fem_mesh_create_node_node_adjacency, 
00052              (int *fem_mesh),  (*fem_mesh) )
00053 
00054 
00055 CLINKAGE void
00056 FEM_Mesh_create_elem_elem_adjacency(int fem_mesh){
00057     const char *caller="FEM_Mesh_create_elem_elem_adjacency"; FEMAPI(caller);
00058     FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,caller);
00059     m->createElemElemAdj();
00060 }
00061 FORTRAN_AS_C(FEM_MESH_CREATE_ELEM_ELEM_ADJACENCY,
00062              FEM_Mesh_create_elem_elem_adjacency,
00063              fem_mesh_create_elem_elem_adjacency, 
00064              (int *fem_mesh),  (*fem_mesh) )
00065 
00066 
00067 
00068 CLINKAGE void
00069 FEM_Mesh_create_elem_node_adjacency(int fem_mesh){
00070   CkPrintf("WARNING: Do Not Call FEM_Mesh_create_elem_node_adjacency(), as the connectivity table should already exist\n");
00071 }
00072 
00073 
00074 CLINKAGE void
00075 FEM_Mesh_get2ElementsOnEdge(int fem_mesh, int n1, int n2, int *e1, int *e2){
00076     const char *caller="FEM_Mesh_get2ElementsOnEdge"; FEMAPI(caller);
00077     FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,caller);
00078     m->get2ElementsOnEdge(n1, n2, e1, e2);
00079 }
00080 FORTRAN_AS_C(FEM_MESH_GET2ELEMENTSONEDGE,
00081              FEM_Mesh_get2ElementsOnEdge,
00082              fem_mesh_get2elementsonedge, 
00083              (int *fem_mesh, int *n1, int *n2, int *e1, int *e2),  
00084              (*fem_mesh,*n1,*n2,e1,e2) )
00085 
00086 
00087 
00088 void FEM_Node::allocateElemAdjacency(){
00089 #ifdef DEBUG
00090   CmiMemoryCheck();
00091 #endif
00092     if(elemAdjacency){
00093         delete elemAdjacency;
00094     }
00095     elemAdjacency = new FEM_VarIndexAttribute(this,FEM_NODE_ELEM_ADJACENCY);
00096     add(elemAdjacency);
00097 }
00098 
00099 void FEM_Node::allocateNodeAdjacency(){
00100 #ifdef DEBUG
00101   CmiMemoryCheck();
00102 #endif
00103     if(nodeAdjacency){
00104         delete nodeAdjacency;
00105     }
00106     nodeAdjacency = new FEM_VarIndexAttribute(this,FEM_NODE_NODE_ADJACENCY);
00107     add(nodeAdjacency);
00108 }
00109 
00110 
00111 //  Fill the node to element adjacency table for both this element and its corresponding ghosts
00112 void FEM_Node::setElemAdjacency(int type, const FEM_Elem &elem){
00113     int nodesPerElem = elem.getNodesPer();
00114     FEM_VarIndexAttribute *adjacencyAttr = elemAdjacency;
00115     CkVec<CkVec<var_id> > &adjacencyTable = elemAdjacency->get();
00116     FEM_VarIndexAttribute *ghostAdjacencyAttr = ((FEM_Node *)getGhost())->elemAdjacency;
00117     CkVec<CkVec<var_id> > &ghostAdjacencyTable = ghostAdjacencyAttr->get();
00118 
00119     // Scan through elements
00120     for(int i=0;i<elem.size();i++){
00121       const int *conn = elem.connFor(i);
00122       for(int j=0;j<nodesPerElem;j++){
00123         int node = conn[j];
00124         if (node!=-1){
00125           if(FEM_Is_ghost_index(node)){
00126             int idx = ghostAdjacencyAttr->findInRow(FEM_To_ghost_index(node),var_id(type,i));
00127             if(idx == -1) {// If not currently in the adjacency list, push onto list
00128               ghostAdjacencyTable[FEM_To_ghost_index(node)].push_back(var_id(type,i));
00129             }
00130           }
00131           else{
00132             int idx = adjacencyAttr->findInRow(node,var_id(type,i));
00133             if(idx == -1) {// If not currently in the adjacency list, push onto list
00134               adjacencyTable[node].push_back(var_id(type,i));
00135             }
00136           }
00137         }
00138       }
00139     }
00140     
00141     // Scan through ghost elements
00142     if(elem.getGhost()){
00143       for(int i=0;i<((FEM_Elem*)elem.getGhost())->size();i++){
00144         const int *conn = ((FEM_Elem*)elem.getGhost())->connFor(i);
00145         for(int j=0;j<nodesPerElem;j++){
00146           int node = conn[j];
00147           if (node!=-1){
00148             if(FEM_Is_ghost_index(node)){
00149               int idx = ghostAdjacencyAttr->findInRow(FEM_To_ghost_index(node),var_id(type,FEM_From_ghost_index(i)));
00150               if(idx == -1){ // If not currently in the adjacency list, push onto list
00151                 ghostAdjacencyTable[FEM_To_ghost_index(node)].push_back(var_id(type,FEM_From_ghost_index(i)));
00152               }
00153             }
00154             
00155             else{
00156               int idx = adjacencyAttr->findInRow(node,var_id(type,FEM_From_ghost_index(i)));
00157               if(idx == -1){// If not currently in the adjacency list, push onto list
00158                 adjacencyTable[node].push_back(var_id(type,FEM_From_ghost_index(i)));
00159               }
00160             }
00161           }
00162         }
00163       }
00164     }
00165 }
00166 
00167 
00168 //  Populate the entire node to node adjacency table
00169 //  Two nodes are considered adjacent if they both are in the connectivity table for a common element.
00170 //  This choice for definition of adjacent nodes does not take into account what are edges 
00171 //  of the element, but it does simplify the computation. It will work fine for 
00172 //  triangles and tetrahedra, but may not make as much sense for more complicated
00173 //  element types where all nodes are not directly connected by edges.
00174 void FEM_Node::setNodeAdjacency(const FEM_Elem &elem){
00175   //CkPrintf("In FEM_Node::setNodeAdjacency()\n");
00176   int nodesPerElem = elem.getNodesPer();
00177   CkVec<CkVec<var_id> > &adjacencyTable = nodeAdjacency->get();
00178   FEM_VarIndexAttribute *ghostAdjacencyAttr = ((FEM_Node *)getGhost())->nodeAdjacency;
00179   CkVec<CkVec<var_id> > &ghostAdjacencyTable = ghostAdjacencyAttr->get();
00180   
00181 #ifdef DEBUG
00182   CmiMemoryCheck();
00183 #endif
00184 
00185   // Add the adjacencies defined by the non-ghost elements
00186   for(int i=0;i<elem.size();i++) {        // for each element of the given type
00187     const int *conn = elem.connFor(i);
00188     for(int j=0;j<nodesPerElem;j++){   // for each node adjacent to the element
00189       const int nodej = conn[j];
00190       if (nodej!=-1){
00191         if(FEM_Is_ghost_index(nodej)) { // A ghost node
00192           for(int k=0;k<nodesPerElem;k++){
00193             const int nodek=conn[k];
00194             if(nodek != nodej){
00195               var_id nodeID = var_id::createNodeID(1,nodek);
00196               int idx = ghostAdjacencyAttr->findInRow(FEM_To_ghost_index(nodej),nodeID);
00197               if(idx == -1){
00198                 //if(nodej==-5|| nodek==-5) CkPrintf("G %d->%d not found adding\n", nodej, nodek);
00199                 ghostAdjacencyTable[FEM_To_ghost_index(nodej)].push_back(nodeID);
00200               }
00201               {
00202                 //if(nodej==-5|| nodek==-5) CkPrintf("G %d->%d found already\n", nodej, nodek);
00203               }
00204             }
00205           }
00206         }
00207         else { // A non-ghost node, almost same as for ghost nodes
00208           for(int k=0;k<nodesPerElem;k++){
00209             const int nodek=conn[k];
00210             if(nodek != nodej){
00211               var_id nodeID = var_id::createNodeID(1,nodek);
00212               int idx = nodeAdjacency->findInRow(nodej,nodeID);
00213               if(idx == -1){
00214                 //if(nodej==-5|| nodek==-5) CkPrintf("NG %d->%d not found--adding\n", nodej, nodek);
00215                 adjacencyTable[nodej].push_back(nodeID);
00216               }
00217               {
00218                 //if(nodej==-5 || nodek==-5) CkPrintf("NG %d->%d found already\n", nodej, nodek);
00219               }
00220             }
00221           }
00222         }
00223       }
00224     }
00225   }
00226 
00227 
00228   for(int i=0;i<((FEM_Elem*)elem.getGhost())->size();i++) {        // for each element of the given type
00229     const int *conn = ((FEM_Elem*)elem.getGhost())->connFor(i);
00230     for(int j=0;j<nodesPerElem;j++){   // for each node adjacent to the element
00231       const int nodej = conn[j];
00232       if (nodej!=-1){
00233         if(FEM_Is_ghost_index(nodej)) { // A ghost node
00234           for(int k=0;k<nodesPerElem;k++){
00235             const int nodek=conn[k];
00236             if(nodek != nodej){
00237               var_id nodeID = var_id::createNodeID(1,nodek);
00238               int idx = ghostAdjacencyAttr->findInRow(FEM_To_ghost_index(nodej),nodeID);
00239               if(idx == -1){
00240                 //if(nodej==-5|| nodek==-5) CkPrintf("G-G %d->%d not found adding\n", nodej, nodek);
00241                 ghostAdjacencyTable[FEM_To_ghost_index(nodej)].push_back(nodeID);
00242               }
00243               {
00244                 //if(nodej==-5|| nodek==-5) CkPrintf("G-G %d->%d found already\n", nodej, nodek);
00245               }
00246             }
00247           }
00248         }
00249         else { // A non-ghost node, almost same as for ghost nodes
00250           for(int k=0;k<nodesPerElem;k++){
00251             const int nodek=conn[k];
00252             if(nodek != nodej){
00253               var_id nodeID = var_id::createNodeID(1,nodek);
00254               int idx = nodeAdjacency->findInRow(nodej,nodeID);
00255               if(idx == -1){
00256                 //if(nodej==-5|| nodek==-5) CkPrintf("G-NG %d->%d not found--adding\n", nodej, nodek);
00257                 adjacencyTable[nodej].push_back(nodeID);
00258               }
00259               {
00260                 //if (nodej==-5 || nodek==-5) CkPrintf("G-NG %d->%d found already\n", nodej, nodek);
00261               }
00262             }
00263           }
00264         }
00265       }
00266     }
00267   }
00268 
00269 
00270 #ifdef DEBUG
00271   CmiMemoryCheck();
00272 #endif
00273 
00274 
00275 }
00276 
00277 
00278 
00279 // Allocate both the FEM_ELEM_ELEM_ADJACENCY and FEM_ELEM_ELEM_ADJ_TYPES attributes
00280 // The values in these two attributes will be generated by 
00281 void FEM_Elem::allocateElemAdjacency(){
00282 
00283 #ifdef DEBUG
00284   CmiMemoryCheck();
00285 #endif
00286 
00287   if(elemAdjacency){
00288     CkPrintf("FEM> WARNING: Deleting previously allocated(?) elemAdjacency array. allocateElemAdjacency() should probably only be called once.\n");
00289     delete elemAdjacency;
00290   }
00291   if(elemAdjacencyTypes){
00292     CkPrintf("FEM> WARNING: Deleting previously allocated(?) elemAdjacencyTypes array. allocateElemAdjacency() should probably only be called once.\n");
00293     delete elemAdjacencyTypes;
00294   }
00295 
00296   elemAdjacency = new FEM_IndexAttribute(this,FEM_ELEM_ELEM_ADJACENCY);
00297   elemAdjacencyTypes = new FEM_IndexAttribute(this,FEM_ELEM_ELEM_ADJ_TYPES);
00298 
00299   elemAdjacency->setLength(size());
00300   elemAdjacency->setWidth(conn->getWidth());
00301   elemAdjacencyTypes->setLength(size());
00302   elemAdjacencyTypes->setWidth(conn->getWidth());
00303     
00304   add(elemAdjacency);
00305   add(elemAdjacencyTypes);
00306 
00307 #ifdef DEBUG
00308   CmiMemoryCheck();
00309 #endif
00310 
00311 }
00312 
00313 
00314 void FEM_Mesh::createNodeElemAdj(){
00315     node.lookup(FEM_NODE_ELEM_ADJACENCY,"FEM_Mesh::createElemNodeAdj");
00316     for(int i=0;i<elem.size();i++){
00317         node.setElemAdjacency(i,elem[i]);
00318     }
00319 #ifdef DEBUG
00320   CmiMemoryCheck();
00321 #endif
00322 }
00323 
00324 void FEM_Mesh::createNodeNodeAdj(){
00325     node.lookup(FEM_NODE_NODE_ADJACENCY,"FEM_Mesh::createNodeNodeAdj");
00326     for(int i=0;i<elem.size();i++){
00327         node.setNodeAdjacency(elem[i]);
00328     }
00329 #ifdef DEBUG
00330   CmiMemoryCheck();
00331 #endif
00332 }
00333 
00334 
00335 /* A created on demand data structure that stores
00336  * the number of tuples(or faces) for each type of element
00337  * along with the number of nodes in each tuple.
00338  * The number of nodes in each tuple depends upon how many
00339  * nodes are in a shared face(an edge will have 2, a triangle 3)
00340  */
00341 FEM_ElemAdj_Layer* FEM_Mesh::getElemAdjLayer(void) {
00342   if (! lastElemAdjLayer) 
00343     lastElemAdjLayer=new FEM_ElemAdj_Layer();
00344   return lastElemAdjLayer;
00345 }
00346 
00347 
00348 /* createElemElemAdj() is similar to splitter::addTuple()
00349  * It will scan through all tuples/faces for each element
00350  * and register them with a tupleTable which will find 
00351  * matching tuples between different elements. The matches  
00352  * are then extracted from the table and the resulting
00353  * adjacencies are marked in the FEM_ELEM_ELEM_ADJACENCY
00354  * and FEM_ELEM_ELEM_ADJ_TYPES attribute fields. The attributes
00355  * store the adjacent element's local id and type respectively.
00356  * The ordering in the attribute tables is by local element id,
00357  * as usual, and is then ordered by the tuple ordering specified
00358  * in the last parameter to FEM_Add_elem2face_tuples, which 
00359  * gets stored in the FEM_ElemAdj_Layer data structure
00360  * created/returned by getElemAdjLayer().
00361  *
00362  * The element id's are NOT indices into the conn array in 
00363  * the case of ghost elements. Ghost elements will have the
00364  * same element type as their corresponding real elements.
00365  * Their id, which gets stored in FEM_ELEM_ELEM_ADJACENCY
00366  * will be a negative number which can be converted to an index
00367  * with FEM_To_ghost_index(). Thus the user MUST use 
00368  * FEM_Is_ghost_index(i) on the values, before accessing 
00369  * them in the conn array, especially since the ghosts
00370  * will have a negative id.
00371  *
00372  * The function assumes that the FEM_ELEM_ELEM_ADJACENCY
00373  * and FEM_ELEM_ELEM_ADJ_TYPES attribute fields already exist.
00374  *
00375  * TODO:
00376  *  
00377  *   Verify the tuple table does not need to be 
00378  *   explicitly deleted.
00379  */
00380 void FEM_Mesh::createElemElemAdj()
00381 {
00382   FEM_ElemAdj_Layer *g = getElemAdjLayer();
00383   if(! g->initialized)
00384     CkAbort("FEM> Cannot call FEM_Mesh_create_elem_elem_adjacency() before registering tuples with FEM_Add_elem2face_tuples()\n");
00385 
00386   const int nodesPerTuple = g->nodesPerTuple;
00387   tupleTable table(nodesPerTuple);
00388 
00389 #ifdef DEBUG
00390   CmiMemoryCheck();
00391 #endif
00392 
00393   // Put tuples into table
00394   for (int t=0;t<elem.size();t++){ // for each element type
00395       if(elem.has(t)) {
00396           const int tuplesPerElem = g->elem[t].tuplesPerElem;
00397           const int numElements = elem[t].size();
00398           // for every element of  type t:
00399           for (int elemIndex=0;elemIndex<numElements;elemIndex++)   {
00400               // insert element into the tuple table
00401               const int *conn=elem[t].connFor(elemIndex);
00402               int tuple[tupleTable::MAX_TUPLE];
00403               FEM_Symmetries_t allSym;
00404               // copy node numbers into tuple
00405               for (int u=0;u<tuplesPerElem;u++) {
00406                   for (int i=0;i<nodesPerTuple;i++) {
00407                     int eidx=g->elem[t].elem2tuple[i+u*g->nodesPerTuple];
00408                     if (eidx==-1)    // "not-there" node
00409                       tuple[i]=-1;   // Don't map via connectivity
00410                     else             // Ordinary node
00411                       tuple[i]=conn[eidx]; 
00412                   }
00413                   // add tuple to table
00414                   table.addTuple(tuple,new elemList(0,elemIndex,t,allSym,u)); 
00415               }
00416           }
00417       
00418           // Put corresponding ghost elements into tuple table
00419           if(elem[t].getGhost() != NULL){
00420               FEM_Elem *ghostElem = (FEM_Elem *)elem[t].getGhost();
00421               const int numElements = ghostElem->size();
00422               // for every element of  type t:
00423               for (int elemIndex=0;elemIndex<numElements;elemIndex++)   {
00424                   // insert element into the tuple table
00425                   const int *conn=ghostElem->connFor(elemIndex);
00426                   int tuple[tupleTable::MAX_TUPLE];
00427                   FEM_Symmetries_t allSym;
00428                   // copy node numbers into tuple
00429                   for (int u=0;u<tuplesPerElem;u++) {
00430                       for (int i=0;i<nodesPerTuple;i++) {
00431                        
00432                           int eidx=g->elem[t].elem2tuple[i+u*g->nodesPerTuple];
00433                           if (eidx==-1) { //"not-there" node--
00434                               tuple[i]=-1; //Don't map via connectivity
00435                           } else { //Ordinary node
00436                               int n=conn[eidx];
00437                               tuple[i]=n; 
00438                           }
00439                       }
00440                       // add tuple to table
00441                       table.addTuple(tuple,new elemList(0,FEM_From_ghost_index(elemIndex),t,allSym,u)); 
00442                   }
00443               }
00444           }
00445           
00446       }
00447   }
00448   
00449 
00450 #ifdef DEBUG
00451   CmiMemoryCheck();
00452 #endif
00453 
00454 
00455   /* Extract adjacencies from table and store into 
00456    * FEM_ELEM_ELEM_ADJACENCY and FEM_ELEM_ELEM_ADJ_TYPES 
00457    * attribute fields
00458    */
00459   for (int t=0;t<elem.size();t++) // for each element type t
00460       if (elem.has(t)) {
00461         elemList *l;
00462         const int tuplesPerElem = g->elem[t].tuplesPerElem;
00463         const int numElements = elem[t].size();
00464         const int numGhostElements = ((FEM_Elem*)(elem[t].getGhost()))->size();
00465         table.beginLookup();
00466         
00467         // directly modify the element adjacency table for element type t
00468         FEM_IndexAttribute *elemAdjTypesAttr = (FEM_IndexAttribute *)elem[t].lookup(FEM_ELEM_ELEM_ADJ_TYPES,"createElemElemAdj");
00469         FEM_IndexAttribute *elemAdjAttr = (FEM_IndexAttribute *)elem[t].lookup(FEM_ELEM_ELEM_ADJACENCY,"createElemElemAdj");
00470         FEM_IndexAttribute *elemAdjTypesAttrGhost = (FEM_IndexAttribute *)elem[t].getGhost()->lookup(FEM_ELEM_ELEM_ADJ_TYPES,"createElemElemAdj");
00471         FEM_IndexAttribute *elemAdjAttrGhost = (FEM_IndexAttribute *)elem[t].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"createElemElemAdj");
00472         CkAssert(elemAdjTypesAttr && elemAdjAttr);
00473         
00474         AllocTable2d<int> &adjTable = elemAdjAttr->get();
00475         int *adjs = adjTable.getData();
00476         AllocTable2d<int> &adjTypesTable = elemAdjTypesAttr->get();
00477         int *adjTypes = adjTypesTable.getData();
00478         
00479         AllocTable2d<int> &adjTableGhost = elemAdjAttrGhost->get();
00480         int *adjsGhost = adjTableGhost.getData();
00481         AllocTable2d<int> &adjTypesTableGhost = elemAdjTypesAttrGhost->get();
00482         int *adjTypesGhost = adjTypesTableGhost.getData();
00483         
00484         // initialize tables
00485         for(int i=0;i<numElements*tuplesPerElem;i++){
00486           adjs[i]=-1;
00487           adjTypes[i]=0;
00488         }
00489         for(int i=0;i<numGhostElements*tuplesPerElem;i++){
00490           adjsGhost[i]=-1;
00491           adjTypesGhost[i]=0;
00492         }
00493         
00494         // look through each elemList that is returned by the tuple table
00495         while (NULL!=(l=table.lookupNext())) {
00496           if (l->next==NULL) { 
00497             // One-entry list: must be a symmetry
00498             // UNHANDLED CASE: not sure exactly what this means
00499           }
00500           else { /* Several elements in list: normal case */
00501                   // for each a,b from the list
00502             for (const elemList *a=l;a!=NULL;a=a->next){
00503               for (const elemList *b=l;b!=NULL;b=b->next){
00504                 // if a and b are different elements
00505                 if((a->localNo != b->localNo) || (a->type != b->type)){
00506                   int j;
00507                   if(a->type == t){ // only update the entries for element type t
00508                     
00509                     if(FEM_Is_ghost_index(a->localNo)){
00510                       j = FEM_To_ghost_index(a->localNo)*tuplesPerElem + a->tupleNo;
00511                       CkAssert(j<numGhostElements*tuplesPerElem);
00512                       adjsGhost[j] = b->localNo;
00513                       adjTypesGhost[j] = b->type;
00514                     }
00515                     else {
00516                       j = a->localNo*tuplesPerElem + a->tupleNo;
00517                       CkAssert(j<numElements*tuplesPerElem);
00518                       adjs[j] = b->localNo;
00519                       adjTypes[j] = b->type;
00520                     }
00521                   }
00522                   
00523                 }
00524               }
00525             }
00526           }
00527         }
00528       }
00529 
00530 #ifdef DEBUG
00531   CmiMemoryCheck();
00532 #endif
00533 
00534   delete g;
00535 }
00536 
00537 
00538 
00539 //  ------- Element-to-element: preserve initial ordering relative to nodes
00542 void FEM_Mesh::e2e_getAll(int e, int *neighbors, int etype) 
00543 {
00544 #ifdef DEBUG
00545   CmiMemoryCheck();
00546 #endif
00547   if (e == -1) return; // non existent element
00548   FEM_IndexAttribute *eAdj;
00549   if(FEM_Is_ghost_index(e)){
00550     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getAll");
00551     AllocTable2d<int> &eAdjs = eAdj->get();
00552     for (int i=0; i<eAdjs.width(); i++) {
00553       neighbors[i] = eAdjs[FEM_To_ghost_index(e)][i];
00554     }
00555   }
00556   else {
00557     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getAll");
00558     AllocTable2d<int> &eAdjs = eAdj->get();
00559     for (int i=0; i<eAdjs.width(); i++) {
00560       neighbors[i] = eAdjs[e][i];
00561     }
00562   }
00563 #ifdef DEBUG
00564   CmiMemoryCheck();
00565 #endif
00566 }
00567 
00569 int FEM_Mesh::e2e_getNbr(int e, short idx, int etype) 
00570 {     
00571 #ifdef DEBUG
00572   CmiMemoryCheck();
00573 #endif
00574   if (e == -1) return -1;
00575   FEM_IndexAttribute *eAdj;
00576   if(FEM_Is_ghost_index(e)){
00577     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getNbr");
00578     AllocTable2d<int> &eAdjs = eAdj->get();
00579     return eAdjs[FEM_To_ghost_index(e)][idx];
00580   }
00581   else {
00582     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getNbr");
00583     AllocTable2d<int> &eAdjs = eAdj->get();
00584     return eAdjs[e][idx];
00585   }
00586 #ifdef DEBUG
00587   CmiMemoryCheck();
00588 #endif
00589 }
00590 
00593 int FEM_Mesh::e2e_getIndex(int e, int nbr, int etype) 
00594 { 
00595 #ifdef DEBUG
00596   CmiMemoryCheck();
00597 #endif
00598   if (e == -1) return -1;
00599   FEM_IndexAttribute *eAdj;
00600   if(FEM_Is_ghost_index(e)){
00601     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getIndex");
00602     AllocTable2d<int> &eAdjs = eAdj->get();
00603     for (int i=0; i<eAdjs.width(); i++) {
00604       if (eAdjs[FEM_To_ghost_index(e)][i] == nbr) {
00605         return i;
00606       }
00607     }
00608   }
00609   else{
00610     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getIndex");
00611     AllocTable2d<int> &eAdjs = eAdj->get();
00612     for (int i=0; i<eAdjs.width(); i++) {
00613       if (eAdjs[e][i] == nbr) {
00614         return i;
00615       }
00616     }
00617   }
00618 #ifdef DEBUG
00619   CmiMemoryCheck();
00620 #endif
00621   return -1;
00622 }
00623 
00626 void FEM_Mesh::e2e_setAll(int e, int *neighbors, int etype) 
00627 {
00628 #ifdef DEBUG
00629   CmiMemoryCheck();
00630 #endif
00631   if (e == -1) return;
00632   FEM_IndexAttribute *eAdj;
00633   if(FEM_Is_ghost_index(e)){
00634     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_setAll");
00635     AllocTable2d<int> &eAdjs = eAdj->get();
00636     //CkPrintf("e2e_setAll: Setting element %d's neighbors to %d,%d,%d\n", e, 
00637     //neighbors[0], neighbors[1], neighbors[2]);
00638     for (int i=0; i<eAdjs.width(); i++) {
00639       eAdjs[FEM_To_ghost_index(e)][i] = neighbors[i];
00640     }
00641   }
00642   else{
00643     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_setAll");
00644     AllocTable2d<int> &eAdjs = eAdj->get();
00645     //CkPrintf("e2e_setAll: Setting element %d's neighbors to %d,%d,%d\n", e, 
00646     //neighbors[0], neighbors[1], neighbors[2]);
00647     for (int i=0; i<eAdjs.width(); i++) {
00648       eAdjs[e][i] = neighbors[i];
00649     }
00650   }
00651 #ifdef DEBUG
00652   CmiMemoryCheck();
00653 #endif
00654   
00655 }
00656 
00658 void FEM_Mesh::e2e_setIndex(int e, short idx, int newElem, int etype) 
00659 {
00660 #ifdef DEBUG
00661   CmiMemoryCheck();
00662 #endif
00663   if (e == -1) return;
00664   FEM_IndexAttribute *eAdj;
00665   if(FEM_Is_ghost_index(e)){
00666     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_setIndex");
00667     AllocTable2d<int> &eAdjs = eAdj->get();
00668     eAdjs[FEM_To_ghost_index(e)][idx] = newElem;
00669   }
00670   else {
00671     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_setIndex");    
00672     AllocTable2d<int> &eAdjs = eAdj->get();
00673     eAdjs[e][idx] = newElem;
00674   }
00675 #ifdef DEBUG
00676   CmiMemoryCheck();
00677 #endif
00678 }
00679 
00681 void FEM_Mesh::e2e_replace(int e, int oldNbr, int newNbr, int etype) 
00682 {
00683 #ifdef DEBUG
00684   CmiMemoryCheck();
00685 #endif
00686   if (e == -1) return;
00687   FEM_IndexAttribute *eAdj;
00688   if(FEM_Is_ghost_index(e)){
00689     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_replace");
00690     AllocTable2d<int> &eAdjs = eAdj->get();
00691     for (int i=0; i<eAdjs.width(); i++) {
00692       if (eAdjs[FEM_To_ghost_index(e)][i] == oldNbr) {
00693         eAdjs[FEM_To_ghost_index(e)][i] = newNbr;
00694         break;
00695       }
00696     }
00697   }
00698   else{
00699     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_replace");
00700     AllocTable2d<int> &eAdjs = eAdj->get();
00701     for (int i=0; i<eAdjs.width(); i++) {
00702       if (eAdjs[e][i] == oldNbr) {
00703         eAdjs[e][i] = newNbr;
00704         break;
00705       }
00706     }
00707   }
00708 #ifdef DEBUG
00709   CmiMemoryCheck();
00710 #endif
00711 }
00712 
00714 void FEM_Mesh::e2e_removeAll(int e, int etype)
00715 {
00716 #ifdef DEBUG
00717   CmiMemoryCheck();
00718 #endif
00719   if (e == -1) return;
00720   FEM_IndexAttribute *eAdj;
00721   if(FEM_Is_ghost_index(e)) {
00722     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_removeAll");
00723     AllocTable2d<int> &eAdjs = eAdj->get();
00724     for (int i=0; i<eAdjs.width(); i++) {
00725       eAdjs[FEM_To_ghost_index(e)][i] = -1;
00726     }
00727   }
00728   else {
00729     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_removeAll");
00730     AllocTable2d<int> &eAdjs = eAdj->get();
00731     for (int i=0; i<eAdjs.width(); i++) {
00732       eAdjs[e][i] = -1;
00733     }
00734   }
00735 #ifdef DEBUG
00736   CmiMemoryCheck();
00737 #endif
00738 }
00739 
00740 //  ------- Element-to-node: preserve initial ordering
00743 void FEM_Mesh::e2n_getAll(int e, int *adjnodes, int etype) 
00744 {
00745 #ifdef DEBUG
00746   CmiMemoryCheck();
00747 #endif
00748   if (e == -1) return;
00749   if(FEM_Is_ghost_index(e)){
00750     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00751     AllocTable2d<int> &conn = eConn->get();
00752     for (int i=0; i<conn.width(); i++) {
00753       adjnodes[i] = conn[FEM_To_ghost_index(e)][i];
00754     }
00755   }
00756   else{
00757     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00758     AllocTable2d<int> &conn = eConn->get();
00759     for (int i=0; i<conn.width(); i++) {
00760       adjnodes[i] = conn[e][i];
00761     }
00762   }
00763 #ifdef DEBUG
00764   CmiMemoryCheck();
00765 #endif
00766 }
00767 
00769 int FEM_Mesh::e2n_getNode(int e, short idx, int etype) 
00770 { 
00771 #ifdef DEBUG
00772   CmiMemoryCheck();
00773 #endif
00774   if (e == -1) return -1;
00775   if(FEM_Is_ghost_index(e)){
00776     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00777     AllocTable2d<int> &conn = eConn->get();
00778     return conn[FEM_To_ghost_index(e)][idx];
00779   }
00780   else{
00781     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00782     AllocTable2d<int> &conn = eConn->get();
00783     return conn[e][idx];
00784   }
00785 #ifdef DEBUG
00786   CmiMemoryCheck();
00787 #endif
00788 }
00789 
00792 short FEM_Mesh::e2n_getIndex(int e, int n, int etype) 
00793 {
00794 #ifdef DEBUG
00795   CmiMemoryCheck();
00796 #endif
00797   if (e == -1) return -1;
00798   if(FEM_Is_ghost_index(e)){
00799     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00800     AllocTable2d<int> &conn = eConn->get();
00801     for (int i=0; i<conn.width(); i++)
00802       if (conn[FEM_To_ghost_index(e)][i] == n) 
00803         return i;
00804   }
00805   else{
00806     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00807     AllocTable2d<int> &conn = eConn->get();
00808     for (int i=0; i<conn.width(); i++)
00809       if (conn[e][i] == n)
00810         return i;
00811   }
00812 #ifdef DEBUG
00813   CmiMemoryCheck();
00814 #endif
00815   return -1;
00816 }
00817 
00820 void FEM_Mesh::e2n_setAll(int e, int *adjnodes, int etype) 
00821 {
00822 #ifdef DEBUG
00823   CmiMemoryCheck();
00824 #endif
00825   if (e == -1) return;
00826   if(FEM_Is_ghost_index(e)){
00827     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00828     AllocTable2d<int> &conn = eConn->get();
00829     for (int i=0; i<conn.width(); i++) {
00830       conn[FEM_To_ghost_index(e)][i] = adjnodes[i];
00831     }
00832   }
00833   else{
00834     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00835     AllocTable2d<int> &conn = eConn->get();
00836     for (int i=0; i<conn.width(); i++) {
00837       conn[e][i] = adjnodes[i];
00838     }
00839   }  
00840 #ifdef DEBUG
00841   CmiMemoryCheck();
00842 #endif
00843 }
00844 
00846 void FEM_Mesh::e2n_setIndex(int e, short idx, int newNode, int etype) 
00847 {
00848 #ifdef DEBUG
00849   CmiMemoryCheck();
00850 #endif
00851   if (e == -1) return;
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     conn[FEM_To_ghost_index(e)][idx] = newNode;
00856   }
00857   else{
00858     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00859     AllocTable2d<int> &conn = eConn->get();
00860     conn[e][idx] = newNode;
00861   }
00862 #ifdef DEBUG
00863   CmiMemoryCheck();
00864 #endif
00865 }
00866 
00868 void FEM_Mesh::e2n_replace(int e, int oldNode, int newNode, int etype) 
00869 {
00870 #ifdef DEBUG
00871   CmiMemoryCheck();
00872 #endif
00873   if (e == -1) return;
00874   if(FEM_Is_ghost_index(e)){
00875     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00876     AllocTable2d<int> &conn = eConn->get();
00877     for (int i=0; i<conn.width(); i++) {
00878       if (conn[FEM_To_ghost_index(e)][i] == oldNode) {
00879         conn[FEM_To_ghost_index(e)][i] = newNode;
00880         break;
00881       }
00882     }
00883   }
00884   else{
00885     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00886     AllocTable2d<int> &conn = eConn->get();
00887     for (int i=0; i<conn.width(); i++) {
00888       if (conn[e][i] == oldNode) {
00889         conn[e][i] = newNode;
00890         break;
00891       }
00892     }
00893   }  
00894 #ifdef DEBUG
00895   CmiMemoryCheck();
00896 #endif
00897 }
00898 
00899 void FEM_Mesh::e2n_removeAll(int e, int etype)
00900 {
00901 #ifdef DEBUG
00902   CmiMemoryCheck();
00903 #endif
00904   if (e == -1) return;
00905   if(FEM_Is_ghost_index(e)){
00906     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00907     AllocTable2d<int> &conn = eConn->get();
00908     for (int i=0; i<conn.width(); i++)
00909       conn[FEM_To_ghost_index(e)][i] = -1;
00910   }
00911   else{
00912     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00913     AllocTable2d<int> &conn = eConn->get();
00914     for (int i=0; i<conn.width(); i++)
00915       conn[e][i] = -1;
00916   }  
00917 #ifdef DEBUG
00918   CmiMemoryCheck();
00919 #endif
00920 }
00921 
00922 
00923 
00924 
00925 //  ------- Node-to-node
00928 void FEM_Mesh::n2n_getAll(int n, int **adjnodes, int *sz) 
00929 {
00930 #ifdef DEBUG
00931   CmiMemoryCheck();
00932 #endif
00933   if (n == -1) {
00934     *sz = 0;
00935     return;
00936   }
00937   if(FEM_Is_ghost_index(n)){
00938     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_getAll");
00939     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &nVec = nAdj->get();
00940     CkVec<FEM_VarIndexAttribute::ID> &nsVec = nVec[FEM_To_ghost_index(n)];
00941     *sz = nsVec.length();
00942     if(*sz != 0) (*adjnodes) = new int[*sz];
00943     for (int i=0; i<(*sz); i++) {
00944       (*adjnodes)[i] = nsVec[i].getSignedId();
00945     }
00946   }
00947   else{
00948     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_getAll");
00949     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &nVec = nAdj->get();
00950     CkVec<FEM_VarIndexAttribute::ID> &nsVec = nVec[n];
00951     *sz = nsVec.length();
00952     if(*sz != 0) (*adjnodes) = new int[*sz];
00953     for (int i=0; i<(*sz); i++) {
00954       (*adjnodes)[i] = nsVec[i].getSignedId();
00955     }
00956   }
00957   
00958 #ifdef DEBUG
00959   CmiMemoryCheck();
00960 #endif
00961 }
00962  
00964 void FEM_Mesh::n2n_add(int n, int newNode) 
00965 {
00966 #ifdef DEBUG
00967   CmiMemoryCheck();
00968 #endif
00969   if (n == -1) return;
00970   if(FEM_Is_ghost_index(n)){
00971     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_add");
00972     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &nVec = nAdj->get();
00973     FEM_VarIndexAttribute::ID nn(0, newNode);
00974     CkVec<FEM_VarIndexAttribute::ID> &nsVec = nVec[FEM_To_ghost_index(n)];
00975     nsVec.push_back(nn);
00976   }
00977   else{
00978     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_add");
00979     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &nVec = nAdj->get();
00980     FEM_VarIndexAttribute::ID nn(0, newNode);
00981     CkVec<FEM_VarIndexAttribute::ID> &nsVec = nVec[n];
00982     nsVec.push_back(nn);
00983   }
00984 #ifdef DEBUG
00985   CmiMemoryCheck();
00986 #endif
00987 }
00988 
00989 
00990 
00991 
00993 void FEM_Mesh::n2n_remove(int n, int oldNode) 
00994 {
00995 #ifdef DEBUG
00996   CmiMemoryCheck();
00997 #endif
00998   if (n == -1) return;
00999   if(FEM_Is_ghost_index(n)){
01000     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_remove");
01001     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &nVec = nAdj->get();
01002     CkVec<FEM_VarIndexAttribute::ID> &nsVec = nVec[FEM_To_ghost_index(n)];
01003     for (int i=0; i<nsVec.length(); i++) {
01004       if (nsVec[i].getSignedId() == oldNode) {
01005         nsVec.remove(i);
01006         break;
01007       }
01008     }
01009   }
01010   else {
01011     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_remove");
01012     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &nVec = nAdj->get();
01013     CkVec<FEM_VarIndexAttribute::ID> &nsVec = nVec[n];
01014     for (int i=0; i<nsVec.length(); i++) {
01015       if (nsVec[i].getSignedId() == oldNode) {
01016         nsVec.remove(i);
01017         break;
01018       }
01019     }
01020   }
01021 #ifdef DEBUG
01022   CmiMemoryCheck();
01023 #endif
01024 }
01025 
01027 int FEM_Mesh::n2n_exists(int n, int queryNode) 
01028 {
01029 #ifdef DEBUG
01030   CmiMemoryCheck();
01031 #endif
01032   if (n == -1) return 0;
01033   if(FEM_Is_ghost_index(n)){
01034     CkAssert(node.getGhost());
01035     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_exists");
01036     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &nVec = nAdj->get();
01037     CkVec<FEM_VarIndexAttribute::ID> &nsVec = nVec[FEM_To_ghost_index(n)];
01038     for (int i=0; i<nsVec.length(); i++)
01039       if (nsVec[i].getSignedId() == queryNode) 
01040         return 1;
01041   }
01042   else {
01043     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_exists");
01044     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &nVec = nAdj->get();
01045     CkVec<FEM_VarIndexAttribute::ID> &nsVec = nVec[n];
01046     for (int i=0; i<nsVec.length(); i++)
01047       if (nsVec[i].getSignedId() == queryNode) 
01048     return 1;
01049   }
01050 #ifdef DEBUG
01051   CmiMemoryCheck();
01052 #endif
01053   return 0;
01054 }
01055 
01057 void FEM_Mesh::n2n_replace(int n, int oldNode, int newNode) 
01058 {
01059 #ifdef DEBUG
01060   CmiMemoryCheck();
01061 #endif
01062   if (n == -1) return;
01063   if(FEM_Is_ghost_index(n)){
01064     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_replace");
01065     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &nVec = nAdj->get();
01066     CkVec<FEM_VarIndexAttribute::ID> &nsVec = nVec[FEM_To_ghost_index(n)];
01067     for (int i=0; i<nsVec.length(); i++) {
01068       if (nsVec[i].getSignedId() == oldNode) {
01069         nsVec[i] = FEM_VarIndexAttribute::ID(0,newNode);
01070         break;
01071       }
01072     }
01073   }
01074   else {
01075     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_replace");
01076     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &nVec = nAdj->get();
01077     CkVec<FEM_VarIndexAttribute::ID> &nsVec = nVec[n];
01078     for (int i=0; i<nsVec.length(); i++) {
01079       if (nsVec[i].getSignedId() == oldNode) {
01080         nsVec[i] = FEM_VarIndexAttribute::ID(0,newNode);
01081         break;
01082       }
01083     }
01084 
01085   }
01086 #ifdef DEBUG
01087   CmiMemoryCheck();
01088 #endif
01089 }
01090 
01092 void FEM_Mesh::n2n_removeAll(int n)
01093 {
01094 #ifdef DEBUG
01095   CmiMemoryCheck();
01096 #endif
01097   if (n == -1) return;
01098   if(FEM_Is_ghost_index(n)){
01099     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_removeAll");  
01100     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &nVec = nAdj->get();
01101     CkVec<FEM_VarIndexAttribute::ID> &nsVec = nVec[FEM_To_ghost_index(n)];
01102     nsVec.free();
01103   }
01104   else{
01105     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_removeAll");  
01106     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &nVec = nAdj->get();
01107     CkVec<FEM_VarIndexAttribute::ID> &nsVec = nVec[n];
01108     nsVec.free();
01109   }
01110 #ifdef DEBUG
01111   CmiMemoryCheck();
01112 #endif
01113 }
01114 
01115 //  ------- Node-to-element
01119 void FEM_Mesh::n2e_getAll(int n, int **adjelements, int *sz) 
01120 {
01121 #ifdef DEBUG
01122   CmiMemoryCheck();
01123 #endif
01124 
01125 
01126   if (n == -1) {
01127     *sz = 0;
01128     return;
01129   }
01130   if(FEM_Is_ghost_index(n)){
01131     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_getAll");  
01132     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &eVec = eAdj->get();
01133     CkVec<FEM_VarIndexAttribute::ID> &nsVec = eVec[FEM_To_ghost_index(n)];
01134     *sz = nsVec.length();
01135     if(*sz !=0) (*adjelements) = new int[*sz];
01136     for (int i=0; i<(*sz); i++) {
01137       (*adjelements)[i] = nsVec[i].getSignedId();
01138     }
01139   }
01140   else {
01141     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_getAll");  
01142     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &eVec = eAdj->get();
01143     CkVec<FEM_VarIndexAttribute::ID> &nsVec = eVec[n];
01144   int len = nsVec.length();
01145     *sz = len;
01146     if(*sz !=0) (*adjelements) = new int[*sz];
01147     for (int i=0; i<(*sz); i++) {
01148       (*adjelements)[i] = nsVec[i].getSignedId();
01149     }
01150   }
01151 
01152 #ifdef DEBUG
01153   CmiMemoryCheck();
01154 #endif
01155 
01156 }
01157  
01159 void FEM_Mesh::n2e_add(int n, int newElem) 
01160 {
01161 #ifdef DEBUG
01162   CmiMemoryCheck();
01163 #endif
01164 
01165   if (n == -1) return;
01166   if(FEM_Is_ghost_index(n)){
01167     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_add");     
01168     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &eVec = eAdj->get();
01169     CkVec<FEM_VarIndexAttribute::ID> &nsVec = eVec[FEM_To_ghost_index(n)];
01170     FEM_VarIndexAttribute::ID ne(0, newElem);
01171     nsVec.push_back(ne);
01172     int *testn2e, testn2ec;
01173     n2e_getAll(n,&testn2e,&testn2ec);
01174     for(int i=0; i<testn2ec; i++) {
01175       if(FEM_Is_ghost_index(testn2e[i]))
01176         CkAssert(elem[0].ghost->is_valid(FEM_From_ghost_index(testn2e[i])));
01177       else 
01178         CkAssert(elem[0].is_valid(testn2e[i]));
01179     }
01180     if(testn2ec!=0) delete[] testn2e;
01181   }
01182   else {
01183     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_add");     
01184     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &eVec = eAdj->get();
01185     CkVec<FEM_VarIndexAttribute::ID> &nsVec = eVec[n];
01186     FEM_VarIndexAttribute::ID ne(0, newElem);
01187     nsVec.push_back(ne);
01188   }
01189 #ifdef DEBUG
01190   CmiMemoryCheck();
01191 #endif
01192 
01193 }
01194 
01196 void FEM_Mesh::n2e_remove(int n, int oldElem) 
01197 {
01198 #ifdef DEBUG
01199   CmiMemoryCheck();
01200 #endif
01201   if (n == -1) return;
01202   if(FEM_Is_ghost_index(n)){
01203     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_remove");
01204     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &eVec = eAdj->get();
01205     CkVec<FEM_VarIndexAttribute::ID> &nsVec = eVec[FEM_To_ghost_index(n)];
01206     for (int i=0; i<nsVec.length(); i++) {
01207       if (nsVec[i].getSignedId() == oldElem) {
01208         nsVec.remove(i);
01209         break;
01210       }
01211     }
01212   }
01213   else {
01214     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_remove");
01215     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &eVec = eAdj->get();
01216     CkVec<FEM_VarIndexAttribute::ID> &nsVec = eVec[n];
01217     for (int i=0; i<nsVec.length(); i++) {
01218       if (nsVec[i].getSignedId() == oldElem) {
01219         nsVec.remove(i);
01220         break;
01221       }
01222     }
01223   }
01224 #ifdef DEBUG
01225   CmiMemoryCheck();
01226 #endif
01227 }
01228  
01230 void FEM_Mesh::n2e_replace(int n, int oldElem, int newElem) 
01231 {
01232 #ifdef DEBUG
01233   CmiMemoryCheck();
01234 #endif
01235   if (n == -1) return;
01236   if(FEM_Is_ghost_index(n)){
01237     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_replace");
01238     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &eVec = eAdj->get();
01239     CkVec<FEM_VarIndexAttribute::ID> &nsVec = eVec[FEM_To_ghost_index(n)];
01240     for (int i=0; i<nsVec.length(); i++) {
01241       if (nsVec[i].getSignedId() == oldElem) {
01242         nsVec[i] = FEM_VarIndexAttribute::ID(0,newElem);
01243         break;
01244       }
01245     }
01246     int *testn2e, testn2ec;
01247     n2e_getAll(n,&testn2e,&testn2ec);
01248     for(int i=0; i<testn2ec; i++) {
01249       if(FEM_Is_ghost_index(testn2e[i]))
01250         CkAssert(elem[0].ghost->is_valid(FEM_From_ghost_index(testn2e[i])));
01251       else 
01252         CkAssert(elem[0].is_valid(testn2e[i]));
01253     }
01254     if(testn2ec!=0) delete[] testn2e;
01255   }
01256   else{
01257     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_replace");
01258     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &eVec = eAdj->get();
01259     CkVec<FEM_VarIndexAttribute::ID> &nsVec = eVec[n];
01260     for (int i=0; i<nsVec.length(); i++) {
01261       if (nsVec[i].getSignedId() == oldElem) {
01262         nsVec[i] = FEM_VarIndexAttribute::ID(0,newElem);
01263         break;
01264       }
01265     }
01266   }
01267 #ifdef DEBUG
01268   CmiMemoryCheck();
01269 #endif
01270 }
01271 
01273 void FEM_Mesh::n2e_removeAll(int n)
01274 {
01275 #ifdef DEBUG
01276   CmiMemoryCheck();
01277 #endif
01278   if (n == -1) return;
01279   if(FEM_Is_ghost_index(n)){
01280     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_removeAll");
01281     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &eVec = eAdj->get();
01282     CkVec<FEM_VarIndexAttribute::ID> &nsVec = eVec[FEM_To_ghost_index(n)];
01283     nsVec.free();
01284   }
01285   else{
01286     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_removeAll");
01287     CkVec<CkVec<FEM_VarIndexAttribute::ID> > &eVec = eAdj->get();
01288     CkVec<FEM_VarIndexAttribute::ID> &nsVec = eVec[n];
01289     nsVec.free();
01290   }
01291 #ifdef DEBUG
01292   CmiMemoryCheck();
01293 #endif
01294 }
01295 
01298 
01299 // I think the break statement may cause some unexpected behaviour here
01300 // The break will break out of the inner for loop back to the outer for loop
01301 // I am not sure exactly what the symantics of this function should be
01302 //   -- Isaac
01303 
01304 int FEM_Mesh::getElementOnEdge(int n1, int n2) 
01305 {
01306   int *n1AdjElems, *n2AdjElems;
01307   int n1NumElems, n2NumElems;
01308   n2e_getAll(n1, &n1AdjElems, &n1NumElems);
01309   n2e_getAll(n2, &n2AdjElems, &n2NumElems);
01310   int ret = -1;
01311   //CkPrintf("%d has %d neighboring elements, %d has %d\n", n1, n1NumElems, n2, n2NumElems);
01312 
01313 #ifdef DEBUG
01314   CmiMemoryCheck();
01315 #endif 
01316 
01317   for (int i=0; i<n1NumElems; i++) {
01318     for (int j=0; j<n2NumElems; j++) {
01319       if (n1AdjElems[i] == n2AdjElems[j]) {
01320         if(n1AdjElems[i] >= 0) {
01321           ret = n1AdjElems[i];
01322           break;
01323         }
01324         else {
01325           ret = n1AdjElems[i];
01326         }
01327       }
01328     }
01329   }
01330   delete[] n1AdjElems;
01331   delete[] n2AdjElems;
01332   return ret; //preferably return a local element, otherwise return a ghost 
01333 }
01334 
01335 
01339 void FEM_Mesh::get2ElementsOnEdge(int n1, int n2, int *result_e1, int *result_e2) 
01340 {
01341   int *n1AdjElems=0, *n2AdjElems=0;
01342   int n1NumElems, n2NumElems;
01343 
01344 #ifdef DEBUG
01345   CmiMemoryCheck();
01346 #endif
01347 
01348   n2e_getAll(n1, &n1AdjElems, &n1NumElems);
01349   n2e_getAll(n2, &n2AdjElems, &n2NumElems);
01350   CkAssert(n1AdjElems!=0);
01351   CkAssert(n2AdjElems!=0);
01352   int found=0;
01353 
01354 #ifdef DEBUG
01355   CmiMemoryCheck();
01356 #endif
01357 
01358   *result_e1=-1;
01359   *result_e2=-1;
01360 
01361   for (int i=0; i<n1NumElems; i++) {
01362     for (int j=0; j<n2NumElems; j++) {
01363       if (n1AdjElems[i] == n2AdjElems[j]) {
01364         if(found==0){
01365           //          CkPrintf("found element1 %d\n", n1AdjElems[i]);
01366           *result_e1 = n1AdjElems[i];
01367           found++;
01368         }
01369         else if(found==1){
01370           //   CkPrintf("found element2 %d\n", n1AdjElems[i]);
01371           *result_e2 = n1AdjElems[i];
01372           found++;
01373         }
01374         else {
01375 #ifdef DEBUG
01376           CmiMemoryCheck();
01377 #endif          
01378           CkPrintf("ERROR: Found a third element(%d) on edge %d,%d \n", n1AdjElems[i], n1, n2);
01379            CkExit();
01380         }
01381       }
01382     }
01383   }
01384 
01385 
01386 #ifdef DEBUG
01387   CmiMemoryCheck();
01388 #endif 
01389   delete[] n1AdjElems;
01390   delete[] n2AdjElems;
01391 }
01392
libs/ck-libs/fem/fem_mesh_adjacency.C
