00001 
00002 
00003 
00004 
00005 
00006 #include "fem_adapt_lock.h"
00007 #include "fem_mesh_modify.h"
00008 
00009 
00010 #define ERVAL -1000000000  //might cause a problem if there are 100million nodes
00011 #define ERVAL1 -1000000001
00012 
00013 int FEM_AdaptL::lockNodes(int *gotlocks, int *lockrnodes, int numRNodes, int *lockwnodes, int numWNodes) {
00014   bool donelocks = false;
00015   int numNodes = numRNodes + numWNodes;
00016   for(int i=0; i<numNodes; i++) gotlocks[i] = 0;
00017   int tryCounts=0;
00018   int ret = 1;
00019 
00020   while(!donelocks) {
00021     for(int i=0; i<numRNodes; i++) {
00022       if(lockrnodes[i]==-1) {
00023     gotlocks[i] = 1;
00024     continue;
00025       }
00026       if(FEM_Is_ghost_index(lockrnodes[i])) {
00027     if(!theMesh->node.ghost->is_valid(FEM_To_ghost_index(lockrnodes[i]))) {
00028       ret = -2;
00029       break;
00030     }
00031       } else {
00032     if(!theMesh->node.is_valid(lockrnodes[i])) {
00033       ret = -2;
00034       break;
00035     }
00036       }
00037       if(gotlocks[i]<=0) gotlocks[i] = FEM_Modify_LockN(theMesh, lockrnodes[i], 1);
00038     }
00039     for(int i=0; i<numWNodes; i++) {
00040       if(lockwnodes[i]==-1) {
00041     gotlocks[numRNodes+i] = 1;
00042     continue;
00043       }
00044       if(FEM_Is_ghost_index(lockwnodes[i])) {
00045     if(!theMesh->node.ghost->is_valid(FEM_To_ghost_index(lockwnodes[i]))) {
00046 #ifdef DEBUG_LOCKS
00047       CkPrintf("[%d] Trying to lock invalid ghost %d\n",theMod->idx,lockwnodes[i]);
00048 #endif
00049       ret = -2;
00050       break;
00051     }
00052       } else {
00053     if(!theMesh->node.is_valid(lockwnodes[i])) {
00054 #ifdef DEBUG_LOCKS
00055       CkPrintf("[%d] Trying to lock invalid node %d\n",theMod->idx,lockwnodes[i]);
00056 #endif
00057       ret = -2;
00058       break;
00059     }
00060       }
00061       if(gotlocks[numRNodes+i]<=0) {
00062     gotlocks[numRNodes+i] = FEM_Modify_LockN(theMesh, lockwnodes[i], 0);
00063       }
00064     }
00065     if(ret==-2) return ret;
00066     bool tmplocks = true;
00067     for(int i=0; i<numNodes; i++) {
00068       tmplocks = tmplocks && (gotlocks[i]>0);
00069     }
00070     if(tmplocks) {
00071       donelocks = true;
00072     }
00073     else {
00074       tryCounts++;
00075       if(tryCounts>=10) return -1;
00076       CthYield(); 
00077     }
00078   }
00079   return 1;
00080 }
00081 
00082 int FEM_AdaptL::unlockNodes(int *gotlocks, int *lockrnodes, int numRNodes, int *lockwnodes, int numWNodes) {
00083   bool donelocks = false;
00084   int numNodes = numRNodes + numWNodes;
00085   int *ungetlocks = (int*)malloc(numNodes*sizeof(int));
00086   
00087   for(int i=0; i<numRNodes; i++) {
00088     if(lockrnodes[i]==-1) {
00089       ungetlocks[i] = 1;
00090       continue;
00091     }
00092     if(FEM_Is_ghost_index(lockrnodes[i])) {
00093       if(!theMesh->node.ghost->is_valid(FEM_To_ghost_index(lockrnodes[i]))) {
00094     gotlocks[i] = -1;
00095     
00096     
00097       }
00098     } else {
00099       if(!theMesh->node.is_valid(lockrnodes[i])) {
00100     gotlocks[i] = -1;
00101     
00102     
00103       }
00104     }
00105     if(gotlocks[i]>0) ungetlocks[i] = FEM_Modify_UnlockN(theMesh, lockrnodes[i], 1);
00106     else ungetlocks[i] = 1;
00107   }
00108   for(int i=0; i<numWNodes; i++) {
00109     if(lockwnodes[i]==-1) {
00110       ungetlocks[numRNodes+i] = 1;
00111       continue;
00112     }
00113     if(FEM_Is_ghost_index(lockwnodes[i])) {
00114       if(!theMesh->node.ghost->is_valid(FEM_To_ghost_index(lockwnodes[i]))) {
00115     gotlocks[i] = -1;
00116 #ifdef DEBUG_LOCKS
00117     CkPrintf("[%d] Trying to unlock invalid ghost %d\n",theMod->idx,lockwnodes[i]);
00118 #endif
00119     
00120     
00121       }
00122     } else {
00123       if(!theMesh->node.is_valid(lockwnodes[i])) {
00124     gotlocks[i] = -1;
00125 #ifdef DEBUG_LOCKS
00126     CkPrintf("[%d] Trying to unlock invalid node %d\n",theMod->idx,lockwnodes[i]);
00127 #endif
00128     
00129     
00130       }
00131     }
00132     if(gotlocks[numRNodes+i]>0) ungetlocks[numRNodes+i] = FEM_Modify_UnlockN(theMesh, lockwnodes[i], 0);
00133     else ungetlocks[numRNodes+i] = 1;
00134   }
00135   bool tmplocks = true;
00136   for(int i=0; i<numNodes; i++) {
00137     tmplocks = tmplocks && (ungetlocks[i]>0);
00138   }
00139   if(tmplocks) donelocks = true;
00140   
00141   
00142   free(ungetlocks);  
00143   return 1;
00144 }
00145 
00146 
00147 int FEM_AdaptL::edge_flip(int n1, int n2) {
00148   int e1, e1_n1, e1_n2, e1_n3, n3;
00149   int e2, e2_n1, e2_n2, e2_n3, n4;
00150   int numNodes = 4;
00151   int *locknodes = (int*)malloc(numNodes*sizeof(int));
00152   int *gotlocks = (int*)malloc(numNodes*sizeof(int));
00153   bool done = false;
00154   int isEdge = 0;
00155   int numtries = 0;
00156   bool warned = false;
00157 
00158   isEdge = findAdjData(n1, n2, &e1, &e2, &e1_n1, &e1_n2, &e1_n3, &e2_n1, &e2_n2, &e2_n3,&n3, &n4);
00159   if(isEdge == -1) {
00160     
00161     free(locknodes);
00162     free(gotlocks);
00163     return -1;
00164   }
00165   locknodes[0] = n1;
00166   locknodes[1] = n2;
00167   locknodes[2] = n3;
00168   locknodes[3] = n4;
00169   while(!done) {
00170     int gotlock = lockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00171     isEdge = findAdjData(n1, n2, &e1, &e2, &e1_n1, &e1_n2, &e1_n3, &e2_n1, &e2_n2, &e2_n3,&n3, &n4);
00172     if(isEdge == -1) {
00173       unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00174       
00175       free(locknodes);
00176       free(gotlocks);
00177       return -1;
00178     }
00179     if(gotlock==1 && locknodes[2]==n3 && locknodes[3]==n4) {
00180       done = true;
00181     }
00182     else {
00183       unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00184       locknodes[2] = n3;
00185       locknodes[3] = n4;
00186       numtries++;
00187       if(numtries>=1000) {
00188     if(!warned) {
00189       
00190       warned = true;
00191     }
00192     numtries = 0;
00193       }
00194       CthYield();
00195     }
00196   }
00197   if ((e1 == -1) || (e2 == -1)) {
00198     unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00199     free(locknodes);
00200     free(gotlocks);
00201     return -1;
00202   }
00203   int ret = edge_flip_help(e1, e2, n1, n2, e1_n1, e1_n2, e1_n3, n3, n4,locknodes);
00204   
00205 
00206 
00207  
00208   unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00209 
00210   free(locknodes);
00211   free(gotlocks);
00212   return ret;
00213 }
00214 
00215 int FEM_AdaptL::edge_bisect(int n1, int n2) {
00216   int e1, e1_n1, e1_n2, e1_n3, n3;
00217   int e2, e2_n1, e2_n2, e2_n3, n4;
00218   int numNodes = 4;
00219   int *locknodes = (int*)malloc(numNodes*sizeof(int));
00220   int *gotlocks = (int*)malloc(numNodes*sizeof(int));
00221   bool done = false;
00222   int isEdge = 0;
00223   int numtries = 0;
00224   bool warned = false;
00225 
00226   isEdge = findAdjData(n1, n2, &e1, &e2, &e1_n1, &e1_n2, &e1_n3, &e2_n1, &e2_n2, &e2_n3,&n3, &n4);
00227   if(isEdge == -1) {
00228     
00229     free(locknodes);
00230     free(gotlocks);
00231     return -1;
00232   }
00233   locknodes[0] = n1;
00234   locknodes[1] = n2;
00235   locknodes[2] = n3;
00236   locknodes[3] = n4;
00237   while(!done) {
00238     int gotlock = lockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00239     isEdge = findAdjData(n1, n2, &e1, &e2, &e1_n1, &e1_n2, &e1_n3, &e2_n1, &e2_n2, &e2_n3,&n3, &n4);
00240     if(isEdge == -1) {
00241       unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00242       
00243       free(locknodes);
00244       free(gotlocks);
00245       return -1;
00246     }
00247     if(gotlock==1 && locknodes[2]==n3 && locknodes[3]==n4) {
00248       done = true;
00249     }
00250     else {
00251       unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00252       locknodes[2] = n3;
00253       locknodes[3] = n4;
00254       numtries++;
00255       if(numtries>=1000) {
00256     if(!warned) {
00257       
00258       warned = true;
00259     }
00260     numtries = 0;
00261       }
00262       CthYield();
00263     }
00264   }
00265   int ret = edge_bisect_help(e1, e2, n1, n2, e1_n1, e1_n2, e1_n3, e2_n1, e2_n2, e2_n3, n3, n4);
00266   unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00267 
00268   free(locknodes);
00269   free(gotlocks);
00270   return ret;
00271 }
00272 
00273 int FEM_AdaptL::vertex_remove(int n1, int n2) {
00274   int e1, e1_n1, e1_n2, e1_n3, n3;
00275   int e2, e2_n1, e2_n2, e2_n3, n4;
00276   int numNodes = 5;
00277   int numElems = 2;
00278   int *locknodes = (int*)malloc(numNodes*sizeof(int));
00279   int *gotlocks = (int*)malloc(numNodes*sizeof(int));
00280   bool done = false;
00281   int isEdge = 0;
00282   int numtries = 0;
00283   bool warned = false;
00284 
00285   isEdge = findAdjData(n1, n2, &e1, &e2, &e1_n1, &e1_n2, &e1_n3, &e2_n1, &e2_n2, &e2_n3,&n3, &n4);
00286   if(isEdge == -1) {
00287     
00288     free(locknodes);
00289     free(gotlocks);
00290     return -1;
00291   }
00292   if (e1 == -1) {
00293     free(locknodes);
00294     free(gotlocks);
00295     return 0;
00296   }
00297   
00298   int *nbrNodes, nnsize, n5;
00299   theMesh->n2n_getAll(n1, &nbrNodes, &nnsize);
00300   if(!(nnsize == 4 || (nnsize==3 && e2==-1))) {
00301     
00302     free(locknodes);
00303     free(gotlocks);
00304     return -1;    
00305   }
00306   for (int i=0; i<nnsize; i++) {
00307     if ((nbrNodes[i] != n2) && (nbrNodes[i] != n3) && (nbrNodes[i] != n4)) {
00308       n5 = nbrNodes[i];
00309       break;
00310     }
00311   }
00312   free(nbrNodes);
00313   locknodes[0] = n1;
00314   locknodes[1] = n2;
00315   locknodes[2] = n3;
00316   locknodes[3] = n4;
00317   locknodes[4] = n5;
00318   while(!done) {
00319     int gotlock = lockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00320     isEdge = findAdjData(n1, n2, &e1, &e2, &e1_n1, &e1_n2, &e1_n3, &e2_n1, &e2_n2, &e2_n3,&n3, &n4);
00321     if(isEdge == -1) {
00322       unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00323       
00324       free(locknodes);
00325       free(gotlocks);
00326       return -1;
00327     }
00328     if (e1 == -1) {
00329       unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00330       free(locknodes);
00331       free(gotlocks);
00332       return 0;
00333     }
00334     
00335     int *nbrNodes, nnsize, n5;
00336     theMesh->n2n_getAll(n1, &nbrNodes, &nnsize);
00337     for (int i=0; i<nnsize; i++) {
00338       if ((nbrNodes[i] != n2) && (nbrNodes[i] != n3) && (nbrNodes[i] != n4)) {
00339     n5 = nbrNodes[i];
00340     break;
00341       }
00342     }
00343     free(nbrNodes);
00344     if(!(nnsize == 4 || (nnsize==3 && e2==-1))) {
00345       unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00346       
00347       free(locknodes);
00348       free(gotlocks);
00349       return -1;    
00350     }
00351     if(gotlock==1 && locknodes[2]==n3 && locknodes[3]==n4 && locknodes[4]==n5) {
00352       done = true;
00353     }
00354     else {
00355       unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00356       locknodes[2] = n3;
00357       locknodes[3] = n4;
00358       locknodes[4] = n5;
00359       numtries++;
00360       if(numtries>=1000) {
00361     if(!warned) {
00362       
00363       warned = true;
00364     }
00365     numtries = 0;
00366       }
00367       CthYield();
00368     }
00369   }
00370   int ret = vertex_remove_help(e1, e2, n1, n2, e1_n1, e1_n2, e1_n3, e2_n1, e2_n2, e2_n3, n3, n4, n5);
00371   unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00372 
00373   free(locknodes);
00374   free(gotlocks);
00375   return ret;
00376 }
00377 
00378 
00379 
00380 
00381 
00382 
00383 
00384 
00385 
00386 
00387 
00388 
00389 
00390 
00391 
00392 
00393 
00394 
00395 
00396 
00397 
00398 
00399 int FEM_AdaptL::edge_contraction(int n1, int n2) {
00400   int e1, e1_n1, e1_n2, e1_n3, n3;
00401   int e2, e2_n1, e2_n2, e2_n3, n4;
00402   int numNodes = 4;
00403   int numElems = 2;
00404   int *locknodes = (int*)malloc(numNodes*sizeof(int));
00405   int *gotlocks = (int*)malloc(numNodes*sizeof(int));
00406   bool done = false;
00407   int isEdge = 0;
00408   int numtries = 0;
00409   int ret = -1;
00410   bool warned = false;
00411   int newe1=-1, newe2=-1;
00412   int acquirecount=0;
00413 
00414   bool invalidcoarsen = false;
00415   if(n1<0 || n2<0) {
00416     invalidcoarsen = true;
00417   }
00418   if(!theMesh->node.is_valid(n1)) {
00419     invalidcoarsen = true;
00420   }
00421   if(!theMesh->node.is_valid(n2)) {
00422     invalidcoarsen = true;
00423   }
00424   if(invalidcoarsen) {
00425     free(locknodes);
00426     free(gotlocks);
00427     CkPrintf("Warning: Trying to coarsen invalid edge %d - %d\n",n1,n2);
00428     return -1; 
00429   }
00430   isEdge = findAdjData(n1, n2, &e1, &e2, &e1_n1, &e1_n2, &e1_n3, &e2_n1, &e2_n2, &e2_n3,&n3, &n4);
00431   if(isEdge == -1) {
00432     CkPrintf("Edge Contract %d->%d not done as it is no longer a valid edge\n",n1,n2);
00433     free(locknodes);
00434     free(gotlocks);
00435     return -1;
00436   }
00437   locknodes[0] = n1;
00438   locknodes[1] = n2;
00439   locknodes[2] = n3;
00440   locknodes[3] = n4;
00441   if (e1 == -1) {
00442     free(locknodes);
00443     free(gotlocks);
00444     return -1;
00445   }
00446   bool locked;
00447   int gotlock = 0;
00448   while(!done) {
00449     if(ret==ERVAL1) {
00450       acquirecount++;
00451       ret=-1; 
00452       
00453       int *newnodes = new int[4];
00454       int newcount=0;
00455       if(newe1!=-1 && newe2!=-1) {
00456     int *e1conns = new int[3];
00457     int *e2conns = new int[3];
00458     theMesh->e2n_getAll(newe1,e1conns,0);
00459     theMesh->e2n_getAll(newe2,e2conns,0);
00460     for(int i=0; i<3; i++) {
00461       for(int j=0; j<3; j++) {
00462         if(e1conns[i] == e2conns[j]) {
00463           newnodes[newcount++] = e1conns[i];
00464           break;
00465         }
00466       }
00467     }
00468     if(newcount!=2) {
00469       
00470       
00471       CkAssert(FEM_Is_ghost_index(newe1) && FEM_Is_ghost_index(newe2));
00472       if(newe1!=e1) {
00473         theMesh->e2n_getAll(newe1,newnodes,0);
00474       }
00475       if(newe2!=e2) {
00476         theMesh->e2n_getAll(newe2,newnodes,0);
00477       }
00478       newcount=3;
00479       
00480       for(int i=0; i<4; i++) {
00481         bool othernodevalid = true;
00482         for(int j=0; j<3; j++) {
00483           if(locknodes[i] == newnodes[j]) othernodevalid = false;
00484         }
00485         if(othernodevalid && FEM_Is_ghost_index(locknodes[i])) {
00486           newnodes[newcount++]=locknodes[i];
00487         }
00488       }
00489     }
00490     else {
00491       for(int i=0; i<3; i++) {
00492         if(e1conns[i]!=newnodes[0] && e1conns[i]!=newnodes[1]) {
00493           newnodes[newcount++] = e1conns[i];
00494           break;
00495         }
00496       }
00497       CkAssert(newcount==3);
00498       for(int i=0; i<3; i++) {
00499         if(e2conns[i]!=newnodes[0] && e2conns[i]!=newnodes[1]) {
00500           newnodes[newcount++] = e2conns[i];
00501           break;
00502         }
00503       }
00504       CkAssert(newcount==4);
00505     }
00506     delete [] e1conns;
00507     delete [] e2conns;
00508       }
00509       else{
00510     if(newe1!=-1) {
00511       theMesh->e2n_getAll(newe1,newnodes,0);
00512       newcount=3;
00513     }
00514     else if(newe2!=-1) {
00515       theMesh->e2n_getAll(newe2,newnodes,0);
00516       newcount=3;
00517     }
00518       }
00519       e1 = newe1;
00520       e2 = newe2;
00521       n1 = newnodes[0];
00522       n2 = newnodes[1];
00523       n3 = newnodes[2];
00524       n4 = newnodes[3];
00525       locknodes[0] = n1;
00526       locknodes[1] = n2;
00527       locknodes[2] = n3;
00528       locknodes[3] = n4;
00529       numNodes = newcount;
00530       delete [] newnodes;
00531       if(numNodes!=0) {
00532     isEdge = findAdjData(n1, n2, &e1, &e2, &e1_n1, &e1_n2, &e1_n3, &e2_n1, &e2_n2, &e2_n3,&n3, &n4);
00533       } else {
00534     isEdge=-1;
00535       }
00536       if(isEdge == -1 || acquirecount>=2) {
00537     if(locked) {
00538       
00539       
00540       unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00541       locked = false;
00542     }
00543     if(acquirecount>2) CkPrintf("Edge contract %d->%d not done as it is causes an acquire livelock\n",n1,n2);
00544     CkPrintf("Edge contract %d->%d not done as it is no longer a valid edge\n",n1,n2);
00545     free(locknodes);
00546     free(gotlocks);
00547     return -1;
00548       }
00549       if (e1==-1 || e2==-1) {
00550     unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00551     free(locknodes);
00552     free(gotlocks);
00553     return -1;
00554       }
00555       if(locknodes[2]==n4 && locknodes[3]==n3) {
00556     n3 = locknodes[2];
00557     n4 = locknodes[3];
00558       }
00559       gotlock = 1;
00560     }
00561     else {
00562       gotlock = lockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00563     }
00564     locked = true;
00565     isEdge = findAdjData(n1, n2, &e1, &e2, &e1_n1, &e1_n2, &e1_n3, &e2_n1, &e2_n2, &e2_n3,&n3, &n4);
00566     if(isEdge == -1) {
00567       if(locked) {
00568     unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00569     locked = false;
00570       }
00571       CkPrintf("Edge contract %d->%d not done as it is no longer a valid edge\n",n1,n2);
00572       free(locknodes);
00573       free(gotlocks);
00574       return -1;
00575     }
00576     
00577 
00578 
00579 
00580 
00581 
00582 
00583 
00584 
00585 
00586 
00587     if(gotlock==1 && locknodes[2]==n3 && locknodes[3]==n4) {
00588       int numtries1=0;
00589       ret = ERVAL;
00590       while(ret==ERVAL && numtries1 < 5) {
00591     newe1=e1; newe2=e2;
00592     ret = edge_contraction_help(&newe1, &newe2, n1, n2, e1_n1, e1_n2, e1_n3, e2_n1, e2_n2, e2_n3, n3, n4);
00593     if(ret == ERVAL1) {
00594       done = false;
00595     }
00596     else if(ret == ERVAL) {
00597       numtries1++;
00598       if(numtries1 >= 5) {
00599         locknodes[2] = n3;
00600         locknodes[3] = n4;
00601         if(locked) {
00602           unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00603           locked = false;
00604         }
00605         numtries+=10;
00606       }
00607       CthYield();
00608     }
00609     else {
00610       done = true;
00611     }
00612       }
00613       if(numtries>=50) {
00614     
00615     
00616     isEdge = findAdjData(n1, n2, &e1, &e2, &e1_n1, &e1_n2, &e1_n3, &e2_n1, &e2_n2, &e2_n3,&n3, &n4);
00617     if(isEdge!=-1) {
00618       locknodes[2] = n3;
00619       locknodes[3] = n4;
00620     }
00621     if(locked) {
00622       unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00623       locked = false;
00624     }
00625     free(locknodes);
00626     free(gotlocks);
00627     return -1;
00628       }
00629     }
00630     else {
00631       if(locked) {
00632     unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00633     locked = false;
00634       }
00635       locknodes[2] = n3;
00636       locknodes[3] = n4;
00637       numtries++;
00638       if(numtries>=50) {
00639     if(!warned) {
00640       
00641       warned = true;
00642     }
00643         
00644     free(locknodes);
00645     free(gotlocks);
00646     return -1;
00647       }
00648       CthYield();
00649     }
00650   }
00651   if(locked) {
00652     int deletednode = -1;
00653     if(ret==n1) deletednode = n2;
00654     else if(ret==n2) deletednode = n1;
00655     if(deletednode!=-1) {
00656       for(int i=0; i<numNodes;i++) {
00657     if(locknodes[i]==deletednode) locknodes[i]=-1;
00658       }
00659     }
00660     unlockNodes(gotlocks, locknodes, 0, locknodes, numNodes);
00661     locked = false;
00662   }
00663   free(locknodes);
00664   free(gotlocks);
00665   return ret;
00666 }
00667 
00668 
00669 int FEM_AdaptL::edge_contraction_help(int *e1P, int *e2P, int n1, int n2, int e1_n1,
00670                      int e1_n2, int e1_n3, int e2_n1, 
00671                      int e2_n2, int e2_n3, int n3, int n4)
00672 {
00673   int e1=*e1P, e2=*e2P;
00674   
00675   int e1chunk=-1, e2chunk=-1;
00676   int index = theMod->getIdx();
00677 
00678   
00679   int n1_shared=0, n2_shared=0;
00680   bool same = true;
00681   n1_shared = theMod->getfmUtil()->isShared(n1);
00682   n2_shared = theMod->getfmUtil()->isShared(n2);
00683   if(n1_shared && n2_shared) {
00684     const IDXL_Rec *irec1 = theMesh->node.shared.getRec(n1);
00685     const IDXL_Rec *irec2 = theMesh->node.shared.getRec(n2);
00686     if(irec1->getShared() == irec2->getShared()) {
00687       for(int i=0; i<irec1->getShared(); i++) {
00688     same = false;
00689     for(int j=0; j<irec2->getShared(); j++) {
00690       if(irec1->getChk(i) == irec2->getChk(j)) {
00691         same = true; break;
00692       }
00693     }
00694     if(!same) break;
00695       }
00696     }
00697     else {
00698       same = false;
00699     }
00700     if(!same) {
00701       return -1;
00702     }
00703   }
00704   int n1_bound, n2_bound;
00705   FEM_Mesh_dataP(theMesh, FEM_NODE, FEM_BOUNDARY, &n1_bound, n1, 1 , FEM_INT, 1);
00706   FEM_Mesh_dataP(theMesh, FEM_NODE, FEM_BOUNDARY, &n2_bound, n2, 1 , FEM_INT, 1);
00707   if((n1_bound != 0) && (n2_bound != 0) && (n1_bound != n2_bound)) {
00708     bool kcorner = isFixedNode(n1);
00709     bool dcorner = isFixedNode(n2);
00710     bool edgeb = isEdgeBoundary(n1,n2);
00711     if((kcorner && !dcorner && edgeb) || (dcorner && !kcorner && edgeb)) {
00712     }
00713     else {
00714       return -1;
00715     }
00716   }
00717 
00718   
00719   int e1new=-1;
00720   int e2new=-1;
00721   if(e1>=0) {
00722     int e1conn[3];
00723     theMesh->e2n_getAll(e1, e1conn, 0);
00724     
00725     
00726     
00727 
00728 
00729 
00730 
00731 
00732 
00733 
00734 
00735 
00736 
00737 
00738 
00739 
00740 
00741 
00742     for(int i=0; i<3; i++) {
00743       
00744     int *e1Elems, esize1;
00745     theMesh->n2e_getAll(e1conn[i], &e1Elems, &esize1);
00746     int e1count=0;
00747     for(int j=0; j<esize1; j++) {
00748       if(e1Elems[j]>=0) e1count++;
00749     }
00750     if(e1count==1) {
00751       int *elemNbrs = new int[3];
00752       int e1ghostelem=-1;
00753       theMesh->e2e_getAll(e1,elemNbrs,0);
00754       for(int k=0; k<3; k++) {
00755         if(elemNbrs[k] < -1) {
00756           e1ghostelem = elemNbrs[k];
00757           break;
00758         }
00759       }
00760       delete [] elemNbrs;
00761       if(e1ghostelem<-1) {
00762         
00763         int e1remoteChk = theMesh->elem[0].ghost->ghostRecv.getRec(FEM_To_ghost_index(e1ghostelem))->getChk(0);
00764         int sharedIdx = theMod->fmUtil->exists_in_IDXL(theMesh,e1,e1remoteChk,3);
00765         CkPrintf("[%d]Edge Contraction, edge %d->%d, chunk %d eating into chunk %d\n",theMod->idx, n1, n2, e1remoteChk, index);
00766         if(FEM_Is_ghost_index(e2)) { 
00767           int e2conn[3];
00768           theMesh->e2n_getAll(e2, e2conn, 0);
00769           int gotlocksn4 = 1, lockn4=-1;
00770           for(int k=0; k<3; k++) {
00771         if(e2conn[k]!=n1 && e2conn[k]!=n2) lockn4=e2conn[k];
00772           }
00773           
00774           
00775           int *n4ns, n4ncount;
00776           bool shouldbeunlocked=true;
00777           theMesh->n2n_getAll(lockn4, &n4ns, &n4ncount);
00778           for(int k=0; k<n4ncount; k++) {
00779         if(n4ns[k]>=0 && n4ns[k]!=n1 && n4ns[k]!=n2) shouldbeunlocked=false; 
00780           }
00781           if(shouldbeunlocked) {
00782         unlockNodes(&gotlocksn4, &lockn4, 0, &lockn4, 1);
00783           }
00784           if(n4ncount!=0) delete[] n4ns;
00785         }
00786         e1new = meshMod[e1remoteChk].eatIntoElement(index,sharedIdx)->i;
00787         if(e1new!=-1) {
00788           e1 = theMod->fmUtil->lookup_in_IDXL(theMesh,e1new,e1remoteChk,4);
00789           
00790           e1 = FEM_To_ghost_index(e1);
00791         }
00792         else e1 = -1;
00793         free(e1Elems);
00794         *e1P = e1;
00795         return ERVAL1;
00796       }
00797     }
00798     
00799     
00800     }
00801   }
00802   if(e2>=0) {
00803     int e2conn[3];
00804     int e2remoteChk=-1;
00805     theMesh->e2n_getAll(e2, e2conn, 0);
00806     
00807     
00808     
00809 
00810 
00811 
00812 
00813 
00814 
00815 
00816 
00817 
00818 
00819 
00820 
00821 
00822 
00823 
00824     for(int i=0; i<3; i++) {
00825       
00826     int *e2Elems, esize2;
00827     theMesh->n2e_getAll(e2conn[i], &e2Elems, &esize2);
00828     int e2count=0;
00829     for(int j=0; j<esize2; j++) {
00830       if(e2Elems[j]>=0) e2count++;
00831     }
00832     if(e2count==1) {
00833       int *elemNbrs = new int[3];
00834       int e2ghostelem=-1;
00835       theMesh->e2e_getAll(e2,elemNbrs,0);
00836       for(int k=0; k<3; k++) {
00837         if(elemNbrs[k] < -1) {
00838           e2ghostelem = elemNbrs[k];
00839           break;
00840         }
00841       }
00842       delete [] elemNbrs;
00843       if(e2ghostelem<-1) {
00844         
00845         int e2remoteChk = theMesh->elem[0].ghost->ghostRecv.getRec(FEM_To_ghost_index(e2ghostelem))->getChk(0);
00846         int sharedIdx = theMod->fmUtil->exists_in_IDXL(theMesh,e2,e2remoteChk,3);
00847         CkPrintf("[%d]Edge Contraction, edge %d->%d, chunk %d eating into chunk %d\n",theMod->idx, n1, n2, e2remoteChk, index);
00848         if(FEM_Is_ghost_index(e1)) {
00849           int e1conn[3];
00850           theMesh->e2n_getAll(e1, e1conn, 0);
00851           int gotlocksn4 = 1, lockn4=-1;
00852           for(int k=0; k<3; k++) {
00853         if(e1conn[k]!=n1 && e1conn[k]!=n2) lockn4=e1conn[k];
00854           }
00855           
00856           
00857           int *n4ns, n4ncount;
00858           bool shouldbeunlocked=true;
00859           theMesh->n2n_getAll(lockn4, &n4ns, &n4ncount);
00860           for(int k=0; k<n4ncount; k++) {
00861         if(n4ns[k]>=0 && n4ns[k]!=n1 && n4ns[k]!=n2) shouldbeunlocked=false; 
00862           }
00863           if(shouldbeunlocked) {
00864         unlockNodes(&gotlocksn4, &lockn4, 0, &lockn4, 1);
00865           }
00866           if(n4ncount!=0) delete[] n4ns;
00867         }
00868         e2new = meshMod[e2remoteChk].eatIntoElement(index,sharedIdx)->i;
00869         if(e2new!=-1) {
00870           e2 = theMod->fmUtil->lookup_in_IDXL(theMesh,e2new,e2remoteChk,4);
00871           
00872           e2 = FEM_To_ghost_index(e2);
00873         }
00874         else e2 = -1;
00875         *e2P = e2;
00876         free(e2Elems);
00877         return ERVAL1;
00878       }
00879       
00880       
00881       }
00882     }
00883   }
00884   if(FEM_Is_ghost_index(e1)) {
00885     
00886     int remChk = theMesh->elem[0].ghost->ghostRecv.getRec(FEM_From_ghost_index(e1))->getChk(0);
00887     int shidx = theMesh->elem[0].ghost->ghostRecv.getRec(FEM_From_ghost_index(e1))->getIdx(0);
00888     
00889     bool shouldacquire = true;
00890     int e1nbrs[3];
00891     theMesh->e2n_getAll(e1,e1nbrs,0);
00892     for(int i=0; i<3;i++) {
00893       if(FEM_Is_ghost_index(e1nbrs[i])) {
00894     shouldacquire=false; break;
00895       }
00896     }
00897     bool shouldacquire1 = meshMod[remChk].willItLose(index,shidx)->b;
00898     if(shouldacquire && shouldacquire1) {
00899       e1new = theMod->fmUtil->eatIntoElement(e1);
00900       *e1P = e1new;
00901       return ERVAL1;
00902     }
00903     else if(shouldacquire1) return -1;
00904   }
00905   if(FEM_Is_ghost_index(e2)) {
00906     
00907     int remChk = theMesh->elem[0].ghost->ghostRecv.getRec(FEM_From_ghost_index(e2))->getChk(0);
00908     int shidx = theMesh->elem[0].ghost->ghostRecv.getRec(FEM_From_ghost_index(e2))->getIdx(0);
00909     
00910     bool shouldacquire = true;
00911     int e2nbrs[3];
00912     theMesh->e2n_getAll(e2,e2nbrs,0);
00913     for(int i=0; i<3;i++) {
00914       if(FEM_Is_ghost_index(e2nbrs[i])) {
00915     shouldacquire=false; break;
00916       }
00917     }
00918     bool shouldacquire1 = meshMod[remChk].willItLose(index,shidx)->b;
00919     if(shouldacquire && shouldacquire1) {
00920       e2new = theMod->fmUtil->eatIntoElement(e2);
00921       *e2P = e2new;
00922       return ERVAL1;
00923     }
00924     else if(shouldacquire1) return -1;
00925   }
00926 
00927 
00928   int *conn = (int*)malloc(3*sizeof(int));
00929   int *adjnodes = (int*)malloc(2*sizeof(int));
00930   adjnodes[0] = n1;
00931   adjnodes[1] = n2;
00932   int *adjelems = (int*)malloc(2*sizeof(int));
00933   adjelems[0] = e1;
00934   adjelems[1] = e2;
00935 
00936   
00937 
00938   
00939   int keepnode=-1, deletenode=-1, shared=0;
00940   
00941   
00942   bool n1fixed = isFixedNode(n1);
00943   bool n2fixed = isFixedNode(n2);
00944   if(n1fixed && n2fixed) { 
00945     free(conn);
00946     free(adjnodes);
00947     free(adjelems);
00948     return -1;
00949   }
00950   if(n1_shared && n2_shared) {
00951     
00952 
00953 
00954 
00955 
00956 
00957 
00958 
00959 
00960 
00961 
00962 
00963 
00964 
00965 
00966 
00967 
00968     if(same) {
00969       if(n1fixed) {
00970     keepnode = n2;
00971     deletenode = n1;
00972       }
00973       else {
00974     keepnode = n1;
00975     deletenode = n2;
00976       }
00977       shared = 2;
00978     } else {
00979       free(conn);
00980       free(adjnodes);
00981       free(adjelems);
00982       return -1; 
00983     }
00984   }
00985   else if(n1_shared && !n2fixed) {
00986     
00987     keepnode = n1;
00988     deletenode = n2;
00989     shared = 1;
00990   } else if(n2_shared && !n1fixed) {
00991     
00992     keepnode = n2;
00993     deletenode = n1;
00994     shared = 1;
00995   } else if(!n1_shared && !n2_shared) {
00996     
00997     if(n2fixed) {
00998       keepnode = n2;
00999       deletenode = n1;
01000     }
01001     else {
01002       keepnode = n1;
01003       deletenode = n2;
01004     }
01005   }
01006   else {
01007     free(conn);
01008     free(adjnodes);
01009     free(adjelems);
01010     return -1;
01011   }
01012   
01013   FEM_Mesh_dataP(theMesh, FEM_NODE, FEM_BOUNDARY, &n1_bound, keepnode, 1 , FEM_INT, 1);
01014   FEM_Mesh_dataP(theMesh, FEM_NODE, FEM_BOUNDARY, &n2_bound, deletenode, 1 , FEM_INT, 1);
01015   
01016   FEM_Interpolate *inp = theMod->getfmInp();
01017   FEM_Interpolate::NodalArgs nm;
01018   if((n1_bound != 0) && (n2_bound != 0) && (n1_bound != n2_bound)) {
01019     bool kcorner;
01020     bool dcorner;
01021     if(keepnode==n1) {
01022       kcorner = n1fixed;
01023       dcorner = n2fixed;
01024     }
01025     else {
01026       kcorner = n2fixed;
01027       dcorner = n1fixed;
01028     }
01029     bool edgeb = isEdgeBoundary(keepnode,deletenode);
01030     if(kcorner && !dcorner && edgeb) {
01031       nm.frac = 1.0;
01032     }
01033     else if(dcorner && !kcorner && edgeb) {
01034       nm.frac = 0.0;
01035     }
01036     else {
01037       nm.frac = 0.5;
01038       free(conn);
01039       free(adjnodes);
01040       free(adjelems);
01041       return -1; 
01042     }
01043   }
01044   else if(n1_bound!=0 && n2_bound!=0) {
01045     nm.frac = 0.5;
01046     
01047     if(isFixedNode(keepnode)) {
01048       nm.frac = 1.0;
01049     } 
01050     else if(isFixedNode(deletenode)) {
01051       nm.frac = 0.0;
01052     }
01053   }
01054   else if(n1_bound != 0) { 
01055     nm.frac = 1.0;
01056   }
01057   else if(n2_bound != 0) {
01058     if(shared==2) {
01059       keepnode = n2;
01060       deletenode = n1;
01061       nm.frac = 1.0;
01062     } else {
01063       nm.frac = 0.0;
01064     }
01065   }
01066   else {
01067     nm.frac = 0.5;
01068   }
01069   nm.nodes[0] = keepnode;
01070   nm.nodes[1] = deletenode;
01071   nm.n = keepnode;
01072   nm.addNode = false;
01073 
01074   
01075   
01076 
01077 
01078 
01079 
01080 
01081 
01082 
01083 
01084 
01085 
01086 
01087 
01088 
01089   int *nbrElems, nesize, echunk;
01090   theMesh->n2e_getAll(deletenode, &nbrElems, &nesize);
01091   bool locked = false;
01092   
01093 
01094 
01095 
01096 
01097 
01098 
01099 
01100 
01101 
01102 
01103 
01104 
01105 
01106 
01107 
01108 
01109 
01110 
01111 
01112 
01113 
01114 
01115 
01116 
01117 
01118 
01119 
01120 
01121 
01122 
01123 
01124 
01125 
01126 
01127 
01128 
01129 
01130 
01131 
01132 
01133 
01134 
01135 
01136 
01137 
01138 
01139 
01140 
01141 
01142 
01143 
01144 
01145 
01146 
01147 
01148 
01149 
01150 
01151 
01152 
01153 
01154 
01155 
01156 
01157 
01158 
01159 
01160 
01161 
01162 
01163 
01164 
01165 
01166 
01167 
01168 
01169 
01170 
01171 
01172 
01173 
01174 
01175 
01176 
01177 
01178 
01179 
01180 
01181 
01182 
01183 
01184 
01185 
01186 
01187 
01188 
01189 
01190 
01191 
01192 
01193 
01194 
01195 
01196 
01197 
01198 
01199 
01200 
01201 
01202 
01203 
01204 
01205 
01206 
01207 
01208 
01209 
01210 
01211 
01212 
01213 
01214 
01215 
01216 
01217 
01218 
01219 
01220 
01221 
01222 
01223 
01224 
01225 
01226 
01227 
01228 
01229 
01230   
01231   CkVec<int> lockedNodes;
01232   int *nnNodes;
01233   int nnsize;
01234   int nncount=0;
01235   int numtries=0;
01236   bool done = false;
01237   while(!done) {
01238     lockedNodes.removeAll();
01239     nncount=0;
01240     theMesh->n2n_getAll(deletenode, &nnNodes, &nnsize);
01241     for(int i=0; i<nnsize; i++) {
01242       if(nnNodes[i]!=n1 && nnNodes[i]!=n2 && nnNodes[i]!=n3 && nnNodes[i]!=n4) {
01243     lockedNodes.push_back(nnNodes[i]);
01244     nncount++;
01245       }
01246     }
01247     int *gotlocks = new int[nncount];
01248     int *lockw = new int[nncount];
01249     for(int i=0; i<nncount; i++) {
01250       gotlocks[i]=-1;
01251       lockw[i] = lockedNodes[i];
01252     }
01253     int gotlock = lockNodes(gotlocks,lockw,0,lockw,nncount);
01254     locked = true;
01255     numtries++;
01256     if(gotlock==1) {
01257       bool isConn = true;
01258       int *nnNodes1;
01259       int nnsize1;
01260       theMesh->n2n_getAll(deletenode, &nnNodes1, &nnsize1);
01261       if(nnsize!=nnsize1) isConn=false;
01262       else {
01263     for(int i=0; i<nnsize; i++) {
01264       isConn = false;
01265       for(int j=0; j<nnsize1; j++) {
01266         if(nnNodes1[i] == nnNodes[j]) {
01267           isConn=true; break;
01268         }
01269       }
01270       if(!isConn) break;
01271     }
01272       }
01273       if(!isConn) { 
01274     unlockNodes(gotlocks,lockw,0,lockw,nncount);
01275     if(numtries>=3) {
01276       if(nnsize1!=0) free(nnNodes1);
01277       if(nnsize!=0) free(nnNodes);
01278       if(nncount!=0) delete [] lockw;
01279       if(nncount!=0) delete [] gotlocks;
01280       free(conn);
01281       free(adjnodes);
01282       free(adjelems);
01283       return ERVAL;
01284     }
01285     CthYield();
01286       }
01287       else done = true;
01288       if(nnsize1!=0) free(nnNodes1);
01289     }
01290     else unlockNodes(gotlocks,lockw,0,lockw,nncount);
01291     if(nnsize!=0) free(nnNodes);
01292     if(nncount!=0) delete [] lockw;
01293     if(nncount!=0) delete [] gotlocks;
01294     if(numtries>=3 && !done) return ERVAL;
01295     free(conn);
01296     free(adjnodes);
01297     free(adjelems);
01298   }
01299 
01300 
01301   
01302   double *n1_coord = new double[2];
01303   double *n2_coord = new double[2];
01304   double *new_coord = new double[2];
01305   FEM_Mesh_dataP(theMesh, FEM_NODE, theMod->fmAdaptAlgs->coord_attr, (void *)n1_coord, nm.nodes[0], 1, FEM_DOUBLE, 2);
01306   FEM_Mesh_dataP(theMesh, FEM_NODE, theMod->fmAdaptAlgs->coord_attr, (void *)n2_coord, nm.nodes[1], 1, FEM_DOUBLE, 2);
01307   new_coord[0] = nm.frac*n1_coord[0] + (1-nm.frac)*n2_coord[0];
01308   new_coord[1] = nm.frac*n1_coord[1] + (1-nm.frac)*n2_coord[1];
01309   int flipSliver = false;
01310   int *nbr1Elems, nesize1;
01311   int *conn1 = new int[3];
01312   theMesh->n2e_getAll(keepnode, &nbr1Elems, &nesize1);
01313   for (int i=0; i<nesize; i++) {
01314     if ((nbrElems[i] != e1) && (nbrElems[i] != e2)) {
01315       theMesh->e2n_getAll(nbrElems[i], conn);
01316       for(int j=0; j<2; j++) {
01317     if (conn[j] == deletenode) {
01318       conn[j] = conn[2];
01319       conn[2] = deletenode;
01320     }
01321       }
01322 
01323       if(theMod->fmAdaptAlgs->didItFlip(conn[0],conn[1],conn[2],new_coord)) {
01324     flipSliver = true;
01325     
01326     break;
01327       }
01328     }
01329   }
01330   if(!flipSliver) {
01331     for (int i=0; i<nesize1; i++) {
01332       if ((nbr1Elems[i] != e1) && (nbr1Elems[i] != e2)) {
01333     theMesh->e2n_getAll(nbr1Elems[i], conn1);
01334     for(int j=0; j<2; j++) {
01335       if (conn1[j] == keepnode) {
01336         conn1[j] = conn1[2];
01337         conn1[2] = keepnode;
01338       }
01339     }
01340     if(theMod->fmAdaptAlgs->didItFlip(conn1[0],conn1[1],conn1[2],new_coord)) {
01341       flipSliver = true;
01342       
01343       break;
01344     }
01345       }
01346     }
01347   }
01348   if(nesize1 != 0) delete [] nbr1Elems;
01349   delete [] conn1;
01350   delete [] n1_coord;
01351   delete [] n2_coord;
01352   delete [] new_coord;
01353   if(flipSliver) {
01354     int size = lockedNodes.size();
01355     int *gotlocks = (int*)malloc(size*sizeof(int));
01356     int *lockw = (int*)malloc(size*sizeof(int));
01357     for(int k=0; k<size; k++) {
01358       gotlocks[k] = 1;
01359       lockw[k] = lockedNodes[k];
01360     }
01361     if(locked) {
01362 #ifdef DEBUG_LOCKS
01363       CkPrintf("[%d]Trying to unlock %d, %d, & %d\n",theMod->idx,conn[0],conn[1],conn[2]);
01364 #endif
01365       unlockNodes(gotlocks, lockw, 0, lockw, size);
01366       locked = false;
01367     }
01368     free(gotlocks);
01369     free(lockw);
01370 
01371     if(nesize!=0) delete[] nbrElems;
01372     free(conn);
01373     free(adjnodes);
01374     free(adjelems);
01375     return -1;
01376   }
01377 
01378   inp->FEM_InterpolateNodeOnEdge(nm); 
01379 
01380 
01381 #ifdef DEBUG_1
01382   CkPrintf("[%d]Edge Contraction, edge %d(%d:%d)->%d(%d:%d) on chunk %d:: deleted %d\n",theMod->idx, n1,n1_bound,n1_shared, n2,n2_bound,n2_shared, theMod->getfmUtil()->getIdx(),deletenode);
01383 #endif
01384   e1chunk = FEM_remove_element(theMesh, e1, 0);
01385   if(e2!=-1) e2chunk = FEM_remove_element(theMesh, e2, 0);
01386   FEM_purge_element(theMesh,e1,0);
01387   if(e2!=-1) FEM_purge_element(theMesh,e2,0);
01388 
01389   for (int i=0; i<nesize; i++) {
01390     if ((nbrElems[i] != e1) && (nbrElems[i] != e2)) {
01391       theMesh->e2n_getAll(nbrElems[i], conn);
01392       for (int j=0; j<3; j++) {
01393     if (conn[j] == deletenode) conn[j] = keepnode;
01394       }
01395       int eTopurge = nbrElems[i];
01396       echunk = FEM_remove_element(theMesh, nbrElems[i], 0);
01397       nbrElems[i] = FEM_add_element(theMesh, conn, 3, 0, echunk); 
01398       theMod->fmUtil->copyElemData(0,eTopurge,nbrElems[i]);
01399       FEM_purge_element(theMesh,eTopurge,0);
01400     }
01401   }
01402 
01403   
01404   int size = lockedNodes.size();
01405   int *gotlocks = (int*)malloc(size*sizeof(int));
01406   int *lockw = (int*)malloc(size*sizeof(int));
01407   for(int k=0; k<size; k++) {
01408     lockw[k] = lockedNodes[k];
01409     gotlocks[k] = 1;
01410   }
01411   if(locked) {
01412 #ifdef DEBUG_LOCKS
01413     CkPrintf("[%d]Done contraction, Trying to unlock %d, %d, & %d\n",theMod->idx,conn[0],conn[1],conn[2]);
01414 #endif
01415     
01416     unlockNodes(gotlocks, lockw, 0, lockw, size);
01417     locked = false;
01418   }
01419   free(gotlocks);
01420   free(lockw);
01421 
01422   FEM_remove_node(theMesh, deletenode);
01423   if(nesize!=0) delete[] nbrElems;
01424   free(conn);
01425   free(adjnodes);
01426   free(adjelems);
01427   return keepnode;
01428 }
01429 
01430 
01431 
01432 #undef DEBUG_1
01433