7601a840d2f32a98cd3c0d37171e0bf7b80d84f9
[charm.git] / src / libs / ck-libs / fem / fem_adapt_algs.C
1 /* File: fem_adapt_algs.C
2  * Authors: Terry Wilmarth, Nilesh Choudhury
3  * 
4  */
5
6
7 // This module implements high level mesh adaptivity algorithms that make use 
8 // of the primitive mesh adaptivity operations provided by fem_adapt(_new).
9 // Ask: TLW
10 #include "fem_adapt_algs.h"
11 #include "fem_mesh_modify.h"
12
13 #define MINAREA 1.0e-12
14 #define MAXAREA 1.0e12
15
16 CtvDeclare(FEM_Adapt_Algs *, _adaptAlgs);
17
18 FEM_Adapt_Algs::FEM_Adapt_Algs(FEM_Mesh *m, femMeshModify *fm, int dimension) 
19
20   theMesh = m; 
21   theMod = fm; 
22   dim = dimension; 
23   //theAdaptor = theMod->fmAdapt;
24   theAdaptor = theMod->fmAdaptL;
25 }
26
27 /* Perform refinements on a mesh.  Tries to maintain/improve element quality as
28    specified by a quality measure qm; if method = 0, refine areas with size 
29    larger than factor down to factor; if method = 1, refine elements down to 
30    sizes specified in sizes array; Negative entries in size array indicate no 
31    refinement. */
32 void FEM_Adapt_Algs::FEM_Refine(int qm, int method, double factor, 
33                                 double *sizes)
34 {
35   numNodes = theMesh->node.size();
36   numElements = theMesh->elem[0].size();
37   (void)Refine(qm, method, factor, sizes);
38 }
39
40 /* Performs refinement; returns number of modifications */
41 int FEM_Adapt_Algs::Refine(int qm, int method, double factor, double *sizes)
42 {
43   // loop through elemsToRefine
44   int elId, mods=0, iter_mods=1;
45   SetMeshSize(method, factor, sizes);
46   refineElements = refineStack = NULL;
47   refineTop = refineHeapSize = 0;
48   while (iter_mods != 0) {
49     iter_mods=0;
50     numNodes = theMesh->node.size();
51     numElements = theMesh->elem[0].size();
52     // sort elements to be refined by quality into elemsToRefine
53     if (refineStack) delete [] refineStack;
54     refineStack = new elemHeap[numElements];
55     if (refineElements) delete [] refineElements;
56     refineElements = new elemHeap[numElements+1];
57     for (int i=0; i<numElements; i++) { 
58       if (theMesh->elem[0].is_valid(i)) {
59         // find maxEdgeLength of i
60         int *eConn = (int*)malloc(3*sizeof(int));
61         double tmpLen, maxEdgeLength;
62         theMesh->e2n_getAll(i, eConn);
63         maxEdgeLength = length(eConn[0], eConn[1]);
64         tmpLen = length(eConn[1], eConn[2]);
65         if (tmpLen > maxEdgeLength) maxEdgeLength = tmpLen;
66         tmpLen = length(eConn[2], eConn[0]);
67         if (tmpLen > maxEdgeLength) maxEdgeLength = tmpLen;
68         double qFactor=getAreaQuality(i);
69         if ((theMesh->elem[0].getMeshSizing(i) > 0.0) &&
70             (maxEdgeLength > (theMesh->elem[0].getMeshSizing(i)*REFINE_TOL))) {
71           Insert(i, qFactor*(1.0/maxEdgeLength), 0);
72         }
73       }
74     }
75     while (refineHeapSize>0 || refineTop > 0) { // loop through the elements
76       if (refineTop>0) {
77         refineTop--;
78         elId=refineStack[refineTop].elID;
79       }
80       else  elId=Delete_Min(0);
81       if ((elId != -1) && (theMesh->elem[0].is_valid(elId))) {
82         int *eConn = (int*)malloc(3*sizeof(int));
83         double tmpLen, maxEdgeLength;
84         theMesh->e2n_getAll(elId, eConn);
85         maxEdgeLength = length(eConn[0], eConn[1]);
86         tmpLen = length(eConn[1], eConn[2]);
87         if (tmpLen > maxEdgeLength) maxEdgeLength = tmpLen;
88         tmpLen = length(eConn[2], eConn[0]);
89         if (tmpLen > maxEdgeLength) maxEdgeLength = tmpLen;
90         if ((theMesh->elem[0].getMeshSizing(elId) > 0.0) &&
91             (maxEdgeLength > (theMesh->elem[0].getMeshSizing(elId)*REFINE_TOL))) {
92           (void)refine_element_leb(elId); // refine the element
93           iter_mods++;
94         }
95       }
96       CthYield(); // give other chunks on the same PE a chance
97     }
98     mods += iter_mods;
99     CkPrintf("ParFUM_Refine: %d modifications in last pass.\n", iter_mods);
100   }
101   CkPrintf("ParFUM_Refine: %d total modifications.\n", mods);
102   return mods;
103 }
104
105 /* Perform coarsening on a mesh.  Tries to maintain/improve element quality as 
106    specified by a quality measure qm; if method = 0, coarsen areas with size 
107    smaller than factor up to factor; if method = 1, coarsen elements up to 
108    sizes specified in sizes array; Negative entries in size array indicate no 
109    coarsening. */
110 void FEM_Adapt_Algs::FEM_Coarsen(int qm, int method, double factor, 
111                                  double *sizes)
112 {
113   CkPrintf("WARNING: ParFUM_Coarsen: Under construction.\n");
114   numNodes = theMesh->node.size();
115   numElements = theMesh->elem[0].size();
116   (void)Coarsen(qm, method, factor, sizes);
117 }
118
119 /* Performs coarsening; returns number of modifications */
120 int FEM_Adapt_Algs::Coarsen(int qm, int method, double factor, double *sizes)
121 {
122   // loop through elemsToRefine
123   int elId, mods=0, iter_mods=1, pass=0;
124   double qFactor;
125   SetMeshSize(method, factor, sizes);
126   coarsenElements = NULL;
127   coarsenHeapSize = 0;
128   while (iter_mods != 0) {
129     iter_mods=0;
130     pass++;
131     numNodes = theMesh->node.size();
132     numElements = theMesh->elem[0].size();
133     // sort elements to be refined by quality into elemsToRefine
134     if (coarsenElements) delete [] coarsenElements;
135     coarsenElements = new elemHeap[numElements+1];
136     coarsenElements[0].len=-2.0;
137     coarsenElements[0].elID=-1;
138     for (int i=0; i<numElements; i++) { 
139       if (theMesh->elem[0].is_valid(i)) {
140         // find minEdgeLength of i
141         int *eConn = (int*)malloc(3*sizeof(int));
142         double tmpLen, minEdgeLength, avgEdgeLength;
143         theMesh->e2n_getAll(i, eConn);
144         avgEdgeLength = minEdgeLength = length(eConn[0], eConn[1]);
145         tmpLen = length(eConn[1], eConn[2]);
146         avgEdgeLength += tmpLen;
147         if (tmpLen < minEdgeLength) minEdgeLength = tmpLen;
148         tmpLen = length(eConn[2], eConn[0]);
149         avgEdgeLength += tmpLen;
150         if (tmpLen < minEdgeLength) minEdgeLength = tmpLen;
151         qFactor=getAreaQuality(i);
152         avgEdgeLength /= 3.0;
153         if (((theMesh->elem[0].getMeshSizing(i) > 0.0) &&
154              (avgEdgeLength < (theMesh->elem[0].getMeshSizing(i)*COARSEN_TOL)))
155             || (qFactor < QUALITY_MIN)) {
156           Insert(i, qFactor*minEdgeLength, 1);
157         }
158       }
159     }
160     while (coarsenHeapSize>0) { // loop through the elements
161       elId=Delete_Min(1);
162       if ((elId != -1) && (theMesh->elem[0].is_valid(elId))) {
163         int *eConn = (int*)malloc(3*sizeof(int));
164         double tmpLen, minEdgeLength, avgEdgeLength;
165         int n1, n2;
166         theMesh->e2n_getAll(elId, eConn);
167         avgEdgeLength = minEdgeLength = length(eConn[0], eConn[1]);
168         n1 = eConn[0]; n2 = eConn[1];
169         tmpLen = length(eConn[1], eConn[2]);
170         avgEdgeLength += tmpLen;
171         if (tmpLen < minEdgeLength) { 
172           minEdgeLength = tmpLen;
173           n1 = eConn[1]; n2 = eConn[2];
174         }
175         tmpLen = length(eConn[2], eConn[0]);
176         avgEdgeLength += tmpLen;
177         if (tmpLen < minEdgeLength) {
178           minEdgeLength = tmpLen;
179           n1 = eConn[2]; n2 = eConn[0];
180         }
181         avgEdgeLength /= 3.0;
182         qFactor=getAreaQuality(elId);
183         // coarsen element's short edge
184         if (((theMesh->elem[0].getMeshSizing(elId) > 0.0) &&
185              (avgEdgeLength < (theMesh->elem[0].getMeshSizing(elId)*COARSEN_TOL)))
186             || (qFactor < QUALITY_MIN)) {
187           if (theAdaptor->edge_contraction(n1, n2) > 0)  iter_mods++;
188         }
189       }
190       CthYield(); // give other chunks on the same PE a chance
191     }
192     mods += iter_mods;
193     CkPrintf("ParFUM_Coarsen: %d modifications in pass %d.\n", iter_mods, pass);
194   }
195   CkPrintf("ParFUM_Coarsen: %d total modifications over %d passes.\n", mods, pass);
196   return mods;
197 }
198
199 /* Smooth the mesh using method according to some quality measure qm */
200 void FEM_Adapt_Algs::FEM_Smooth(int qm, int method)
201 {
202   CkPrintf("WARNING: ParFUM_Smooth: Not yet implemented.\n");
203 }
204
205 /* Repair the mesh according to some quality measure qm */
206 void FEM_Adapt_Algs::FEM_Repair(int qm)
207 {
208   CkPrintf("WARNING: ParFUM_Repair: Not yet implemented.\n");
209 }
210
211 /* Remesh entire mesh according to quality measure qm. If method = 0, set 
212    entire mesh size to factor; if method = 1, use sizes array; if method = 2, 
213    uses existing regional sizes and scale by factor*/
214 void FEM_Adapt_Algs::FEM_Remesh(int qm, int method, double factor, 
215                                 double *sizes)
216 {
217   CkPrintf("WARNING: ParFUM_Remesh: Under construction.\n");
218 }
219
220 /* Set sizes on elements throughout the mesh; note: size is edge length */
221 void FEM_Adapt_Algs::SetMeshSize(int method, double factor, double *sizes)
222 {
223   if (method == 0) {
224     for (int i=0; i<numElements; i++) {
225       theMesh->elem[0].setMeshSizing(i, factor);
226     }
227   }
228   else if (method == 1) {
229     theMesh->elem[0].setMeshSizing(sizes);
230   }
231   else if (method == 2) {
232     double avgEdgeLength = 0.0;
233     for (int i=0; i<numElements; i++) {
234       int eConn[4];
235       theMesh->e2n_getAll(i, eConn);
236       avgEdgeLength = length(eConn[0], eConn[1]) + length(eConn[1], eConn[2]) +
237         length(eConn[2], eConn[0]);
238       if (dim == 3) {
239         avgEdgeLength += length(eConn[0], eConn[3]) + length(eConn[1], eConn[3]) +
240           length(eConn[2], eConn[3]);
241         avgEdgeLength /= 6.0;
242       }
243       else {
244         avgEdgeLength /= 3.0;
245       }
246       theMesh->elem[0].setMeshSizing(i, factor*avgEdgeLength);
247     }
248   }
249   else if (method == 3) { // scale existing mesh size
250     for (int i=0; i<numElements; i++) {
251       theMesh->elem[0].setMeshSizing(i, factor*theMesh->elem[0].getMeshSizing(i));
252     }
253   }
254   else if (method == 4) { // mesh sizing has been set independently; use as is
255     return;
256   }
257 }
258
259 int FEM_Adapt_Algs::simple_refine(double targetA, double xmin, double ymin, double xmax, double ymax) {
260   int noEle = theMesh->elem[0].size();
261   int *con = (int*)malloc(3*sizeof(int));
262   double *areas = (double*)malloc(noEle*sizeof(double));
263   int *map1 = (int*)malloc(noEle*sizeof(int));
264   double *n1_coord = (double*)malloc(2*sizeof(double));
265   double *n2_coord = (double*)malloc(2*sizeof(double));
266   double *n3_coord = (double*)malloc(2*sizeof(double));
267
268   for(int i=0; i<noEle; i++) {
269     if(theMesh->elem[0].is_valid(i)) {
270       theMesh->e2n_getAll(i,con,0);
271       getCoord(con[0], n1_coord);
272       getCoord(con[1], n2_coord);
273       getCoord(con[2], n3_coord);
274       //do a refinement only if it has any node within the refinement box
275       if((n1_coord[0]<xmax && n1_coord[0]>xmin && n1_coord[1]<ymax && n1_coord[1]>ymin) || (n2_coord[0]<xmax && n2_coord[0]>xmin && n2_coord[1]<ymax && n2_coord[1]>ymin) || (n3_coord[0]<xmax && n3_coord[0]>xmin && n3_coord[1]<ymax && n3_coord[1]>ymin)) {
276         areas[i] = getArea(n1_coord, n2_coord, n3_coord);
277         } 
278       else {
279         areas[i] = MINAREA; //make it believe that this triangle does not need refinement
280       }
281     } else {
282       areas[i] = MINAREA;
283     }
284     map1[i] = i;
285   }
286
287   for(int i=0; i<noEle; i++) {
288     for(int j=i+1; j<noEle; j++) {
289       if(areas[j] > areas[i]) {
290         double tmp = areas[j];
291         areas[j] = areas[i];
292         areas[i] = tmp;
293         int t = map1[j];
294         map1[j] = map1[i];
295         map1[i] = t;
296       }
297     }
298   }
299
300   for(int i=0; i<noEle; i++) {
301     if(theMesh->elem[0].is_valid(map1[i])) {
302       if(areas[i] > targetA) {
303         refine_element_leb(map1[i]);
304       }
305     }
306   }
307   free(con);
308   free(areas);
309   free(map1);
310   free(n1_coord);
311   free(n2_coord);
312   free(n3_coord);
313
314   return 1;
315 }
316
317 int FEM_Adapt_Algs::simple_coarsen(double targetA, double xmin, double ymin, double xmax, double ymax) {
318   int noEle = theMesh->elem[0].size();
319   int *con = (int*)malloc(3*sizeof(int));
320   double *areas = (double*)malloc(noEle*sizeof(double));
321   int *map1 = (int*)malloc(noEle*sizeof(int));
322   double *n1_coord = (double*)malloc(2*sizeof(double));
323   double *n2_coord = (double*)malloc(2*sizeof(double));
324   double *n3_coord = (double*)malloc(2*sizeof(double));
325   int *shortestEdge = (int *)malloc(2*sizeof(int));
326   bool adapted = true;
327
328   while(adapted) {
329     adapted = false;
330     for(int i=0; i<noEle; i++) {
331       if(theMesh->elem[0].is_valid(i)) {
332         theMesh->e2n_getAll(i,con,0);
333         getCoord(con[0], n1_coord);
334         getCoord(con[1], n2_coord);
335         getCoord(con[2], n3_coord);
336         //do a coarsening only if it has any node within the coarsen box
337         if((n1_coord[0]<xmax && n1_coord[0]>xmin && n1_coord[1]<ymax && n1_coord[1]>ymin) || (n2_coord[0]<xmax && n2_coord[0]>xmin && n2_coord[1]<ymax && n2_coord[1]>ymin) || (n3_coord[0]<xmax && n3_coord[0]>xmin && n3_coord[1]<ymax && n3_coord[1]>ymin)) {
338           areas[i] = getArea(n1_coord, n2_coord, n3_coord);
339         } 
340         else {
341           areas[i] = MAXAREA; //make it believe that this triangle is big enough
342         }
343       } else {
344         areas[i] = MAXAREA;
345       }
346       map1[i] = i;
347     }
348     
349     for(int i=0; i<noEle; i++) {
350       for(int j=i+1; j<noEle; j++) {
351         if(areas[j] < areas[i]) {
352           double tmp = areas[j];
353           areas[j] = areas[i];
354           areas[i] = tmp;
355           int t = map1[j];
356           map1[j] = map1[i];
357           map1[i] = t;
358         }
359       }
360     }
361     
362     for(int i=0; i<noEle; i++) {
363       if(theMesh->elem[0].is_valid(map1[i])) {
364         if(areas[i] < targetA) {
365           //find the nodes along the smallest edge & coarsen the edge
366           theMesh->e2n_getAll(map1[i],con,0);
367           getShortestEdge(con[0], con[1], con[2], shortestEdge);
368           int ret = theAdaptor->edge_contraction(shortestEdge[0], shortestEdge[1]);
369           if(ret != -1) adapted = true;
370         }
371       }
372       //if(adapted) break;
373     }
374     //if(adapted) break;
375   }
376
377   free(con);
378   free(areas);
379   free(map1);
380   free(n1_coord);
381   free(n2_coord);
382   free(n3_coord);
383   free(shortestEdge);
384
385   return 1;
386 }
387
388 void FEM_Adapt_Algs::tests() {
389   //test the mesh for slivered triangles
390
391   theMod->fmUtil->StructureTest(theMesh);
392   theMod->fmUtil->IdxlListTest(theMesh);
393   return;
394 }
395
396 // =====================  BEGIN refine_element_leb ========================= 
397 /* Given an element e, if e's longest edge f is also the longest edge
398    of e's neighbor across f, g, split f by adding a new node in the 
399    center of f, and splitting both e and g into two elements.  If g
400    does not have f as it's longest edge, recursively call refine_element_leb 
401    on g, and start over. */ 
402 int FEM_Adapt_Algs::refine_element_leb(int e) {
403   int *eConn = (int*)malloc(3*sizeof(int));
404   int fixNode, otherNode, opNode, longEdge, nbr; 
405   double eLens[3], longEdgeLen = 0.0;
406
407   if(e==-1) {
408     free(eConn);
409     return -1;
410   }
411
412   theMesh->e2n_getAll(e, eConn);
413   eLens[0] = length(eConn[0], eConn[1]);
414   eLens[1] = length(eConn[1], eConn[2]);
415   eLens[2] = length(eConn[2], eConn[0]);
416   for (int i=0; i<3; i++) {
417     if (eLens[i] > longEdgeLen) {
418       longEdgeLen = eLens[i];
419       longEdge = i;
420       fixNode = eConn[i];
421       otherNode = eConn[(i+1)%3];
422       opNode = eConn[(i+2)%3];
423     }
424   }
425   free(eConn);
426
427   nbr = theMesh->e2e_getNbr(e, longEdge);
428   if (nbr == -1) // e's longEdge is on physical boundary
429     return theAdaptor->edge_bisect(fixNode, otherNode);
430   int nbrOpNode = theAdaptor->e2n_getNot(nbr, fixNode, otherNode);
431   double fixEdgeLen = length(fixNode, nbrOpNode);
432   double otherEdgeLen = length(otherNode, nbrOpNode);
433   if ((fixEdgeLen > longEdgeLen) || (otherEdgeLen > longEdgeLen)) { 
434     // longEdge is not nbr's longest edge
435     int newNode = theAdaptor->edge_bisect(fixNode, otherNode);
436     if(newNode==-1) return -1;
437     int propElem, propNode; // get the element to propagate on
438     if (fixEdgeLen > otherEdgeLen) {
439       propElem = theAdaptor->findElementWithNodes(newNode, fixNode, nbrOpNode);
440       propNode = fixNode;
441     }
442     else {
443       propElem = theAdaptor->findElementWithNodes(newNode,otherNode,nbrOpNode);
444       propNode = otherNode;
445     }
446
447     //if propElem is ghost, then it's propagating to neighbor, otherwise not
448     if(!FEM_Is_ghost_index(propElem)) {
449       refine_flip_element_leb(propElem,propNode,newNode,nbrOpNode,longEdgeLen);
450     }
451     else {
452       int localChk, nbrChk;
453       localChk = theMod->getfmUtil()->getIdx();
454       nbrChk = theMod->getfmUtil()->getRemoteIdx(theMesh,propElem,0);
455       int propNodeT = theAdaptor->getSharedNodeIdxl(propNode, nbrChk);
456       int newNodeT = theAdaptor->getSharedNodeIdxl(newNode, nbrChk);
457       int nbrghost = (nbrOpNode>=0)?0:1;
458       int nbrOpNodeT = (nbrOpNode>=0)?(theAdaptor->getSharedNodeIdxl(nbrOpNode, nbrChk)):(theAdaptor->getGhostNodeIdxl(nbrOpNode, nbrChk));
459       int propElemT = theAdaptor->getGhostElementIdxl(propElem, nbrChk);
460       meshMod[nbrChk].refine_flip_element_leb(localChk, propElemT, propNodeT, newNodeT,nbrOpNodeT,nbrghost,longEdgeLen);
461     }
462     return newNode;
463   }
464   else return theAdaptor->edge_bisect(fixNode, otherNode); // longEdge on nbr
465 }
466 void FEM_Adapt_Algs::refine_flip_element_leb(int e, int p, int n1, int n2, 
467                                              double le) 
468 {
469   int newNode = refine_element_leb(e);
470   if(newNode == -1) return;
471   (void) theAdaptor->edge_flip(n1, n2);
472   if (length(p, newNode) > le) {
473     int localChk = theMod->getfmUtil()->getIdx();
474     int newElem = theAdaptor->findElementWithNodes(newNode, n1, p);
475     refine_flip_element_leb(newElem, p, n1, newNode, le);
476   }
477 }
478 // ========================  END refine_element_leb ========================
479
480 //HELPER functions
481
482 double FEM_Adapt_Algs::length(int n1, int n2)
483 {
484   double *n1_coord = (double*)malloc(dim*sizeof(double));
485   double *n2_coord = (double*)malloc(dim*sizeof(double));
486
487   getCoord(n1, n1_coord);
488   getCoord(n2, n2_coord);
489
490   double ret = length(n1_coord, n2_coord);
491
492   free(n1_coord);
493   free(n2_coord);
494   return ret;
495 }
496
497 double FEM_Adapt_Algs::length(double *n1_coord, double *n2_coord) { 
498   double d, ds_sum=0.0;
499
500   for (int i=0; i<dim; i++) {
501     d = n1_coord[i] - n2_coord[i];
502     ds_sum += d*d;
503   }
504   return (sqrt(ds_sum));
505 }
506
507
508 double FEM_Adapt_Algs::getArea(int n1, int n2, int n3)
509 {
510   double *n1_coord = (double*)malloc(dim*sizeof(double));
511   double *n2_coord = (double*)malloc(dim*sizeof(double));
512   double *n3_coord = (double*)malloc(dim*sizeof(double));
513
514   getCoord(n1, n1_coord);
515   getCoord(n2, n2_coord);
516   getCoord(n3, n3_coord);
517
518   double ret = getArea(n1_coord, n2_coord, n3_coord);
519
520   free(n1_coord);
521   free(n2_coord);
522   free(n3_coord);
523   return ret;
524 }
525
526 double FEM_Adapt_Algs::getArea(double *n1_coord, double *n2_coord, double *n3_coord) {
527   double area=0.0;
528   double aLen, bLen, cLen, sLen, d, ds_sum;
529
530   ds_sum = 0.0;
531   for (int i=0; i<2; i++) {
532     d = n1_coord[i] - n2_coord[i];
533     ds_sum += d*d;
534   }
535   aLen = sqrt(ds_sum);
536   ds_sum = 0.0;
537   for (int i=0; i<2; i++) {
538     d = n2_coord[i] - n3_coord[i];
539     ds_sum += d*d;
540   }
541   bLen = sqrt(ds_sum);
542   ds_sum = 0.0;
543   for (int i=0; i<2; i++) {
544     d = n3_coord[i] - n1_coord[i];
545     ds_sum += d*d;
546   }
547   cLen = sqrt(ds_sum);
548   sLen = (aLen+bLen+cLen)/2;
549   if(sLen-aLen < 0) return 0.0;
550   else if(sLen-bLen < 0) return 0.0;
551   else if(sLen-cLen < 0) return 0.0; //area too small to note
552   return (sqrt(sLen*(sLen-aLen)*(sLen-bLen)*(sLen-cLen)));
553 }
554
555 bool FEM_Adapt_Algs::didItFlip(int n1, int n2, int n3, double *n4_coord)
556 {
557   //n3 is the node to be deleted, n4 is the new node to be added
558   double *n1_coord = (double*)malloc(dim*sizeof(double));
559   double *n2_coord = (double*)malloc(dim*sizeof(double));
560   double *n3_coord = (double*)malloc(dim*sizeof(double));
561
562   getCoord(n1, n1_coord);
563   getCoord(n2, n2_coord);
564   getCoord(n3, n3_coord);
565
566   double ret_old = getSignedArea(n1_coord, n2_coord, n3_coord);
567   double ret_new = getSignedArea(n1_coord, n2_coord, n4_coord);
568
569   free(n1_coord);
570   free(n2_coord);
571   free(n3_coord);
572
573   if(ret_old > SLIVERAREA && ret_new < -SLIVERAREA) return true; //it is a flip
574   else if(ret_old < -SLIVERAREA && ret_new > SLIVERAREA) return true; //it is a flip
575   else if(fabs(ret_new) < SLIVERAREA) return true; // it is a sliver
576   else return false;
577 }
578
579
580 bool FEM_Adapt_Algs::didItFlip(double *n1_coord, double *n2_coord, double *n3_coord, double *n4_coord)
581 {
582   double ret_old = getSignedArea(n1_coord, n2_coord, n3_coord);
583   double ret_new = getSignedArea(n1_coord, n2_coord, n4_coord);
584   if(ret_old > MINAREA && ret_new < -MINAREA) return true;
585   else if(ret_old < -MINAREA && ret_new > MINAREA) return true;
586   else return false;
587 }
588
589 double FEM_Adapt_Algs::getSignedArea(double *n1_coord, double *n2_coord, double *n3_coord) {
590   double area=0.0;
591   double vec1_x, vec1_y, vec2_x, vec2_y;
592
593   vec1_x = n1_coord[0] - n2_coord[0];
594   vec1_y = n1_coord[1] - n2_coord[1];
595   vec2_x = n3_coord[0] - n2_coord[0];
596   vec2_y = n3_coord[1] - n2_coord[1];
597
598   area = vec1_x*vec2_y - vec2_x*vec1_y;
599   return area;
600 }
601
602 int FEM_Adapt_Algs::getCoord(int n1, double *crds) {
603   if(!FEM_Is_ghost_index(n1)) {
604     FEM_Mesh_dataP(theMesh, FEM_NODE, coord_attr, (void *)crds, n1, 1, FEM_DOUBLE, dim);
605   }
606   else {
607     int numchunks;
608     IDXL_Share **chunks1;
609     theMod->fmUtil->getChunkNos(0,n1,&numchunks,&chunks1);
610     int index = theMod->idx;
611     for(int j=0; j<numchunks; j++) {
612       int chk = chunks1[j]->chk;
613       if(chk==index) continue;
614       int ghostidx = theMod->fmUtil->exists_in_IDXL(theMesh,n1,chk,2);
615       double2Msg *d = meshMod[chk].getRemoteCoord(index,ghostidx);
616       crds[0] = d->i;
617       crds[1] = d->j;
618       for(int j=0; j<numchunks; j++) {
619         delete chunks1[j];
620       }
621       if(numchunks != 0) free(chunks1);
622       break;
623     }
624   }
625   return 1;
626 }
627
628 int FEM_Adapt_Algs::getShortestEdge(int n1, int n2, int n3, int* shortestEdge) {
629   double *n1_coord = (double*)malloc(dim*sizeof(double));
630   double *n2_coord = (double*)malloc(dim*sizeof(double));
631   double *n3_coord = (double*)malloc(dim*sizeof(double));
632
633   getCoord(n1, n1_coord);
634   getCoord(n2, n2_coord);
635   getCoord(n3, n3_coord);
636
637   double aLen = length(n1_coord, n2_coord);
638   int shortest = 0;
639
640   double bLen = length(n2_coord, n3_coord);
641   if(bLen < aLen) shortest = 1;
642
643   double cLen = length(n3_coord, n1_coord);
644   if((cLen < aLen) && (cLen < bLen)) shortest = 2;
645
646   if(shortest==0) {
647     shortestEdge[0] = n1;
648     shortestEdge[1] = n2;
649   }
650   else if (shortest==1) {
651     shortestEdge[0] = n2;
652     shortestEdge[1] = n3;
653   }
654   else {
655     shortestEdge[0] = n3;
656     shortestEdge[1] = n1;
657   }
658   free(n1_coord);
659   free(n2_coord);
660   free(n3_coord);
661   return 1;
662 }
663
664
665 void FEM_Adapt_Algs::Insert(int eIdx, double len, int cFlag)
666 {
667   int i;
668   if (cFlag) {
669     i = ++coarsenHeapSize; 
670     while ((coarsenElements[i/2].len>=len) && (i != 1)) {
671       coarsenElements[i].len=coarsenElements[i/2].len;
672       coarsenElements[i].elID=coarsenElements[i/2].elID;
673       i/=2;                     
674     }
675     coarsenElements[i].elID=eIdx;
676     coarsenElements[i].len=len; 
677   }
678   else {
679     i = ++refineHeapSize; 
680     while ((refineElements[i/2].len>=len) && (i != 1)) {
681       refineElements[i].len=refineElements[i/2].len;
682       refineElements[i].elID=refineElements[i/2].elID;
683       i/=2;                     
684     }
685     refineElements[i].elID=eIdx;
686     refineElements[i].len=len; 
687   }
688 }
689
690 // removes and returns the minimum element of the heap 
691 int FEM_Adapt_Algs::Delete_Min(int cflag)
692 {
693   int Child, i, Min_ID; 
694   if (cflag) {
695     Min_ID=coarsenElements[1].elID;
696     for (i=1; i*2 <= coarsenHeapSize-1; i=Child) { // Find smaller child
697       Child = i*2; // child is left child  
698       if (Child != coarsenHeapSize)  // right child exists
699         if (coarsenElements[Child+1].len < coarsenElements[Child].len)
700           Child++; 
701       // Percolate one level
702       if (coarsenElements[coarsenHeapSize].len >= coarsenElements[Child].len) {
703         coarsenElements[i].elID = coarsenElements[Child].elID;
704         coarsenElements[i].len = coarsenElements[Child].len;
705       }
706       else break; 
707     }
708     coarsenElements[i].elID = coarsenElements[coarsenHeapSize].elID;
709     coarsenElements[i].len = coarsenElements[coarsenHeapSize].len; 
710     coarsenHeapSize--;
711     return Min_ID; 
712   }
713   else {
714     Min_ID=refineElements[1].elID;
715     for (i=1; i*2 <= refineHeapSize-1; i=Child) { // Find smaller child
716       Child = i*2;       // child is left child  
717       if (Child !=refineHeapSize)  // right child exists
718         if (refineElements[Child+1].len < refineElements[Child].len)
719           Child++; 
720       // Percolate one level
721       if (refineElements[refineHeapSize].len >= refineElements[Child].len){  
722         refineElements[i].elID = refineElements[Child].elID;   
723         refineElements[i].len = refineElements[Child].len;
724       }
725       else break; 
726     }
727     refineElements[i].elID = refineElements[refineHeapSize].elID;
728     refineElements[i].len = refineElements[refineHeapSize].len; 
729     refineHeapSize--;
730     return Min_ID; 
731   }
732 }
733
734 double FEM_Adapt_Algs::getAreaQuality(int elem)
735 {
736   double f, q, len[3];
737   int n[3];
738   double currentArea;
739   double *n1_coord = (double*)malloc(dim*sizeof(double));
740   double *n2_coord = (double*)malloc(dim*sizeof(double));
741   double *n3_coord = (double*)malloc(dim*sizeof(double));
742   theMesh->e2n_getAll(elem, n);
743   getCoord(n[0], n1_coord);
744   getCoord(n[1], n2_coord);
745   getCoord(n[2], n3_coord);
746
747   currentArea = getArea(n1_coord, n2_coord, n3_coord);
748
749   len[0] = length(n1_coord, n2_coord);
750   len[1] = length(n2_coord, n3_coord);
751   len[2] = length(n3_coord, n1_coord);
752   f = 4.0*sqrt(3.0); //proportionality constant
753   q = (f*currentArea)/(len[0]*len[0]+len[1]*len[1]+len[2]*len[2]);
754   return q;
755 }
756
757 // FEM_Mesh_mooth
758 //  Inputs  : meshP - a pointer to the FEM_Mesh object to smooth
759 //          : nodes - an array of local node numbers to be smoothed.  Send
760 //                NULL pointer to smooth all nodes.
761 //          : nNodes - the number of nodes to be smoothed.  This must 
762 //                    be the total number of nodes on this chunk if the nodes
763 //                    array is NULL
764 //          : attrNo - the attribute number where the coords are registered
765 //  Shifts nodes around to improve mesh quality.  FEM_BOUNDARY attribute
766 //  and interpolator function must be registered by user to maintain 
767 //  boundary information.
768 void  FEM_Adapt_Algs::FEM_mesh_smooth(FEM_Mesh *meshP, int *nodes, int nNodes, int attrNo) {
769   vector2d newPos, *coords, *ghostCoords;
770   int idx, nNod, nGn, gIdxN, *boundVals, nodesInChunk, mesh;
771   int *adjnodes;
772
773   mesh=FEM_Mesh_default_read();
774   nodesInChunk = FEM_Mesh_get_length(mesh,FEM_NODE);
775   nGn = FEM_Mesh_get_length(mesh, FEM_GHOST + FEM_NODE);
776   
777   boundVals = new int[nodesInChunk];
778   coords = new vector2d[nodesInChunk+nGn];
779
780   FEM_Mesh_data(mesh, FEM_NODE, FEM_BOUNDARY, (int*) boundVals, 0, nodesInChunk, FEM_INT, 1);    
781
782   FEM_Mesh_data(mesh, FEM_NODE, attrNo, (double*)coords, 0, nodesInChunk, FEM_DOUBLE, 2);
783
784   IDXL_Layout_t coord_layout = IDXL_Layout_create(IDXL_DOUBLE, 2);
785   FEM_Update_ghost_field(coord_layout,-1, coords); 
786   ghostCoords = &(coords[nodesInChunk]);
787   for (int i=0; i<nNodes; i++)
788   {
789     if (nodes==NULL) idx=i;
790     else idx=nodes[i];
791     newPos.x=0;
792     newPos.y=0;
793     CkAssert(idx<nodesInChunk);  
794     if (FEM_is_valid(mesh, FEM_NODE, idx) && boundVals[idx]>-1) //node must be internal
795     {
796       meshP->n2n_getAll(idx, &adjnodes, &nNod);
797       for (int j=0; j<nNod; j++) { //for all adjacent nodes, find coords
798         if (adjnodes[j]<-1) {
799           gIdxN = FEM_From_ghost_index(adjnodes[j]);
800           newPos.x += ghostCoords[gIdxN].x;
801           newPos.y += ghostCoords[gIdxN].y;
802         }
803         else {
804           newPos.x += coords[adjnodes[j]].x;
805           newPos.y += coords[adjnodes[j]].y;
806         }     
807       }
808       newPos.x/=nNod;
809       newPos.y/=nNod;
810       FEM_set_entity_coord2(mesh, FEM_NODE, idx, newPos.x, newPos.y);
811       delete [] adjnodes;
812     }
813   }
814
815   delete [] coords;
816   delete [] boundVals;
817 }
818