Charm++: Charm Source Code Documentation

00001 /* File: fem_adapt_algs.C
00002  * Authors: Terry Wilmarth, Nilesh Choudhury
00003  * 
00004  */
00005 
00006 
00011 #include "ParFUM.h"
00012 #include "ParFUM_internals.h"
00013 
00014 
00015 #define MINAREA 1.0e-18
00016 #define MAXAREA 1.0e12
00017 
00018 #define GRADATION 1.2
00019 //#define ADAPT_VERBOSE
00020 
00021 CtvDeclare(FEM_Adapt_Algs *, _adaptAlgs);
00022 
00023 FEM_Adapt_Algs::FEM_Adapt_Algs(FEM_Mesh *m, femMeshModify *fm) 
00024 { 
00025   theMesh = m; 
00026   theMod = fm; 
00027   //theAdaptor = theMod->fmAdapt;
00028   theAdaptor = theMod->fmAdaptL;
00029 }
00030 FEM_Adapt_Algs::FEM_Adapt_Algs(femMeshModify *fm) {
00031   theMesh = NULL;
00032   theMod = fm;
00033   theAdaptor = theMod->fmAdaptL;
00034 }
00035 
00036 FEM_Adapt_Algs::~FEM_Adapt_Algs() {
00037 }
00038 
00039 
00040 
00047 void FEM_Adapt_Algs::FEM_Refine(int qm, int method, double factor, 
00048                 double *sizes)
00049 {
00050   SetMeshSize(method, factor, sizes);
00051   GradateMesh(GRADATION);
00052   (void)Refine(qm, method, factor, sizes);
00053 }
00054 
00057 int FEM_Adapt_Algs::Refine(int qm, int method, double factor, double *sizes)
00058 {
00059   // loop through elemsToRefine
00060   int elId, mods=0, iter_mods=1;
00061   int elemWidth = theMesh->elem[0].getConn().width();
00062   refineElements = refineStack = NULL;
00063   refineTop = refineHeapSize = 0;
00064   int elemConn[3];
00065   int count = 0;
00066   while (iter_mods != 0 && count<20) {
00067     count++;
00068     iter_mods=0;
00069     numNodes = theMesh->node.size();
00070     numElements = theMesh->elem[0].size();
00071     // sort elements to be refined by quality into elemsToRefine
00072     if (refineStack) delete [] refineStack;
00073     refineStack = new elemHeap[numElements];
00074     if (refineElements) delete [] refineElements;
00075     refineElements = new elemHeap[numElements+1];
00076     for (int i=0; i<numElements; i++) { 
00077       if (theMesh->elem[0].is_valid(i)) {
00078     // find maxEdgeLength of i
00079     double tmpLen, avgEdgeLength=0.0, maxEdgeLength = 0.0;
00080     theMesh->e2n_getAll(i, elemConn);
00081     bool notclear = false;
00082     for(int j=0; j<elemWidth; j++) {
00083       int nd = elemConn[j];
00084       if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
00085       //else cannot look up because connectivity might be screwed up,
00086       //I am hoping that at least one of the nodes will be local
00087     }
00088     if(notclear) continue;
00089     for (int j=0; j<elemWidth-1; j++) {
00090       for (int k=j+1; k<elemWidth; k++) {
00091         tmpLen = length(elemConn[j], elemConn[k]);
00092         if(tmpLen < -1.0) {notclear = true;}
00093         avgEdgeLength += tmpLen;
00094         if (tmpLen > maxEdgeLength) maxEdgeLength = tmpLen;
00095       }
00096     }
00097     if(notclear) continue;
00098     avgEdgeLength /= 3.0;
00099     double qFactor=getAreaQuality(i);
00100     if (theMesh->elem[0].getMeshSizing(i) <= 0.0) 
00101       CkPrintf("WARNING: mesh element %d has no sizing!\n", i);
00102     if ((theMesh->elem[0].getMeshSizing(i) > 0.0) &&
00103         (avgEdgeLength > (theMesh->elem[0].getMeshSizing(i)*REFINE_TOL))){
00104         //|| (qFactor < QUALITY_MIN)) {
00105       Insert(i, qFactor*(1.0/maxEdgeLength), 0);
00106     }
00107       }
00108     }
00109     while (refineHeapSize>0 || refineTop > 0) { // loop through the elements
00110       int n1, n2;
00111       if (refineTop>0) {
00112     refineTop--;
00113     elId=refineStack[refineTop].elID;
00114       }
00115       else  elId=Delete_Min(0);
00116       if ((elId != -1) && (theMesh->elem[0].is_valid(elId))) {
00117     int n1, n2;
00118     double tmpLen, avgEdgeLength=0.0, maxEdgeLength = 0.0;
00119     int maxEdgeIdx=0;
00120     theMesh->e2n_getAll(elId, elemConn);
00121     bool notclear = false;
00122     for(int j=0; j<elemWidth; j++) {
00123       int nd = elemConn[j];
00124       if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
00125       //else cannot look up because connectivity might be screwed up,
00126       //I am hoping that at least one of the nodes will be local
00127     }
00128     if(notclear) continue;
00129     for (int j=0; j<elemWidth-1; j++) {
00130       for (int k=j+1; k<elemWidth; k++) {
00131         tmpLen = length(elemConn[j], elemConn[k]);
00132         if(tmpLen < -1.0) {notclear = true;}
00133         avgEdgeLength += tmpLen;
00134         if (tmpLen > maxEdgeLength) { 
00135           maxEdgeLength = tmpLen;
00136           if(k==elemWidth-1) {
00137         maxEdgeIdx = j+1;
00138           }
00139           else maxEdgeIdx = j;
00140           n1 = elemConn[j]; n2 = elemConn[k];
00141         }
00142       }
00143     }
00144     if(notclear) continue;
00145     avgEdgeLength /= 3.0;
00146     //double qFactor=getAreaQuality(elId);
00147     if ((theMesh->elem[0].getMeshSizing(elId) > 0.0) &&
00148         (avgEdgeLength>(theMesh->elem[0].getMeshSizing(elId)*REFINE_TOL))){
00149       //|| (qFactor < QUALITY_MIN)) {
00150       //decide if we should do a flip or bisect
00151       bool flag = flipOrBisect(elId,n1,n2,maxEdgeIdx,maxEdgeLength);
00152       if(flag) {
00153 #ifdef DEBUG_FLIP
00154         if(theAdaptor->edge_flip(n1, n2) > 0)  iter_mods++;
00155 #endif
00156       }
00157       else {
00158         if(theAdaptor->edge_bisect(n1, n2) > 0)  iter_mods++;
00159       }
00160     }
00161       }
00162       CthYield(); // give other chunks on the same PE a chance
00163     }
00164     mods += iter_mods;
00165 #ifdef ADAPT_VERBOSE
00166     CkPrintf("[%d]ParFUM_Refine: %d modifications in last pass.\n",theMod->idx,iter_mods);
00167 #endif
00168 #ifdef DEBUG_QUALITY
00169     tests(false);
00170 #endif
00171   }
00172 #ifdef ADAPT_VERBOSE
00173   CkPrintf("[%d]ParFUM_Refine: %d total modifications.\n",theMod->idx,mods);
00174 #endif
00175   delete[] refineStack;
00176   delete[] refineElements;
00177   return mods;
00178 }
00179 
00186 void FEM_Adapt_Algs::FEM_Coarsen(int qm, int method, double factor, 
00187                  double *sizes)
00188 {
00189   SetMeshSize(method, factor, sizes);
00190   GradateMesh(GRADATION);
00191   (void)Coarsen(qm, method, factor, sizes);
00192 }
00193 
00196 int FEM_Adapt_Algs::Coarsen(int qm, int method, double factor, double *sizes)
00197 {
00198   // loop through elemsToRefine
00199   int elId, mods=0, iter_mods=1, pass=0;
00200   int elemWidth = theMesh->elem[0].getConn().width();
00201   double qFactor;
00202   coarsenElements = NULL;
00203   coarsenHeapSize = 0;
00204   int elemConn[3];
00205   while (iter_mods != 0) {
00206     iter_mods=0;
00207     pass++;
00208     numNodes = theMesh->node.size();
00209     numElements = theMesh->elem[0].size();
00210     // sort elements to be refined by quality into elemsToRefine
00211     if (coarsenElements) delete [] coarsenElements;
00212     coarsenElements = new elemHeap[numElements+1];
00213     coarsenElements[0].len=-2.0;
00214     coarsenElements[0].elID=-1;
00215     for (int i=0; i<numElements; i++) { 
00216       if (theMesh->elem[0].is_valid(i)) {
00217     // find minEdgeLength of i
00218     theMesh->e2n_getAll(i, elemConn);
00219     bool notclear = false;
00220     for(int j=0; j<elemWidth; j++) {
00221       int nd = elemConn[j];
00222       if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
00223       //else cannot look up because connectivity might be screwed up,
00224       //I am hoping that at least one of the nodes will be local
00225     }
00226     if(notclear) continue;
00227     double tmpLen, avgEdgeLength=0.0, 
00228       minEdgeLength = length(elemConn[0], elemConn[1]);
00229     for (int j=0; j<elemWidth-1; j++) {
00230       for (int k=j+1; k<elemWidth; k++) {
00231         tmpLen = length(elemConn[j], elemConn[k]);
00232         if(tmpLen < -1.0) {notclear = true;}
00233         avgEdgeLength += tmpLen;
00234         if (tmpLen < minEdgeLength) minEdgeLength = tmpLen;
00235       }
00236     }
00237     if(notclear) continue;
00238     avgEdgeLength /= 3.0;
00239     qFactor=getAreaQuality(i);
00240     if (((theMesh->elem[0].getMeshSizing(i) > 0.0) &&
00241          (avgEdgeLength < (theMesh->elem[0].getMeshSizing(i)*COARSEN_TOL)))
00242          || (qFactor < QUALITY_MIN)) {
00243         //){
00244       //CkPrintf("Marking elem %d for coarsening\n", i);
00245       Insert(i, qFactor*minEdgeLength, 1);
00246     }
00247       }
00248     }
00249     while (coarsenHeapSize>0) { // loop through the elements
00250       elId=Delete_Min(1);
00251       if ((elId != -1) && (theMesh->elem[0].is_valid(elId))) {
00252     theMesh->e2n_getAll(elId, elemConn);
00253     bool notclear = false;
00254     for(int j=0; j<elemWidth; j++) {
00255       int nd = elemConn[j];
00256       if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
00257       //else cannot look up because connectivity might be screwed up,
00258       //I am hoping that at least one of the nodes will be local
00259     }
00260     if(notclear) continue;
00261     int n1=elemConn[0], n2=elemConn[1];
00262     double tmpLen, avgEdgeLength=0.0, 
00263       minEdgeLength = length(n1, n2);
00264     for (int j=0; j<elemWidth-1; j++) {
00265       for (int k=j+1; k<elemWidth; k++) {
00266         tmpLen = length(elemConn[j], elemConn[k]);
00267         if(tmpLen < -1.0) {notclear = true;}
00268         avgEdgeLength += tmpLen;
00269         if (tmpLen < minEdgeLength) {
00270           minEdgeLength = tmpLen;
00271           n1 = elemConn[j]; n2 = elemConn[k];
00272         }
00273       }
00274     }
00275     if(notclear) continue;
00276     CkAssert(n1!=-1 && n2!=-1);
00277     avgEdgeLength /= 3.0;
00278     qFactor=getAreaQuality(elId);
00279     // coarsen element's short edge
00280     if (((theMesh->elem[0].getMeshSizing(elId) > 0.0) &&
00281          (avgEdgeLength < (theMesh->elem[0].getMeshSizing(elId)*COARSEN_TOL)))
00282           || (qFactor < QUALITY_MIN)) {
00283       //){
00284 
00285       int eNbr = theMesh->e2e_getNbr(elId, theMesh->e2n_getIndex(elId,n1)+
00286                      theMesh->e2n_getIndex(elId,n2)-1, 0);
00287       // determine if eNbr should also be coarsened
00288       if ((eNbr >= 0) && (theMesh->elem[0].is_valid(eNbr))) {
00289         theMesh->e2n_getAll(eNbr, elemConn);
00290         avgEdgeLength=0.0;
00291         for (int j=0; j<elemWidth-1; j++) {
00292           for (int k=j+1; k<elemWidth; k++) {
00293         avgEdgeLength += length(elemConn[j], elemConn[k]);
00294           }
00295         }
00296         avgEdgeLength /= 3.0;
00297         qFactor=getAreaQuality(eNbr);
00298         if (((theMesh->elem[0].getMeshSizing(eNbr) > 0.0) &&
00299          (avgEdgeLength < (theMesh->elem[0].getMeshSizing(eNbr))))
00300         || (qFactor < QUALITY_MIN)) {
00301           //CkPrintf("Coarsening elem %d which has desired edge length %6.6e and average edge length %6.6e.\n", elId, theMesh->elem[0].getMeshSizing(elId), avgEdgeLength);
00302           if (theAdaptor->edge_contraction(n1, n2) > 0)  iter_mods++;
00303         }
00304       }
00305       else {
00306         //CkPrintf("Coarsening elem %d which has desired edge length %6.6e and average edge length %6.6e.\n", elId, theMesh->elem[0].getMeshSizing(elId), avgEdgeLength);
00307         if (theAdaptor->edge_contraction(n1, n2) > 0)  iter_mods++;
00308       }
00309     }
00310       }
00311       CthYield(); // give other chunks on the same PE a chance
00312     }
00313     mods += iter_mods;
00314 #ifdef ADAPT_VERBOSE
00315     CkPrintf("[%d]ParFUM_Coarsen: %d modifications in pass %d.\n",theMod->idx,iter_mods,pass);
00316 #endif
00317 #ifdef DEBUG_QUALITY
00318     tests(false);
00319 #endif
00320   }
00321 #ifdef ADAPT_VERBOSE
00322   CkPrintf("[%d]ParFUM_Coarsen: %d total modifications over %d passes.\n",theMod->idx,mods,pass);
00323 #endif
00324   delete[] coarsenElements;
00325   return mods;
00326 }
00327 
00331 void FEM_Adapt_Algs::FEM_AdaptMesh(int qm, int method, double factor, 
00332                    double *sizes)
00333 {
00334   MPI_Comm comm=(MPI_Comm)FEM_chunk::get("FEM_Update_mesh")->defaultComm;
00335   MPI_Barrier(comm);
00336 #ifdef ADAPT_VERBOSE
00337   CkPrintf("[%d]BEGIN: FEM_AdaptMesh...\n",theMod->idx);
00338 #endif
00339 #ifdef DEBUG_QUALITY
00340   tests(true);
00341 #endif
00342   SetMeshSize(method, factor, sizes);
00343   MPI_Barrier(comm);
00344   GradateMesh(GRADATION);
00345   MPI_Barrier(comm);
00346   (void)Refine(qm, method, factor, sizes);
00347   MPI_Barrier(comm);
00348   GradateMesh(GRADATION);
00349   MPI_Barrier(comm);
00350   (void)Coarsen(qm, method, factor, sizes);
00351 #ifdef DEBUG_QUALITY
00352   MPI_Barrier(comm);
00353 #endif
00354   FEM_Repair(qm);
00355   MPI_Barrier(comm);
00356 #ifdef DEBUG_QUALITY
00357   tests(true);
00358   MPI_Barrier(comm);
00359 #endif
00360 #ifdef ADAPT_VERBOSE
00361   CkPrintf("[%d]...END: FEM_AdaptMesh.\n",theMod->idx);
00362 #endif
00363 }
00364 
00367 void FEM_Adapt_Algs::FEM_Smooth(int qm, int method)
00368 {
00369   CkPrintf("WARNING: ParFUM_Smooth: Not yet implemented.\n");
00370 }
00371 
00382 void  FEM_Adapt_Algs::FEM_mesh_smooth(FEM_Mesh *meshP, int *nodes, int nNodes, int attrNo) {
00383   vector2d newPos, *coords, *ghostCoords;
00384   int idx, nNod, nGn, gIdxN, *boundVals, nodesInChunk, mesh;
00385   int *adjnodes;
00386 
00387   mesh=FEM_Mesh_default_read();
00388   nodesInChunk = FEM_Mesh_get_length(mesh,FEM_NODE);
00389   nGn = FEM_Mesh_get_length(mesh, FEM_GHOST + FEM_NODE);
00390   
00391   boundVals = new int[nodesInChunk];
00392   coords = new vector2d[nodesInChunk+nGn];
00393 
00394   FEM_Mesh_data(mesh, FEM_NODE, FEM_BOUNDARY, (int*) boundVals, 0, nodesInChunk, FEM_INT, 1);    
00395 
00396   FEM_Mesh_data(mesh, FEM_NODE, attrNo, (double*)coords, 0, nodesInChunk, FEM_DOUBLE, 2);
00397 
00398   IDXL_Layout_t coord_layout = IDXL_Layout_create(IDXL_DOUBLE, 2);
00399   FEM_Update_ghost_field(coord_layout,-1, coords); 
00400   ghostCoords = &(coords[nodesInChunk]);
00401   for (int i=0; i<nNodes; i++)
00402   {
00403     if (nodes==NULL) idx=i;
00404     else idx=nodes[i];
00405     newPos.x=0;
00406     newPos.y=0;
00407     CkAssert(idx<nodesInChunk);  
00408     if (FEM_is_valid(mesh, FEM_NODE, idx) && boundVals[idx]==0) //node must be internal
00409     {
00410       meshP->n2n_getAll(idx, adjnodes, nNod);
00411       for (int j=0; j<nNod; j++) { //for all adjacent nodes, find coords
00412     if (adjnodes[j]<-1) {
00413       gIdxN = FEM_From_ghost_index(adjnodes[j]);
00414       newPos.x += ghostCoords[gIdxN].x;
00415       newPos.y += ghostCoords[gIdxN].y;
00416     }
00417     else {
00418       newPos.x += coords[adjnodes[j]].x;
00419       newPos.y += coords[adjnodes[j]].y;
00420     }     
00421       }
00422       newPos.x/=nNod;
00423       newPos.y/=nNod;
00424       FEM_set_entity_coord2(mesh, FEM_NODE, idx, newPos.x, newPos.y);
00425       delete [] adjnodes;
00426     }
00427   }
00428 
00429   delete [] coords;
00430   delete [] boundVals;
00431 }
00432 
00436 void FEM_Adapt_Algs::FEM_Repair(int qm)
00437 {
00438   double avgQual = 0.0, minQual = getAreaQuality(0);
00439   int numBadElems = 0;
00440   int elemConn[3];
00441 #ifdef ADAPT_VERBOSE
00442   CkPrintf("[%d]WARNING: ParFUM_Repair: Under construction.\n",theMod->idx);
00443   numElements = theMesh->elem[0].size();
00444   for (int i=0; i<numElements; i++) { 
00445     if (theMesh->elem[0].is_valid(i)) {
00446       double qFactor=getAreaQuality(i);
00447       avgQual += qFactor;
00448       if (qFactor <  QUALITY_MIN) {
00449     numBadElems++;
00450     if (qFactor < minQual) minQual = qFactor;
00451       }
00452     }
00453   }
00454   avgQual /= numElements;
00455   CkPrintf("BEFORE FEM_Repair: Average Element Quality = %2.6f, Min = %2.6f (1.0 is perfect)\n", avgQual, minQual);
00456   //CkPrintf("BEFORE FEM_Repair: Average Element Quality = %2.6f, Min = %2.6f (1.0 is perfect)\n  %d out of %d elements were below the minimum quality tolerance of %2.6f\n", avgQual, minQual, numBadElems, numElements, QUALITY_MIN);
00457 #endif
00458 
00459   int elemWidth = theMesh->elem[0].getConn().width();
00460   int changes=1, totalChanges=0;
00461   int count=0;
00462   while (changes!=0 && count<4) {
00463     count++;
00464     changes = 0;
00465     numElements = theMesh->elem[0].size();
00466     for (int i=0; i<numElements; i++) { 
00467       if (theMesh->elem[0].is_valid(i)) {
00468     double qFactor=getAreaQuality(i);
00469     if (qFactor <  0.75*QUALITY_MIN) {
00470       int elId = i;
00471       theMesh->e2n_getAll(elId, elemConn);
00472       bool notclear = false;
00473       for(int j=0; j<elemWidth; j++) {
00474         int nd = elemConn[j];
00475         if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
00476         //else cannot look up because connectivity might be screwed up,
00477         //I am hoping that at least one of the nodes will be local
00478       }
00479       if(notclear) continue;
00480       int n1=elemConn[0], n2=elemConn[1];
00481       //too bad.. should not decide without locking!!
00482       //these values might change by the time lock is acquired
00483       double len1 = length(elemConn[0],elemConn[1]);
00484       double len2 = length(elemConn[1],elemConn[2]);
00485       double len3 = length(elemConn[2],elemConn[0]);
00486       if(len1<-1.0 || len2<-1.0 || len3<-1.0) continue;
00487       double avglen=(len1+len2+len3)/3.0;
00488       int maxn1=0, maxn2=1;
00489       int minn1=0, minn2=1;
00490       double maxlen=len1;
00491       int maxed=0;
00492       if(len2>maxlen) {
00493         maxlen = len2;
00494         maxn1 = 1;
00495         maxn2 = 2;
00496         maxed = 1;
00497       }
00498       if(len3>maxlen) {
00499         maxlen = len3;
00500         maxn1 = 2;
00501         maxn2 = 0;
00502         maxed = 2;
00503       }
00504       double minlen = len1;
00505       if(len2<minlen) {
00506         minlen = len2;
00507         minn1 = 1;
00508         minn2 = 2;
00509       }
00510       if(len3<minlen) {
00511         minlen = len3;
00512         minn1 = 2;
00513         minn2 = 0;
00514       }
00515       double otherlen = 3*avglen - maxlen - minlen;
00516       if (maxlen > 0.95*(minlen+otherlen)) { // refine
00517         //decide if this should be a bisect or flip, 
00518         //depends on if the longest edge on this element is also the longest on the
00519         //element sharing this edge, if not, flip it
00520         int nbrEl = theMesh->e2e_getNbr(i,maxed);
00521         double len4=0.0, len5=0.0;
00522         if(nbrEl!=-1) {
00523           int con1[3];
00524           theMesh->e2n_getAll(nbrEl,con1);
00525           int nbrnode=-1;
00526           for(int j=0; j<3; j++) {
00527         if(con1[j]!=elemConn[maxn1] && con1[j]!=elemConn[maxn2]) {
00528           nbrnode = con1[j];
00529           break;
00530         }
00531           }
00532           len4 = length(elemConn[maxn1],nbrnode);
00533           len5 = length(elemConn[maxn2],nbrnode);
00534         }
00535         int success = -1;
00536         if(len4>maxlen || len5>maxlen) {
00537 #ifdef DEBUG_FLIP
00538           success = theAdaptor->edge_flip(elemConn[maxn1], elemConn[maxn2]);
00539 #endif
00540         }
00541         else {
00542           success = theAdaptor->edge_bisect(elemConn[maxn1], elemConn[maxn2]);
00543         }
00544         if (success >= 0) {
00545           //CkPrintf("Refined bad element!\n");
00546           changes++;
00547         }
00548       }
00549       else if (minlen < 0.10*(maxlen+otherlen)) { // coarsen
00550         int success = theAdaptor->edge_contraction(elemConn[minn1], elemConn[minn2]);
00551         if (success >= 0) { 
00552           //CkPrintf("Coarsened bad element!\n");
00553           changes++;
00554         }
00555       }
00556       else {
00557         //CkPrintf("Leaving one bad element alone...\n");
00558       }
00559     }
00560       }
00561     }
00562     totalChanges += changes;
00563   }
00564 
00565 #ifdef ADAPT_VERBOSE
00566   numElements = theMesh->elem[0].size();
00567   numBadElems = 0;
00568   avgQual = 0.0;
00569   minQual = getAreaQuality(0);
00570   for (int i=0; i<numElements; i++) { 
00571     if (theMesh->elem[0].is_valid(i)) {
00572       double qFactor=getAreaQuality(i);
00573       avgQual += qFactor;
00574       if (qFactor <  QUALITY_MIN) {
00575     numBadElems++;
00576     if (qFactor < minQual) minQual = qFactor;
00577       }
00578     }
00579   }
00580   avgQual /= numElements;
00581   CkPrintf("[%d]AFTER FEM_Repair: Average Element Quality = %2.6f, Min = %2.6f (1.0 is perfect) No. of repairs %d\n",theMod->idx,avgQual, minQual,totalChanges);
00582 #ifdef DEBUG_QUALITY
00583     tests(false);
00584 #endif
00585   //  CkPrintf("AFTER FEM_Repair: Average Element Quality = %2.6f, Min = %2.6f (1.0 is perfect)\n  %d out of %d elements were below the minimum quality tolerance of %2.6f\n", avgQual, minQual, numBadElems, numElements, QUALITY_MIN);
00586 #endif
00587 }
00588 
00594 void FEM_Adapt_Algs::FEM_Remesh(int qm, int method, double factor, 
00595                 double *sizes)
00596 {
00597   CkPrintf("WARNING: ParFUM_Remesh: Under construction.\n");
00598 }
00599 
00600 
00601 
00604 void FEM_Adapt_Algs::SetReferenceMesh()
00605 {
00606   // TODO: do we need to run this loop for element types other than 0?
00607   double avgLength = 0.0;
00608   int width = theMesh->elem[0].getConn().width();
00609   int numElements = theMesh->elem[0].size();
00610   int elemConn[3];
00611   
00612   for (int i=0; i<numElements; ++i, avgLength=0) {
00613     if(theMesh->elem[0].is_valid(i)) {
00614       theMesh->e2n_getAll(i, elemConn);
00615       for (int j=0; j<width-1; ++j) {
00616     avgLength += length(elemConn[j], elemConn[j+1]);
00617       }
00618       avgLength += length(elemConn[0], elemConn[width-1]);
00619       avgLength /= width;
00620       theMesh->elem[0].setMeshSizing(i, avgLength);      
00621     }
00622   }
00623 }
00624 
00630 void FEM_Adapt_Algs::GradateMesh(double smoothness)
00631 {
00632     const double beta = smoothness;
00633 
00634     double maxShock, minShock;
00635     int iteration = 0, updates = 0;
00636 
00637     int* adjNodes, *boundNodes;
00638     int nadjNodes, nnodes;
00639     int meshNum = FEM_Mesh_default_read();
00640     //if (smoothness < 1.0) {
00641       //printf("");
00642     //}
00643 
00644     nnodes = theMesh->node.size();
00645     boundNodes = new int[nnodes];
00646     FEM_Mesh_data(meshNum, FEM_NODE, FEM_BOUNDARY, 
00647             boundNodes, 0, nnodes, FEM_INT, 1);
00648 
00649 
00650     //printf("Running h-shock mesh gradation with beta=%.3f\n", beta);
00651     fflush(NULL);
00652 
00653 #ifndef GRADATION_ITER_LIMIT
00654 #define GRADATION_ITER_LIMIT    10
00655 #endif
00656     
00657     do {
00658         maxShock = 0;
00659         minShock = 1e10;
00660         
00661         for (int node=0; node<nnodes; ++node) {
00662             if (boundNodes[node]!= 0 || !FEM_is_valid(meshNum, FEM_NODE, node))
00663                continue;
00664             //if (!FEM_is_valid(meshNum, FEM_NODE, node))
00665             //    continue;
00666             
00667             theMesh->n2n_getAll(node, adjNodes, nadjNodes);
00668             for (int adjNode=0; adjNode<nadjNodes; ++adjNode) {
00669                 double edgelen = length(node, adjNodes[adjNode]);
00670                 
00671                 // get adjacent elemnents and their sizes
00672                 int e1, e2;
00673                 theMesh->get2ElementsOnEdge(node, adjNodes[adjNode], &e1, &e2);
00674                 if (e1 <= -1 || e2 <= -1) continue; //this will not smooth across boundaries
00675                 
00676                 double s1, s2;
00677                 s1 = theMesh->elem[0].getMeshSizing(e1);
00678                 s2 = theMesh->elem[0].getMeshSizing(e2);
00679                 if (s1 <= 0 || s2 <= 0) continue;
00680                 
00681                 // h-shock=max(size ratio)^(1/edge length)
00682                 CkAssert(s1 >= 0 && s2 >= 0 && "Bad size");
00683                 //CkAssert(edgelen > 1e-6 && "Length 0 edge");
00684                 CkAssert(edgelen == edgelen && "Length inf edge");
00685 
00686                 double ratio = (s1 > s2) ? s1/s2 : s2/s1;
00687                 CkAssert (ratio >= 1.0 && ratio == ratio && "Bad ratio");
00688                 
00689                 // WARNING WARNING WARNING
00690                 // TEST ONLY, THIS IS NOT CORRECT
00691                 if (ratio > beta) {
00692                     if (s1 > s2) {
00693                         theMesh->elem[0].setMeshSizing(e1, s1 - (s1-s2)/3);
00694                         theMesh->elem[0].setMeshSizing(e2, s2 + (s1-s2)/3);
00695                     } else {
00696                         theMesh->elem[0].setMeshSizing(e2, s2 - (s2-s1)/3);
00697                         theMesh->elem[0].setMeshSizing(e1, s1 + (s2-s1)/3);
00698                     }
00699                     updates++;
00700                 }
00701                 if (ratio > maxShock) maxShock = ratio;
00702                 if (ratio < minShock) minShock = ratio;
00703                 
00704                 
00706                 //double hs = pow(ratio, 1.0/edgelen);
00707                 //
00708                 //if (hs > maxShock) maxShock = hs;
00709                 //if (hs < minShock) minShock = hs;
00710 
00715                 //
00716                 //if (hs > beta) {
00717                 //    double etasq = pow(beta, edgelen) / ratio;
00718                 //    etasq *= etasq;
00719 
00720                 //    //if (hs > 100) {
00721                 //    //    printf("hs: %8.5f\ns1: %8.5f\ns2: %8.5f\nedgelen: %8.5f\nratio: %8.5f\netasq: %8.5f", hs, s1, s2, edgelen, ratio, etasq);
00722                 //    //    abort();
00723                 //    //}
00724                 //    
00725                 //    if (s1 > s2) {
00726                 //        theMesh->elem[0].setMeshSizing(e1, s1 / etasq);
00727                 //    } else {
00728                 //        theMesh->elem[0].setMeshSizing(e2, s2 / etasq);
00729                 //    }
00730                 //    updates++;
00731                 //}
00732                 
00733             }
00734             free(adjNodes);
00735         } 
00736         
00737     //printf("Finished iteration %d\n", iteration);
00738         //printf("Max shock:%8.3f\n", maxShock);
00739         //printf("Min shock:%8.3f\n", minShock);
00740         //printf("Target:%8.3f\n", beta);
00741         
00742     } while (maxShock > beta && ++iteration < GRADATION_ITER_LIMIT);
00743 
00744     //printf("%d total updates in %d iterations in GradateMesh\n", updates, iteration);
00745     fflush(NULL);
00746 
00747     delete[] boundNodes;
00748 
00749     return;
00750 }
00751 
00753 void FEM_Adapt_Algs::SetMeshSize(int method, double factor, double *sizes)
00754 {
00755   numNodes = theMesh->node.size();
00756   numElements = theMesh->elem[0].size();
00757   int elemConn[3];
00758 
00759   if (method == 0) { // set uniform sizing specified in factor
00760 
00761     for (int i=0; i<numElements; i++) {
00762       theMesh->elem[0].setMeshSizing(i, factor);
00763     }
00764     //CkPrintf("ParFUM_SetMeshSize: UNIFORM %4.6e\n", factor);
00765   }
00766   else if (method == 1) { // copy sizing from array
00767     for (int i=0; i<numElements; i++) {
00768       if (sizes[i] > 0.0) {
00769     theMesh->elem[0].setMeshSizing(i, sizes[i]);
00770       }
00771     }
00772     //CkPrintf("ParFUM_SetMeshSize: SIZES input\n");
00773   }
00774   else if (method == 2) { // calculate current sizing and scale by factor
00775     double avgEdgeLength = 0.0;
00776     int width = theMesh->elem[0].getConn().width();
00777     int numEdges=3;
00778     if (dim==3) numEdges=6;
00779     for (int i=0; i<numElements; i++) {
00780       if(theMesh->elem[0].is_valid(i)) {
00781     theMesh->e2n_getAll(i, elemConn);
00782     for (int j=0; j<width-1; j++) {
00783       for (int k=j+1; k<width; k++) {
00784         avgEdgeLength += length(elemConn[j], elemConn[k]);
00785       }
00786     }
00787     avgEdgeLength += length(elemConn[0], elemConn[width-1]);
00788     avgEdgeLength /= (double)numEdges;
00789     theMesh->elem[0].setMeshSizing(i, factor*avgEdgeLength);
00790       }
00791     }
00792     //CkPrintf("ParFUM_SetMeshSize: CALCULATED & SCALED \n");
00793   }
00794   else if (method == 3) { // scale existing sizes by array sizes
00795     for (int i=0; i<numElements; i++) {
00796       if (sizes[i] > 0.0) {
00797     theMesh->elem[0].setMeshSizing(i, sizes[i]*theMesh->elem[0].getMeshSizing(i));
00798       }
00799     }
00800   }
00801   else if (method == 4) { // scale existing sizes by factor
00802     for (int i=0; i<numElements; i++) {
00803       theMesh->elem[0].setMeshSizing(i, factor*theMesh->elem[0].getMeshSizing(i));
00804     }
00805   }
00806   else if (method == 5) { // mesh sizing has been set independently; use as is
00807     //CkPrintf("ParFUM_SetMeshSize: USE EXISTING SIZES \n");
00808   }
00809   //  CkPrintf("Current mesh sizing: ");
00810   //for (int i=0; i<numElements; i++) {
00811   //CkPrintf("%4.6e ", theMesh->elem[0].getMeshSizing(i));
00812   //}
00813 }
00814 
00817 void FEM_Adapt_Algs::Insert(int eIdx, double len, int cFlag)
00818 {
00819   int i;
00820   if (cFlag) {
00821     i = ++coarsenHeapSize; 
00822     while ((coarsenElements[i/2].len>=len) && (i != 1)) {
00823       coarsenElements[i].len=coarsenElements[i/2].len;
00824       coarsenElements[i].elID=coarsenElements[i/2].elID;
00825       i/=2;
00826     }
00827     coarsenElements[i].elID=eIdx;
00828     coarsenElements[i].len=len; 
00829   }
00830   else {
00831     i = ++refineHeapSize; 
00832     while ((refineElements[i/2].len>=len) && (i != 1)) {
00833       refineElements[i].len=refineElements[i/2].len;
00834       refineElements[i].elID=refineElements[i/2].elID;
00835       i/=2;
00836     }
00837     refineElements[i].elID=eIdx;
00838     refineElements[i].len=len; 
00839   }
00840 }
00841 
00844 int FEM_Adapt_Algs::Delete_Min(int cflag)
00845 {
00846   int Child, i, Min_ID; 
00847   if (cflag) {
00848     Min_ID=coarsenElements[1].elID;
00849     for (i=1; i*2 <= coarsenHeapSize-1; i=Child) { // Find smaller child
00850       Child = i*2; // child is left child  
00851       if (Child != coarsenHeapSize)  // right child exists
00852     if (coarsenElements[Child+1].len < coarsenElements[Child].len)
00853       Child++; 
00854       // Percolate one level
00855       if (coarsenElements[coarsenHeapSize].len >= coarsenElements[Child].len) {
00856     coarsenElements[i].elID = coarsenElements[Child].elID;
00857     coarsenElements[i].len = coarsenElements[Child].len;
00858       }
00859       else break; 
00860     }
00861     coarsenElements[i].elID = coarsenElements[coarsenHeapSize].elID;
00862     coarsenElements[i].len = coarsenElements[coarsenHeapSize].len; 
00863     coarsenHeapSize--;
00864     return Min_ID; 
00865   }
00866   else {
00867     Min_ID=refineElements[1].elID;
00868     for (i=1; i*2 <= refineHeapSize-1; i=Child) { // Find smaller child
00869       Child = i*2;       // child is left child  
00870       if (Child !=refineHeapSize)  // right child exists
00871     if (refineElements[Child+1].len < refineElements[Child].len)
00872       Child++; 
00873       // Percolate one level
00874       if (refineElements[refineHeapSize].len >= refineElements[Child].len){  
00875     refineElements[i].elID = refineElements[Child].elID;   
00876     refineElements[i].len = refineElements[Child].len;
00877       }
00878       else break; 
00879     }
00880     refineElements[i].elID = refineElements[refineHeapSize].elID;
00881     refineElements[i].len = refineElements[refineHeapSize].len; 
00882     refineHeapSize--;
00883     return Min_ID; 
00884   }
00885 }
00886 
00887 
00888 
00889 // =====================  BEGIN refine_element_leb ========================= 
00899 int FEM_Adapt_Algs::refine_element_leb(int e) {
00900   int fixNode, otherNode, opNode, longEdge, nbr; 
00901   double eLens[3], longEdgeLen = 0.0;
00902   int elemConn[3];
00903 
00904   if(e==-1) {
00905     return -1;
00906   }
00907 
00908   theMesh->e2n_getAll(e, elemConn);
00909   eLens[0] = length(elemConn[0], elemConn[1]);
00910   eLens[1] = length(elemConn[1], elemConn[2]);
00911   eLens[2] = length(elemConn[2], elemConn[0]);
00912   for (int i=0; i<3; i++) {
00913     if (eLens[i] > longEdgeLen) {
00914       longEdgeLen = eLens[i];
00915       longEdge = i;
00916       fixNode = elemConn[i];
00917       otherNode = elemConn[(i+1)%3];
00918       opNode = elemConn[(i+2)%3];
00919     }
00920   }
00921 
00922   nbr = theMesh->e2e_getNbr(e, longEdge);
00923   if (nbr == -1) // e's longEdge is on physical boundary
00924     return theAdaptor->edge_bisect(fixNode, otherNode);
00925   int nbrOpNode = theAdaptor->e2n_getNot(nbr, fixNode, otherNode);
00926   double fixEdgeLen = length(fixNode, nbrOpNode);
00927   double otherEdgeLen = length(otherNode, nbrOpNode);
00928   if ((fixEdgeLen > longEdgeLen) || (otherEdgeLen > longEdgeLen)) { 
00929     // longEdge is not nbr's longest edge
00930     int newNode = theAdaptor->edge_bisect(fixNode, otherNode);
00931     if(newNode==-1) return -1;
00932     int propElem, propNode; // get the element to propagate on
00933     if (fixEdgeLen > otherEdgeLen) {
00934       propElem = theAdaptor->findElementWithNodes(newNode, fixNode, nbrOpNode);
00935       propNode = fixNode;
00936     }
00937     else {
00938       propElem = theAdaptor->findElementWithNodes(newNode,otherNode,nbrOpNode);
00939       propNode = otherNode;
00940     }
00941 
00942     //if propElem is ghost, then it's propagating to neighbor, otherwise not
00943     if(!FEM_Is_ghost_index(propElem)) {
00944       refine_flip_element_leb(propElem,propNode,newNode,nbrOpNode,longEdgeLen);
00945     }
00946     else {
00947       int localChk, nbrChk;
00948       localChk = theMod->getfmUtil()->getIdx();
00949       nbrChk = theMod->getfmUtil()->getRemoteIdx(theMesh,propElem,0);
00950       int propNodeT = theAdaptor->getSharedNodeIdxl(propNode, nbrChk);
00951       int newNodeT = theAdaptor->getSharedNodeIdxl(newNode, nbrChk);
00952       int nbrghost = (nbrOpNode>=0)?0:1;
00953       int nbrOpNodeT = (nbrOpNode>=0)?(theAdaptor->getSharedNodeIdxl(nbrOpNode, nbrChk)):(theAdaptor->getGhostNodeIdxl(nbrOpNode, nbrChk));
00954       int propElemT = theAdaptor->getGhostElementIdxl(propElem, nbrChk);
00955       meshMod[nbrChk].refine_flip_element_leb(localChk, propElemT, propNodeT, newNodeT,nbrOpNodeT,nbrghost,longEdgeLen);
00956     }
00957     return newNode;
00958   }
00959   else return theAdaptor->edge_bisect(fixNode, otherNode); // longEdge on nbr
00960 }
00961 
00966 void FEM_Adapt_Algs::refine_flip_element_leb(int e, int p, int n1, int n2, 
00967                          double le) 
00968 {
00969   int newNode = refine_element_leb(e);
00970   if(newNode == -1) return;
00971   (void) theAdaptor->edge_flip(n1, n2);
00972   if (length(p, newNode) > le) {
00973     int localChk = theMod->getfmUtil()->getIdx();
00974     int newElem = theAdaptor->findElementWithNodes(newNode, n1, p);
00975     refine_flip_element_leb(newElem, p, n1, newNode, le);
00976   }
00977 }
00978 // ========================  END refine_element_leb ========================
00979 
00980 
00981 
00984 int FEM_Adapt_Algs::simple_refine(double targetA, double xmin, double ymin, double xmax, double ymax) {
00985   bool adapted = true;
00986   refineElements = refineStack = NULL;
00987   refineTop = refineHeapSize = 0;
00988   int elemConn[3];
00989   double coordsn1[2], coordsn2[2], coordsn3[2];
00990   int elemWidth=3;
00991 
00992   while(adapted) {
00993     adapted = false;
00994     int noEle = theMesh->elem[0].size();
00995     refineElements = new elemHeap[noEle+1];
00996     for(int i=0; i<noEle; i++) {
00997       if(theMesh->elem[0].is_valid(i)) {
00998     theMesh->e2n_getAll(i,elemConn,0);
00999     bool notclear = false;
01000     for(int j=0; j<elemWidth; j++) {
01001       int nd = elemConn[j];
01002       if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
01003       //else cannot look up because connectivity might be screwed up,
01004       //I am hoping that at least one of the nodes will be local
01005     }
01006     if(notclear) continue;
01007     getCoord(elemConn[0], coordsn1);
01008     getCoord(elemConn[1], coordsn2);
01009     getCoord(elemConn[2], coordsn3);
01010     if(coordsn1[0]<-1.0 ||coordsn1[1]<-1.0 || coordsn2[0]<-1.0 ||coordsn2[1]<-1.0 || coordsn3[0]<-1.0 ||coordsn3[1]<-1.0) continue;
01011     //do a refinement only if it has any node within the refinement box
01012     if((coordsn1[0]<xmax && coordsn1[0]>xmin && coordsn1[1]<ymax && coordsn1[1]>ymin) || (coordsn2[0]<xmax && coordsn2[0]>xmin && coordsn2[1]<ymax && coordsn2[1]>ymin) || (coordsn3[0]<xmax && coordsn3[0]>xmin && coordsn3[1]<ymax && coordsn3[1]>ymin)) {
01013       Insert(i,-getArea(coordsn1, coordsn2, coordsn3),0); 
01014       //negative area to trick the minheap to become a maxheap!
01015     }
01016     else {
01017       Insert(i,-MINAREA,0); //make it believe that this triangle does not need refinement
01018     }
01019       } else {
01020     Insert(i,-MINAREA,0);
01021       }
01022     }
01023 
01024     //bad bubblesort
01025     /*for(int i=0; i<noEle; i++) {
01026       for(int j=i+1; j<noEle; j++) {
01027     if(areas[j] > areas[i]) {
01028       double tmp = areas[j];
01029       areas[j] = areas[i];
01030       areas[i] = tmp;
01031       int t = map1[j];
01032       map1[j] = map1[i];
01033       map1[i] = t;
01034     }
01035       }
01036       }*/
01037     
01038     /*for(int i=0; i<noEle; i++) {
01039       if(theMesh->elem[0].is_valid(map1[i])) {
01040     if(areas[i] > targetA) {
01041       int ret = refine_element_leb(map1[i]);
01042       if(ret!=-1) {
01043         adapted = true;
01044       }
01045     }
01046     }
01047       //if(adapted) break;
01048       }*/
01049     while(refineHeapSize>0) {
01050       int tmp = Delete_Min(0);
01051       if(tmp!=-1 && theMesh->elem[0].is_valid(tmp)) {
01052     //since we are reusing element indices, so do not trust earlier areas
01053     theMesh->e2n_getAll(tmp,elemConn,0);
01054     bool notclear = false;
01055     for(int j=0; j<3; j++) {
01056       int nd = elemConn[j];
01057       if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
01058       //else cannot look up because connectivity might be screwed up,
01059       //I am hoping that at least one of the nodes will be local
01060     }
01061     if(notclear) continue;
01062     getCoord(elemConn[0], coordsn1);
01063     getCoord(elemConn[1], coordsn2);
01064     getCoord(elemConn[2], coordsn3);
01065     if(coordsn1[0]<-1.0 ||coordsn1[1]<-1.0 || coordsn2[0]<-1.0 ||coordsn2[1]<-1.0 || coordsn3[0]<-1.0 ||coordsn3[1]<-1.0) continue;
01066     if(getArea(coordsn1, coordsn2, coordsn3) > targetA) {
01067       double len1 = length(coordsn1,coordsn2);
01068       double len2 = length(coordsn2,coordsn3);
01069       double len3 = length(coordsn3,coordsn1);
01070       int n1=elemConn[0], n2=elemConn[1];
01071       double max=len1;
01072       if(len2>max) {
01073         max = len2; n1 = elemConn[1]; n2 = elemConn[2];
01074       }
01075       if(len3>max) {
01076         max = len3; n1 = elemConn[2]; n2 = elemConn[0];
01077       }
01078       //int ret = refine_element_leb(tmp);
01079       int ret = theAdaptor->edge_bisect(n1,n2);
01080       if(ret!=-1) adapted = true;
01081     }
01082       }
01083       //if(adapted) break;
01084     }
01085     delete[] refineElements;
01086     //if(adapted) break;
01087   }
01088 
01089   if(adapted) return -1;
01090   else return 1;
01091 }
01092 
01095 int FEM_Adapt_Algs::simple_coarsen(double targetA, double xmin, double ymin, double xmax, double ymax) {
01096   int noEle = theMesh->elem[0].size();
01097   int *shortestEdge = (int *)malloc(2*sizeof(int));
01098   bool adapted = true;
01099   coarsenElements = NULL;
01100   coarsenHeapSize = 0;
01101   int elemConn[3];
01102   double coordsn1[2], coordsn2[2], coordsn3[2];
01103   int elemWidth=3;
01104 
01105   while(adapted) {
01106     adapted = false;
01107     coarsenElements = new elemHeap[noEle+1];
01108     coarsenElements[0].len=-2.0;
01109     coarsenElements[0].elID=-1;
01110     for(int i=0; i<noEle; i++) {
01111       if(theMesh->elem[0].is_valid(i)) {
01112     theMesh->e2n_getAll(i,elemConn,0);
01113     bool notclear = false;
01114     for(int j=0; j<elemWidth; j++) {
01115       int nd = elemConn[j];
01116       if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
01117       //else cannot look up because connectivity might be screwed up,
01118       //I am hoping that at least one of the nodes will be local
01119     }
01120     if(notclear) continue;
01121     getCoord(elemConn[0], coordsn1);
01122     getCoord(elemConn[1], coordsn2);
01123     getCoord(elemConn[2], coordsn3);
01124     if(coordsn1[0]<-1.0 ||coordsn1[1]<-1.0 || coordsn2[0]<-1.0 ||coordsn2[1]<-1.0 || coordsn3[0]<-1.0 ||coordsn3[1]<-1.0) continue;
01125     //do a coarsening only if it has any node within the coarsen box
01126     if((coordsn1[0]<xmax && coordsn1[0]>xmin && coordsn1[1]<ymax && coordsn1[1]>ymin) || (coordsn2[0]<xmax && coordsn2[0]>xmin && coordsn2[1]<ymax && coordsn2[1]>ymin) || (coordsn3[0]<xmax && coordsn3[0]>xmin && coordsn3[1]<ymax && coordsn3[1]>ymin)) {
01127       Insert(i,getArea(coordsn1, coordsn2, coordsn3),1);
01128     } 
01129     else {
01130       Insert(i,MAXAREA,1); //make it believe that this triangle is big enough
01131     }
01132       } else {
01133     Insert(i,MAXAREA,1);
01134       }
01135     }
01136     
01137     /*for(int i=0; i<noEle; i++) {
01138       for(int j=i+1; j<noEle; j++) {
01139     if(areas[j] < areas[i]) {
01140       double tmp = areas[j];
01141       areas[j] = areas[i];
01142       areas[i] = tmp;
01143       int t = map1[j];
01144       map1[j] = map1[i];
01145       map1[i] = t;
01146     }
01147       }
01148       }*/
01149     
01150     /*for(int i=0; i<noEle; i++) {
01151       if(theMesh->elem[0].is_valid(map1[i])) {
01152     if(areas[i] < targetA) {
01153       //find the nodes along the smallest edge & coarsen the edge
01154       theMesh->e2n_getAll(map1[i],elemConn,0);
01155       getShortestEdge(elemConn[0], elemConn[1], elemConn[2], shortestEdge);
01156       int ret = theAdaptor->edge_contraction(shortestEdge[0], shortestEdge[1]);
01157       if(ret != -1) adapted = true;
01158     }
01159       }
01160       //if(adapted) break;
01161       }*/
01162     while(coarsenHeapSize>0) {
01163       int tmp = Delete_Min(1);
01164       if(tmp!=-1 && theMesh->elem[0].is_valid(tmp)) {
01165     theMesh->e2n_getAll(tmp,elemConn,0);
01166     bool notclear = false;
01167     for(int j=0; j<3; j++) {
01168       int nd = elemConn[j];
01169       if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
01170       //else cannot look up because connectivity might be screwed up,
01171       //I am hoping that at least one of the nodes will be local
01172     }
01173     if(notclear) continue;
01174     getCoord(elemConn[0], coordsn1);
01175     getCoord(elemConn[1], coordsn2);
01176     getCoord(elemConn[2], coordsn3);
01177     if(coordsn1[0]<-1.0 ||coordsn1[1]<-1.0 || coordsn2[0]<-1.0 ||coordsn2[1]<-1.0 || coordsn3[0]<-1.0 ||coordsn3[1]<-1.0) continue;
01178     if(getArea(coordsn1, coordsn2, coordsn3) < targetA) {
01179       getShortestEdge(elemConn[0], elemConn[1], elemConn[2], shortestEdge);
01180       int ret = theAdaptor->edge_contraction(shortestEdge[0], shortestEdge[1]);
01181       if(ret != -1) adapted = true;
01182     }
01183       }
01184       //if(adapted) break;
01185     }
01186     delete[] coarsenElements;
01187     //if(adapted) break;
01188   }
01189   free(shortestEdge);
01190 
01191   if(adapted) return -1;
01192   else return 1;
01193 }
01194 
01195 
01196 
01197 //HELPER functions
01198 double FEM_Adapt_Algs::length(int n1, int n2)
01199 {
01200   double coordsn1[2], coordsn2[2];
01201   getCoord(n1, coordsn1);
01202   getCoord(n2, coordsn2);
01203 
01204   double ret = length(coordsn1, coordsn2);
01205   return ret;
01206 }
01207 
01208 double FEM_Adapt_Algs::length(double *n1_coord, double *n2_coord) { 
01209   double d, ds_sum=0.0;
01210 
01211   for (int i=0; i<dim; i++) {
01212     if(n1_coord[i]<-1.0 || n2_coord[i]<-1.0) return -2.0;
01213     d = n1_coord[i] - n2_coord[i];
01214     ds_sum += d*d;
01215   }
01216   return (sqrt(ds_sum));
01217 }
01218 
01219 double FEM_Adapt_Algs::getArea(int n1, int n2, int n3)
01220 {
01221   double coordsn1[2], coordsn2[2], coordsn3[2];
01222   getCoord(n1, coordsn1);
01223   getCoord(n2, coordsn2);
01224   getCoord(n3, coordsn3);
01225 
01226   double ret = getArea(coordsn1, coordsn2, coordsn3);
01227   return ret;
01228 }
01229 
01230 double FEM_Adapt_Algs::getArea(double *n1_coord, double *n2_coord, double *n3_coord) {
01231   double area=0.0;
01232   double aLen, bLen, cLen, sLen, d, ds_sum;
01233 
01234   ds_sum = 0.0;
01235   for (int i=0; i<2; i++) {
01236     d = n1_coord[i] - n2_coord[i];
01237     ds_sum += d*d;
01238   }
01239   aLen = sqrt(ds_sum);
01240   ds_sum = 0.0;
01241   for (int i=0; i<2; i++) {
01242     d = n2_coord[i] - n3_coord[i];
01243     ds_sum += d*d;
01244   }
01245   bLen = sqrt(ds_sum);
01246   ds_sum = 0.0;
01247   for (int i=0; i<2; i++) {
01248     d = n3_coord[i] - n1_coord[i];
01249     ds_sum += d*d;
01250   }
01251   cLen = sqrt(ds_sum);
01252   sLen = (aLen+bLen+cLen)/2;
01253   if(sLen-aLen < 0) return 0.0;
01254   else if(sLen-bLen < 0) return 0.0;
01255   else if(sLen-cLen < 0) return 0.0; //area too small to note
01256   return (sqrt(sLen*(sLen-aLen)*(sLen-bLen)*(sLen-cLen)));
01257 }
01258 
01259 double FEM_Adapt_Algs::getSignedArea(int n1, int n2, int n3) {
01260   double coordsn1[2], coordsn2[2], coordsn3[2];
01261   getCoord(n1, coordsn1);
01262   getCoord(n2, coordsn2);
01263   getCoord(n3, coordsn3);
01264 
01265   double ret = getSignedArea(coordsn1, coordsn2, coordsn3);
01266   return ret;
01267 }
01268 
01269 double FEM_Adapt_Algs::getSignedArea(double *n1_coord, double *n2_coord, double *n3_coord) {
01270   double area=0.0;
01271   double vec1_x, vec1_y, vec2_x, vec2_y;
01272 
01273   vec1_x = n1_coord[0] - n2_coord[0];
01274   vec1_y = n1_coord[1] - n2_coord[1];
01275   vec2_x = n3_coord[0] - n2_coord[0];
01276   vec2_y = n3_coord[1] - n2_coord[1];
01277 
01278   area = vec1_x*vec2_y - vec2_x*vec1_y;
01279   return area;
01280 }
01281 
01282 int FEM_Adapt_Algs::getCoord(int n1, double *crds) {
01283   if(!FEM_Is_ghost_index(n1)) {
01284     FEM_Mesh_dataP(theMesh, FEM_NODE, coord_attr, (void *)crds, n1, 1, FEM_DOUBLE, dim);
01285   }
01286   else {
01287     int numchunks;
01288     IDXL_Share **chunks1;
01289     theMod->fmUtil->getChunkNos(0,n1,&numchunks,&chunks1);
01290     int index = theMod->idx;
01291     for(int j=0; j<numchunks; j++) {
01292       int chk = chunks1[j]->chk;
01293       if(chk==index) continue;
01294       int ghostidx = theMod->fmUtil->exists_in_IDXL(theMesh,n1,chk,2);
01295       double2Msg *d = meshMod[chk].getRemoteCoord(index,ghostidx);
01296       crds[0] = d->i;
01297       crds[1] = d->j;
01298       for(int j=0; j<numchunks; j++) {
01299     delete chunks1[j];
01300       }
01301       if(numchunks != 0) free(chunks1);
01302       delete d;
01303       break;
01304     }
01305   }
01306   return 1;
01307 }
01308 
01309 int FEM_Adapt_Algs::getShortestEdge(int n1, int n2, int n3, int* shortestEdge) {
01310   //note that getCoord might be a remote call, which means 
01311   //it might not have the same value in the mem
01312   //location if the memory is reused by someone else meanwhile!
01313   double coordsn1[2], coordsn2[2], coordsn3[2];
01314   getCoord(n1, coordsn1);
01315   getCoord(n2, coordsn2);
01316   getCoord(n3, coordsn3);
01317 
01318   double aLen = length(coordsn1, coordsn2);
01319   int shortest = 0;
01320 
01321   double bLen = length(coordsn2, coordsn3);
01322   if(bLen < aLen) shortest = 1;
01323 
01324   double cLen = length(coordsn3, coordsn1);
01325   if((cLen < aLen) && (cLen < bLen)) shortest = 2;
01326 
01327   if(shortest==0) {
01328     shortestEdge[0] = n1;
01329     shortestEdge[1] = n2;
01330   }
01331   else if (shortest==1) {
01332     shortestEdge[0] = n2;
01333     shortestEdge[1] = n3;
01334   }
01335   else {
01336     shortestEdge[0] = n3;
01337     shortestEdge[1] = n1;
01338   }
01339   return 1;
01340 }
01341 
01342 
01343 
01347 double FEM_Adapt_Algs::getAreaQuality(int elem)
01348 {
01349   double f, q, len[3];
01350   int n[3];
01351   double currentArea;
01352   double coordsn1[2], coordsn2[2], coordsn3[2];
01353 
01354   theMesh->e2n_getAll(elem, n);
01355   getCoord(n[0], coordsn1);
01356   getCoord(n[1], coordsn2);
01357   getCoord(n[2], coordsn3);
01358 
01359   currentArea = getArea(coordsn1, coordsn2, coordsn3);
01360 
01361   len[0] = length(coordsn1, coordsn2);
01362   len[1] = length(coordsn2, coordsn3);
01363   len[2] = length(coordsn3, coordsn1);
01364   f = 4.0*sqrt(3.0); //proportionality constant
01365   q = (f*currentArea)/(len[0]*len[0]+len[1]*len[1]+len[2]*len[2]);
01366   return q;
01367 }
01368 
01373 void FEM_Adapt_Algs::ensureQuality(int n1, int n2, int n3) {
01374   double coordsn1[2], coordsn2[2], coordsn3[2];
01375   getCoord(n1, coordsn1);
01376   getCoord(n2, coordsn2);
01377   getCoord(n3, coordsn3);
01378   double area = getSignedArea(coordsn1, coordsn2, coordsn3);
01379   double len1 = length(coordsn1, coordsn2);
01380   double len2 = length(coordsn2, coordsn3);
01381   double len3 = length(coordsn3, coordsn1);
01382   double max = len1;
01383   if(len2>max) max = len2;
01384   if(len3>max) max = len3;
01385   //shortest edge
01386   double min = len1;
01387   if(len2<min) min = len2;
01388   if(len3<min) min = len3;
01389   double shortest_al = fabs(area/max);
01390   double largestR = max/shortest_al;
01391   CkAssert(largestR<=100.0 && -area > SLIVERAREA);
01392 }
01393 
01396 bool FEM_Adapt_Algs::controlQualityF(int n1, int n2, int n3, int n4) {
01397   //n1 or n2 will be replaced by n4
01398   double coordsn1[2], coordsn2[2], coordsn3[2], coordsn4[2];
01399   if(n4==-1) return false;
01400   getCoord(n1, coordsn1);
01401   getCoord(n2, coordsn2);
01402   getCoord(n3, coordsn3);
01403   getCoord(n4, coordsn4);
01404   bool flag = false;
01405   if(!flag) flag = controlQualityC(coordsn3,coordsn1,coordsn2,coordsn4);
01406   if(!flag) flag = controlQualityC(coordsn4,coordsn2,coordsn1,coordsn3);
01407   return flag;
01408 }
01409 
01412 bool FEM_Adapt_Algs::controlQualityR(int n1, int n2, int n3, int n4) {
01413   //n1 or n2 will be replaced by n4
01414   double coordsn1[2], coordsn2[2], coordsn3[2], coordsn4[2], coordsn5[2];
01415   if(n4==-1) return false;
01416   getCoord(n1, coordsn1);
01417   getCoord(n2, coordsn2);
01418   getCoord(n3, coordsn3);
01419   getCoord(n4, coordsn4);
01420   coordsn5[0] = (coordsn1[0]+coordsn2[0])*0.5;
01421   coordsn5[1] = (coordsn1[1]+coordsn2[1])*0.5;
01422   bool flag = false;
01423   if(!flag) flag = theMod->fmAdaptAlgs->controlQualityR(coordsn1,coordsn3,coordsn5);
01424   if(!flag) flag = theMod->fmAdaptAlgs->controlQualityR(coordsn2,coordsn3,coordsn5);
01425   if(!flag) flag = theMod->fmAdaptAlgs->controlQualityR(coordsn1,coordsn4,coordsn5);
01426   if(!flag) flag = theMod->fmAdaptAlgs->controlQualityR(coordsn2,coordsn4,coordsn5);
01427   return flag;
01428 }
01429 
01430 bool FEM_Adapt_Algs::controlQualityR(double *coordsn1, double *coordsn2, double *coordsn3) {
01431   double area = getArea(coordsn1, coordsn2, coordsn3);
01432   //do some quality preservation
01433   double len1 = length(coordsn1, coordsn2);
01434   double len2 = length(coordsn2, coordsn3);
01435   double len3 = length(coordsn3, coordsn1);
01436   //longest edge
01437   double max = len1;
01438   if(len2>max) max = len2;
01439   if(len3>max) max = len3;
01440   double shortest_al = area/max;
01441   double largestR = max/shortest_al;
01442   if(largestR>50.0) return true;
01443   else return false;
01444 }
01445 
01450 bool FEM_Adapt_Algs::controlQualityC(int n1, int n2, int n3, double *n4_coord) {
01451   //n3 is the node to be deleted, n4 is the new node to be added
01452   double coordsn1[2], coordsn2[2], coordsn3[2];
01453   getCoord(n1, coordsn1);
01454   getCoord(n2, coordsn2);
01455   getCoord(n3, coordsn3);
01456   return controlQualityC(coordsn1,coordsn2,coordsn3,n4_coord);
01457 }
01458 
01459 bool FEM_Adapt_Algs::controlQualityC(double *coordsn1, double *coordsn2, double *coordsn3, double *n4_coord) {
01460   double ret_old = getSignedArea(coordsn1, coordsn2, coordsn3);
01461   double ret_new = getSignedArea(coordsn1, coordsn2, n4_coord);
01462   //do some quality preservation
01463   double len1 = length(coordsn1, coordsn2);
01464   double len2 = length(coordsn2, n4_coord);
01465   double len3 = length(n4_coord, coordsn1);
01466   //longest edge
01467   double max = len1;
01468   if(len2>max) max = len2;
01469   if(len3>max) max = len3;
01470   //shortest edge
01471   double min = len1;
01472   if(len2<min) min = len2;
01473   if(len3<min) min = len3;
01474   double shortest_al = ret_new/max;
01475   double largestR = max/shortest_al;
01476   if(ret_old > SLIVERAREA && ret_new < -SLIVERAREA) return true; //it is a flip
01477   else if(ret_old < -SLIVERAREA && ret_new > SLIVERAREA) return true; //it is a flip
01478   else if(fabs(ret_new) < SLIVERAREA) return true; // it is a sliver
01479   //else if(fabs(shortest_al) < 1e-5) return true; //shortest altitude is too small
01480   //else if(fabs(min) < 1e-5) return true; //shortest edge is too small
01481   else if(fabs(largestR)>50.0) return true;
01482   else return false;
01483 }
01484 
01485 
01486 
01492 bool FEM_Adapt_Algs::flipOrBisect(int elId, int n1, int n2, int maxEdgeIdx, double maxlen) {
01493   //return true if it should flip
01494   int nbrEl = theMesh->e2e_getNbr(elId,maxEdgeIdx);
01495   double len4=0.0, len5=0.0;
01496   if(nbrEl!=-1) {
01497     int con1[3];
01498     theMesh->e2n_getAll(nbrEl,con1);
01499     int nbrnode=-1;
01500     for(int j=0; j<3; j++) {
01501       if(con1[j]!=n1 && con1[j]!=n2) {
01502     nbrnode = con1[j];
01503     break;
01504       }
01505     }
01506     len4 = length(n1,nbrnode);
01507     len5 = length(n2,nbrnode);
01508   }
01509   if(len4>1.2*maxlen || len5>1.2*maxlen) {
01510     return true;
01511   }
01512   else return false;
01513 }
01514 
01515 void FEM_Adapt_Algs::tests(bool b=true) {
01516   //test the mesh for slivered triangles
01517   if(!b) {
01518     theMod->fmUtil->AreaTest(theMesh);
01519     return;
01520   }
01521   theMod->fmUtil->AreaTest(theMesh);
01522   theMod->fmUtil->StructureTest(theMesh);
01523   theMod->fmUtil->IdxlListTest(theMesh);
01524   theMod->fmUtil->residualLockTest(theMesh);
01525   /*for(int i=0; i<theMesh->node.size(); i++) {
01526     if(theMesh->node.is_valid(i)) CkPrintf("Valid -- ");
01527     else  CkPrintf("Invalid -- ");
01528     theMod->fmUtil->FEM_Print_coords(theMesh,i);
01529     }*/
01530   return;
01531 }
libs/ck-libs/ParFUM/adapt_algs.C
