Charm++: Charm Source Code Documentation

00001 /* File: bulk_adapt.C
00002  * Authors: Terry Wilmarth
00003  */
00004 
00008 #include "ParFUM.h"
00009 #include "ParFUM_internals.h"
00010 #include "ParFUM_SA.h"
00011 #include "bulk_adapt.h"
00012 #include "bulk_adapt_ops.h"
00013 
00014 #define MINAREA 1.0e-18
00015 #define MAXAREA 1.0e12
00016 #define GRADATION 1.2
00017 #define ADAPT_VERBOSE
00018 
00019 CtvDeclare(Bulk_Adapt *, _bulkAdapt);
00020 
00027 void Bulk_Adapt::ParFUM_Refine(int qm, int method, double factor, double *sizes)
00028 {
00029   SetMeshSize(method, factor, sizes);
00030   ParFUM_GradateMesh(GRADATION);
00031   (void)Refine_h(qm, method, factor, sizes);
00032 }
00033 
00035 int Bulk_Adapt::Refine_h(int qm, int method, double factor, double *sizes)
00036 {
00037   int numNodes, numElements;
00038   int elId, mods, iter_mods;
00039   int elemWidth = theMesh->elem[0].getConn().width();
00040   Element_Bucket elemsToRefine;
00041   int elemConn[3];
00042   int count = 0;
00043   int myId = FEM_My_partition();
00044 
00045   mods=0;
00046   iter_mods=1;
00047 
00048   while ((iter_mods != 0) || (count < 20)) {
00049     if (iter_mods == 0)
00050       count++;
00051     else 
00052       count = 0;
00053     iter_mods=0;
00054     numNodes = theMesh->node.size();
00055     numElements = theMesh->elem[0].size();
00056     // sort elements to be refined by quality into elemsToRefine
00057     elemsToRefine.Clear();
00058     elemsToRefine.Alloc(numElements);
00059     for (int i=0; i<numElements; i++) { 
00060       if (theMesh->elem[0].is_valid(i)) {
00061     // find avg and max edge length of element i
00062     double avgEdgeLength, maxEdgeLength, minEdgeLength;
00063     int maxEdge, minEdge;
00064     findEdgeLengths(i, &avgEdgeLength, &maxEdgeLength, &maxEdge, 
00065             &minEdgeLength, &minEdge);
00066     //double qFactor=getAreaQuality(i);
00067     if (theMesh->elem[0].getMeshSizing(i) <= 0.0) 
00068       CkPrintf("WARNING: mesh element %d has no sizing!\n", i);
00069     if ((theMesh->elem[0].getMeshSizing(i) > 0.0) &&
00070       (avgEdgeLength > (theMesh->elem[0].getMeshSizing(i)*REFINE_TOL))){
00071         //|| (qFactor < QUALITY_MIN)) {
00072       //elemsToRefine.Insert(i, qFactor*(1.0/(avgEdgeLength+maxEdgeLength)), 0);
00073       elemsToRefine.Insert(i, 1.0/maxEdgeLength);
00074 #ifdef ADAPT_VERBOSE
00075       CkPrintf("[%d]ParFUM_Refine: Added element %d to refine list: targetSize=%6.4f maxEdgeLength=%6.4f maxEdge=%d\n", myId, i, theMesh->elem[0].getMeshSizing(i), maxEdgeLength, maxEdge);
00076 #endif    
00077     }
00078       }
00079     }
00080 
00081     while (!elemsToRefine.IsBucketEmpty()) { // loop through the elements
00082       int n1, n2;
00083       double len;
00084 
00085       elId = elemsToRefine.Remove(&len);
00086 
00087       if ((elId != -1) && (theMesh->elem[0].is_valid(elId))) {
00088     int n1, n2;
00089     double avgEdgeLength, maxEdgeLength, minEdgeLength;
00090     int maxEdge, minEdge;
00091     findEdgeLengths(elId, &avgEdgeLength, &maxEdgeLength, &maxEdge, 
00092             &minEdgeLength, &minEdge);
00093     //double qFactor=getAreaQuality(elId);
00094 
00095     //if (len == qFactor*(1.0/(avgEdgeLength+maxEdgeLength))) { // this elem does not appear to be modified; refine it
00096     if (len == (1.0/maxEdgeLength)) { // this elem does not appear to be modified; refine it
00097       if ((theMesh->elem[0].getMeshSizing(elId) > 0.0) &&
00098           (avgEdgeLength>(theMesh->elem[0].getMeshSizing(elId)*REFINE_TOL))){
00099         //|| (qFactor < QUALITY_MIN)) {
00100         //decide if we should do a flip or bisect
00101 #ifdef ADAPT_VERBOSE
00102         CkPrintf("[%d]ParFUM_Refine: Refining element %d: targetSize=%6.4f maxEdgeLength=%6.4f maxEdge=%d\n", myId, elId, theMesh->elem[0].getMeshSizing(elId), maxEdgeLength, maxEdge);
00103 #endif    
00104         RegionID lockRegionID;
00105         int success;
00106         lockRegionID.localID = -1;
00107         success = theMesh->parfumSA->bulkAdapt->lock_3D_region(elId, 0, maxEdge, maxEdgeLength, 
00108                                    &lockRegionID);
00109         if (success == 2) { // lock obtained straight away
00110 #ifdef ADAPT_VERBOSE
00111           CkPrintf("[%d]ParFUM_Refine: Refining element %d: GOT THE LOCK\n", myId, elId);
00112 #endif    
00113           iter_mods = iter_mods + 1;
00114           CkAssert(theMesh->parfumSA->holdingLock == lockRegionID);
00115           (void) theMesh->parfumSA->bulkAdapt->edge_bisect(elId, 0, maxEdge, dim, lockRegionID);
00116           theMesh->parfumSA->bulkAdapt->unlock_3D_region(lockRegionID);
00117         }
00118         else if (success == 1) { // lock is pending
00119 #ifdef ADAPT_VERBOSE
00120           CkPrintf("[%d]ParFUM_Refine: Refining element %d: LOCK PENDING\n", myId, elId);
00121 #endif    
00122           while (success == 1) {
00123         CthYield();
00124         double DavgEdgeLength, DmaxEdgeLength, DminEdgeLength;
00125         int DmaxEdge, DminEdge;
00126         findEdgeLengths(elId, &DavgEdgeLength, &DmaxEdgeLength, &DmaxEdge, 
00127                 &DminEdgeLength, &DminEdge);
00128         if (len == (1.0/DmaxEdgeLength)) { // element is unchanged: i.e. still needs refinement
00129           success = theMesh->parfumSA->bulkAdapt->lock_3D_region(elId, 0, DmaxEdge, DmaxEdgeLength, 
00130                                      &lockRegionID);
00131         }
00132         else {
00133           theMesh->parfumSA->bulkAdapt->unpend_3D_region(lockRegionID);
00134 #ifdef ADAPT_VERBOSE
00135           CkPrintf("[%d]ParFUM_Refine: Refining element %d: FAILED TO GET LOCK\n", myId, elId);
00136 #endif    
00137           elemsToRefine.Insert(elId, 1.0/DmaxEdgeLength);     
00138           break; // element was modified by a different operation
00139         }
00140           }
00141           if (success==2) { // broke out of loop with a successful lock
00142 #ifdef ADAPT_VERBOSE
00143         CkPrintf("[%d]ParFUM_Refine: Refining element %d: GOT THE LOCK\n", myId, elId);
00144 #endif    
00145         iter_mods = iter_mods + 1;
00146         CkAssert(theMesh->parfumSA->holdingLock == lockRegionID);
00147         (void) theMesh->parfumSA->bulkAdapt->edge_bisect(elId, 0, maxEdge, dim, lockRegionID);
00148         theMesh->parfumSA->bulkAdapt->unlock_3D_region(lockRegionID);
00149           }
00150           else if (success==0) {
00151         theMesh->parfumSA->bulkAdapt->unpend_3D_region(lockRegionID);
00152 #ifdef ADAPT_VERBOSE
00153         CkPrintf("[%d]ParFUM_Refine: Refining element %d: FAILED TO GET LOCK\n", myId, elId);
00154 #endif    
00155         elemsToRefine.Insert(elId, 1.0/maxEdgeLength);    
00156           }
00157         }
00158         else if (success == 0) { // lock failed immediately
00159 #ifdef ADAPT_VERBOSE
00160           CkPrintf("[%d]ParFUM_Refine: Refining element %d: FAILED TO GET LOCK\n", myId, elId);
00161 #endif    
00162           elemsToRefine.Insert(elId, 1.0/maxEdgeLength);      
00163         }
00164       }
00165     }
00166     else { // elem was modified; return to bucket
00167       elemsToRefine.Insert(elId, 1.0/maxEdgeLength);      
00168     }
00169       }
00170       CthYield(); // give other chunks on the same PE a chance
00171     }
00172     mods += iter_mods;
00173 #ifdef ADAPT_VERBOSE
00174     CkPrintf("[%d]ParFUM_Refine: %d modifications in last pass.\n",myId,iter_mods);
00175 #endif
00176   }
00177 #ifdef ADAPT_VERBOSE
00178   CkPrintf("[%d]ParFUM_Refine: %d total modifications.\n",myId,mods);
00179 #endif
00180   elemsToRefine.Clear();
00181   return mods;
00182 }
00183 
00190 void Bulk_Adapt::ParFUM_Coarsen(int qm, int method, double factor, 
00191                  double *sizes)
00192 {
00193   SetMeshSize(method, factor, sizes);
00194   ParFUM_GradateMesh(GRADATION);
00195   //  (void)Coarsen_h(qm, method, factor, sizes);
00196 }
00197 
00200 int Bulk_Adapt::Coarsen_h(int qm, int method, double factor, double *sizes)
00201 {
00202   return 0;
00203   /*
00204   // loop through elemsToRefine
00205   int elId, mods=0, iter_mods=1, pass=0;
00206   int elemWidth = theMesh->elem[0].getConn().width();
00207   double qFactor;
00208   coarsenElements = NULL;
00209   coarsenHeapSize = 0;
00210   int elemConn[3];
00211   while (iter_mods != 0) {
00212     iter_mods=0;
00213     pass++;
00214     numNodes = theMesh->node.size();
00215     numElements = theMesh->elem[0].size();
00216     // sort elements to be refined by quality into elemsToRefine
00217     if (coarsenElements) delete [] coarsenElements;
00218     coarsenElements = new elemHeap[numElements+1];
00219     coarsenElements[0].len=-2.0;
00220     coarsenElements[0].elID=-1;
00221     for (int i=0; i<numElements; i++) { 
00222       if (theMesh->elem[0].is_valid(i)) {
00223     // find minEdgeLength of i
00224     theMesh->e2n_getAll(i, elemConn);
00225     bool notclear = false;
00226     for(int j=0; j<elemWidth; j++) {
00227       int nd = elemConn[j];
00228       if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
00229       //else cannot look up because connectivity might be screwed up,
00230       //I am hoping that at least one of the nodes will be local
00231     }
00232     if(notclear) continue;
00233     double tmpLen, avgEdgeLength=0.0, 
00234       minEdgeLength = length(elemConn[0], elemConn[1]);
00235     for (int j=0; j<elemWidth-1; j++) {
00236       for (int k=j+1; k<elemWidth; k++) {
00237         tmpLen = length(elemConn[j], elemConn[k]);
00238         if(tmpLen < -1.0) {notclear = true;}
00239         avgEdgeLength += tmpLen;
00240         if (tmpLen < minEdgeLength) minEdgeLength = tmpLen;
00241       }
00242     }
00243     if(notclear) continue;
00244     avgEdgeLength /= 3.0;
00245     qFactor=getAreaQuality(i);
00246     if (((theMesh->elem[0].getMeshSizing(i) > 0.0) &&
00247          (avgEdgeLength < (theMesh->elem[0].getMeshSizing(i)*COARSEN_TOL)))
00248          || (qFactor < QUALITY_MIN)) {
00249         //){
00250       //CkPrintf("Marking elem %d for coarsening\n", i);
00251       Insert(i, qFactor*minEdgeLength, 1);
00252     }
00253       }
00254     }
00255     while (coarsenHeapSize>0) { // loop through the elements
00256       elId=Delete_Min(1);
00257       if ((elId != -1) && (theMesh->elem[0].is_valid(elId))) {
00258     theMesh->e2n_getAll(elId, elemConn);
00259     bool notclear = false;
00260     for(int j=0; j<elemWidth; j++) {
00261       int nd = elemConn[j];
00262       if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
00263       //else cannot look up because connectivity might be screwed up,
00264       //I am hoping that at least one of the nodes will be local
00265     }
00266     if(notclear) continue;
00267     int n1=elemConn[0], n2=elemConn[1];
00268     double tmpLen, avgEdgeLength=0.0, 
00269       minEdgeLength = length(n1, n2);
00270     for (int j=0; j<elemWidth-1; j++) {
00271       for (int k=j+1; k<elemWidth; k++) {
00272         tmpLen = length(elemConn[j], elemConn[k]);
00273         if(tmpLen < -1.0) {notclear = true;}
00274         avgEdgeLength += tmpLen;
00275         if (tmpLen < minEdgeLength) {
00276           minEdgeLength = tmpLen;
00277           n1 = elemConn[j]; n2 = elemConn[k];
00278         }
00279       }
00280     }
00281     if(notclear) continue;
00282     CkAssert(n1!=-1 && n2!=-1);
00283     avgEdgeLength /= 3.0;
00284     qFactor=getAreaQuality(elId);
00285     // coarsen element's short edge
00286     if (((theMesh->elem[0].getMeshSizing(elId) > 0.0) &&
00287          (avgEdgeLength < (theMesh->elem[0].getMeshSizing(elId)*COARSEN_TOL)))
00288           || (qFactor < QUALITY_MIN)) {
00289       //){
00290 
00291       int eNbr = theMesh->e2e_getNbr(elId, theMesh->e2n_getIndex(elId,n1)+
00292                      theMesh->e2n_getIndex(elId,n2)-1, 0);
00293       // determine if eNbr should also be coarsened
00294       if ((eNbr >= 0) && (theMesh->elem[0].is_valid(eNbr))) {
00295         theMesh->e2n_getAll(eNbr, elemConn);
00296         avgEdgeLength=0.0;
00297         for (int j=0; j<elemWidth-1; j++) {
00298           for (int k=j+1; k<elemWidth; k++) {
00299         avgEdgeLength += length(elemConn[j], elemConn[k]);
00300           }
00301         }
00302         avgEdgeLength /= 3.0;
00303         qFactor=getAreaQuality(eNbr);
00304         if (((theMesh->elem[0].getMeshSizing(eNbr) > 0.0) &&
00305          (avgEdgeLength < (theMesh->elem[0].getMeshSizing(eNbr))))
00306         || (qFactor < QUALITY_MIN)) {
00307           //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);
00308           if (theAdaptor->edge_contraction(n1, n2) > 0)  iter_mods++;
00309         }
00310       }
00311       else {
00312         //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);
00313         if (theAdaptor->edge_contraction(n1, n2) > 0)  iter_mods++;
00314       }
00315     }
00316       }
00317       CthYield(); // give other chunks on the same PE a chance
00318     }
00319     mods += iter_mods;
00320 #ifdef ADAPT_VERBOSE
00321     CkPrintf("[%d]ParFUM_Coarsen: %d modifications in pass %d.\n",theMod->idx,iter_mods,pass);
00322 #endif
00323 #ifdef DEBUG_QUALITY
00324     tests(false);
00325 #endif
00326   }
00327 #ifdef ADAPT_VERBOSE
00328   CkPrintf("[%d]ParFUM_Coarsen: %d total modifications over %d passes.\n",theMod->idx,mods,pass);
00329 #endif
00330   delete[] coarsenElements;
00331   return mods;
00332   */
00333 }
00334 
00338 void Bulk_Adapt::ParFUM_AdaptMesh(int qm, int method, double factor, 
00339                    double *sizes)
00340 {
00341   /*
00342   MPI_Comm comm=(MPI_Comm)FEM_chunk::get("FEM_Update_mesh")->defaultComm;
00343   MPI_Barrier(comm);
00344 #ifdef ADAPT_VERBOSE
00345   CkPrintf("[%d]BEGIN: FEM_AdaptMesh...\n",theMod->idx);
00346 #endif
00347 #ifdef DEBUG_QUALITY
00348   tests(true);
00349 #endif
00350   SetMeshSize(method, factor, sizes);
00351   MPI_Barrier(comm);
00352   GradateMesh(GRADATION);
00353   MPI_Barrier(comm);
00354   (void)Refine(qm, method, factor, sizes);
00355   MPI_Barrier(comm);
00356   GradateMesh(GRADATION);
00357   MPI_Barrier(comm);
00358   (void)Coarsen(qm, method, factor, sizes);
00359 #ifdef DEBUG_QUALITY
00360   MPI_Barrier(comm);
00361 #endif
00362   FEM_Repair(qm);
00363   MPI_Barrier(comm);
00364 #ifdef DEBUG_QUALITY
00365   tests(true);
00366   MPI_Barrier(comm);
00367 #endif
00368 #ifdef ADAPT_VERBOSE
00369   CkPrintf("[%d]...END: FEM_AdaptMesh.\n",theMod->idx);
00370 #endif
00371   */
00372 }
00373 
00376 void Bulk_Adapt::ParFUM_Smooth(int qm, int method)
00377 {
00378   CkPrintf("WARNING: Bulk_Adapt::Smooth: Not yet implemented.\n");
00379 }
00380 
00391 /*
00392 void  FEM_Adapt_Algs::FEM_mesh_smooth(FEM_Mesh *meshP, int *nodes, int nNodes, int attrNo) {
00393   vector2d newPos, *coords, *ghostCoords;
00394   int idx, nNod, nGn, gIdxN, *boundVals, nodesInChunk, mesh;
00395   int *adjnodes;
00396 
00397   mesh=FEM_Mesh_default_read();
00398   nodesInChunk = FEM_Mesh_get_length(mesh,FEM_NODE);
00399   nGn = FEM_Mesh_get_length(mesh, FEM_GHOST + FEM_NODE);
00400   
00401   boundVals = new int[nodesInChunk];
00402   coords = new vector2d[nodesInChunk+nGn];
00403 
00404   FEM_Mesh_data(mesh, FEM_NODE, FEM_BOUNDARY, (int*) boundVals, 0, nodesInChunk, FEM_INT, 1);    
00405 
00406   FEM_Mesh_data(mesh, FEM_NODE, attrNo, (double*)coords, 0, nodesInChunk, FEM_DOUBLE, 2);
00407 
00408   IDXL_Layout_t coord_layout = IDXL_Layout_create(IDXL_DOUBLE, 2);
00409   FEM_Update_ghost_field(coord_layout,-1, coords); 
00410   ghostCoords = &(coords[nodesInChunk]);
00411   for (int i=0; i<nNodes; i++)
00412   {
00413     if (nodes==NULL) idx=i;
00414     else idx=nodes[i];
00415     newPos.x=0;
00416     newPos.y=0;
00417     CkAssert(idx<nodesInChunk);  
00418     if (FEM_is_valid(mesh, FEM_NODE, idx) && boundVals[idx]==0) //node must be internal
00419     {
00420       meshP->n2n_getAll(idx, adjnodes, nNod);
00421       for (int j=0; j<nNod; j++) { //for all adjacent nodes, find coords
00422     if (adjnodes[j]<-1) {
00423       gIdxN = FEM_From_ghost_index(adjnodes[j]);
00424       newPos.x += ghostCoords[gIdxN].x;
00425       newPos.y += ghostCoords[gIdxN].y;
00426     }
00427     else {
00428       newPos.x += coords[adjnodes[j]].x;
00429       newPos.y += coords[adjnodes[j]].y;
00430     }     
00431       }
00432       newPos.x/=nNod;
00433       newPos.y/=nNod;
00434       FEM_set_entity_coord2(mesh, FEM_NODE, idx, newPos.x, newPos.y);
00435       delete [] adjnodes;
00436     }
00437   }
00438 
00439   delete [] coords;
00440   delete [] boundVals;
00441 }
00442 */
00443 
00447 void Bulk_Adapt::ParFUM_Repair(int qm)
00448 {
00449   /*
00450   double avgQual = 0.0, minQual = getAreaQuality(0);
00451   int numBadElems = 0;
00452   int elemConn[3];
00453 #ifdef ADAPT_VERBOSE
00454   CkPrintf("[%d]WARNING: ParFUM_Repair: Under construction.\n",theMod->idx);
00455   numElements = theMesh->elem[0].size();
00456   for (int i=0; i<numElements; i++) { 
00457     if (theMesh->elem[0].is_valid(i)) {
00458       double qFactor=getAreaQuality(i);
00459       avgQual += qFactor;
00460       if (qFactor <  QUALITY_MIN) {
00461     numBadElems++;
00462     if (qFactor < minQual) minQual = qFactor;
00463       }
00464     }
00465   }
00466   avgQual /= numElements;
00467   CkPrintf("BEFORE FEM_Repair: Average Element Quality = %2.6f, Min = %2.6f (1.0 is perfect)\n", avgQual, minQual);
00468   //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);
00469 #endif
00470 
00471   int elemWidth = theMesh->elem[0].getConn().width();
00472   int changes=1, totalChanges=0;
00473   int count=0;
00474   while (changes!=0 && count<4) {
00475     count++;
00476     changes = 0;
00477     numElements = theMesh->elem[0].size();
00478     for (int i=0; i<numElements; i++) { 
00479       if (theMesh->elem[0].is_valid(i)) {
00480     double qFactor=getAreaQuality(i);
00481     if (qFactor <  0.75*QUALITY_MIN) {
00482       int elId = i;
00483       theMesh->e2n_getAll(elId, elemConn);
00484       bool notclear = false;
00485       for(int j=0; j<elemWidth; j++) {
00486         int nd = elemConn[j];
00487         if(nd>=0) notclear = theMod->fmLockN[nd].haslocks();
00488         //else cannot look up because connectivity might be screwed up,
00489         //I am hoping that at least one of the nodes will be local
00490       }
00491       if(notclear) continue;
00492       int n1=elemConn[0], n2=elemConn[1];
00493       //too bad.. should not decide without locking!!
00494       //these values might change by the time lock is acquired
00495       double len1 = length(elemConn[0],elemConn[1]);
00496       double len2 = length(elemConn[1],elemConn[2]);
00497       double len3 = length(elemConn[2],elemConn[0]);
00498       if(len1<-1.0 || len2<-1.0 || len3<-1.0) continue;
00499       double avglen=(len1+len2+len3)/3.0;
00500       int maxn1=0, maxn2=1;
00501       int minn1=0, minn2=1;
00502       double maxlen=len1;
00503       int maxed=0;
00504       if(len2>maxlen) {
00505         maxlen = len2;
00506         maxn1 = 1;
00507         maxn2 = 2;
00508         maxed = 1;
00509       }
00510       if(len3>maxlen) {
00511         maxlen = len3;
00512         maxn1 = 2;
00513         maxn2 = 0;
00514         maxed = 2;
00515       }
00516       double minlen = len1;
00517       if(len2<minlen) {
00518         minlen = len2;
00519         minn1 = 1;
00520         minn2 = 2;
00521       }
00522       if(len3<minlen) {
00523         minlen = len3;
00524         minn1 = 2;
00525         minn2 = 0;
00526       }
00527       double otherlen = 3*avglen - maxlen - minlen;
00528       if (maxlen > 0.95*(minlen+otherlen)) { // refine
00529         //decide if this should be a bisect or flip, 
00530         //depends on if the longest edge on this element is also the longest on the
00531         //element sharing this edge, if not, flip it
00532         int nbrEl = theMesh->e2e_getNbr(i,maxed);
00533         double len4=0.0, len5=0.0;
00534         if(nbrEl!=-1) {
00535           int con1[3];
00536           theMesh->e2n_getAll(nbrEl,con1);
00537           int nbrnode=-1;
00538           for(int j=0; j<3; j++) {
00539         if(con1[j]!=elemConn[maxn1] && con1[j]!=elemConn[maxn2]) {
00540           nbrnode = con1[j];
00541           break;
00542         }
00543           }
00544           len4 = length(elemConn[maxn1],nbrnode);
00545           len5 = length(elemConn[maxn2],nbrnode);
00546         }
00547         int success = -1;
00548         if(len4>maxlen || len5>maxlen) {
00549 #ifdef DEBUG_FLIP
00550           success = theAdaptor->edge_flip(elemConn[maxn1], elemConn[maxn2]);
00551 #endif
00552         }
00553         else {
00554           success = theAdaptor->edge_bisect(elemConn[maxn1], elemConn[maxn2]);
00555         }
00556         if (success >= 0) {
00557           //CkPrintf("Refined bad element!\n");
00558           changes++;
00559         }
00560       }
00561       else if (minlen < 0.10*(maxlen+otherlen)) { // coarsen
00562         int success = theAdaptor->edge_contraction(elemConn[minn1], elemConn[minn2]);
00563         if (success >= 0) { 
00564           //CkPrintf("Coarsened bad element!\n");
00565           changes++;
00566         }
00567       }
00568       else {
00569         //CkPrintf("Leaving one bad element alone...\n");
00570       }
00571     }
00572       }
00573     }
00574     totalChanges += changes;
00575   }
00576 
00577 #ifdef ADAPT_VERBOSE
00578   numElements = theMesh->elem[0].size();
00579   numBadElems = 0;
00580   avgQual = 0.0;
00581   minQual = getAreaQuality(0);
00582   for (int i=0; i<numElements; i++) { 
00583     if (theMesh->elem[0].is_valid(i)) {
00584       double qFactor=getAreaQuality(i);
00585       avgQual += qFactor;
00586       if (qFactor <  QUALITY_MIN) {
00587     numBadElems++;
00588     if (qFactor < minQual) minQual = qFactor;
00589       }
00590     }
00591   }
00592   avgQual /= numElements;
00593   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);
00594 #ifdef DEBUG_QUALITY
00595     tests(false);
00596 #endif
00597   //  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);
00598 #endif
00599   */
00600 }
00601 
00607 void Bulk_Adapt::ParFUM_Remesh(int qm, int method, double factor, double *sizes)
00608 {
00609   CkPrintf("WARNING: ParFUM_Remesh: Under construction.\n");
00610 }
00611 
00612 
00615 void Bulk_Adapt::ParFUM_SetReferenceMesh()
00616 {
00617   // TODO: do we need to run this loop for element types other than 0?
00618   double avgLength = 0.0;
00619   int width = theMesh->elem[0].getConn().width();
00620   int numElements = theMesh->elem[0].size();
00621   int elemConn[3];
00622   
00623   for (int i=0; i<numElements; ++i, avgLength=0) {
00624     if(theMesh->elem[0].is_valid(i)) {
00625       theMesh->e2n_getAll(i, elemConn);
00626       for (int j=0; j<width-1; ++j) {
00627     avgLength += length(elemConn[j], elemConn[j+1]);
00628       }
00629       avgLength += length(elemConn[0], elemConn[width-1]);
00630       avgLength /= width;
00631       theMesh->elem[0].setMeshSizing(i, avgLength);      
00632     }
00633   }
00634 }
00635 
00641 void Bulk_Adapt::ParFUM_GradateMesh(double smoothness)
00642 {
00643   const double beta = smoothness;
00644   double maxShock, minShock;
00645   int iteration = 0, updates = 0;
00646   int *adjNodes, *boundNodes;
00647   int nadjNodes, nnodes;
00648   int meshNum = FEM_Mesh_default_read();
00649   //if (smoothness < 1.0) {
00650   //printf("");
00651   //}
00652   nnodes = theMesh->node.size();
00653   boundNodes = new int[nnodes];
00654   FEM_Mesh_data(meshNum, FEM_NODE, FEM_BOUNDARY, 
00655         boundNodes, 0, nnodes, FEM_INT, 1);
00656   
00657   //printf("Running h-shock mesh gradation with beta=%.3f\n", beta);
00658   fflush(NULL);
00659   
00660 #ifndef GRADATION_ITER_LIMIT
00661 #define GRADATION_ITER_LIMIT    10
00662 #endif
00663   
00664   do {
00665     maxShock = 0;
00666     minShock = 1e10;
00667     
00668     for (int node=0; node<nnodes; ++node) {
00669       if (boundNodes[node]!= 0 || !FEM_is_valid(meshNum, FEM_NODE, node))
00670     continue;
00671       //if (!FEM_is_valid(meshNum, FEM_NODE, node))
00672       //    continue;
00673       
00674       theMesh->n2n_getAll(node, adjNodes, nadjNodes);
00675       for (int adjNode=0; adjNode<nadjNodes; ++adjNode) {
00676     double edgelen = length(node, adjNodes[adjNode]);
00677     
00678     // get adjacent elemnents and their sizes
00679     int e1, e2;
00680     theMesh->get2ElementsOnEdge(node, adjNodes[adjNode], &e1, &e2);
00681     if (e1 <= -1 || e2 <= -1) continue; //this will not smooth across boundaries
00682     
00683     double s1, s2;
00684     s1 = theMesh->elem[0].getMeshSizing(e1);
00685     s2 = theMesh->elem[0].getMeshSizing(e2);
00686     if (s1 <= 0 || s2 <= 0) continue;
00687     
00688     // h-shock=max(size ratio)^(1/edge length)
00689     CkAssert(s1 >= 0 && s2 >= 0 && "Bad size");
00690     //CkAssert(edgelen > 1e-6 && "Length 0 edge");
00691     CkAssert(edgelen == edgelen && "Length inf edge");
00692     
00693     double ratio = (s1 > s2) ? s1/s2 : s2/s1;
00694     CkAssert (ratio >= 1.0 && ratio == ratio && "Bad ratio");
00695     
00696     // WARNING WARNING WARNING
00697     // TEST ONLY, THIS IS NOT CORRECT
00698     if (ratio > beta) {
00699       if (s1 > s2) {
00700         theMesh->elem[0].setMeshSizing(e1, s1 - (s1-s2)/3);
00701         theMesh->elem[0].setMeshSizing(e2, s2 + (s1-s2)/3);
00702       } else {
00703         theMesh->elem[0].setMeshSizing(e2, s2 - (s2-s1)/3);
00704         theMesh->elem[0].setMeshSizing(e1, s1 + (s2-s1)/3);
00705       }
00706       updates++;
00707     }
00708     if (ratio > maxShock) maxShock = ratio;
00709     if (ratio < minShock) minShock = ratio;
00710     
00711     
00713     //double hs = pow(ratio, 1.0/edgelen);
00714     //
00715     //if (hs > maxShock) maxShock = hs;
00716     //if (hs < minShock) minShock = hs;
00717     
00722     //
00723     //if (hs > beta) {
00724     //    double etasq = pow(beta, edgelen) / ratio;
00725     //    etasq *= etasq;
00726     
00727     //    //if (hs > 100) {
00728     //    //    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);
00729     //    //    abort();
00730     //    //}
00731     //    
00732     //    if (s1 > s2) {
00733     //        theMesh->elem[0].setMeshSizing(e1, s1 / etasq);
00734     //    } else {
00735     //        theMesh->elem[0].setMeshSizing(e2, s2 / etasq);
00736     //    }
00737     //    updates++;
00738     //}
00739     
00740       }
00741       if (nadjNodes > 0)
00742     delete[] adjNodes;
00743     } 
00744     
00745     //printf("Finished iteration %d\n", iteration);
00746     //printf("Max shock:%8.3f\n", maxShock);
00747     //printf("Min shock:%8.3f\n", minShock);
00748     //printf("Target:%8.3f\n", beta);
00749     
00750   } while (maxShock > beta && ++iteration < GRADATION_ITER_LIMIT);
00751   
00752   //printf("%d total updates in %d iterations in GradateMesh\n", updates, iteration);
00753   fflush(NULL);
00754   delete[] boundNodes;
00755   return;
00756 }
00757 
00759 void Bulk_Adapt::SetMeshSize(int method, double factor, double *sizes)
00760 {
00761   int numNodes = theMesh->node.size();
00762   int numElements = theMesh->elem[0].size();
00763   int elemConn[3];
00764 
00765   if (method == 0) { // set uniform sizing specified in factor
00766 
00767     for (int i=0; i<numElements; i++) {
00768       theMesh->elem[0].setMeshSizing(i, factor);
00769     }
00770     //CkPrintf("ParFUM_SetMeshSize: UNIFORM %4.6e\n", factor);
00771   }
00772   else if (method == 1) { // copy sizing from array
00773     for (int i=0; i<numElements; i++) {
00774       if (sizes[i] > 0.0) {
00775     theMesh->elem[0].setMeshSizing(i, sizes[i]);
00776       }
00777     }
00778     //CkPrintf("ParFUM_SetMeshSize: SIZES input\n");
00779   }
00780   else if (method == 2) { // calculate current sizing and scale by factor
00781     double avgEdgeLength = 0.0;
00782     int width = theMesh->elem[0].getConn().width();
00783     int numEdges=3;
00784     if (dim==3) numEdges=6;
00785     for (int i=0; i<numElements; i++) {
00786       if(theMesh->elem[0].is_valid(i)) {
00787     theMesh->e2n_getAll(i, elemConn);
00788     for (int j=0; j<width-1; j++) {
00789       for (int k=j+1; k<width; k++) {
00790         avgEdgeLength += length(elemConn[j], elemConn[k]);
00791       }
00792     }
00793     avgEdgeLength += length(elemConn[0], elemConn[width-1]);
00794     avgEdgeLength /= (double)numEdges;
00795     theMesh->elem[0].setMeshSizing(i, factor*avgEdgeLength);
00796       }
00797     }
00798     //CkPrintf("ParFUM_SetMeshSize: CALCULATED & SCALED \n");
00799   }
00800   else if (method == 3) { // scale existing sizes by array sizes
00801     for (int i=0; i<numElements; i++) {
00802       if (sizes[i] > 0.0) {
00803     theMesh->elem[0].setMeshSizing(i, sizes[i]*theMesh->elem[0].getMeshSizing(i));
00804       }
00805     }
00806   }
00807   else if (method == 4) { // scale existing sizes by factor
00808     for (int i=0; i<numElements; i++) {
00809       theMesh->elem[0].setMeshSizing(i, factor*theMesh->elem[0].getMeshSizing(i));
00810     }
00811   }
00812   else if (method == 5) { // mesh sizing has been set independently; use as is
00813     //CkPrintf("ParFUM_SetMeshSize: USE EXISTING SIZES \n");
00814   }
00815   //  CkPrintf("Current mesh sizing: ");
00816   //for (int i=0; i<numElements; i++) {
00817   //CkPrintf("%4.6e ", theMesh->elem[0].getMeshSizing(i));
00818   //}
00819 }
00820 
00821 
00822 //HELPER functions
00823 double Bulk_Adapt::length(int n1, int n2)
00824 {
00825   double coordsn1[2], coordsn2[2];
00826   getCoord(n1, coordsn1);
00827   getCoord(n2, coordsn2);
00828 
00829   double ret = length(coordsn1, coordsn2);
00830   return ret;
00831 }
00832 
00833 double Bulk_Adapt::length3D(int n1, int n2) {
00834   double coordsn1[3], coordsn2[3];
00835   getCoord(n1, coordsn1);
00836   getCoord(n2, coordsn2);
00837 
00838   double ret = length3D(coordsn1, coordsn2);
00839   return ret;
00840 }
00841 
00842 double Bulk_Adapt::length(double *n1_coord, double *n2_coord) { 
00843   double d, ds_sum=0.0;
00844 
00845   for (int i=0; i<dim; i++) {
00846     if(n1_coord[i]<-1.0 || n2_coord[i]<-1.0) return -2.0;
00847     d = n1_coord[i] - n2_coord[i];
00848     ds_sum += d*d;
00849   }
00850   return (sqrt(ds_sum));
00851 }
00852 
00853 double Bulk_Adapt::length3D(double *n1_coord, double *n2_coord) { 
00854   double d, ds_sum=0.0;
00855 
00856   for (int i=0; i<3; i++) {
00857     if(n1_coord[i]<-1.0 || n2_coord[i]<-1.0) return -2.0;
00858     d = n1_coord[i] - n2_coord[i];
00859     ds_sum += d*d;
00860   }
00861   return (sqrt(ds_sum));
00862 }
00863 
00864 void Bulk_Adapt::findEdgeLengths(int elemID, double *avgEdgeLength, double *maxEdgeLength, int *maxEdge, 
00865                  double *minEdgeLength, int *minEdge)
00866 {
00867   FEM_Elem &elem = theMesh->elem[elemType]; // elem is local elements
00868   int *conn = elem.connFor(elemID);
00869   //CkPrintf("conn[%d]:[%d,%d,%d,%d]\n", elemID, conn[0], conn[1], conn[2], conn[3]);
00870   FEM_DataAttribute *coord = theMesh->node.getCoord(); // all local coords
00871   double *n0co = (coord->getDouble()).getRow(conn[0]);
00872   double *n1co = (coord->getDouble()).getRow(conn[1]);
00873   double *n2co = (coord->getDouble()).getRow(conn[2]);
00874   double *n3co = (coord->getDouble()).getRow(conn[3]);
00875   //CkPrintf(" node %d = (%6.4f, %6.4f, %6.4f)\n", conn[0], n0co[0], n0co[1], n0co[2]);
00876   //CkPrintf(" node %d = (%6.4f, %6.4f, %6.4f)\n", conn[1], n1co[0], n1co[1], n1co[2]);
00877   //CkPrintf(" node %d = (%6.4f, %6.4f, %6.4f)\n", conn[2], n2co[0], n2co[1], n2co[2]);
00878   //CkPrintf(" node %d = (%6.4f, %6.4f, %6.4f)\n", conn[3], n3co[0], n3co[1], n3co[2]);
00879 
00880   // assuming tets for now
00881   double edgeLengths[6];
00882   edgeLengths[0] = length3D(conn[0], conn[1]);
00883   edgeLengths[1] = length3D(conn[0], conn[2]);
00884   edgeLengths[2] = length3D(conn[0], conn[3]);
00885   edgeLengths[3] = length3D(conn[1], conn[2]);
00886   edgeLengths[4] = length3D(conn[1], conn[3]);
00887   edgeLengths[5] = length3D(conn[2], conn[3]);
00888   //CkPrintf("edge 0 has length %6.4f\n", edgeLengths[0]);
00889   (*maxEdgeLength) = (*minEdgeLength) = (*avgEdgeLength) = edgeLengths[0];
00890   (*maxEdge) = (*minEdge) = 0;
00891   for (int i=1; i<6; i++) {
00892     //CkPrintf("edge %d has length %6.4f\n", i, edgeLengths[i]);
00893     (*avgEdgeLength) += edgeLengths[i];
00894     if (edgeLengths[i] > (*maxEdgeLength)) {
00895       (*maxEdgeLength) = edgeLengths[i]; 
00896       (*maxEdge) = i;
00897     }
00898     else if (edgeLengths[i] < (*minEdgeLength)) {
00899       (*minEdgeLength) = edgeLengths[i];
00900       (*minEdge) = i;
00901     }
00902   }
00903   (*avgEdgeLength) /= 6;
00904   CkAssert((*minEdgeLength) > 0.0);
00905 }
00906 
00907 
00908 double Bulk_Adapt::getArea(int n1, int n2, int n3)
00909 {
00910   double coordsn1[2], coordsn2[2], coordsn3[2];
00911   getCoord(n1, coordsn1);
00912   getCoord(n2, coordsn2);
00913   getCoord(n3, coordsn3);
00914 
00915   double ret = getArea(coordsn1, coordsn2, coordsn3);
00916   return ret;
00917 }
00918 
00919 double Bulk_Adapt::getArea(double *n1_coord, double *n2_coord, double *n3_coord) {
00920   double area=0.0;
00921   double aLen, bLen, cLen, sLen, d, ds_sum;
00922 
00923   ds_sum = 0.0;
00924   for (int i=0; i<2; i++) {
00925     d = n1_coord[i] - n2_coord[i];
00926     ds_sum += d*d;
00927   }
00928   aLen = sqrt(ds_sum);
00929   ds_sum = 0.0;
00930   for (int i=0; i<2; i++) {
00931     d = n2_coord[i] - n3_coord[i];
00932     ds_sum += d*d;
00933   }
00934   bLen = sqrt(ds_sum);
00935   ds_sum = 0.0;
00936   for (int i=0; i<2; i++) {
00937     d = n3_coord[i] - n1_coord[i];
00938     ds_sum += d*d;
00939   }
00940   cLen = sqrt(ds_sum);
00941   sLen = (aLen+bLen+cLen)/2;
00942   if(sLen-aLen < 0) return 0.0;
00943   else if(sLen-bLen < 0) return 0.0;
00944   else if(sLen-cLen < 0) return 0.0; //area too small to note
00945   return (sqrt(sLen*(sLen-aLen)*(sLen-bLen)*(sLen-cLen)));
00946 }
00947 
00948 double Bulk_Adapt::getSignedArea(int n1, int n2, int n3) {
00949   double coordsn1[2], coordsn2[2], coordsn3[2];
00950   getCoord(n1, coordsn1);
00951   getCoord(n2, coordsn2);
00952   getCoord(n3, coordsn3);
00953 
00954   double ret = getSignedArea(coordsn1, coordsn2, coordsn3);
00955   return ret;
00956 }
00957 
00958 double Bulk_Adapt::getSignedArea(double *n1_coord, double *n2_coord, double *n3_coord) {
00959   double area=0.0;
00960   double vec1_x, vec1_y, vec2_x, vec2_y;
00961 
00962   vec1_x = n1_coord[0] - n2_coord[0];
00963   vec1_y = n1_coord[1] - n2_coord[1];
00964   vec2_x = n3_coord[0] - n2_coord[0];
00965   vec2_y = n3_coord[1] - n2_coord[1];
00966 
00967   area = vec1_x*vec2_y - vec2_x*vec1_y;
00968   return area;
00969 }
00970 
00971 int Bulk_Adapt::getCoord(int n1, double *crds) {
00972   FEM_DataAttribute *coord = theMesh->node.getCoord(); // entire local coords
00973   double *nodeCoords = (coord->getDouble()).getRow(n1); // ptrs to ACTUAL coords!
00974   crds[0] = nodeCoords[0];
00975   crds[1] = nodeCoords[1];
00976   if (dim == 3) 
00977     crds[2] = nodeCoords[2];
00978   return 1;
00979 }
00980 
00981 int Bulk_Adapt::getShortestEdge(int n1, int n2, int n3, int* shortestEdge) {
00982   //note that getCoord might be a remote call, which means 
00983   //it might not have the same value in the mem
00984   //location if the memory is reused by someone else meanwhile!
00985   double coordsn1[2], coordsn2[2], coordsn3[2];
00986   getCoord(n1, coordsn1);
00987   getCoord(n2, coordsn2);
00988   getCoord(n3, coordsn3);
00989 
00990   double aLen = length(coordsn1, coordsn2);
00991   int shortest = 0;
00992 
00993   double bLen = length(coordsn2, coordsn3);
00994   if(bLen < aLen) shortest = 1;
00995 
00996   double cLen = length(coordsn3, coordsn1);
00997   if((cLen < aLen) && (cLen < bLen)) shortest = 2;
00998 
00999   if(shortest==0) {
01000     shortestEdge[0] = n1;
01001     shortestEdge[1] = n2;
01002   }
01003   else if (shortest==1) {
01004     shortestEdge[0] = n2;
01005     shortestEdge[1] = n3;
01006   }
01007   else {
01008     shortestEdge[0] = n3;
01009     shortestEdge[1] = n1;
01010   }
01011   return 1;
01012 }
01013 
01014 
01015 
01019 double Bulk_Adapt::getAreaQuality(int elem)
01020 {
01021   double f, q, len[3];
01022   int n[3];
01023   double currentArea;
01024   double coordsn1[2], coordsn2[2], coordsn3[2];
01025 
01026   theMesh->e2n_getAll(elem, n);
01027   getCoord(n[0], coordsn1);
01028   getCoord(n[1], coordsn2);
01029   getCoord(n[2], coordsn3);
01030 
01031   currentArea = getArea(coordsn1, coordsn2, coordsn3);
01032 
01033   len[0] = length(coordsn1, coordsn2);
01034   len[1] = length(coordsn2, coordsn3);
01035   len[2] = length(coordsn3, coordsn1);
01036   f = 4.0*sqrt(3.0); //proportionality constant
01037   q = (f*currentArea)/(len[0]*len[0]+len[1]*len[1]+len[2]*len[2]);
01038   return q;
01039 }
01040 
01045 void Bulk_Adapt::ensureQuality(int n1, int n2, int n3) {
01046   double coordsn1[2], coordsn2[2], coordsn3[2];
01047   getCoord(n1, coordsn1);
01048   getCoord(n2, coordsn2);
01049   getCoord(n3, coordsn3);
01050   double area = getSignedArea(coordsn1, coordsn2, coordsn3);
01051   double len1 = length(coordsn1, coordsn2);
01052   double len2 = length(coordsn2, coordsn3);
01053   double len3 = length(coordsn3, coordsn1);
01054   double max = len1;
01055   if(len2>max) max = len2;
01056   if(len3>max) max = len3;
01057   //shortest edge
01058   double min = len1;
01059   if(len2<min) min = len2;
01060   if(len3<min) min = len3;
01061   double shortest_al = fabs(area/max);
01062   double largestR = max/shortest_al;
01063   CkAssert(largestR<=100.0 && -area > SLIVERAREA);
01064 }
01065 
01068 bool Bulk_Adapt::controlQualityF(int n1, int n2, int n3, int n4) {
01069   //n1 or n2 will be replaced by n4
01070   double coordsn1[2], coordsn2[2], coordsn3[2], coordsn4[2];
01071   if(n4==-1) return false;
01072   getCoord(n1, coordsn1);
01073   getCoord(n2, coordsn2);
01074   getCoord(n3, coordsn3);
01075   getCoord(n4, coordsn4);
01076   bool flag = false;
01077   if(!flag) flag = controlQualityC(coordsn3,coordsn1,coordsn2,coordsn4);
01078   if(!flag) flag = controlQualityC(coordsn4,coordsn2,coordsn1,coordsn3);
01079   return flag;
01080 }
01081 
01082 bool Bulk_Adapt::controlQualityR(double *coordsn1, double *coordsn2, double *coordsn3) {
01083   double area = getArea(coordsn1, coordsn2, coordsn3);
01084   //do some quality preservation
01085   double len1 = length(coordsn1, coordsn2);
01086   double len2 = length(coordsn2, coordsn3);
01087   double len3 = length(coordsn3, coordsn1);
01088   //longest edge
01089   double max = len1;
01090   if(len2>max) max = len2;
01091   if(len3>max) max = len3;
01092   double shortest_al = area/max;
01093   double largestR = max/shortest_al;
01094   if(largestR>50.0) return true;
01095   else return false;
01096 }
01097 
01102 bool Bulk_Adapt::controlQualityC(int n1, int n2, int n3, double *n4_coord) {
01103   //n3 is the node to be deleted, n4 is the new node to be added
01104   double coordsn1[2], coordsn2[2], coordsn3[2];
01105   getCoord(n1, coordsn1);
01106   getCoord(n2, coordsn2);
01107   getCoord(n3, coordsn3);
01108   return controlQualityC(coordsn1,coordsn2,coordsn3,n4_coord);
01109 }
01110 
01111 bool Bulk_Adapt::controlQualityC(double *coordsn1, double *coordsn2, double *coordsn3, double *n4_coord) {
01112   double ret_old = getSignedArea(coordsn1, coordsn2, coordsn3);
01113   double ret_new = getSignedArea(coordsn1, coordsn2, n4_coord);
01114   //do some quality preservation
01115   double len1 = length(coordsn1, coordsn2);
01116   double len2 = length(coordsn2, n4_coord);
01117   double len3 = length(n4_coord, coordsn1);
01118   //longest edge
01119   double max = len1;
01120   if(len2>max) max = len2;
01121   if(len3>max) max = len3;
01122   //shortest edge
01123   double min = len1;
01124   if(len2<min) min = len2;
01125   if(len3<min) min = len3;
01126   double shortest_al = ret_new/max;
01127   double largestR = max/shortest_al;
01128   if(ret_old > SLIVERAREA && ret_new < -SLIVERAREA) return true; //it is a flip
01129   else if(ret_old < -SLIVERAREA && ret_new > SLIVERAREA) return true; //it is a flip
01130   else if(fabs(ret_new) < SLIVERAREA) return true; // it is a sliver
01131   //else if(fabs(shortest_al) < 1e-5) return true; //shortest altitude is too small
01132   //else if(fabs(min) < 1e-5) return true; //shortest edge is too small
01133   else if(fabs(largestR)>50.0) return true;
01134   else return false;
01135 }
01136 
01137 
01138 
01144 bool Bulk_Adapt::flipOrBisect(int elId, int n1, int n2, int maxEdgeIdx, double maxlen) {
01145   //return true if it should flip
01146   int nbrEl = theMesh->e2e_getNbr(elId,maxEdgeIdx);
01147   double len4=0.0, len5=0.0;
01148   if(nbrEl!=-1) {
01149     int con1[3];
01150     theMesh->e2n_getAll(nbrEl,con1);
01151     int nbrnode=-1;
01152     for(int j=0; j<3; j++) {
01153       if(con1[j]!=n1 && con1[j]!=n2) {
01154     nbrnode = con1[j];
01155     break;
01156       }
01157     }
01158     len4 = length(n1,nbrnode);
01159     len5 = length(n2,nbrnode);
01160   }
01161   if(len4>1.2*maxlen || len5>1.2*maxlen) {
01162     return true;
01163   }
01164   else return false;
01165 }
01166
libs/ck-libs/ParFUM/bulk_adapt.C
