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 CDECL 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 CDECL 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 CDECL 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 CDECL 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 CDECL 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 CDECL 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 }
00513 
00514 
00515 
00520 void FEM_Mesh::e2e_getAll(int e, int *neighbors, int etype) 
00521 {
00522   if (e == -1) return; // non existent element
00523   FEM_IndexAttribute *eAdj;
00524   if(FEM_Is_ghost_index(e)){
00525     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getAll");
00526     AllocTable2d<int> &eAdjs = eAdj->get();
00527     for (int i=0; i<eAdjs.width(); i++) {
00528       neighbors[i] = eAdjs[FEM_To_ghost_index(e)][i];
00529     }
00530   }
00531   else {
00532     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getAll");
00533     AllocTable2d<int> &eAdjs = eAdj->get();
00534     for (int i=0; i<eAdjs.width(); i++) {
00535       neighbors[i] = eAdjs[e][i];
00536     }
00537   }
00538 }
00539 
00542 int FEM_Mesh::e2e_getNbr(int e, short idx, int etype) 
00543 {     
00544   if (e == -1) return -1;
00545   FEM_IndexAttribute *eAdj;
00546   if(FEM_Is_ghost_index(e)){
00547     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getNbr");
00548     AllocTable2d<int> &eAdjs = eAdj->get();
00549     return eAdjs[FEM_To_ghost_index(e)][idx];
00550   }
00551   else {
00552     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getNbr");
00553     AllocTable2d<int> &eAdjs = eAdj->get();
00554     return eAdjs[e][idx];
00555   }
00556 }
00557 
00558 
00559 ElemID FEM_Mesh::e2e_getElem(ElemID elem, int idx){
00560         return e2e_getElem(elem.id, idx, elem.type);
00561 }
00562 
00563 ElemID FEM_Mesh::e2e_getElem(ElemID *elem, int idx){
00564         return e2e_getElem(elem->id, idx, elem->type);
00565 }
00566 
00567 
00568 
00570 ElemID FEM_Mesh::e2e_getElem(int e, int idx, int etype){
00571   if (e == -1){
00572     ElemID ele(-1,-1);    
00573     return ele;
00574   }
00575   
00576   if(FEM_Is_ghost_index(e)){
00577 
00578           FEM_IndexAttribute *eAdj;
00579           FEM_IndexAttribute *eAdjType;
00580 
00581           eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getElem");
00582           AllocTable2d<int> &eAdjs = eAdj->get();
00583           eAdjType = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJ_TYPES,"e2e_getElem");
00584           AllocTable2d<int> &eAdjTypes = eAdjType->get();
00585 
00586           int t =  eAdjTypes[FEM_To_ghost_index(e)][idx];
00587           int id =  eAdjs[FEM_To_ghost_index(e)][idx];
00588           ElemID ele(t,id);
00589 
00590           return ele;
00591   }
00592   else {
00593 
00594           CkAssert(elem.has(etype) && e < elem[etype].size());
00595 
00596           FEM_IndexAttribute *eAdj;
00597           FEM_IndexAttribute *eAdjType;
00598 
00599           eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getElem");
00600           AllocTable2d<int> &eAdjs = eAdj->get();
00601           eAdjType = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJ_TYPES,"e2e_getElem");
00602           AllocTable2d<int> &eAdjTypes = eAdjType->get();
00603 
00604           int t =  eAdjTypes[e][idx];
00605           int id =  eAdjs[e][idx];
00606           ElemID ele(t,id);
00607 
00608           return ele;
00609   }
00610   
00611 }
00612 
00613 
00617 int FEM_Mesh::e2e_getIndex(int e, int nbr, int etype) 
00618 { 
00619   if (e == -1) return -1;
00620   FEM_IndexAttribute *eAdj;
00621   if(FEM_Is_ghost_index(e)){
00622     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getIndex");
00623     AllocTable2d<int> &eAdjs = eAdj->get();
00624     for (int i=0; i<eAdjs.width(); i++) {
00625       if (eAdjs[FEM_To_ghost_index(e)][i] == nbr) {
00626         return i;
00627       }
00628     }
00629   }
00630   else{
00631     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_getIndex");
00632     AllocTable2d<int> &eAdjs = eAdj->get();
00633     for (int i=0; i<eAdjs.width(); i++) {
00634       if (eAdjs[e][i] == nbr) {
00635         return i;
00636       }
00637     }
00638   }
00639   return -1;
00640 }
00641 
00645 void FEM_Mesh::e2e_setAll(int e, int *neighbors, int etype) 
00646 {
00647   if (e == -1) return;
00648   FEM_IndexAttribute *eAdj;
00649   if(FEM_Is_ghost_index(e)){
00650     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_setAll");
00651     AllocTable2d<int> &eAdjs = eAdj->get();
00652     //CkPrintf("e2e_setAll: Setting element %d's neighbors to %d,%d,%d\n", e, 
00653     //neighbors[0], neighbors[1], neighbors[2]);
00654     for (int i=0; i<eAdjs.width(); i++) {
00655       eAdjs[FEM_To_ghost_index(e)][i] = neighbors[i];
00656     }
00657   }
00658   else{
00659     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_setAll");
00660     AllocTable2d<int> &eAdjs = eAdj->get();
00661     //CkPrintf("e2e_setAll: Setting element %d's neighbors to %d,%d,%d\n", e, 
00662     //neighbors[0], neighbors[1], neighbors[2]);
00663     for (int i=0; i<eAdjs.width(); i++) {
00664       eAdjs[e][i] = neighbors[i];
00665     }
00666   }
00667 }
00668 
00671 void FEM_Mesh::e2e_setIndex(int e, short idx, int newElem, int etype) 
00672 {
00673   if (e == -1) return;
00674   FEM_IndexAttribute *eAdj;
00675   if(FEM_Is_ghost_index(e)){
00676     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_setIndex");
00677     AllocTable2d<int> &eAdjs = eAdj->get();
00678     eAdjs[FEM_To_ghost_index(e)][idx] = newElem;
00679   }
00680   else {
00681     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_setIndex");    
00682     AllocTable2d<int> &eAdjs = eAdj->get();
00683     eAdjs[e][idx] = newElem;
00684   }
00685 }
00686 
00689 void FEM_Mesh::e2e_replace(int e, int oldNbr, int newNbr, int etype) 
00690 {
00691   if (e == -1) return;
00692   FEM_IndexAttribute *eAdj;
00693   if(FEM_Is_ghost_index(e)){
00694     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_replace");
00695     AllocTable2d<int> &eAdjs = eAdj->get();
00696     for (int i=0; i<eAdjs.width(); i++) {
00697       if (eAdjs[FEM_To_ghost_index(e)][i] == oldNbr) {
00698         eAdjs[FEM_To_ghost_index(e)][i] = newNbr;
00699         break;
00700       }
00701     }
00702   }
00703   else{
00704     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_replace");
00705     AllocTable2d<int> &eAdjs = eAdj->get();
00706     for (int i=0; i<eAdjs.width(); i++) {
00707       if (eAdjs[e][i] == oldNbr) {
00708         eAdjs[e][i] = newNbr;
00709         break;
00710       }
00711     }
00712   }
00713 }
00714 
00715 
00716 
00719 void FEM_Mesh::e2e_replace(ElemID e, ElemID oldNbr, ElemID newNbr) 
00720 {
00721   if (e.id == -1) return;
00722 
00723 //      CkPrintf("replacing %d,%d with %d,%d\n", oldNbr.type, oldNbr.id, newNbr.type, newNbr.id);
00724 
00725   
00726   FEM_IndexAttribute *eAdj;
00727   if(FEM_Is_ghost_index(e.id)){
00728           FEM_IndexAttribute *elemAdjTypesAttrGhost = (FEM_IndexAttribute *)elem[e.type].getGhost()->lookup(FEM_ELEM_ELEM_ADJ_TYPES,"e2e_replace");
00729           FEM_IndexAttribute *elemAdjAttrGhost = (FEM_IndexAttribute *)elem[e.type].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_replace");
00730    
00731           AllocTable2d<int> &eAdjs = elemAdjAttrGhost->get();
00732           AllocTable2d<int> &eAdjTypes = elemAdjTypesAttrGhost->get();
00733           
00734           for (int i=0; i<eAdjs.width(); i++) {
00735                   if (eAdjs[FEM_To_ghost_index(e.id)][i] == oldNbr.id && eAdjTypes[FEM_To_ghost_index(e.id)][i] == oldNbr.type ) {
00736                           eAdjs[FEM_To_ghost_index(e.id)][i] = newNbr.id;
00737                           eAdjTypes[FEM_To_ghost_index(e.id)][i] = newNbr.type;
00738                           break;
00739                   }
00740           }
00741   }
00742   else{
00743     eAdj = (FEM_IndexAttribute *)elem[e.type].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_replace");
00744   
00745     FEM_IndexAttribute *elemAdjTypesAttr = (FEM_IndexAttribute *)elem[e.type].lookup(FEM_ELEM_ELEM_ADJ_TYPES,"e2e_replace");
00746     FEM_IndexAttribute *elemAdjAttr = (FEM_IndexAttribute *)elem[e.type].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_replace");
00747     
00748     AllocTable2d<int> &eAdjs = elemAdjAttr->get();
00749     AllocTable2d<int> &eAdjTypes = elemAdjTypesAttr->get();
00750 
00751     for (int i=0; i<eAdjs.width(); i++) {
00752         if (eAdjs[e.id][i] == oldNbr.id && eAdjTypes[e.id][i] == oldNbr.type ) {
00753 //              CkPrintf("eAdjTypes[e.id][i] was %d eAdjs[e.id][i] was %d\n", eAdjTypes[e.id][i], eAdjs[e.id][i]);
00754                 eAdjs[e.id][i] = newNbr.id;
00755                 eAdjTypes[e.id][i] = newNbr.type;
00756 //              CkPrintf("newNbr.type = %d\n", newNbr.type);
00757 //              CkPrintf("eAdjTypes[e.id][i] is %d eAdjs[e.id][i] is %d\n", eAdjTypes[e.id][i], eAdjs[e.id][i]);
00758                 break;
00759         }
00760     }
00761     
00762   }
00763 }
00764 
00765 
00768 void FEM_Mesh::e2e_removeAll(int e, int etype)
00769 {
00770   if (e == -1) return;
00771   FEM_IndexAttribute *eAdj;
00772   if(FEM_Is_ghost_index(e)) {
00773     eAdj = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_removeAll");
00774     AllocTable2d<int> &eAdjs = eAdj->get();
00775     for (int i=0; i<eAdjs.width(); i++) {
00776       eAdjs[FEM_To_ghost_index(e)][i] = -1;
00777     }
00778   }
00779   else {
00780     eAdj = (FEM_IndexAttribute *)elem[etype].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_removeAll");
00781     AllocTable2d<int> &eAdjs = eAdj->get();
00782     for (int i=0; i<eAdjs.width(); i++) {
00783       eAdjs[e][i] = -1;
00784     }
00785   }
00786 }
00787 
00788 
00791 void FEM_Mesh::e2e_printAll(ElemID e) 
00792 {
00793         if (e.id == -1) return; // non existent element
00794         
00795         FEM_IndexAttribute *eAdj;
00796         FEM_IndexAttribute *eAdjType;
00797         
00798         if(FEM_Is_ghost_index(e.id)){
00799                 eAdj = (FEM_IndexAttribute *)elem[e.type].getGhost()->lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_printAll");
00800                 eAdjType = (FEM_IndexAttribute *)elem[e.type].getGhost()->lookup(FEM_ELEM_ELEM_ADJ_TYPES,"e2e_printAll");
00801         }
00802         else {
00803                 eAdj = (FEM_IndexAttribute *)elem[e.type].lookup(FEM_ELEM_ELEM_ADJACENCY,"e2e_printAll");
00804                 eAdjType = (FEM_IndexAttribute *)elem[e.type].lookup(FEM_ELEM_ELEM_ADJ_TYPES,"e2e_printAll");
00805         }
00806         
00807         
00808         AllocTable2d<int> &eAdjs = eAdj->get();
00809         AllocTable2d<int> &eAdjTypes = eAdjType->get();
00810         CkAssert(eAdjs.width() == eAdjTypes.width());
00811 
00812         CkAssert(e.getSignedId()>=0);
00813 
00814         for (int i=0; i<eAdjs.width(); i++) {
00815                 CkPrintf("Element %d,%d is adjacent to %d,%d\n", e.type, e.id, eAdjTypes[e.getSignedId()][i], eAdjs[e.getSignedId()][i]);
00816 
00817         }
00818 
00819 }
00820 
00821 
00822 
00827 void FEM_Mesh::e2n_getAll(int e, int *adjnodes, int etype) 
00828 {
00829   if (e == -1) return;
00830   if(FEM_Is_ghost_index(e)){
00831     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00832     AllocTable2d<int> &conn = eConn->get();
00833     for (int i=0; i<conn.width(); i++) {
00834       adjnodes[i] = conn[FEM_To_ghost_index(e)][i];
00835     }
00836   }
00837   else{
00838     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00839     AllocTable2d<int> &conn = eConn->get();
00840     for (int i=0; i<conn.width(); i++) {
00841       adjnodes[i] = conn[e][i];
00842     }
00843   }
00844 }
00845 
00848 int FEM_Mesh::e2n_getNode(int e, short idx, int etype) 
00849 { 
00850   if (e == -1) return -1;
00851   if(FEM_Is_ghost_index(e)){
00852     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00853     AllocTable2d<int> &conn = eConn->get();
00854     return conn[FEM_To_ghost_index(e)][idx];
00855   }
00856   else{
00857     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00858     AllocTable2d<int> &conn = eConn->get();
00859     return conn[e][idx];
00860   }
00861 }
00862 
00866 short FEM_Mesh::e2n_getIndex(int e, int n, int etype) 
00867 {
00868   if (e == -1) return -1;
00869   if(FEM_Is_ghost_index(e)){
00870     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00871     AllocTable2d<int> &conn = eConn->get();
00872     for (int i=0; i<conn.width(); i++)
00873       if (conn[FEM_To_ghost_index(e)][i] == n) 
00874         return i;
00875   }
00876   else{
00877     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00878     AllocTable2d<int> &conn = eConn->get();
00879     for (int i=0; i<conn.width(); i++)
00880       if (conn[e][i] == n)
00881         return i;
00882   }
00883   return -1;
00884 }
00885 
00889 void FEM_Mesh::e2n_setAll(int e, int *adjnodes, int etype) 
00890 {
00891   if (e == -1) return;
00892   if(FEM_Is_ghost_index(e)){
00893     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00894     AllocTable2d<int> &conn = eConn->get();
00895     for (int i=0; i<conn.width(); i++) {
00896       conn[FEM_To_ghost_index(e)][i] = adjnodes[i];
00897     }
00898   }
00899   else{
00900     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00901     AllocTable2d<int> &conn = eConn->get();
00902     for (int i=0; i<conn.width(); i++) {
00903       conn[e][i] = adjnodes[i];
00904     }
00905   }  
00906 }
00907 
00910 void FEM_Mesh::e2n_setIndex(int e, short idx, int newNode, int etype) 
00911 {
00912   if (e == -1) return;
00913   if(FEM_Is_ghost_index(e)){
00914     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00915     AllocTable2d<int> &conn = eConn->get();
00916     conn[FEM_To_ghost_index(e)][idx] = newNode;
00917   }
00918   else{
00919     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00920     AllocTable2d<int> &conn = eConn->get();
00921     conn[e][idx] = newNode;
00922   }
00923 }
00924 
00927 void FEM_Mesh::e2n_replace(int e, int oldNode, int newNode, int etype) 
00928 {
00929   if (e == -1) return;
00930   if(FEM_Is_ghost_index(e)){
00931     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00932     AllocTable2d<int> &conn = eConn->get();
00933     for (int i=0; i<conn.width(); i++) {
00934       if (conn[FEM_To_ghost_index(e)][i] == oldNode) {
00935         conn[FEM_To_ghost_index(e)][i] = newNode;
00936         break;
00937       }
00938     }
00939   }
00940   else{
00941     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00942     AllocTable2d<int> &conn = eConn->get();
00943     for (int i=0; i<conn.width(); i++) {
00944       if (conn[e][i] == oldNode) {
00945         conn[e][i] = newNode;
00946         break;
00947       }
00948     }
00949   }  
00950 }
00951 
00952 void FEM_Mesh::e2n_removeAll(int e, int etype)
00953 {
00954   if (e == -1) return;
00955   if(FEM_Is_ghost_index(e)){
00956     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].getGhost()->lookup(FEM_CONN,"e2n_getAll");
00957     AllocTable2d<int> &conn = eConn->get();
00958     for (int i=0; i<conn.width(); i++)
00959           conn[FEM_To_ghost_index(e)][i] = -1;
00960   }
00961   else{
00962     FEM_IndexAttribute *eConn = (FEM_IndexAttribute *)elem[etype].lookup(FEM_CONN,"e2n_getAll");
00963     AllocTable2d<int> &conn = eConn->get();
00964     for (int i=0; i<conn.width(); i++)
00965           conn[e][i] = -1;
00966   }  
00967 }
00968 
00969 
00970 
00971 
00973 int FEM_Mesh::n2n_getLength(int n) {
00974   if (n == -1) {
00975     return 0;
00976   }
00977   FEM_VarIndexAttribute *eAdj;
00978   if(FEM_Is_ghost_index(n)){
00979     eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_getLength");
00980     n = FEM_To_ghost_index(n);
00981   }
00982   else {
00983     eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_getLength");
00984   }
00985   CkVec<CkVec<ElemID> > &eVec = eAdj->get();
00986   CkVec<ElemID> &nsVec = eVec[n];
00987   return nsVec.length();
00988 }
00989 
00993 void FEM_Mesh::n2n_getAll(int n, int *&adjnodes, int &sz) 
00994 {
00995   if (n == -1) {
00996     sz = 0;
00997     return;
00998   }
00999   if(FEM_Is_ghost_index(n)){
01000         FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_getAll");
01001         CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01002         CkVec<ElemID> &nsVec = nVec[FEM_To_ghost_index(n)];
01003         sz = nsVec.length();
01004         if(sz > 0) adjnodes = new int[sz];
01005         for (int i=0; i<sz; i++) {
01006           adjnodes[i] = nsVec[i].getSignedId();
01007         }
01008   }
01009   else{
01010         FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_getAll");
01011         CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01012         CkVec<ElemID> &nsVec = nVec[n];
01013         sz = nsVec.length();
01014         if(sz > 0) adjnodes = new int[sz];
01015         for (int i=0; i<sz; i++) {
01016           adjnodes[i] = nsVec[i].getSignedId();
01017         }
01018   }
01019 }
01020  
01023 void FEM_Mesh::n2n_add(int n, int newNode) 
01024 {
01025   if (n == -1) return;
01026   if(FEM_Is_ghost_index(n)){
01027         FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_add");
01028         CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01029         ElemID nn(0, newNode);
01030         CkVec<ElemID> &nsVec = nVec[FEM_To_ghost_index(n)];
01031         nsVec.push_back(nn);
01032   }
01033   else{
01034         FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_add");
01035         CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01036         ElemID nn(0, newNode);
01037         CkVec<ElemID> &nsVec = nVec[n];
01038         nsVec.push_back(nn);
01039   }
01040 }
01041 
01042 
01043 
01044 
01047 void FEM_Mesh::n2n_remove(int n, int oldNode) 
01048 {
01049   if (n == -1) return;
01050   if(FEM_Is_ghost_index(n)){
01051         FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_remove");
01052         CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01053         CkVec<ElemID> &nsVec = nVec[FEM_To_ghost_index(n)];
01054         for (int i=0; i<nsVec.length(); i++) {
01055           if (nsVec[i].getSignedId() == oldNode) {
01056                 nsVec.remove(i);
01057                 break;
01058           }
01059         }
01060   }
01061   else {
01062         FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_remove");
01063         CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01064         CkVec<ElemID> &nsVec = nVec[n];
01065         for (int i=0; i<nsVec.length(); i++) {
01066           if (nsVec[i].getSignedId() == oldNode) {
01067                 nsVec.remove(i);
01068                 break;
01069           }
01070         }
01071   }
01072 }
01073 
01076 int FEM_Mesh::n2n_exists(int n, int queryNode) 
01077 {
01078   if (n == -1) return 0;
01079   if(FEM_Is_ghost_index(n)){
01080         CkAssert(node.getGhost());
01081         FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_exists");
01082         CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01083         CkVec<ElemID> &nsVec = nVec[FEM_To_ghost_index(n)];
01084         for (int i=0; i<nsVec.length(); i++)
01085           if (nsVec[i].getSignedId() == queryNode) 
01086                 return 1;
01087   }
01088   else {
01089     FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_exists");
01090     CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01091     CkVec<ElemID> &nsVec = nVec[n];
01092     for (int i=0; i<nsVec.length(); i++)
01093       if (nsVec[i].getSignedId() == queryNode) 
01094         return 1;
01095   }
01096   return 0;
01097 }
01098 
01101 void FEM_Mesh::n2n_replace(int n, int oldNode, int newNode) 
01102 {
01103   if (n == -1) return;
01104   if(FEM_Is_ghost_index(n)){
01105         FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_replace");
01106         CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01107         CkVec<ElemID> &nsVec = nVec[FEM_To_ghost_index(n)];
01108         for (int i=0; i<nsVec.length(); i++) {
01109           if (nsVec[i].getSignedId() == oldNode) {
01110                 nsVec[i] = ElemID(0,newNode);
01111                 break;
01112           }
01113         }
01114   }
01115   else {
01116         FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_replace");
01117         CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01118         CkVec<ElemID> &nsVec = nVec[n];
01119         for (int i=0; i<nsVec.length(); i++) {
01120           if (nsVec[i].getSignedId() == oldNode) {
01121                 nsVec[i] = ElemID(0,newNode);
01122                 break;
01123           }
01124         }
01125 
01126   }
01127 }
01128 
01131 void FEM_Mesh::n2n_removeAll(int n)
01132 {
01133   if (n == -1) return;
01134   if(FEM_Is_ghost_index(n)){
01135         FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_NODE_ADJACENCY,"n2n_removeAll");  
01136         CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01137         CkVec<ElemID> &nsVec = nVec[FEM_To_ghost_index(n)];
01138         nsVec.free();
01139   }
01140   else{
01141         FEM_VarIndexAttribute *nAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_NODE_ADJACENCY,"n2n_removeAll");  
01142         CkVec<CkVec<ElemID> > &nVec = nAdj->get();
01143         CkVec<ElemID> &nsVec = nVec[n];
01144         nsVec.free();
01145   }
01146 }
01147 
01149 int FEM_Mesh::n2e_getLength(int n) {
01150   if (n == -1) {
01151     return 0;
01152   }
01153   FEM_VarIndexAttribute *eAdj;
01154   if(FEM_Is_ghost_index(n)){
01155     eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_getLength");
01156     n = FEM_To_ghost_index(n);
01157   }
01158   else {
01159     eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_getLength");
01160   }
01161   CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01162   CkVec<ElemID> &nsVec = eVec[n];
01163   return nsVec.length();
01164 }
01165 
01169 ElemID FEM_Mesh::n2e_getElem(int n, int whichIdx){
01170 
01171   if(FEM_Is_ghost_index(n)){
01172     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_getElem");
01173     CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01174     CkVec<ElemID> &nsVec = eVec[FEM_To_ghost_index(n)];
01175     assert(whichIdx < nsVec.length());
01176     return  nsVec[whichIdx];
01177     
01178   }
01179   else {
01180     FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_getAll");
01181     CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01182     CkVec<ElemID> &nsVec = eVec[n];
01183     assert(whichIdx < nsVec.length());
01184     return  nsVec[whichIdx];
01185   }
01186   
01187 }
01188 
01189 
01190 
01195 void FEM_Mesh::n2e_getAll(int n, int *&adjelements, int &sz) 
01196 {
01197   if (n == -1) {
01198     sz = 0;
01199     return;
01200   }
01201   if(FEM_Is_ghost_index(n)){
01202         FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_getAll");  
01203         CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01204         CkVec<ElemID> &nsVec = eVec[FEM_To_ghost_index(n)];
01205         sz = nsVec.length();
01206         if(sz > 0) adjelements = new int[sz];
01207         for (int i=0; i<sz; i++) {
01208           adjelements[i] = nsVec[i].getSignedId();
01209         }
01210   }
01211   else {
01212         FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_getAll");  
01213         CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01214         CkVec<ElemID> &nsVec = eVec[n];
01215         int len = nsVec.length();
01216         sz = len;
01217         if(sz > 0) adjelements = new int[sz];
01218         for (int i=0; i<sz; i++) {
01219           adjelements[i] = nsVec[i].getSignedId();
01220         }
01221   }
01222 }
01223  
01224 
01231 const CkVec<ElemID> &  FEM_Mesh::n2e_getAll(int n) 
01232 {
01233         assert(n!=-1); 
01234 
01235         if(FEM_Is_ghost_index(n)){
01236                 FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_getAll");  
01237                 CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01238                 return eVec[FEM_To_ghost_index(n)];
01239         }
01240         else {
01241                 FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_getAll");  
01242                 CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01243                 return eVec[n];
01244         }
01245 }
01246  
01247 
01248 
01251 void FEM_Mesh::n2e_add(int n, int newElem) 
01252 {
01253   if (n == -1) return;
01254   if(FEM_Is_ghost_index(n)){
01255         FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_add");     
01256         CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01257         CkVec<ElemID> &nsVec = eVec[FEM_To_ghost_index(n)];
01258         ElemID ne(0, newElem);
01259         nsVec.push_back(ne);
01260         int *testn2e, testn2ec;
01261         n2e_getAll(n,testn2e,testn2ec);
01262         for(int i=0; i<testn2ec; i++) {
01263           if(FEM_Is_ghost_index(testn2e[i]))
01264             CkAssert(elem[0].ghost->is_valid(FEM_From_ghost_index(testn2e[i])));
01265           else 
01266             CkAssert(elem[0].is_valid(testn2e[i]));
01267         }
01268         if(testn2ec!=0) delete[] testn2e;
01269   }
01270   else {
01271         FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_add");     
01272         CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01273         CkVec<ElemID> &nsVec = eVec[n];
01274         ElemID ne(0, newElem);
01275         nsVec.push_back(ne);
01276   }
01277 }
01278 
01281 void FEM_Mesh::n2e_remove(int n, int oldElem) 
01282 {
01283   if (n == -1) return;
01284   if(FEM_Is_ghost_index(n)){
01285         FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_remove");
01286         CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01287         CkVec<ElemID> &nsVec = eVec[FEM_To_ghost_index(n)];
01288         for (int i=0; i<nsVec.length(); i++) {
01289           if (nsVec[i].getSignedId() == oldElem) {
01290                 nsVec.remove(i);
01291                 break;
01292           }
01293         }
01294   }
01295   else {
01296         FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_remove");
01297         CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01298         CkVec<ElemID> &nsVec = eVec[n];
01299         for (int i=0; i<nsVec.length(); i++) {
01300           if (nsVec[i].getSignedId() == oldElem) {
01301                 nsVec.remove(i);
01302                 break;
01303           }
01304         }
01305   }
01306 }
01307  
01310 void FEM_Mesh::n2e_replace(int n, int oldElem, int newElem) 
01311 {
01312   if (n == -1) return;
01313   if(FEM_Is_ghost_index(n)){
01314         FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_replace");
01315         CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01316         CkVec<ElemID> &nsVec = eVec[FEM_To_ghost_index(n)];
01317         for (int i=0; i<nsVec.length(); i++) {
01318           if (nsVec[i].getSignedId() == oldElem) {
01319                 nsVec[i] = ElemID(0,newElem);
01320                 break;
01321           }
01322         }
01323         int *testn2e, testn2ec;
01324         n2e_getAll(n,testn2e,testn2ec);
01325         for(int i=0; i<testn2ec; i++) {
01326           if(FEM_Is_ghost_index(testn2e[i]))
01327             CkAssert(elem[0].ghost->is_valid(FEM_From_ghost_index(testn2e[i])));
01328           else 
01329             CkAssert(elem[0].is_valid(testn2e[i]));
01330         }
01331         if(testn2ec!=0) delete[] testn2e;
01332   }
01333   else{
01334         FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_replace");
01335         CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01336         CkVec<ElemID> &nsVec = eVec[n];
01337         for (int i=0; i<nsVec.length(); i++) {
01338           if (nsVec[i].getSignedId() == oldElem) {
01339                 nsVec[i] = ElemID(0,newElem);
01340                 break;
01341           }
01342         }
01343   }
01344 }
01345 
01348 void FEM_Mesh::n2e_removeAll(int n)
01349 {
01350   if (n == -1) return;
01351   if(FEM_Is_ghost_index(n)){
01352         FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.getGhost()->lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_removeAll");
01353         CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01354         CkVec<ElemID> &nsVec = eVec[FEM_To_ghost_index(n)];
01355         nsVec.free();
01356   }
01357   else{
01358         FEM_VarIndexAttribute *eAdj = (FEM_VarIndexAttribute *)node.lookup(FEM_NODE_ELEM_ADJACENCY,"n2e_removeAll");
01359         CkVec<CkVec<ElemID> > &eVec = eAdj->get();
01360         CkVec<ElemID> &nsVec = eVec[n];
01361         nsVec.free();
01362   }
01363 }
01364 
01373 int FEM_Mesh::getElementOnEdge(int n1, int n2) 
01374 {
01375   int *n1AdjElems, *n2AdjElems;
01376   int n1NumElems, n2NumElems;
01377   n2e_getAll(n1, n1AdjElems, n1NumElems);
01378   n2e_getAll(n2, n2AdjElems, n2NumElems);
01379   int ret = -1;
01380   //CkPrintf("%d has %d neighboring elements, %d has %d\n", n1, n1NumElems, n2, n2NumElems);
01381 
01382   bool flag = false;
01383   for (int i=0; i<n1NumElems; i++) {
01384     for (int j=0; j<n2NumElems; j++) {
01385       if (n1AdjElems[i] == n2AdjElems[j]) {
01386         if(n1AdjElems[i] >= 0) {
01387           ret = n1AdjElems[i];
01388           flag = true;
01389           break;
01390         }
01391         else {
01392           ret = n1AdjElems[i];
01393         }
01394       }
01395     }
01396     if(flag) break;
01397   }
01398   delete[] n1AdjElems;
01399   delete[] n2AdjElems;
01400   return ret; //preferably return a local element, otherwise return a ghost 
01401 }
01402 
01403 
01408 void FEM_Mesh::get2ElementsOnEdge(int n1, int n2, int *result_e1, int *result_e2) 
01409 {
01410   int *n1AdjElems=0, *n2AdjElems=0;
01411   int n1NumElems, n2NumElems;
01412 
01413   if(n1==n2){
01414     CkPrintf("ERROR: You called get2ElementsOnEdge() with two identical nodes %d, and %d \n", n1, n2);
01415     CkExit();
01416   }
01417 
01418   n2e_getAll(n1, n1AdjElems, n1NumElems);
01419   n2e_getAll(n2, n2AdjElems, n2NumElems);
01420   CkAssert(n1AdjElems!=0);
01421   CkAssert(n2AdjElems!=0);
01422   int found=0;
01423 
01424   *result_e1=-1;
01425   *result_e2=-1;
01426 
01427   for (int i=0; i<n1NumElems; i++) {
01428     for (int j=0; j<n2NumElems; j++) {
01429       if (n1AdjElems[i] == n2AdjElems[j]) {
01430         if(found==0){
01431           //          CkPrintf("found element1 %d\n", n1AdjElems[i]);
01432           *result_e1 = n1AdjElems[i];
01433           found++;
01434         }
01435         else if(found==1){
01436           //   CkPrintf("found element2 %d\n", n1AdjElems[i]);
01437           *result_e2 = n1AdjElems[i];
01438           found++;
01439         }
01440         else {
01441           CkPrintf("ERROR: Found a third element(%d) on edge %d,%d \n", n1AdjElems[i], n1, n2);
01442            CkExit();
01443         }
01444       }
01445     }
01446   }
01447   delete[] n1AdjElems;
01448   delete[] n2AdjElems;
01449 }
01450 
01454 int getID(FEM_Mesh *mesh, int elem){
01455     int id = 0;
01456 
01457     if(elem<0){
01458       AllocTable2d<int> &idTable = ((FEM_DataAttribute *)mesh->elem[0].getGhost()->lookup(FEM_DATA+3,"getID"))->getInt();
01459       id = idTable(0,FEM_To_ghost_index(elem));
01460     } else {
01461       AllocTable2d<int> &idTable = ((FEM_DataAttribute *)mesh->elem[0].lookup(FEM_DATA+3,"getID"))->getInt();
01462       id = idTable(0,elem);
01463     }
01464     
01465   return id;
01466 }
01467 
01468 
01476 void FEM_Mesh::get2ElementsOnEdgeSorted(int n1, int n2, int *result_e1, int *result_e2)
01477 {
01478   int *n1AdjElems=0, *n2AdjElems=0;
01479   int n1NumElems, n2NumElems;
01480   
01481   get2ElementsOnEdge(n1, n2, result_e1, result_e2);
01482 
01483   //   Swap the two elements if the first has greater id than the second 
01484     if(getID(this,*result_e1) > getID(this,*result_e2)){
01485     int temp = *result_e1;
01486     *result_e1 = *result_e2;
01487     *result_e2 = temp;
01488   }
01489 }
01490 
01491 
01492 
01493 
01495 int FEM_Mesh::countElementsOnEdge(int n1, int n2) {
01496   if (n1==n2) {
01497     CkPrintf("ERROR: You called countElementsOnEdge() with two identical nodes %d, and %d \n", n1, n2);
01498     CkExit();
01499   }
01500   
01501   int *n1AdjElems=0, *n2AdjElems=0;
01502   int n1NumElems, n2NumElems;
01503 
01504   CkAssert(node.is_valid_any_idx(n1));
01505   CkAssert(node.is_valid_any_idx(n2));
01506 
01507   n2e_getAll(n1, n1AdjElems, n1NumElems);
01508   n2e_getAll(n2, n2AdjElems, n2NumElems);
01509   CkAssert(n1AdjElems!=0);
01510   CkAssert(n2AdjElems!=0);
01511   int count=0;
01512 
01513 
01514   for (int i=0; i<n1NumElems; i++) {
01515     for (int j=0; j<n2NumElems; j++) {
01516       if (n1AdjElems[i] == n2AdjElems[j]) {
01517         count++;
01518       }
01519     }
01520   }
01521   delete[] n1AdjElems;
01522   delete[] n2AdjElems;
01523 
01524   return count;
01525 }
01526 
01527
libs/ck-libs/ParFUM/mesh_adjacency.C
