First cut of GradateMesh, moved SetReferenceMesh into fem_adapt_algs
[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) { // mesh sizing has been set independently; use as is
250     return;
251   }
252 }
253
254
255 void FEM_Adapt_Algs::SetReferenceMesh()
256 {
257     // for each element, set its size to its average edge length
258     // TODO: do we need to run this loop for element types other than 0?
259     
260     double avgLength;
261     int width = theMesh->elem[0].getConn().width();
262     int* eConn = (int*)malloc(width*sizeof(int));
263     int numElements = theMesh->elem[0].size();
264     
265     for (int i=0; i<numElements; ++i, avgLength=0) {
266         theMesh->e2n_getAll(i, eConn);
267         
268         for (int j=0; j<width-1; ++j) {
269             avgLength += length(eConn[j], eConn[j+1]);
270         }
271         avgLength += length(eConn[0], eConn[width-1]);
272         avgLength /= width;
273         theMesh->elem[0].setMeshSizing(i, avgLength);      
274     }
275     free(eConn);
276 }
277
278
279 void FEM_Adapt_Algs::GradateMesh(double smoothness)
280 {
281     // Resize mesh elements to avoid jumps in element size
282     // Algorithm based on h-shock correction, described in
283     // Mesh Gradation Control, Borouchaki et al
284     // IJNME43 1998 www.ann.jussieu.fr/~frey/publications/ijnme4398.pdf
285
286     const double beta = smoothness;
287
288     double maxShock = 0;
289     double minShock = 1e10;
290     int iteration = 0;
291
292     int* adjNodes, *boundNodes;
293     int nadjNodes, nnodes;
294     int meshNum = FEM_Mesh_default_read();
295
296     nnodes = theMesh->node.size();
297     boundNodes = new int[nnodes];
298     FEM_Mesh_data(meshNum, FEM_NODE, FEM_BOUNDARY, 
299             boundNodes, 0, nnodes, FEM_INT, 1);
300
301
302     printf("Running h-shock mesh gradation with beta=%.3f\n", beta);
303     while (maxShock > beta) {
304         
305         for (int node=0; node<nnodes; ++node) {
306             // only do internal nodes (?)
307             if (boundNodes[node] < 0 || !FEM_is_valid(meshNum, FEM_NODE, node))
308                 continue;
309             
310             theMesh->n2n_getAll(node, &adjNodes, &nadjNodes);
311             for (int adjNode=0; adjNode<nadjNodes; ++adjNode) {
312                 // get length of edge
313                 // (do we need to worry about ghost/non-ghost?)
314                 double edgelen = length(node, adjNode);
315                 
316                 // get adjacent elemnents and their sizes
317                 int e1, e2;
318                 theMesh->get2ElementsOnEdge(node, adjNode, &e1, &e2);
319                 double s1, s2;
320                 s1 = theMesh->elem[0].getMeshSizing(e1);
321                 s2 = theMesh->elem[0].getMeshSizing(e2);
322                 
323                 // h-shock=max(size ratio)^(1/edge length)
324                 double ratio = (s1 > s2) ? s1/s2 : s2/s1;
325                 double hs = pow(ratio, 1.0/edgelen);
326                 if (hs > maxShock) maxShock = hs;
327                 if (hs < minShock) minShock = hs;
328                 
329                 // if hs > beta, resize the larger elt:
330                 // new size = old size / eta^2
331                 // eta = (beta / h-shock)^(edge length)
332                 //     = beta^(edge length) / size ratio
333                 if (hs > beta) {
334                     double etasq = pow(pow(beta, edgelen) / ratio, 2.0);
335                     if (s1 > s2) {
336                         theMesh->elem[0].setMeshSizing(e1, etasq);
337                     } else {
338                         theMesh->elem[0].setMeshSizing(e2, etasq);
339                     }
340                 }
341             }
342         }
343         printf("Finished iteration %d\n", iteration++);
344         printf("Max shock:%8.3f\n", maxShock);
345         printf("Min shock:%8.3f\n", minShock);
346         printf("Target:%8.3f\n", beta);
347     }
348 }
349
350
351 int FEM_Adapt_Algs::simple_refine(double targetA, double xmin, double ymin, double xmax, double ymax) {
352   int noEle = theMesh->elem[0].size();
353   int *con = (int*)malloc(3*sizeof(int));
354   double *areas = (double*)malloc(noEle*sizeof(double));
355   int *map1 = (int*)malloc(noEle*sizeof(int));
356   double *n1_coord = (double*)malloc(2*sizeof(double));
357   double *n2_coord = (double*)malloc(2*sizeof(double));
358   double *n3_coord = (double*)malloc(2*sizeof(double));
359
360   for(int i=0; i<noEle; i++) {
361     if(theMesh->elem[0].is_valid(i)) {
362       theMesh->e2n_getAll(i,con,0);
363       getCoord(con[0], n1_coord);
364       getCoord(con[1], n2_coord);
365       getCoord(con[2], n3_coord);
366       //do a refinement only if it has any node within the refinement box
367       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)) {
368         areas[i] = getArea(n1_coord, n2_coord, n3_coord);
369         } 
370       else {
371         areas[i] = MINAREA; //make it believe that this triangle does not need refinement
372       }
373     } else {
374       areas[i] = MINAREA;
375     }
376     map1[i] = i;
377   }
378
379   for(int i=0; i<noEle; i++) {
380     for(int j=i+1; j<noEle; j++) {
381       if(areas[j] > areas[i]) {
382         double tmp = areas[j];
383         areas[j] = areas[i];
384         areas[i] = tmp;
385         int t = map1[j];
386         map1[j] = map1[i];
387         map1[i] = t;
388       }
389     }
390   }
391
392   for(int i=0; i<noEle; i++) {
393     if(theMesh->elem[0].is_valid(map1[i])) {
394       if(areas[i] > targetA) {
395         refine_element_leb(map1[i]);
396       }
397     }
398   }
399   free(con);
400   free(areas);
401   free(map1);
402   free(n1_coord);
403   free(n2_coord);
404   free(n3_coord);
405
406   return 1;
407 }
408
409 int FEM_Adapt_Algs::simple_coarsen(double targetA, double xmin, double ymin, double xmax, double ymax) {
410   int noEle = theMesh->elem[0].size();
411   int *con = (int*)malloc(3*sizeof(int));
412   double *areas = (double*)malloc(noEle*sizeof(double));
413   int *map1 = (int*)malloc(noEle*sizeof(int));
414   double *n1_coord = (double*)malloc(2*sizeof(double));
415   double *n2_coord = (double*)malloc(2*sizeof(double));
416   double *n3_coord = (double*)malloc(2*sizeof(double));
417   int *shortestEdge = (int *)malloc(2*sizeof(int));
418   bool adapted = true;
419
420   while(adapted) {
421     adapted = false;
422     for(int i=0; i<noEle; i++) {
423       if(theMesh->elem[0].is_valid(i)) {
424         theMesh->e2n_getAll(i,con,0);
425         getCoord(con[0], n1_coord);
426         getCoord(con[1], n2_coord);
427         getCoord(con[2], n3_coord);
428         //do a coarsening only if it has any node within the coarsen box
429         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)) {
430           areas[i] = getArea(n1_coord, n2_coord, n3_coord);
431         } 
432         else {
433           areas[i] = MAXAREA; //make it believe that this triangle is big enough
434         }
435       } else {
436         areas[i] = MAXAREA;
437       }
438       map1[i] = i;
439     }
440     
441     for(int i=0; i<noEle; i++) {
442       for(int j=i+1; j<noEle; j++) {
443         if(areas[j] < areas[i]) {
444           double tmp = areas[j];
445           areas[j] = areas[i];
446           areas[i] = tmp;
447           int t = map1[j];
448           map1[j] = map1[i];
449           map1[i] = t;
450         }
451       }
452     }
453     
454     for(int i=0; i<noEle; i++) {
455       if(theMesh->elem[0].is_valid(map1[i])) {
456         if(areas[i] < targetA) {
457           //find the nodes along the smallest edge & coarsen the edge
458           theMesh->e2n_getAll(map1[i],con,0);
459           getShortestEdge(con[0], con[1], con[2], shortestEdge);
460           int ret = theAdaptor->edge_contraction(shortestEdge[0], shortestEdge[1]);
461           if(ret != -1) adapted = true;
462         }
463       }
464       //if(adapted) break;
465     }
466     //if(adapted) break;
467   }
468
469   free(con);
470   free(areas);
471   free(map1);
472   free(n1_coord);
473   free(n2_coord);
474   free(n3_coord);
475   free(shortestEdge);
476
477   return 1;
478 }
479
480 void FEM_Adapt_Algs::tests() {
481   //test the mesh for slivered triangles
482
483   theMod->fmUtil->StructureTest(theMesh);
484   theMod->fmUtil->IdxlListTest(theMesh);
485   return;
486 }
487
488 // =====================  BEGIN refine_element_leb ========================= 
489 /* Given an element e, if e's longest edge f is also the longest edge
490    of e's neighbor across f, g, split f by adding a new node in the 
491    center of f, and splitting both e and g into two elements.  If g
492    does not have f as it's longest edge, recursively call refine_element_leb 
493    on g, and start over. */ 
494 int FEM_Adapt_Algs::refine_element_leb(int e) {
495   int *eConn = (int*)malloc(3*sizeof(int));
496   int fixNode, otherNode, opNode, longEdge, nbr; 
497   double eLens[3], longEdgeLen = 0.0;
498
499   if(e==-1) {
500     free(eConn);
501     return -1;
502   }
503
504   theMesh->e2n_getAll(e, eConn);
505   eLens[0] = length(eConn[0], eConn[1]);
506   eLens[1] = length(eConn[1], eConn[2]);
507   eLens[2] = length(eConn[2], eConn[0]);
508   for (int i=0; i<3; i++) {
509     if (eLens[i] > longEdgeLen) {
510       longEdgeLen = eLens[i];
511       longEdge = i;
512       fixNode = eConn[i];
513       otherNode = eConn[(i+1)%3];
514       opNode = eConn[(i+2)%3];
515     }
516   }
517   free(eConn);
518
519   nbr = theMesh->e2e_getNbr(e, longEdge);
520   if (nbr == -1) // e's longEdge is on physical boundary
521     return theAdaptor->edge_bisect(fixNode, otherNode);
522   int nbrOpNode = theAdaptor->e2n_getNot(nbr, fixNode, otherNode);
523   double fixEdgeLen = length(fixNode, nbrOpNode);
524   double otherEdgeLen = length(otherNode, nbrOpNode);
525   if ((fixEdgeLen > longEdgeLen) || (otherEdgeLen > longEdgeLen)) { 
526     // longEdge is not nbr's longest edge
527     int newNode = theAdaptor->edge_bisect(fixNode, otherNode);
528     if(newNode==-1) return -1;
529     int propElem, propNode; // get the element to propagate on
530     if (fixEdgeLen > otherEdgeLen) {
531       propElem = theAdaptor->findElementWithNodes(newNode, fixNode, nbrOpNode);
532       propNode = fixNode;
533     }
534     else {
535       propElem = theAdaptor->findElementWithNodes(newNode,otherNode,nbrOpNode);
536       propNode = otherNode;
537     }
538
539     //if propElem is ghost, then it's propagating to neighbor, otherwise not
540     if(!FEM_Is_ghost_index(propElem)) {
541       refine_flip_element_leb(propElem,propNode,newNode,nbrOpNode,longEdgeLen);
542     }
543     else {
544       int localChk, nbrChk;
545       localChk = theMod->getfmUtil()->getIdx();
546       nbrChk = theMod->getfmUtil()->getRemoteIdx(theMesh,propElem,0);
547       int propNodeT = theAdaptor->getSharedNodeIdxl(propNode, nbrChk);
548       int newNodeT = theAdaptor->getSharedNodeIdxl(newNode, nbrChk);
549       int nbrghost = (nbrOpNode>=0)?0:1;
550       int nbrOpNodeT = (nbrOpNode>=0)?(theAdaptor->getSharedNodeIdxl(nbrOpNode, nbrChk)):(theAdaptor->getGhostNodeIdxl(nbrOpNode, nbrChk));
551       int propElemT = theAdaptor->getGhostElementIdxl(propElem, nbrChk);
552       meshMod[nbrChk].refine_flip_element_leb(localChk, propElemT, propNodeT, newNodeT,nbrOpNodeT,nbrghost,longEdgeLen);
553     }
554     return newNode;
555   }
556   else return theAdaptor->edge_bisect(fixNode, otherNode); // longEdge on nbr
557 }
558 void FEM_Adapt_Algs::refine_flip_element_leb(int e, int p, int n1, int n2, 
559                                              double le) 
560 {
561   int newNode = refine_element_leb(e);
562   if(newNode == -1) return;
563   (void) theAdaptor->edge_flip(n1, n2);
564   if (length(p, newNode) > le) {
565     int localChk = theMod->getfmUtil()->getIdx();
566     int newElem = theAdaptor->findElementWithNodes(newNode, n1, p);
567     refine_flip_element_leb(newElem, p, n1, newNode, le);
568   }
569 }
570 // ========================  END refine_element_leb ========================
571
572 //HELPER functions
573
574 double FEM_Adapt_Algs::length(int n1, int n2)
575 {
576   double *n1_coord = (double*)malloc(dim*sizeof(double));
577   double *n2_coord = (double*)malloc(dim*sizeof(double));
578
579   getCoord(n1, n1_coord);
580   getCoord(n2, n2_coord);
581
582   double ret = length(n1_coord, n2_coord);
583
584   free(n1_coord);
585   free(n2_coord);
586   return ret;
587 }
588
589 double FEM_Adapt_Algs::length(double *n1_coord, double *n2_coord) { 
590   double d, ds_sum=0.0;
591
592   for (int i=0; i<dim; i++) {
593     d = n1_coord[i] - n2_coord[i];
594     ds_sum += d*d;
595   }
596   return (sqrt(ds_sum));
597 }
598
599
600 double FEM_Adapt_Algs::getArea(int n1, int n2, int n3)
601 {
602   double *n1_coord = (double*)malloc(dim*sizeof(double));
603   double *n2_coord = (double*)malloc(dim*sizeof(double));
604   double *n3_coord = (double*)malloc(dim*sizeof(double));
605
606   getCoord(n1, n1_coord);
607   getCoord(n2, n2_coord);
608   getCoord(n3, n3_coord);
609
610   double ret = getArea(n1_coord, n2_coord, n3_coord);
611
612   free(n1_coord);
613   free(n2_coord);
614   free(n3_coord);
615   return ret;
616 }
617
618 double FEM_Adapt_Algs::getArea(double *n1_coord, double *n2_coord, double *n3_coord) {
619   double area=0.0;
620   double aLen, bLen, cLen, sLen, d, ds_sum;
621
622   ds_sum = 0.0;
623   for (int i=0; i<2; i++) {
624     d = n1_coord[i] - n2_coord[i];
625     ds_sum += d*d;
626   }
627   aLen = sqrt(ds_sum);
628   ds_sum = 0.0;
629   for (int i=0; i<2; i++) {
630     d = n2_coord[i] - n3_coord[i];
631     ds_sum += d*d;
632   }
633   bLen = sqrt(ds_sum);
634   ds_sum = 0.0;
635   for (int i=0; i<2; i++) {
636     d = n3_coord[i] - n1_coord[i];
637     ds_sum += d*d;
638   }
639   cLen = sqrt(ds_sum);
640   sLen = (aLen+bLen+cLen)/2;
641   if(sLen-aLen < 0) return 0.0;
642   else if(sLen-bLen < 0) return 0.0;
643   else if(sLen-cLen < 0) return 0.0; //area too small to note
644   return (sqrt(sLen*(sLen-aLen)*(sLen-bLen)*(sLen-cLen)));
645 }
646
647 bool FEM_Adapt_Algs::didItFlip(int n1, int n2, int n3, double *n4_coord)
648 {
649   //n3 is the node to be deleted, n4 is the new node to be added
650   double *n1_coord = (double*)malloc(dim*sizeof(double));
651   double *n2_coord = (double*)malloc(dim*sizeof(double));
652   double *n3_coord = (double*)malloc(dim*sizeof(double));
653
654   getCoord(n1, n1_coord);
655   getCoord(n2, n2_coord);
656   getCoord(n3, n3_coord);
657
658   double ret_old = getSignedArea(n1_coord, n2_coord, n3_coord);
659   double ret_new = getSignedArea(n1_coord, n2_coord, n4_coord);
660
661   free(n1_coord);
662   free(n2_coord);
663   free(n3_coord);
664
665   if(ret_old > SLIVERAREA && ret_new < -SLIVERAREA) return true; //it is a flip
666   else if(ret_old < -SLIVERAREA && ret_new > SLIVERAREA) return true; //it is a flip
667   else if(fabs(ret_new) < SLIVERAREA) return true; // it is a sliver
668   else return false;
669 }
670
671
672 bool FEM_Adapt_Algs::didItFlip(double *n1_coord, double *n2_coord, double *n3_coord, double *n4_coord)
673 {
674   double ret_old = getSignedArea(n1_coord, n2_coord, n3_coord);
675   double ret_new = getSignedArea(n1_coord, n2_coord, n4_coord);
676   if(ret_old > MINAREA && ret_new < -MINAREA) return true;
677   else if(ret_old < -MINAREA && ret_new > MINAREA) return true;
678   else return false;
679 }
680
681 double FEM_Adapt_Algs::getSignedArea(double *n1_coord, double *n2_coord, double *n3_coord) {
682   double area=0.0;
683   double vec1_x, vec1_y, vec2_x, vec2_y;
684
685   vec1_x = n1_coord[0] - n2_coord[0];
686   vec1_y = n1_coord[1] - n2_coord[1];
687   vec2_x = n3_coord[0] - n2_coord[0];
688   vec2_y = n3_coord[1] - n2_coord[1];
689
690   area = vec1_x*vec2_y - vec2_x*vec1_y;
691   return area;
692 }
693
694 int FEM_Adapt_Algs::getCoord(int n1, double *crds) {
695   if(!FEM_Is_ghost_index(n1)) {
696     FEM_Mesh_dataP(theMesh, FEM_NODE, coord_attr, (void *)crds, n1, 1, FEM_DOUBLE, dim);
697   }
698   else {
699     int numchunks;
700     IDXL_Share **chunks1;
701     theMod->fmUtil->getChunkNos(0,n1,&numchunks,&chunks1);
702     int index = theMod->idx;
703     for(int j=0; j<numchunks; j++) {
704       int chk = chunks1[j]->chk;
705       if(chk==index) continue;
706       int ghostidx = theMod->fmUtil->exists_in_IDXL(theMesh,n1,chk,2);
707       double2Msg *d = meshMod[chk].getRemoteCoord(index,ghostidx);
708       crds[0] = d->i;
709       crds[1] = d->j;
710       for(int j=0; j<numchunks; j++) {
711         delete chunks1[j];
712       }
713       if(numchunks != 0) free(chunks1);
714       break;
715     }
716   }
717   return 1;
718 }
719
720 int FEM_Adapt_Algs::getShortestEdge(int n1, int n2, int n3, int* shortestEdge) {
721   double *n1_coord = (double*)malloc(dim*sizeof(double));
722   double *n2_coord = (double*)malloc(dim*sizeof(double));
723   double *n3_coord = (double*)malloc(dim*sizeof(double));
724
725   getCoord(n1, n1_coord);
726   getCoord(n2, n2_coord);
727   getCoord(n3, n3_coord);
728
729   double aLen = length(n1_coord, n2_coord);
730   int shortest = 0;
731
732   double bLen = length(n2_coord, n3_coord);
733   if(bLen < aLen) shortest = 1;
734
735   double cLen = length(n3_coord, n1_coord);
736   if((cLen < aLen) && (cLen < bLen)) shortest = 2;
737
738   if(shortest==0) {
739     shortestEdge[0] = n1;
740     shortestEdge[1] = n2;
741   }
742   else if (shortest==1) {
743     shortestEdge[0] = n2;
744     shortestEdge[1] = n3;
745   }
746   else {
747     shortestEdge[0] = n3;
748     shortestEdge[1] = n1;
749   }
750   free(n1_coord);
751   free(n2_coord);
752   free(n3_coord);
753   return 1;
754 }
755
756
757 void FEM_Adapt_Algs::Insert(int eIdx, double len, int cFlag)
758 {
759   int i;
760   if (cFlag) {
761     i = ++coarsenHeapSize; 
762     while ((coarsenElements[i/2].len>=len) && (i != 1)) {
763       coarsenElements[i].len=coarsenElements[i/2].len;
764       coarsenElements[i].elID=coarsenElements[i/2].elID;
765       i/=2;                     
766     }
767     coarsenElements[i].elID=eIdx;
768     coarsenElements[i].len=len; 
769   }
770   else {
771     i = ++refineHeapSize; 
772     while ((refineElements[i/2].len>=len) && (i != 1)) {
773       refineElements[i].len=refineElements[i/2].len;
774       refineElements[i].elID=refineElements[i/2].elID;
775       i/=2;                     
776     }
777     refineElements[i].elID=eIdx;
778     refineElements[i].len=len; 
779   }
780 }
781
782 // removes and returns the minimum element of the heap 
783 int FEM_Adapt_Algs::Delete_Min(int cflag)
784 {
785   int Child, i, Min_ID; 
786   if (cflag) {
787     Min_ID=coarsenElements[1].elID;
788     for (i=1; i*2 <= coarsenHeapSize-1; i=Child) { // Find smaller child
789       Child = i*2; // child is left child  
790       if (Child != coarsenHeapSize)  // right child exists
791         if (coarsenElements[Child+1].len < coarsenElements[Child].len)
792           Child++; 
793       // Percolate one level
794       if (coarsenElements[coarsenHeapSize].len >= coarsenElements[Child].len) {
795         coarsenElements[i].elID = coarsenElements[Child].elID;
796         coarsenElements[i].len = coarsenElements[Child].len;
797       }
798       else break; 
799     }
800     coarsenElements[i].elID = coarsenElements[coarsenHeapSize].elID;
801     coarsenElements[i].len = coarsenElements[coarsenHeapSize].len; 
802     coarsenHeapSize--;
803     return Min_ID; 
804   }
805   else {
806     Min_ID=refineElements[1].elID;
807     for (i=1; i*2 <= refineHeapSize-1; i=Child) { // Find smaller child
808       Child = i*2;       // child is left child  
809       if (Child !=refineHeapSize)  // right child exists
810         if (refineElements[Child+1].len < refineElements[Child].len)
811           Child++; 
812       // Percolate one level
813       if (refineElements[refineHeapSize].len >= refineElements[Child].len){  
814         refineElements[i].elID = refineElements[Child].elID;   
815         refineElements[i].len = refineElements[Child].len;
816       }
817       else break; 
818     }
819     refineElements[i].elID = refineElements[refineHeapSize].elID;
820     refineElements[i].len = refineElements[refineHeapSize].len; 
821     refineHeapSize--;
822     return Min_ID; 
823   }
824 }
825
826 double FEM_Adapt_Algs::getAreaQuality(int elem)
827 {
828   double f, q, len[3];
829   int n[3];
830   double currentArea;
831   double *n1_coord = (double*)malloc(dim*sizeof(double));
832   double *n2_coord = (double*)malloc(dim*sizeof(double));
833   double *n3_coord = (double*)malloc(dim*sizeof(double));
834   theMesh->e2n_getAll(elem, n);
835   getCoord(n[0], n1_coord);
836   getCoord(n[1], n2_coord);
837   getCoord(n[2], n3_coord);
838
839   currentArea = getArea(n1_coord, n2_coord, n3_coord);
840
841   len[0] = length(n1_coord, n2_coord);
842   len[1] = length(n2_coord, n3_coord);
843   len[2] = length(n3_coord, n1_coord);
844   f = 4.0*sqrt(3.0); //proportionality constant
845   q = (f*currentArea)/(len[0]*len[0]+len[1]*len[1]+len[2]*len[2]);
846   return q;
847 }
848
849 // FEM_Mesh_mooth
850 //  Inputs  : meshP - a pointer to the FEM_Mesh object to smooth
851 //          : nodes - an array of local node numbers to be smoothed.  Send
852 //                NULL pointer to smooth all nodes.
853 //          : nNodes - the number of nodes to be smoothed.  This must 
854 //                    be the total number of nodes on this chunk if the nodes
855 //                    array is NULL
856 //          : attrNo - the attribute number where the coords are registered
857 //  Shifts nodes around to improve mesh quality.  FEM_BOUNDARY attribute
858 //  and interpolator function must be registered by user to maintain 
859 //  boundary information.
860 void  FEM_Adapt_Algs::FEM_mesh_smooth(FEM_Mesh *meshP, int *nodes, int nNodes, int attrNo) {
861   vector2d newPos, *coords, *ghostCoords;
862   int idx, nNod, nGn, gIdxN, *boundVals, nodesInChunk, mesh;
863   int *adjnodes;
864
865   mesh=FEM_Mesh_default_read();
866   nodesInChunk = FEM_Mesh_get_length(mesh,FEM_NODE);
867   nGn = FEM_Mesh_get_length(mesh, FEM_GHOST + FEM_NODE);
868   
869   boundVals = new int[nodesInChunk];
870   coords = new vector2d[nodesInChunk+nGn];
871
872   FEM_Mesh_data(mesh, FEM_NODE, FEM_BOUNDARY, (int*) boundVals, 0, nodesInChunk, FEM_INT, 1);    
873
874   FEM_Mesh_data(mesh, FEM_NODE, attrNo, (double*)coords, 0, nodesInChunk, FEM_DOUBLE, 2);
875
876   IDXL_Layout_t coord_layout = IDXL_Layout_create(IDXL_DOUBLE, 2);
877   FEM_Update_ghost_field(coord_layout,-1, coords); 
878   ghostCoords = &(coords[nodesInChunk]);
879   for (int i=0; i<nNodes; i++)
880   {
881     if (nodes==NULL) idx=i;
882     else idx=nodes[i];
883     newPos.x=0;
884     newPos.y=0;
885     CkAssert(idx<nodesInChunk);  
886     if (FEM_is_valid(mesh, FEM_NODE, idx) && boundVals[idx]>-1) //node must be internal
887     {
888       meshP->n2n_getAll(idx, &adjnodes, &nNod);
889       for (int j=0; j<nNod; j++) { //for all adjacent nodes, find coords
890         if (adjnodes[j]<-1) {
891           gIdxN = FEM_From_ghost_index(adjnodes[j]);
892           newPos.x += ghostCoords[gIdxN].x;
893           newPos.y += ghostCoords[gIdxN].y;
894         }
895         else {
896           newPos.x += coords[adjnodes[j]].x;
897           newPos.y += coords[adjnodes[j]].y;
898         }     
899       }
900       newPos.x/=nNod;
901       newPos.y/=nNod;
902       FEM_set_entity_coord2(mesh, FEM_NODE, idx, newPos.x, newPos.y);
903       delete [] adjnodes;
904     }
905   }
906
907   delete [] coords;
908   delete [] boundVals;
909 }
910