examples: add diagnostic printout
[charm.git] / examples / ParFUM / 2Dexample / pgm.C
1 #include "pgm.h"
2 #include "mpi.h"
3
4
5 #define SUMMARY_ON
6
7 extern void _registerParFUM(void);
8
9 double getArea(double *n1_coord, double *n2_coord, double *n3_coord) {
10   double area=0.0;
11   double aLen, bLen, cLen, sLen, d, ds_sum;
12
13   ds_sum = 0.0;
14   for (int i=0; i<2; i++) {
15     d = n1_coord[i] - n2_coord[i];
16     ds_sum += d*d;
17   }
18   aLen = sqrt(ds_sum);
19   ds_sum = 0.0;
20   for (int i=0; i<2; i++) {
21     d = n2_coord[i] - n3_coord[i];
22     ds_sum += d*d;
23   }
24   bLen = sqrt(ds_sum);
25   ds_sum = 0.0;
26   for (int i=0; i<2; i++) {
27     d = n3_coord[i] - n1_coord[i];
28     ds_sum += d*d;
29   }
30   cLen = sqrt(ds_sum);
31   sLen = (aLen+bLen+cLen)/2;
32   if(sLen-aLen < 0) return 0.0;
33   else if(sLen-bLen < 0) return 0.0;
34   else if(sLen-cLen < 0) return 0.0; //area too small to note
35   return (sqrt(sLen*(sLen-aLen)*(sLen-bLen)*(sLen-cLen)));
36 }
37
38 void publish_data_netfem(int i,  myGlobals g, MPI_Comm comm) {
39   MPI_Barrier(comm);
40   if (1) { //Publish data to the net
41     int mesh=FEM_Mesh_default_read(); // Tell framework we are reading data from the mesh
42     int rank = FEM_My_partition();
43     g.nnodes=FEM_Mesh_get_length(mesh,FEM_NODE); // Get number of nodes
44     g.coord=new vector2d[g.nnodes];
45     int count = 0;
46     vector2d *coord = new vector2d[g.nnodes];
47     int *maptovalid = new int[g.nnodes];
48     double *nodeid = new double[g.nnodes];
49     // Get node positions
50     FEM_Mesh_data(mesh, FEM_NODE, FEM_DATA+0, (double*)g.coord, 0, g.nnodes, FEM_DOUBLE, 2);
51     for(int j=0; j<g.nnodes; j++) {
52       if(FEM_is_valid(mesh,FEM_NODE+0,j)) {
53         coord[count].x = g.coord[j].x;
54         coord[count].y = g.coord[j].y;
55         maptovalid[j] = count;
56         nodeid[count] = j;
57         count++;
58       }
59     }
60     NetFEM n=NetFEM_Begin(rank,i,2,NetFEM_WRITE);
61     NetFEM_Nodes(n,count,(double *)coord,"Position (m)");
62     //NetFEM_Nodes(n,g.nnodes,(double *)g.coord,"Position (m)");
63     NetFEM_Scalar(n,nodeid,1,"Node ID");
64
65     // Get element data
66     g.nelems=FEM_Mesh_get_length(mesh,FEM_ELEM+0); // Get number of elements
67     g.conn=new connRec[g.nelems];
68     connRec *conn = new connRec[g.nelems];
69     double *elid = new double[g.nelems];
70     count = 0;
71     // Get connectivity for elements
72     FEM_Mesh_data(mesh, FEM_ELEM+0, FEM_CONN, (int *)g.conn, 0, g.nelems, FEM_INDEX_0, 3);
73     double totalArea = 0.0;
74     for(int j=0; j<g.nelems; j++) {
75       if(FEM_is_valid(mesh,FEM_ELEM+0,j)) {
76         conn[count][0] = maptovalid[g.conn[j][0]];
77         conn[count][1] = maptovalid[g.conn[j][1]];
78         conn[count][2] = maptovalid[g.conn[j][2]];
79         elid[count] = j;
80         totalArea += getArea(coord[conn[count][0]],coord[conn[count][1]],coord[conn[count][2]]);
81         if(totalArea != totalArea) {
82           CkPrintf("NAN\n");
83         }
84         count++;
85       }
86     }
87     NetFEM_Elements(n,count,3,(int *)conn,"Triangles");
88     //NetFEM_Elements(n,g.nelems,3,(int *)g.conn,"Triangles");
89     NetFEM_Scalar(n,elid,1,"Element ID");
90     NetFEM_End(n);
91
92     double finalArea;
93     CkPrintf("Chunk[%d]: local area: %.12f\n",rank,totalArea);
94     MPI_Reduce((void*)&totalArea,(void*)&finalArea,1,MPI_DOUBLE,MPI_SUM,0,comm);
95     if(rank == 0) CkPrintf("Chunk[%d]: total area: %.12f\n",rank,finalArea);
96
97     delete [] g.coord;
98     delete [] g.conn;
99     delete [] coord;
100     delete [] conn;
101     delete [] maptovalid;
102     delete [] nodeid;
103     delete [] elid;
104   }
105 }
106
107
108 extern "C" void
109 init(void)
110 {
111   CkPrintf("init started\n");
112   double startTime=CmiWallTimer();
113   const char *eleName="mesh1.tri";//"adpmm/xxx.1.ele";//*/"88mesh/mesh1.tri";
114   const char *nodeName="mesh1.node";//"adpmm/xxx.1.node";//*/"88mesh/mesh1.node";
115   int nPts=0; //Number of nodes
116   vector2d *pts=0; //Node coordinates
117   int *bounds;
118
119   CkPrintf("Reading node coordinates from %s\n",nodeName);
120   //Open and read the node coordinate file
121   {
122     char line[1024];
123     FILE *f=fopen(nodeName,"r");
124     if (f==NULL) die("Can't open node file!");
125     fgets(line,1024,f);
126     if (1!=sscanf(line,"%d",&nPts)) die("Can't read number of points!");
127     pts=new vector2d[nPts];
128     bounds = new int[nPts];
129     for (int i=0;i<nPts;i++) {
130       int ptNo;
131       if (NULL==fgets(line,1024,f)) die("Can't read node input line!");
132       if (4!=sscanf(line,"%d%lf%lf%d",&ptNo,&pts[i].x,&pts[i].y,&bounds[i])) 
133         die("Can't parse node input line!");
134     }
135     fclose(f);
136   }
137  
138   int nEle=0;
139   connRec *ele=NULL;
140   CkPrintf("Reading elements from %s\n",eleName);
141   //Open and read the element connectivity file
142   {
143     char line[1024];
144     FILE *f=fopen(eleName,"r");
145     if (f==NULL) die("Can't open element file!");
146     fgets(line,1024,f);
147     if (1!=sscanf(line,"%d",&nEle)) die("Can't read number of elements!");
148     ele=new connRec[nEle];
149     for (int i=0;i<nEle;i++) {
150       int elNo;
151       if (NULL==fgets(line,1024,f)) die("Can't read element input line!");
152       if (4!=sscanf(line,"%d%d%d%d",&elNo,&ele[i][0],&ele[i][1],&ele[i][2])) 
153         die("Can't parse element input line!");  
154       ele[i][0]--; //Fortran to C indexing
155       ele[i][1]--; //Fortran to C indexing
156       ele[i][2]--; //Fortran to C indexing
157     }
158     fclose(f);
159   }
160
161
162   int fem_mesh=FEM_Mesh_default_write(); // Tell framework we are writing to the mesh
163
164   CkPrintf("Passing node coords to framework\n");
165
166
167   FEM_Mesh_data(fem_mesh,        // Add nodes to the current mesh
168                 FEM_NODE,        // We are registering nodes
169                 FEM_DATA+0,      // Register the point locations which are normally 
170                                  // the first data elements for an FEM_NODE
171                 (double *)pts,   // The array of point locations
172                 0,               // 0 based indexing
173                 nPts,            // The number of points
174                 FEM_DOUBLE,      // Coordinates are doubles
175                 2);              // Points have dimension 2 (x,y)
176  
177  
178   CkPrintf("Passing node bounds to framework\n");
179
180
181   FEM_Mesh_data(fem_mesh,        // Add nodes to the current mesh
182                 FEM_NODE,        // We are registering nodes
183                 FEM_BOUNDARY,      // Register the point bound info 
184                                  // the first data elements for an FEM_NODE
185                 (int *)bounds,   // The array of point bound info
186                 0,               // 0 based indexing
187                 nPts,            // The number of points
188                 FEM_INT,         // bounds are ints
189                 1);              // Points have dimension 1
190
191   CkPrintf("Passing elements to framework\n");
192
193   FEM_Mesh_data(fem_mesh,      // Add nodes to the current mesh
194                 FEM_ELEM+0,      // We are registering elements with type 0
195                                  // The next type of element could be registered with FEM_ELEM+1
196                 FEM_CONN,        // Register the connectivity table for this
197                                  // data elements for this type of FEM entity
198                 (int *)ele,      // The array of point locations
199                 0,               // 0 based indexing
200                 nEle,            // The number of elements
201                 FEM_INDEX_0,     // We use zero based node numbering
202                 3);              // Elements have degree 3, since triangles are defined 
203                                  // by three nodes
204
205
206
207   // Register values to the elements so we can keep track of them after partitioning
208   double *values = (double*)malloc(nEle*sizeof(double));
209   for(int i=0;i<nEle;i++)values[i]=i;
210
211   // triangles
212   FEM_Mesh_data(fem_mesh,      // Add nodes to the current mesh
213                 FEM_ELEM,      // We are registering elements with type 1
214                 FEM_DATA,   
215                 (int *)values,   // The array of point locations
216                 0,               // 0 based indexing
217                 nEle,            // The number of elements
218                 FEM_DOUBLE,         // We use zero based node numbering
219                 1);
220  
221
222   //boundary conditions
223   FEM_Mesh_data(fem_mesh,      // Add nodes to the current mesh
224                 FEM_NODE,      // We are registering elements with type 1
225                 FEM_BOUNDARY,   
226                 (int *)bounds,   // The array of point locations
227                 0,               // 0 based indexing
228                 nPts,            // The number of elements
229                 FEM_INT,         // We use zero based node numbering
230                 1);
231
232
233
234   delete[] ele;
235   delete[] pts;
236   delete[] bounds;
237   free(values);
238   
239   double *sizingData = new double[nEle];
240   for (int i=0; i<nEle; i++) sizingData[i]=-1.0;
241   FEM_Mesh_data(fem_mesh, FEM_ELEM+0, FEM_MESH_SIZING, sizingData, 0, nEle,
242                   FEM_DOUBLE, 1);
243   delete [] sizingData;
244
245   // add ghost layers
246
247   const int triangleFaces[6] = {0,1,2};
248   CkPrintf("Adding Ghost layers\n");
249   FEM_Add_ghost_layer(1,1);
250   FEM_Add_ghost_elem(0,3,triangleFaces);
251
252   CkPrintf("Finished with init (Reading took %.f s)\n",CmiWallTimer()-startTime);
253
254 }
255
256
257 // A driver() function 
258 // driver() is required in all FEM programs
259 extern "C" void
260 driver(void)
261 {
262   int nnodes,nelems,nelems2,ignored;
263   int i, myId=FEM_My_partition();
264   myGlobals g;
265   FEM_Register(&g,(FEM_PupFn)pup_myGlobals);
266   
267   _registerParFUM();
268
269   printf("partition %d is in driver\n", myId);
270
271   int mesh=FEM_Mesh_default_read(); // Tell framework we are reading data from the mesh
272   
273   // Get node data
274   nnodes=FEM_Mesh_get_length(mesh,FEM_NODE); // Get number of nodes
275   g.nnodes=nnodes;
276   g.coord=new vector2d[nnodes];
277   // Get node positions
278   FEM_Mesh_data(mesh, FEM_NODE, FEM_DATA+0, (double*)g.coord, 0, nnodes, FEM_DOUBLE, 2);
279   // Get element data
280   nelems=FEM_Mesh_get_length(mesh,FEM_ELEM+0); // Get number of elements
281   g.nelems=nelems;
282   g.conn=new connRec[nelems];
283   // Get connectivity for elements
284   FEM_Mesh_data(mesh, FEM_ELEM+0, FEM_CONN, (int *)g.conn, 0, nelems, FEM_INDEX_0, 3);
285
286
287   double* tridata =new double[nelems];
288   FEM_Mesh_data(mesh, FEM_ELEM+0, FEM_DATA, (int *)tridata, 0, nelems, FEM_DOUBLE, 1);  
289
290   int nelemsghost=FEM_Mesh_get_length(mesh,FEM_ELEM+0+FEM_GHOST); 
291   double* trighostdata =new double[nelemsghost];
292   FEM_Mesh_data(mesh, FEM_ELEM+0+FEM_GHOST, FEM_DATA, (int *)trighostdata, 0, nelemsghost, FEM_DOUBLE, 1);  
293   int nnodesghost=FEM_Mesh_get_length(mesh,FEM_NODE+0+FEM_GHOST); 
294   double* nodeghostdata =new double[2*nnodesghost];
295   FEM_Mesh_data(mesh, FEM_NODE+0+FEM_GHOST, FEM_DATA, (int *)nodeghostdata, 0, nnodesghost, FEM_DOUBLE, 2);  
296
297
298   {
299         const int triangleFaces[6] = {0,1,1,2,2,0};
300         FEM_Add_elem2face_tuples(mesh, 0, 2, 3, triangleFaces);
301         FEM_Mesh_create_elem_elem_adjacency(mesh);
302         FEM_Mesh_allocate_valid_attr(mesh, FEM_ELEM+0);
303         FEM_Mesh_allocate_valid_attr(mesh, FEM_NODE);
304         FEM_Mesh_create_node_elem_adjacency(mesh);
305         FEM_Mesh_create_node_node_adjacency(mesh);
306
307         int netIndex = 0;
308         int rank = 0;
309         MPI_Comm comm=MPI_COMM_WORLD;
310         //MPI_Group iwgroup,commgroup;
311         //MPI_Comm_group(MPI_COMM_WORLD, &commgroup);
312         //MPI_Comm_create(MPI_COMM_WORLD,commgroup,&comm);
313         MPI_Comm_rank(comm,&rank);
314         
315         publish_data_netfem(netIndex,g,comm); netIndex++;
316 #ifdef SUMMARY_ON
317         FEM_Print_Mesh_Summary(mesh);
318 #endif
319
320         MPI_Barrier(comm);
321         FEM_REF_INIT(mesh,2);
322         
323         FEM_Mesh *meshP = FEM_Mesh_lookup(FEM_Mesh_default_read(),"driver");
324         FEM_AdaptL *ada = meshP->getfmMM()->getfmAdaptL();
325         int ret_op = -1;
326
327         FEM_Adapt_Algs *adaptAlgs= meshP->getfmMM()->getfmAdaptAlgs();
328         adaptAlgs->FEM_Adapt_Algs_Init(FEM_DATA+0,FEM_DATA+4);
329         FEM_Interpolate *interp = meshP->getfmMM()->getfmInp();
330         //interp->FEM_SetInterpolateNodeEdgeFnPtr(interpolate);
331
332         MPI_Barrier(comm);
333
334         //CkPrintf("Shadow arrays have been bound\n");
335         /*
336 #ifdef SUMMARY_ON
337         FEM_Print_Mesh_Summary(mesh);
338 #endif
339                 
340         CkPrintf("Marking 5 nodes and one element as invalid\n");
341         FEM_set_entity_invalid(mesh, FEM_NODE, 5);
342         FEM_set_entity_invalid(mesh, FEM_NODE, 6);
343         FEM_set_entity_invalid(mesh, FEM_NODE, 7);
344         FEM_set_entity_invalid(mesh, FEM_NODE, 8);
345         FEM_set_entity_invalid(mesh, FEM_NODE, 9);      
346         FEM_set_entity_invalid(mesh, FEM_ELEM, 9);
347
348 #ifdef SUMMARY_ON
349         FEM_Print_Mesh_Summary(mesh);
350 #endif
351         
352         CkPrintf("Marking 5 nodes and one element as valid again\n");
353         FEM_set_entity_valid(mesh, FEM_NODE, 5);
354         FEM_set_entity_valid(mesh, FEM_NODE, 6);
355         FEM_set_entity_valid(mesh, FEM_NODE, 7);
356         FEM_set_entity_valid(mesh, FEM_NODE, 8);
357         FEM_set_entity_valid(mesh, FEM_NODE, 9);        
358         FEM_set_entity_valid(mesh, FEM_ELEM, 9);
359
360 #ifdef SUMMARY_ON
361         FEM_Print_Mesh_Summary(mesh);
362 #endif
363
364   g.S11=new double[nelems];
365   g.S22=new double[nelems];
366   g.S12=new double[nelems];
367   // Get connectivity for elements
368   FEM_Mesh_data(mesh, FEM_ELEM+0, FEM_CONN, (int *)g.conn, 0, nelems, FEM_INDEX_0, 3);  
369
370
371   //Initialize associated data
372   g.R_net=new vector2d[nnodes]; //Net force
373   g.d=new vector2d[nnodes];//Node displacement
374   g.v=new vector2d[nnodes];//Node velocity
375   g.a=new vector2d[nnodes];//Node acceleration
376   for (i=0;i<nnodes;i++)
377     g.R_net[i]=g.d[i]=g.v[i]=g.a[i]=vector2d(0.0);
378
379 //Apply a small initial perturbation to positions
380   for (i=0;i<nnodes;i++) {
381           const double max=1.0e-10/15.0; //Tiny perturbation
382           g.d[i].x+=max*(i&15);
383           g.d[i].y+=max*((i+5)&15);
384   }
385
386   int fid=FEM_Create_simple_field(FEM_DOUBLE,2);
387
388   //Timeloop
389   if (myId==0)
390     CkPrintf("Entering timeloop\n");
391   int tSteps=5000;
392   double startTime, totalStart;
393   startTime=totalStart=CkWallTimer();
394   for (int t=0;t<tSteps;t++) {
395     if (1) { //Structural mechanics
396
397         int adjs[3];
398         int elemid;
399         if(rank == 0) {
400           adjs[0] = 15;
401           adjs[1] = 16;
402           adjs[2] = 21; // -5;
403           elemid = 28;
404         } else if(rank == 1) {
405           adjs[0] = 19;
406           adjs[1] = 5;
407           adjs[2] = 7;
408           elemid = 21;
409         } else if(rank == 2) {
410           adjs[0] = 8;
411           adjs[1] = 11;
412           adjs[2] = 6;
413           elemid = 7;
414         } else {
415           adjs[0] = 0;
416           adjs[1] = 1;
417           adjs[2] = 2;
418           elemid = 0;
419         }
420         int newel1 = 0;
421 #ifdef SUMMARY_ON
422         FEM_Print_Mesh_Summary(mesh);
423         //FEM_Print_n2e(mesh,adjs[0]);
424         //FEM_Print_n2e(mesh,adjs[1]);
425         //FEM_Print_n2e(mesh,adjs[2]);
426         //FEM_Print_n2n(mesh,adjs[0]);
427         //FEM_Print_n2n(mesh,adjs[1]);
428         //FEM_Print_n2n(mesh,adjs[2]);
429         //FEM_Print_e2n(mesh,newel1);
430         //FEM_Print_e2e(mesh,newel1);
431 #endif
432
433         //FEM_Modify_Lock(mesh, adjs, 3, adjs, 0);
434         if(rank == 0) {
435           FEM_remove_element(mesh, elemid, 0, 1);
436         }
437         //FEM_Modify_Unlock(mesh);
438         MPI_Barrier(comm);
439 #ifdef SUMMARY_ON
440         FEM_Print_Mesh_Summary(mesh);
441 #endif
442
443         //FEM_Modify_Lock(mesh, adjs, 3, adjs, 0);
444         if(rank == 0) {
445           newel1 = FEM_add_element(mesh,adjs,3,0,0);
446           CkPrintf("New Element\n");
447         }
448         //FEM_Modify_Unlock(mesh);
449         publish_data_netfem(netIndex,g,comm); netIndex++;
450 #ifdef SUMMARY_ON
451         FEM_Print_Mesh_Summary(mesh);
452         //FEM_Print_n2e(mesh,adjs[0]);
453         //FEM_Print_n2e(mesh,adjs[1]);
454         //FEM_Print_n2e(mesh,adjs[2]);
455         //FEM_Print_n2n(mesh,adjs[0]);
456         //FEM_Print_n2n(mesh,adjs[1]);
457         //FEM_Print_n2n(mesh,adjs[2]);
458         //FEM_Print_e2n(mesh,newel1);
459         //FEM_Print_e2e(mesh,newel1);
460 #endif
461
462         if(rank==0){
463           FEM_Print_Mesh_Summary(mesh);
464           CkPrintf("%d: Removing element \n", rank);
465           
466           int nelemsghost   =FEM_Mesh_get_length(mesh,FEM_ELEM+0+FEM_GHOST); 
467           int numvalidghost =FEM_count_valid(mesh,FEM_ELEM+0+FEM_GHOST);
468           CkPrintf("nelemsghost=%d numvalidghost=%d\n", nelemsghost, numvalidghost);
469         
470           for(int i=1;i<20;i++){
471                 if(FEM_is_valid(mesh, FEM_ELEM+FEM_GHOST, i)){
472                   double data[1];
473                   FEM_Mesh_data(mesh, FEM_ELEM+FEM_GHOST, FEM_DATA, (int *)data, i, 1, FEM_DOUBLE, 1);  
474
475                   CkPrintf("%d: Eating ghost element %d with value %f\n", rank, i, data[1]);
476                   int conn[3];
477                   
478                   FEM_Mesh_data(mesh, FEM_ELEM+FEM_GHOST, FEM_CONN, (int *)conn, i, 1, FEM_INDEX_0, 3);
479                   CkPrintf("conn for element is: %d %d %d\n", conn[0], conn[1], conn[2]);
480                   FEM_Modify_Lock(mesh, conn, 3, conn, 0);
481                   FEM_remove_element(mesh, FEM_From_ghost_index(i), 0, 1);
482                   FEM_Modify_Unlock(mesh);
483
484                   MPI_Barrier(comm);
485                   FEM_Print_Mesh_Summary(mesh);
486
487                   FEM_Modify_Lock(mesh, conn, 3, conn, 0);
488                   FEM_add_element(mesh, conn, 3, 0, rank); // add locally
489                   FEM_Modify_Unlock(mesh);
490                   CkPrintf("New conn for element is: %d %d %d\n", conn[0], conn[1], conn[2]);
491                   
492                   publish_data_netfem(netIndex,g,comm); netIndex++;
493                   FEM_Print_Mesh_Summary(mesh);
494                 }
495                 else{
496                   //  CkPrintf("invalid element %d\n", i);
497                 }
498           }
499         }
500         else {
501           CkPrintf("Rank %d\n", rank);
502           for(int i=1;i<20;i++){
503             MPI_Barrier(comm);
504             FEM_Print_Mesh_Summary(mesh);
505
506             publish_data_netfem(netIndex,g,comm); netIndex++;
507             FEM_Print_Mesh_Summary(mesh);
508           }
509         }
510         
511         publish_data_netfem(netIndex,g,comm); netIndex++;
512         */      
513         /*
514         CkPrintf("Starting Local edge flips on individual chunks\n");
515         int flip[4];
516         if(rank == 0) {
517           flip[0] = 20;
518           flip[1] = 21;
519           flip[2] = 0;
520           flip[3] = 3;
521         }
522         else if(rank == 1) {
523           flip[0] = 9;
524           flip[1] = 10;
525           flip[2] = 0;
526           flip[3] = 4;
527         }
528         else if(rank == 2) {
529           flip[0] = 1;
530           flip[1] = 2;
531           flip[2] = 6;
532           flip[3] = -5;
533         }
534         else {
535           flip[0] = 0;
536           flip[1] = 1;
537           flip[2] = 2;
538           flip[3] = 3;
539         }
540 #ifdef SUMMARY_ON
541         FEM_Print_Mesh_Summary(mesh);
542 #endif
543         ret_op = ada->edge_flip(flip[0],flip[1]);
544         publish_data_netfem(netIndex,g,comm); netIndex++;
545         adaptAlgs->tests();
546         MPI_Barrier(comm);
547         //ret_op = ada->edge_flip(flip[2],flip[3]);
548         //publish_data_netfem(netIndex,g,comm); netIndex++;
549 #ifdef SUMMARY_ON
550         FEM_Print_Mesh_Summary(mesh);
551 #endif
552
553         CkPrintf("Starting shared edge flips on individual chunks\n");
554         int sflip[4];
555         if(rank == 0) {
556           sflip[0] = 19;
557           sflip[1] = 18;
558           sflip[2] = 1;
559           sflip[3] = -4;
560         }
561         else if(rank == 1) {
562           sflip[0] = 5;
563           sflip[1] = 6;
564           sflip[2] = 7;
565           sflip[3] = -5;
566         }
567         else if(rank == 2) {
568           sflip[0] = 11;
569           sflip[1] = 2;
570           sflip[2] = 0;
571           sflip[3] = -2;
572         }
573         else {
574           sflip[0] = 0;
575           sflip[1] = 1;
576           sflip[2] = 2;
577           sflip[3] = 3;
578         }
579 #ifdef SUMMARY_ON
580         FEM_Print_Mesh_Summary(mesh);
581 #endif
582         ret_op = ada->edge_flip(sflip[0],sflip[1]);
583         publish_data_netfem(netIndex,g,comm); netIndex++;
584         if(ret_op > 0) {
585           if(sflip[2]<0) sflip[2] = ret_op;
586           else if(sflip[3]<0) sflip[3] = ret_op;
587         }
588         ret_op = ada->edge_flip(sflip[2],sflip[3]);
589         publish_data_netfem(netIndex,g,comm); netIndex++;
590 #ifdef SUMMARY_ON
591         FEM_Print_Mesh_Summary(mesh);
592 #endif
593         
594         CkPrintf("Starting Local edge bisect on individual chunks\n");
595         int bisect[2];
596         if(rank == 0) {
597           bisect[0] = 16;
598           bisect[1] = 21;
599         }
600         else if(rank == 1) {
601           bisect[0] = 5;
602           bisect[1] = 6;
603         }
604         else if(rank == 2) {
605           bisect[0] = 8;
606           bisect[1] = 11;
607         }
608         else {
609           bisect[0] = 0;
610           bisect[1] = 1;
611         }
612         if(rank==2) ret_op = ada->edge_bisect(bisect[0],bisect[1]);
613         publish_data_netfem(netIndex,g,comm); netIndex++;
614         adaptAlgs->tests();
615         MPI_Barrier(comm);
616         
617 #ifdef SUMMARY_ON
618         FEM_Print_Mesh_Summary(mesh);
619 #endif
620         
621         CkPrintf("Starting Local vertex remove on individual chunks\n");
622         int vr[2];
623         if(rank == 0) {
624           vr[0] = ret_op;
625           vr[1] = 6;
626         }
627         else if(rank == 1) {
628           vr[0] = ret_op;
629           vr[1] = 13;
630         }
631         else if(rank == 2) {
632           vr[0] = ret_op;
633           vr[1] = 21;
634         }
635         else {
636           vr[0] = ret_op;
637           vr[1] = 1;
638         }
639 #ifdef SUMMARY_ON
640         FEM_Print_Mesh_Summary(mesh);
641 #endif
642         ret_op = ada->vertex_remove(vr[0],vr[1]);
643         publish_data_netfem(netIndex,g,comm); netIndex++;
644 #ifdef SUMMARY_ON
645         FEM_Print_Mesh_Summary(mesh);
646 #endif
647         
648         
649         CkPrintf("Starting shared edge bisect on individual chunks\n");
650         int sbisect[2];
651         if(rank == 0) {
652           sbisect[0] = 1;
653           sbisect[1] = 19;
654         }
655         else if(rank == 1) {
656           sbisect[0] = 0;
657           sbisect[1] = 21;
658         }
659         else if(rank == 2) {
660           sbisect[0] = 1;
661           sbisect[1] = 9;
662         }
663         else {
664           sbisect[0] = 0;
665           sbisect[1] = 1;
666         }
667         
668 #ifdef SUMMARY_ON
669         FEM_Print_Mesh_Summary(mesh);
670 #endif
671         ret_op = ada->edge_bisect(sbisect[0],sbisect[1]);
672         publish_data_netfem(netIndex,g,comm); netIndex++;
673         adaptAlgs->tests();
674         MPI_Barrier(comm);
675 #ifdef SUMMARY_ON
676         FEM_Print_Mesh_Summary(mesh);
677 #endif
678
679         
680         CkPrintf("Starting shared vertex remove on individual chunks\n");
681
682         int svr[2];
683         if(rank == 0) {
684           svr[0] = ret_op;
685           svr[1] = 19;
686         }
687         else if(rank == 1) {
688           svr[0] = ret_op;
689           svr[1] = 21;
690         }
691         else if(rank == 2) {
692           svr[0] = ret_op;
693           svr[1] = 20;
694         }
695         else {
696           svr[0] = ret_op;
697           svr[1] = 1;
698         }
699 #ifdef SUMMARY_ON
700         FEM_Print_Mesh_Summary(mesh);
701 #endif
702         ret_op = ada->vertex_remove(svr[0],svr[1]);
703         publish_data_netfem(netIndex,g,comm); netIndex++;
704 #ifdef SUMMARY_ON
705         FEM_Print_Mesh_Summary(mesh);
706 #endif
707
708         CkPrintf("Starting Local edge contract on individual chunks\n");
709         int contract[2];
710         if(rank == 0) {
711           contract[0] = 28;
712           contract[1] = 30;
713         }
714         else if(rank == 1) {
715           contract[0] = 10;
716           contract[1] = 9;
717         }
718         else if(rank == 2) {
719           contract[0] = 1;
720           contract[1] = 2;
721         }
722         else {
723           contract[0] = 0;
724           contract[1] = 1;
725         }
726         
727 #ifdef SUMMARY_ON
728         FEM_Print_Mesh_Summary(mesh);
729 #endif
730         ret_op = ada->edge_contraction(contract[0],contract[1]);
731         publish_data_netfem(netIndex,g,comm); netIndex++;
732         adaptAlgs->tests();
733         MPI_Barrier(comm);
734
735         //if(rank==2) adaptAlgs->simple_coarsen(0.00004);
736         if(rank==0) ret_op = ada->edge_contraction(27,28);
737         publish_data_netfem(netIndex,g,comm); netIndex++;
738         adaptAlgs->tests();
739         MPI_Barrier(comm);
740
741         CkPrintf("Starting Local edge bisect on individual chunks\n");
742         int bisect[2];
743         if(rank == 0) {
744           bisect[0] = 16;
745           bisect[1] = 21;
746         }
747         else if(rank == 1) {
748           bisect[0] = 10;
749           bisect[1] = 9;
750         }
751         else if(rank == 2) {
752           bisect[0] = 5;
753           bisect[1] = 3;
754         }
755         else {
756           bisect[0] = 0;
757           bisect[1] = 1;
758         }
759         ret_op = ada->edge_bisect(bisect[0],bisect[1]);
760         publish_data_netfem(netIndex,g,comm); netIndex++;
761         adaptAlgs->tests();
762         MPI_Barrier(comm);
763         
764 #ifdef SUMMARY_ON
765         FEM_Print_Mesh_Summary(mesh);
766 #endif
767
768         //CkPrintf("Starting Local vertex split on individual chunks\n");
769         /*
770         int vs[3];
771         if(rank == 0) {
772           vs[0] = ret_op;
773           vs[1] = 9;
774           vs[2] = 13;
775         }
776         else if(rank == 1) {
777           vs[0] = ret_op;
778           vs[1] = 8;
779           vs[2] = 7;
780         }
781         else if(rank == 2) {
782           vs[0] = ret_op;
783           vs[1] = 14;
784           vs[2] = 23;
785         }
786         else {
787           vs[0] = ret_op;
788           vs[1] = 2;
789           vs[2] = 3;
790         }
791 #ifdef SUMMARY_ON
792         //FEM_Print_Mesh_Summary(mesh);
793 #endif
794         //ret_op = ada->vertex_split(vs[0],vs[1],vs[2]);
795         //publish_data_netfem(netIndex,g,comm); netIndex++;
796 #ifdef SUMMARY_ON
797         FEM_Print_Mesh_Summary(mesh);
798 #endif
799         */
800         /*
801         CkPrintf("Starting shared edge contract on individual chunks\n");
802         int scontract[2];
803         if(rank == 0) {
804           scontract[0] = 9;
805           scontract[1] = 10;
806         }
807         else if(rank == 1) {
808           scontract[0] = 5;
809           scontract[1] = 6;
810         }
811         else if(rank == 2) {
812           scontract[0] = 11;
813           scontract[1] = 2;
814         }
815         else {
816           scontract[0] = 0;
817           scontract[1] = 1;
818         }
819         
820 #ifdef SUMMARY_ON
821         FEM_Print_Mesh_Summary(mesh);
822 #endif
823         ret_op = ada->edge_contraction(scontract[0],scontract[1]);
824         publish_data_netfem(netIndex,g,comm); netIndex++;
825         /*
826         //CkPrintf("Starting shared vertex split on individual chunks\n");
827         int svs[3];
828         if(rank == 0) {
829           svs[0] = ret_op;
830           svs[1] = 1;
831           svs[2] = -6;
832         }
833         else if(rank == 1) {
834           svs[0] = ret_op;
835           svs[1] = 7;
836           svs[2] = 7;
837         }
838         else if(rank == 2) {
839           svs[0] = ret_op;
840           svs[1] = 0;
841           svs[2] = -2;
842         }
843         else {
844           svs[0] = ret_op;
845           svs[1] = 2;
846           svs[2] = 3;
847         }
848 #ifdef SUMMARY_ON
849         //FEM_Print_Mesh_Summary(mesh);
850 #endif
851         //ret_op = ada->vertex_split(svs[0],svs[1],svs[2]);
852         //publish_data_netfem(netIndex,g,comm); netIndex++;
853 #ifdef SUMMARY_ON
854         FEM_Print_Mesh_Summary(mesh);
855 #endif
856         */
857
858         /*
859         CkPrintf("Starting LEB on individual chunks\n");
860         int *leb_elem = (int*)malloc(1*sizeof(int));
861         if(rank==0) {
862           leb_elem[0] = 2;
863         }
864         else if(rank==1) {
865           leb_elem[0] = 13; //4;
866         }
867         else if(rank==2) {
868           leb_elem[0] = 20; //26;
869         }
870         else if (rank == 3){
871           leb_elem[0] = 14;
872         }
873         else {
874           leb_elem[0] = 0;
875         }
876
877         adaptAlgs->refine_element_leb(leb_elem[0]);
878         publish_data_netfem(netIndex,g,comm); netIndex++;
879         */
880         /*
881           int nEle;
882           //for(int tstep = 0; tstep < 2; tstep++) {
883           nEle = FEM_Mesh_get_length(mesh, FEM_ELEM);   
884           for (int i=0; i<nEle; i++)
885           if (FEM_is_valid(mesh, FEM_ELEM, i))
886           adaptAlgs->refine_element_leb(i);
887           publish_data_netfem(netIndex,g,comm); netIndex++;
888           FEM_Print_Mesh_Summary(mesh);
889           //}
890           */
891
892       
893       double targetArea = 0.0001;
894       
895       for(int tstep = 0; tstep < 0; tstep++) {
896         int ret = -1;
897         //for(int tstep1=0; tstep1<60; tstep1++) {
898         while(ret==-1) {
899           ret = adaptAlgs->simple_refine(targetArea);
900           publish_data_netfem(netIndex,g,comm); netIndex++;
901           adaptAlgs->tests();
902           MPI_Barrier(comm);
903         }
904         //int *nodes = new int[g.nnodes];
905         //for (int i=0; i<g.nnodes; i++) nodes[i]=i;    
906         //FEM_mesh_smooth(mesh, nodes, g.nnodes, FEM_DATA+0);
907         //publish_data_netfem(netIndex,g,comm); netIndex++;
908         //delete [] nodes;
909
910 #ifdef SUMMARY_ON
911         FEM_Print_Mesh_Summary(mesh);
912 #endif
913         ret = -1;
914         //for(int tstep1=0; tstep1<60; tstep1++) {
915         while(ret==-1) {
916           ret = adaptAlgs->simple_coarsen(targetArea);
917           publish_data_netfem(netIndex,g,comm); netIndex++;
918           adaptAlgs->tests();
919           MPI_Barrier(comm);
920         }
921         //int *nodes = new int[g.nnodes];
922         //for (int i=0; i<g.nnodes; i++) nodes[i]=i;
923         //FEM_mesh_smooth(mesh, nodes, g.nnodes, FEM_DATA+0);
924         //publish_data_netfem(netIndex,g,comm); netIndex++;
925         //delete [] nodes;
926
927 #ifdef SUMMARY_ON
928         FEM_Print_Mesh_Summary(mesh);
929 #endif
930       }
931
932       targetArea *= 0.5;
933       for(int tstep = 0; tstep < 0; tstep++) {
934         int ret = -1;
935         //for(int tstep1=0; tstep1<60; tstep1++) {
936         while(ret==-1) {
937           ret = adaptAlgs->simple_refine(targetArea);
938           publish_data_netfem(netIndex,g,comm); netIndex++;
939           MPI_Barrier(comm);
940           adaptAlgs->tests();
941           MPI_Barrier(comm);
942         }
943         if(rank==0) CkPrintf("[%d] Iteration No. %d",rank,tstep);
944         //int *nodes = new int[g.nnodes];
945         //for (int i=0; i<g.nnodes; i++) nodes[i]=i;    
946         //FEM_mesh_smooth(mesh, nodes, g.nnodes, FEM_DATA+0);
947         //publish_data_netfem(netIndex,g,comm); netIndex++;
948         //delete [] nodes;
949
950
951         targetArea *= 0.5;
952 #ifdef SUMMARY_ON
953         FEM_Print_Mesh_Summary(mesh);
954 #endif
955       }
956       targetArea *= 2.0;
957       
958       for(int tstep = 0; tstep < 2; tstep++) {
959         int ret = -1;
960         //for(int tstep1=0; tstep1<60;tstep1++) {
961         while(ret==-1) {
962           ret = adaptAlgs->simple_coarsen(targetArea);
963           //MPI_Barrier(comm);
964           publish_data_netfem(netIndex,g,comm); netIndex++;
965           adaptAlgs->tests();
966           MPI_Barrier(comm);
967         }
968         //int *nodes = new int[g.nnodes];
969         //for (int i=0; i<g.nnodes; i++) nodes[i]=i;
970         //FEM_mesh_smooth(mesh, nodes, g.nnodes, FEM_DATA+0);
971         //publish_data_netfem(netIndex,g,comm); netIndex++;
972         //delete [] nodes;
973
974         targetArea *= 2.0;
975 #ifdef SUMMARY_ON
976         FEM_Print_Mesh_Summary(mesh);
977 #endif
978       }
979
980       //wave propagation on a bar
981       targetArea = 0.00004;
982       double xmin = 0.00;
983       double xmax = 0.1;
984       double ymin = 0.00;
985       double ymax = 0.01;
986       for(int tstep = 0; tstep < 0; tstep++) {
987         targetArea = 0.000002;
988         adaptAlgs->simple_refine(targetArea, xmin, ymin, xmax, ymax);
989         publish_data_netfem(netIndex,g,comm); netIndex++;
990 #ifdef SUMMARY_ON
991         FEM_Print_Mesh_Summary(mesh);
992 #endif
993         targetArea = 0.0000014;
994         adaptAlgs->simple_coarsen(targetArea, xmin, ymin, xmax, ymax);
995         publish_data_netfem(netIndex,g,comm); netIndex++;
996 #ifdef SUMMARY_ON
997         FEM_Print_Mesh_Summary(mesh);
998 #endif
999         ymin += 0.01;
1000         ymax += 0.01;
1001       }
1002
1003       //crack propagation on a block
1004       targetArea = 0.00004;
1005       xmin = 0.00;
1006       xmax = 0.2;
1007       double xcrackmin = 0.09;
1008       double xcrackmax = 0.10;
1009       ymin = 0.00;
1010       ymax = 0.02;
1011       for(int tstep = 0; tstep < 0; tstep++) {
1012         targetArea = 0.000025;
1013         adaptAlgs->simple_refine(targetArea, xmin, ymin, xmax, ymax);
1014         publish_data_netfem(netIndex,g,comm); netIndex++;
1015 #ifdef SUMMARY_ON
1016         FEM_Print_Mesh_Summary(mesh);
1017 #endif
1018         targetArea = 0.00005;
1019         adaptAlgs->simple_coarsen(targetArea, xmin, ymin, xmax, ymax);
1020         //publish_data_netfem(netIndex,g,comm); netIndex++;
1021 #ifdef SUMMARY_ON
1022         FEM_Print_Mesh_Summary(mesh);
1023 #endif
1024         /*if(tstep > 2) {
1025           targetArea = 0.000025;
1026           adaptAlgs->simple_refine(targetArea, xcrackmin, ymin, xcrackmax, ymax);
1027           //publish_data_netfem(netIndex,g,comm); netIndex++;
1028 #ifdef SUMMARY_ON
1029           FEM_Print_Mesh_Summary(mesh);
1030 #endif
1031           xcrackmin -= 0.004;
1032           xcrackmax += 0.004;
1033         }
1034         */
1035
1036         ymin += 0.02;
1037         ymax += 0.02;
1038       }
1039
1040       bool adapted = true;
1041       for(int j=0; j<0; j++) {
1042         int i=0;
1043         adapted = false;
1044         if(rank==2) {
1045           for(i=0; i<meshP->elem[0].ghost->size(); i++) {
1046             if(meshP->elem[0].ghost->is_valid(i)) {
1047               adapted = true;
1048               break;
1049             }
1050           }
1051         }
1052         if(adapted) {
1053           //lock the nodes
1054           int conn[3];
1055           bool done = false;
1056           int *gotlocks = (int*)malloc(3*sizeof(int));
1057           bool bailout = false;
1058           while(!done) {
1059             if(meshP->elem[0].ghost->is_valid(i)) {
1060               meshP->e2n_getAll(FEM_To_ghost_index(i),conn,0);
1061               int gotlock = ada->lockNodes(gotlocks, conn, 0, conn, 3);
1062               if(gotlock==1) done = true;
1063             }
1064             else {
1065               bailout = true;
1066               break;
1067             }
1068           }
1069           if(!bailout) {
1070             int newEl = meshP->getfmMM()->fmUtil->eatIntoElement(FEM_To_ghost_index(i));
1071             meshP->e2n_getAll(newEl,conn,0);
1072             //FEM_Modify_correctLockN(meshP, conn[0]);
1073             //FEM_Modify_correctLockN(meshP, conn[1]);
1074             //FEM_Modify_correctLockN(meshP, conn[2]);
1075             ada->unlockNodes(gotlocks, conn, 0, conn, 3);
1076             free(gotlocks);
1077           }
1078         }
1079         publish_data_netfem(netIndex,g,comm); netIndex++;
1080         adaptAlgs->tests();
1081         MPI_Barrier(comm);
1082       }
1083
1084       /*
1085     //Debugging/perf. output
1086     double curTime=CkWallTimer();
1087     double total=curTime-startTime;
1088     startTime=curTime;
1089     if (myId==0 && (t%64==0)) {
1090             CkPrintf("%d %.6f sec for loop %d \n",CkNumPes(),total,t);
1091             if (0) {
1092               CkPrintf("    Triangle 0:\n");
1093               for (int j=0;j<3;j++) {
1094                     int n=g.conn[0][j];
1095                     CkPrintf("    Node %d: coord=(%.4f,%.4f)  d=(%.4g,%.4g)\n",
1096                              n,g.coord[n].x,g.coord[n].y,g.d[n].x,g.d[n].y);
1097               }
1098             }
1099     }
1100     // perform migration-based load balancing 
1101     if (t%1024==0)
1102       FEM_Migrate();
1103     
1104     if (t%1024==0) { //Publish data to the net
1105             NetFEM n=NetFEM_Begin(FEM_My_partition(),t,2,NetFEM_POINTAT);
1106             
1107             NetFEM_Nodes(n,nnodes,(double *)g.coord,"Position (m)");
1108             NetFEM_Vector(n,(double *)g.d,"Displacement (m)");
1109             NetFEM_Vector(n,(double *)g.v,"Velocity (m/s)");
1110             
1111             NetFEM_Elements(n,nelems,3,(int *)g.conn,"Triangles");
1112                 NetFEM_Scalar(n,g.S11,1,"X Stress (pure)");
1113                 NetFEM_Scalar(n,g.S22,1,"Y Stress (pure)");
1114                 NetFEM_Scalar(n,g.S12,1,"Shear Stress (pure)");
1115             
1116             NetFEM_End(n);
1117     }
1118   }
1119 */
1120
1121       CkPrintf("chunk %d Waiting for Synchronization\n",rank);
1122       MPI_Barrier(comm);
1123       CkPrintf("Synchronized\n");
1124 #ifdef SUMMARY_ON
1125       FEM_Print_Mesh_Summary(mesh);
1126 #endif
1127       publish_data_netfem(netIndex,g,comm); netIndex++;
1128       
1129       CkExit();
1130   }
1131   
1132 }
1133
1134
1135 // A PUP function to allow for migration and load balancing of mesh partitions.
1136 // The PUP function is not needed if no migration or load balancing is desired.
1137 void pup_myGlobals(pup_er p,myGlobals *g) 
1138 {
1139   FEM_Print("-------- called pup routine -------");
1140   pup_int(p,&g->nnodes);
1141   pup_int(p,&g->nelems);
1142   int nnodes=g->nnodes, nelems=g->nelems;
1143   if (pup_isUnpacking(p)) {
1144     g->coord=new vector2d[nnodes];
1145     g->conn=new connRec[nelems];
1146     //g->R_net=new vector2d[nnodes]; //Net force
1147     //g->d=new vector2d[nnodes];//Node displacement
1148     //g->v=new vector2d[nnodes];//Node velocity
1149     //g->a=new vector2d[nnodes];
1150     //g->S11=new double[nelems];
1151     //g->S22=new double[nelems];
1152     //g->S12=new double[nelems];
1153   }
1154   pup_doubles(p,(double *)g->coord,2*nnodes);
1155   pup_ints(p,(int *)g->conn,3*nelems);
1156   //pup_doubles(p,(double *)g->R_net,2*nnodes);
1157   //pup_doubles(p,(double *)g->d,2*nnodes);
1158   //pup_doubles(p,(double *)g->v,2*nnodes);
1159   //pup_doubles(p,(double *)g->a,2*nnodes);
1160   //pup_doubles(p,(double *)g->S11,nelems);
1161   //pup_doubles(p,(double *)g->S22,nelems);
1162   //pup_doubles(p,(double *)g->S12,nelems);
1163   if (pup_isDeleting(p)) {
1164     delete[] g->coord;
1165     delete[] g->conn;
1166     //delete[] g->R_net;
1167     //delete[] g->d;
1168     //delete[] g->v;
1169     //delete[] g->a;
1170     //delete[] g->S11;
1171     //delete[] g->S22;
1172     //delete[] g->S12;
1173   }
1174 }
1175
1176
1177 void FEM_mesh_smooth(int mesh, int *nodes, int nNodes, int attrNo)
1178 {
1179   vector2d *centroids, newPos, *coords, *ghostCoords, *vGcoords;
1180   int nEle, nGn, *boundVals, nodesInChunk, nVg;
1181   int neighbors[3], *adjelems;
1182   int gIdxN;
1183   int j=0;
1184   double x[3], y[3];
1185   FEM_Mesh *meshP = FEM_Mesh_lookup(mesh, "driver");
1186
1187   nodesInChunk = FEM_Mesh_get_length(mesh,FEM_NODE);
1188   boundVals = new int[nodesInChunk];
1189   nGn = FEM_Mesh_get_length(mesh, FEM_GHOST + FEM_NODE);
1190   coords = new vector2d[nodesInChunk+nGn];
1191
1192   FEM_Mesh_data(mesh, FEM_NODE, FEM_BOUNDARY, (int*) boundVals, 0, nodesInChunk, FEM_INT, 1);    
1193
1194   FEM_Mesh_data(mesh, FEM_NODE, attrNo, (double*)coords, 0, nodesInChunk, FEM_DOUBLE, 2);
1195   for (int i=0; i<(nodesInChunk); i++) {
1196     //CkPrintf(" coords[%d]: (%f, %f)\n", i, coords[i].x, coords[i].y);
1197   }
1198   IDXL_Layout_t coord_layout = IDXL_Layout_create(IDXL_DOUBLE, 2);
1199   FEM_Update_ghost_field(coord_layout,-1, coords); 
1200   ghostCoords = &(coords[nodesInChunk]);
1201   /*
1202   for (int i=0; i<nGn;i++) {
1203     if (FEM_is_valid(mesh, FEM_GHOST+FEM_NODE, i)) {
1204       CkPrintf("ghost %d is valid \n", i);        
1205       // vGcoords[j]=ghostCoords[i];
1206       //j++;
1207     }
1208     else
1209       CkPrintf("ghost %d is invalid \n", i);
1210   }
1211   */
1212   for (int i=0; i<(nodesInChunk+nGn); i++) {
1213     //CkPrintf(" coords[%d]: (%f, %f)\n", i, coords[i].x, coords[i].y);
1214   }
1215 //  FEM_Mesh_data(FEM_Mesh_default_write(), FEM_GHOST+FEM_NODE, attrNo, (double*)ghostCoords, 0, nGn, FEM_DOUBLE, 2);
1216  
1217   for (int i=0; i<nNodes; i++)
1218   {
1219     newPos.x=0;
1220     newPos.y=0;
1221     CkAssert(nodes[i]<nodesInChunk);    
1222     if (FEM_is_valid(mesh, FEM_NODE, i) && boundVals[i]>-1) //node must be internal
1223     {
1224       meshP->n2e_getAll(i, adjelems, nEle);
1225       centroids = new vector2d[nEle];
1226       
1227       for (int j=0; j<nEle; j++) { //for all adjacent elements, find centroids
1228         meshP->e2n_getAll(adjelems[j], neighbors);
1229         for (int k=0; k<3; k++) {
1230           if (neighbors[k]<-1) {
1231             gIdxN = FEM_From_ghost_index(neighbors[k]);
1232             x[k] = ghostCoords[gIdxN].x;
1233             y[k] = ghostCoords[gIdxN].y;
1234           }
1235           else {
1236             x[k] = coords[neighbors[k]].x;
1237             y[k] = coords[neighbors[k]].y;
1238           }
1239         }     
1240         centroids[j].x=(x[0]+x[1]+x[2])/3.0;
1241         centroids[j].y=(y[0]+y[1]+y[2])/3.0;
1242         newPos.x += centroids[j].x;
1243         newPos.y += centroids[j].y;
1244       }
1245       newPos.x/=nEle;
1246       newPos.y/=nEle;
1247       FEM_set_entity_coord2(mesh, FEM_NODE, nodes[i], newPos.x, newPos.y);
1248       delete [] centroids;
1249       delete [] adjelems;
1250     }
1251   }
1252   delete [] coords;
1253   delete [] boundVals;
1254 }
1255
1256 void interpolate(FEM_Interpolate::NodalArgs args, FEM_Mesh *meshP)
1257 {
1258   //CkPrintf("INTERPOLATOR!!!!!!!!!!!\n");
1259   int length = meshP->node.realsize();
1260   int *boundVals= new int[length];
1261
1262   FEM_Mesh_dataP(meshP, FEM_NODE, FEM_BOUNDARY, (int*) boundVals, 0, length , FEM_INT, 1);   
1263   CkVec<FEM_Attribute *>*attrs = (meshP->node).getAttrVec();
1264   for (int i=0; i<attrs->size(); i++) {
1265     FEM_Attribute *a = (FEM_Attribute *)(*attrs)[i];
1266     if (a->getAttr() < FEM_ATTRIB_TAG_MAX || a->getAttr()==FEM_BOUNDARY) {
1267       if (a->getAttr()==FEM_BOUNDARY) {
1268         
1269         int n1_bound =boundVals[args.nodes[0]]; 
1270         int n2_bound =boundVals[args.nodes[1]]; 
1271         if (n1_bound == n2_bound && n1_bound < 0) {
1272           a->copyEntity(args.n, *a, args.nodes[0]);
1273         } else if (n1_bound != n2_bound && n1_bound<0 && n2_bound < 0){
1274           a->copyEntity(args.n, *a, args.nodes[0]); //a node which is not on the boundary
1275         } else  if(n1_bound < 0){
1276           a->copyEntity(args.n, *a, args.nodes[1]); 
1277         } else {
1278           a->copyEntity(args.n, *a, args.nodes[0]); 
1279         }
1280       }
1281       else {
1282         FEM_DataAttribute *d = (FEM_DataAttribute *)a;
1283         d->interpolate(args.nodes[0], args.nodes[1], args.n, args.frac);
1284       }
1285     }
1286   }
1287   delete boundVals;
1288 }
1289
1290
1291