src/libs/ck-libs/TMRC2D/pgm-coarsen.C

   1 /*
   2 Triangle Mesh Refinement (TMR) demo code.
   3
   4 Reads in a Triangle mesh.
   5 Passes it to the REFINE2D subsystem.
   6 Asks for refinements occasionally.
   7 */
   8 #include <stdlib.h>
   9 #include <unistd.h>
  10 #include <stdio.h>
  11 #include <math.h>
  12 #include "charm++.h"
  13 #include "fem.h"
  14 /*My HEADERS*/
  15 #include "ckvector3d.h"
  16 #include "charm-api.h"
  17 #include "fem_mesh.h"
  18 #include "netfem.h"
  19 #include "refine.h"
  20 #include "femrefine.h"
  21 #include "pgm.h"
  22
  23 //The material constants c, as computed by fortran mat_const
  24 // I think the units here are Pascals (N/m^2)
  25 const double matConst[4]={3.692e9,  1.292e9,  3.692e9,  1.200e9 };
  26
  27 //Material density, Kg/m^3
  28 const double density=5.0*1000.0;
  29 //Plate thickness, meters
  30 const double thickness=0.0001;
  31
  32 //The timestep, in seconds
  33 // This aught to be adjusted when the mesh changes!
  34 const double dt=1.0e-12;
  35
  36 static void die(const char *str) {
  37   CkError("Fatal error: %s\n",str);
  38   CkExit();
  39 }
  40
  41 /**schak*/
  42 void resize_nodes(void *data,int *len,int *max);
  43 void resize_elems(void *data,int *len,int *max);
  44
  45 #define NANCHECK 1 /*Check for NaNs at each timestep*/
  46
  47 extern "C" void
  48 init(void)
  49 {
  50   CkPrintf("init started\n");
  51
  52   const char *eleName="out.1024.ele";
  53   const char *nodeName="out.1024.node";
  54   int nPts=0; //Number of nodes
  55   vector2d *pts=0; //Node coordinates
  56
  57   CkPrintf("Reading node coordinates from %s\n",nodeName);
  58   //Open and read the node coordinate file
  59   {
  60     char line[1024];
  61     FILE *f=fopen(nodeName,"r");
  62     if (f==NULL) die("Can't open node file!");
  63     fgets(line,1024,f);
  64     if (1!=sscanf(line,"%d",&nPts)) die("Can't read number of points!");
  65     pts=new vector2d[nPts];
  66     for (int i=0;i<nPts;i++) {
  67       int ptNo;
  68       if (NULL==fgets(line,1024,f)) die("Can't read node input line!");
  69       if (3!=sscanf(line,"%d%lf%lf",&ptNo,&pts[i].x,&pts[i].y))
  70         die("Can't parse node input line!");
  71     }
  72     fclose(f);
  73   }
  74   CkPrintf("Passing node coords to framework\n");
  75
  76   // Register the node entity and its data arrays that will be used later. This
  77   // needs to be done so that the width of the data segments are set correctly
  78   // and can be used later
  79   FEM_Register_entity(FEM_Mesh_default_write(),FEM_NODE,NULL,nPts,nPts,resize_nodes);
  80   for(int k=0;k<=4;k++){
  81     if(k != 0){
  82       vector2d *t = new vector2d[nPts];
  83       FEM_Register_array(FEM_Mesh_default_write(),FEM_NODE,FEM_DATA+k,t,FEM_DOUBLE,2);
  84     }else{
  85       FEM_Register_array(FEM_Mesh_default_write(),FEM_NODE,FEM_DATA+k,pts,FEM_DOUBLE,2);
  86     }
  87   }
  88   double *td = new double[nPts];
  89   FEM_Register_array(FEM_Mesh_default_write(),FEM_NODE,FEM_DATA+5,td,FEM_DOUBLE,1);
  90   int *validNodes = new int[nPts];
  91   for(int ii=0;ii<nPts;ii++){
  92     validNodes[ii]=1;
  93   }
  94   FEM_Register_array(FEM_Mesh_default_write(),FEM_NODE,FEM_VALID,validNodes,FEM_INT,1);
  95
  96   int nEle=0;
  97   int *ele=NULL;
  98   CkPrintf("Reading elements from %s\n",eleName);
  99   //Open and read the element connectivity file
 100   {
 101     char line[1024];
 102     FILE *f=fopen(eleName,"r");
 103     if (f==NULL) die("Can't open element file!");
 104     fgets(line,1024,f);
 105     if (1!=sscanf(line,"%d",&nEle)) die("Can't read number of elements!");
 106     ele=new int[3*nEle];
 107     for (int i=0;i<nEle;i++) {
 108       int elNo;
 109       if (NULL==fgets(line,1024,f)) die("Can't read element input line!");
 110       if (4!=sscanf(line,"%d%d%d%d",&elNo,&ele[3*i+0],&ele[3*i+1],&ele[3*i+2]))
 111         die("Can't parse element input line!");
 112       ele[3*i+0]--; //Fortran to C indexing
 113       ele[3*i+1]--; //Fortran to C indexing
 114       ele[3*i+2]--; //Fortran to C indexing
 115     }
 116     fclose(f);
 117   }
 118
 119   CkPrintf("Passing elements to framework\n");
 120   // Register the Element entity and its connectivity array. Register the
 121   // data arrays to set up the widths correctly at the beginning
 122   FEM_Register_entity(FEM_Mesh_default_write(),FEM_ELEM,NULL,nEle,nEle,resize_elems);
 123   FEM_Register_array(FEM_Mesh_default_write(),FEM_ELEM,FEM_CONN,ele,FEM_INDEX_0,3);
 124
 125   for(int k=0;k<3;k++){
 126     void *t = new double[nEle];
 127     FEM_Register_array(FEM_Mesh_default_write(),FEM_ELEM,FEM_DATA+k,t,FEM_DOUBLE,1);
 128   }
 129   int *validElem = new int[nEle];
 130   for(int ii=0;ii<nEle;ii++){
 131     validElem[ii]=1;
 132   }
 133   FEM_Register_array(FEM_Mesh_default_write(),FEM_ELEM,FEM_VALID,validElem,FEM_INT,1);
 134
 135   /*Build the ghost layer for refinement border*/
 136   FEM_Add_ghost_layer(2,0); /*2 nodes/tuple, do not add ghost nodes*/
 137   //FEM_Add_ghost_layer(2,1); /*2 nodes/tuple, do not add ghost nodes*/
 138   const static int tri2edge[6]={0,1, 1,2, 2,0};
 139   FEM_Add_ghost_elem(0,3,tri2edge);
 140   CkPrintf("Finished with init\n");
 141 }
 142
 143 struct myGlobals {
 144   int nnodes,maxnodes;
 145   int nelems,maxelems;
 146   int *conn; //Element connectivity table
 147   vector2d *coord; //Undeformed coordinates of each node
 148   vector2d *R_net, *d, *v, *a; //Physical fields of each node
 149   double *m_i; //Inverse of mass at each node
 150   int m_i_fid; //Field ID for m_i
 151   int *validNode,*validElem;
 152   double *S11, *S22, *S12; //Stresses for each element
 153 };
 154
 155 void pup_myGlobals(pup_er p,myGlobals *g)
 156 {
 157   FEM_Print("-------- called pup routine -------");
 158   pup_int(p,&g->nnodes);
 159   pup_int(p,&g->nelems);
 160   pup_int(p,&g->maxelems);
 161   pup_int(p,&g->maxnodes);
 162   int nnodes=g->nnodes, nelems=g->nelems;
 163   if (pup_isUnpacking(p)) {
 164     g->coord=new vector2d[g->maxnodes];
 165     g->conn=new int[3*g->maxelems];
 166     g->R_net=new vector2d[g->maxnodes]; //Net force
 167     g->d=new vector2d[g->maxnodes];//Node displacement
 168     g->v=new vector2d[g->maxnodes];//Node velocity
 169     g->a=new vector2d[g->maxnodes];
 170     g->m_i=new double[g->maxnodes];
 171     g->S11=new double[g->maxelems];
 172     g->S22=new double[g->maxelems];
 173     g->S12=new double[g->maxelems];
 174     g->validNode = new int[g->maxnodes];
 175     g->validElem = new int[g->maxelems];
 176   }
 177   pup_doubles(p,(double *)g->coord,2*nnodes);
 178   pup_ints(p,(int *)g->conn,3*nelems);
 179   pup_doubles(p,(double *)g->R_net,2*nnodes);
 180   pup_doubles(p,(double *)g->d,2*nnodes);
 181   pup_doubles(p,(double *)g->v,2*nnodes);
 182   pup_doubles(p,(double *)g->a,2*nnodes);
 183   pup_doubles(p,(double *)g->m_i,nnodes);
 184   pup_doubles(p,(double *)g->S11,nelems);
 185   pup_doubles(p,(double *)g->S22,nelems);
 186   pup_doubles(p,(double *)g->S12,nelems);
 187   pup_ints(p,(int *)g->validNode,nnodes);
 188   pup_ints(p,(int *)g->validElem,nelems);
 189   if (pup_isDeleting(p)) {
 190     delete[] g->coord;
 191     delete[] g->conn;
 192     delete[] g->R_net;
 193     delete[] g->d;
 194     delete[] g->v;
 195     delete[] g->a;
 196     delete[] g->m_i;
 197     delete[] g->S11;
 198     delete[] g->S22;
 199     delete[] g->S12;
 200     delete[] g->validNode;
 201     delete[] g->validElem;
 202   }
 203 }
 204
 205 //Return the signed area of triangle i
 206 double calcArea(myGlobals &g, int i)
 207 {
 208   int n1=g.conn[3*i+0];
 209   int n2=g.conn[3*i+1];
 210   int n3=g.conn[3*i+2];
 211   vector2d a=g.coord[n1];
 212   vector2d b=g.coord[n2];
 213   vector2d c=g.coord[n3];
 214   c-=a; b-=a;
 215   double area=0.5*(b.x*c.y-c.x*b.y);
 216   return area;
 217 }
 218
 219 // Check the quality of triangle i
 220 void checkTriangle(myGlobals &g, int i)
 221 {
 222   double area=calcArea(g,i);
 223   if (area<0) {
 224     CkError("Triangle %d of chunk %d is inverted! (area=%g)\n",
 225             i,FEM_My_partition(),area);
 226     CkAbort("Inverted triangle");
 227   }
 228   if (area<1.0e-15) {
 229     CkError("Triangle %d of chunk %d is a sliver!\n",i,FEM_My_partition());
 230     CkAbort("Sliver triangle");
 231   }
 232 }
 233
 234 //Compute forces on constant-strain triangles:
 235 void CST_NL(const vector2d *coor,const int *lm,vector2d *R_net,
 236             const vector2d *d,const double *c,
 237             int numnp,int numel,
 238             double *S11o,double *S22o,double *S12o);
 239
 240 //Update node position, velocity, accelleration based on net force.
 241 void advanceNodes(const double dt,int nnodes,const vector2d *coord,
 242                   vector2d *R_net,vector2d *a,vector2d *v,vector2d *d,
 243                   const double *m_i,bool dampen)
 244 {
 245   const vector2d z(0,0);
 246   const double shearForce=1.0e-11/(dt*dt);
 247   bool someNaNs=false;
 248   int i;
 249   for (i=0;i<nnodes;i++) {
 250     vector2d R_n=R_net[i];
 251 #if NANCHECK
 252     if (((R_n.x-R_n.x)!=0)) {
 253       CkPrintf("R_net[%d]=NaN at (%.4f,%.4f)   ",i,coord[i].x,coord[i].y);
 254       CmiAbort("nan node");
 255       someNaNs=true;
 256     }
 257     if (fabs(d[i].x)>1.0) {
 258       CkPrintf("d[%d] %f large at (%.4f,%.4f)   ",i,d[i].x,coord[i].x,coord[i].y);
 259       someNaNs=true;
 260     }
 261 #endif
 262     R_net[i]=z;
 263     //Apply boundary conditions (HACK: hardcoded!)
 264     if (1) {
 265       if (coord[i].x<0.00001)
 266         R_n.y+=shearForce/m_i[i]; //Bottom edge pushed hard down
 267       if (coord[i].y>0.02-0.00001)
 268         R_n=z; //Top edge held in place
 269     }
 270     //Update displacement and velocity
 271     vector2d aNew=R_n*m_i[i];
 272     v[i]+=(dt*0.5)*(aNew+a[i]);
 273     d[i]+=dt*v[i]+(dt*dt*0.5)*aNew;
 274     a[i]=aNew;
 275     //if (coord[i].y>0.02-0.00001) d[i].y=0.0; //Top edge in horizontal slot
 276   }
 277   if (dampen)
 278     for (i=0;i<nnodes;i++)
 279       v[i]*=0.9; //Dampen velocity slightly (prevents eventual blowup)
 280
 281   if (someNaNs) {
 282     CkPrintf("Nodes all NaN!\n");
 283     CkAbort("Node forces NaN!");
 284   }
 285 }
 286
 287 //Fill out the m_i array with the inverse of the node masses
 288 void calcMasses(myGlobals &g) {
 289   int i;
 290   double *m_i=g.m_i;
 291   //Zero out node masses
 292   for (i=0;i<g.nnodes;i++) m_i[i]=0.0;
 293   //Add mass from surrounding triangles:
 294   for (i=0;i<g.nelems;i++) {
 295     if(g.validElem[i]){
 296       int n1=g.conn[3*i+0];
 297       int n2=g.conn[3*i+1];
 298       int n3=g.conn[3*i+2];
 299       double area=calcArea(g,i);
 300       //                if (1 || i%100==0) CkPrintf("Triangle %d (%d %d %d) has area %.3g\n",i,n1,n2,n3,area);
 301       double mass=0.333*density*(thickness*area);
 302       m_i[n1]+=mass;
 303       m_i[n2]+=mass;
 304       m_i[n3]+=mass;
 305     }
 306   }
 307   //Include mass from other processors
 308   FEM_Update_field(g.m_i_fid,m_i);
 309   //Invert masses to get m_i
 310   for (i=0;i<g.nnodes;i++) {
 311     double mass=m_i[i];
 312     if (mass<1.0e-10) m_i[i]=1.0; //Disconnected node (!)
 313     else m_i[i]=1.0/mass;
 314   }
 315 }
 316
 317 void init_myGlobal(myGlobals *g){
 318   g->coord = g->R_net = g->d = g->v = g->a = NULL;
 319   g->m_i = NULL;
 320   g->conn = NULL;
 321   g->S11 = g->S22 = g->S12 = NULL;
 322 }
 323
 324
 325 void resize_nodes(void *data,int *len,int *max){
 326   printf("[%d] resize nodes called len %d max %d\n",FEM_My_partition(),*len,*max);
 327   FEM_Register_entity(FEM_Mesh_default_read(),FEM_NODE,data,*len,*max,resize_nodes);
 328   myGlobals *g = (myGlobals *)data;
 329   vector2d *coord=g->coord,*R_net=g->R_net,*d=g->d,*v=g->v,*a=g->a;
 330   double *m_i=g->m_i;
 331   int *validNode = g->validNode;
 332
 333   g->coord=new vector2d[*max];
 334   g->coord[0].x = 0.9;
 335   g->coord[0].y = 0.8;
 336   g->maxnodes = *max;
 337   g->R_net=new vector2d[g->maxnodes]; //Net force
 338   g->d=new vector2d[g->maxnodes];//Node displacement
 339   g->v=new vector2d[g->maxnodes];//Node velocity
 340   g->a=new vector2d[g->maxnodes];//Node accelleration
 341   g->m_i=new double[g->maxnodes];//Node mass
 342   g->validNode = new int[g->maxnodes]; //is the node valid
 343
 344   if(coord != NULL){
 345     for(int k=0;k<*len;k++){
 346       printf("before resize node %d ( %.6f %.6f ) \n",k,coord[k].x,coord[k].y);
 347     }
 348   }
 349
 350   FEM_Register_array(FEM_Mesh_default_read(),FEM_NODE,FEM_DATA,(void *)g->coord,FEM_DOUBLE,2);
 351   FEM_Register_array(FEM_Mesh_default_read(),FEM_NODE,FEM_DATA+1,(void *)g->R_net,FEM_DOUBLE,2);
 352   FEM_Register_array(FEM_Mesh_default_read(),FEM_NODE,FEM_DATA+2,(void *)g->d,FEM_DOUBLE,2);
 353   FEM_Register_array(FEM_Mesh_default_read(),FEM_NODE,FEM_DATA+3,(void *)g->v,FEM_DOUBLE,2);
 354   FEM_Register_array(FEM_Mesh_default_read(),FEM_NODE,FEM_DATA+4,(void *)g->a,FEM_DOUBLE,2);
 355   FEM_Register_array_layout(FEM_Mesh_default_read(),FEM_NODE,FEM_DATA+5,(void *)g->m_i,g->m_i_fid);
 356   FEM_Register_array(FEM_Mesh_default_read(),FEM_NODE,FEM_VALID,(void *)g->validNode,FEM_INT,1);
 357
 358   for(int k=0;k<*len;k++){
 359     printf("after resize node %d ( %.6f %.6f )\n",k,g->coord[k].x,g->coord[k].y);
 360   }
 361
 362   if(coord != NULL){
 363     delete [] coord;
 364     delete [] R_net;
 365     delete [] d;
 366     delete [] v;
 367     delete [] a;
 368     delete [] m_i;
 369     delete [] validNode;
 370   }
 371 };
 372
 373 void resize_elems(void *data,int *len,int *max){
 374   printf("[%d] resize elems called len %d max %d\n",FEM_My_partition(),*len,*max);
 375   FEM_Register_entity(FEM_Mesh_default_read(),FEM_ELEM,data,*len,*max,resize_elems);
 376   myGlobals *g = (myGlobals *)data;
 377   int *conn=g->conn;
 378   double *S11 = g->S11,*S22 = g->S22,*S12 = g->S12;
 379   int *validElem = g->validElem;
 380
 381   g->conn = new int[3*(*max)];
 382   g->maxelems = *max;
 383   g->S11=new double[g->maxelems];
 384   g->S22=new double[g->maxelems];
 385   g->S12=new double[g->maxelems];
 386   g->validElem = new int[g->maxelems];
 387
 388   FEM_Register_array(FEM_Mesh_default_read(),FEM_ELEM,FEM_CONN,(void *)g->conn,FEM_INDEX_0,3);
 389   CkPrintf("Connectivity array starts at %p \n",g->conn);
 390   FEM_Register_array(FEM_Mesh_default_read(),FEM_ELEM,FEM_DATA,(void *)g->S11,FEM_DOUBLE,1);
 391   FEM_Register_array(FEM_Mesh_default_read(),FEM_ELEM,FEM_DATA+1,(void *)g->S22,FEM_DOUBLE,1);
 392   FEM_Register_array(FEM_Mesh_default_read(),FEM_ELEM,FEM_DATA+2,(void *)g->S12,FEM_DOUBLE,1);
 393   FEM_Register_array(FEM_Mesh_default_read(),FEM_ELEM,FEM_VALID,(void *)g->validElem,FEM_INT,1);
 394
 395   if(conn != NULL){
 396     delete [] conn;
 397     delete [] S11;
 398     delete [] S22;
 399     delete [] S12;
 400     delete [] validElem;
 401   }
 402 };
 403
 404 void repeat_after_split(void *data){
 405   myGlobals *g = (myGlobals *)data;
 406   g->nelems = FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_ELEM);
 407   g->nnodes = FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_NODE);
 408   for(int k=0;k<g->nnodes;k++){
 409     if(g->validNode[k]){
 410       printf(" node %d ( %.6f %.6f )\n",k,g->coord[k].x,g->coord[k].y);
 411     }
 412   }
 413   calcMasses(*g);
 414 };
 415
 416
 417 extern "C" void
 418 driver(void)
 419 {
 420   int ignored;
 421   int i;
 422   int myChunk=FEM_My_partition();
 423
 424   /*Add a refinement object to FEM array*/
 425   CkPrintf("[%d] begin init\n",myChunk);
 426   FEM_REFINE2D_Init();
 427   CkPrintf("[%d] end init\n",myChunk);
 428
 429   myGlobals g;
 430   FEM_Register(&g,(FEM_PupFn)pup_myGlobals);
 431   init_myGlobal(&g);
 432
 433   g.nnodes = FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_NODE);
 434   int maxNodes = g.nnodes;
 435   g.maxnodes=2*maxNodes;
 436   g.m_i_fid=FEM_Create_field(FEM_DOUBLE,1,0,sizeof(double));
 437   resize_nodes((void *)&g,&g.nnodes,&maxNodes);
 438   int nghost=0;
 439   g.nelems=FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_ELEM);
 440   g.maxelems=g.nelems;
 441   resize_elems((void *)&g,&g.nelems,&g.maxelems);
 442
 443   FEM_REFINE2D_Newmesh(FEM_Mesh_default_read(),FEM_NODE,FEM_ELEM);
 444
 445   //Initialize associated data
 446   for (i=0;i<g.maxnodes;i++) {
 447     g.R_net[i]=g.d[i]=g.v[i]=g.a[i]=vector2d(0.0);
 448   }
 449
 450   //Apply a small initial perturbation to positions
 451   for (i=0;i<g.nnodes;i++) {
 452     const double max=1.0e-15/15.0; //Tiny perturbation
 453     g.d[i].x+=max*(i&15);
 454     g.d[i].y+=max*((i+5)&15);
 455   }
 456
 457   int fid=FEM_Create_field(FEM_DOUBLE,2,0,sizeof(vector2d));
 458
 459   for (i=0;i<g.nelems;i++){
 460     checkTriangle(g,i);
 461   }
 462   sleep(5);
 463   //Timeloop
 464   if (CkMyPe()==0){
 465     CkPrintf("Entering timeloop\n");
 466   }
 467   //  int tSteps=0x70FF00FF;
 468   int tSteps=10;
 469   int z=13;
 470   calcMasses(g);
 471   double startTime=CkWallTimer();
 472   double curArea=2.5e-5/1024;
 473   int t = 0;
 474   if (1) { //Publish data to the net
 475     NetFEM n=NetFEM_Begin(myChunk,t,2,NetFEM_WRITE);
 476     int count=0;
 477     double *vcoord = new double[2*g.nnodes];
 478     double *vnodeid = new double[g.nnodes];
 479     int *maptovalid = new int[g.nnodes];
 480     for(int i=0;i<g.nnodes;i++){
 481       if(g.validNode[i]){
 482         vcoord[2*count] = ((double *)g.coord)[2*i];
 483         vcoord[2*count+1] = ((double *)g.coord)[2*i+1];
 484         maptovalid[i] = count;
 485         printf("~~~~~~~ %d %d %.6lf %.6lf \n",count,i,vcoord[2*count],vcoord[2*count+1]);
 486         vnodeid[count] = i;
 487         count++;
 488       }
 489     }
 490     NetFEM_Nodes(n,count,(double *)vcoord,"Position (m)");
 491     NetFEM_Scalar(n,vnodeid,1,"Node ID");
 492     /*    NetFEM_Vector(n,(double *)g.d,"Displacement (m)");
 493           NetFEM_Vector(n,(double *)g.v,"Velocity (m/s)");*/
 494     count=0;
 495     int *vconn = new int[3*g.nelems];
 496     double *vid = new double[3*g.nelems];
 497     for(int i=0;i<g.nelems;i++){
 498       if(g.validElem[i]){
 499         vconn[3*count] = maptovalid[g.conn[3*i]];
 500         vconn[3*count+1] = maptovalid[g.conn[3*i+1]];
 501         vconn[3*count+2] = maptovalid[g.conn[3*i+2]];
 502         printf("~~~~~~~ %d %d < %d,%d %d,%d %d,%d >\n",count,i,vconn[3*count],g.conn[3*i],vconn[3*count+1],g.conn[3*i+1],vconn[3*count+2],g.conn[3*i+2]);
 503         vid[count]=count;
 504         count++;
 505       }
 506     }
 507     NetFEM_Elements(n,count,3,(int *)vconn,"Triangles");
 508     NetFEM_Scalar(n,vid,1,"Element ID");
 509     /*  NetFEM_Scalar(n,g.S22,1,"Y Stress (pure)");
 510         NetFEM_Scalar(n,g.S12,1,"Shear Stress (pure)");*/
 511     NetFEM_End(n);
 512     delete [] vcoord;
 513     delete [] vconn;
 514     delete [] maptovalid;
 515     delete [] vid;
 516     delete [] vnodeid;
 517   }
 518   for (t=1;t<=tSteps;t++) {
 519     /*    if (1) { //Structural mechanics
 520     //Compute forces on nodes exerted by elements
 521     CST_NL(g.coord,g.conn,g.R_net,g.d,matConst,g.nnodes,g.nelems,g.S11,g.S22,g.S12);
 522     //Communicate net force on shared nodes
 523     FEM_Update_field(fid,g.R_net);
 524     //Advance node positions
 525     advanceNodes(dt,g.nnodes,g.coord,g.R_net,g.a,g.v,g.d,g.m_i,(t%4)==0);
 526     }*/
 527
 528     //Debugging/perf. output
 529     double curTime=CkWallTimer();
 530     double total=curTime-startTime;
 531     startTime=curTime;
 532     /*    if (CkMyPe()==0 && (t%64==0))
 533           CkPrintf("%d %.6f sec for loop %d \n",CkNumPes(),total,t);*/
 534     /*   if (0 && t%16==0) {
 535          CkPrintf("    Triangle 0:\n");
 536          for (int j=0;j<3;j++) {
 537          int n=g.conn[0][j];
 538          CkPrintf("    Node %d: coord=(%.4f,%.4f)  d=(%.4g,%.4g)\n",
 539          n,g.coord[n].x,g.coord[n].y,g.d[n].x,g.d[n].y);
 540          }
 541          }*/
 542     //    if (t%512==0)
 543     //      FEM_Migrate();
 544
 545     double *areas=new double[g.nelems];
 546     for (i=0;i<g.nelems;i++) {
 547       areas[i]=1.5*calcArea(g,i);
 548 //      areas[i]=calcArea(g,i);
 549     }
 550
 551     //coarsen all steps
 552 //    areas[z] *= 2.0;
 553     z += 3;
 554     CkPrintf("[%d] Starting coarsening step: %d nodes, %d elements to %.3g\n", myChunk,g.nnodes,g.nelems,curArea);
 555     FEM_REFINE2D_Coarsen(FEM_Mesh_default_read(),FEM_NODE,(double *)g.coord,FEM_ELEM,areas);
 556     repeat_after_split((void *)&g);
 557     g.nelems = FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_ELEM);
 558     g.nnodes = FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_NODE);
 559     CkPrintf("[%d] Done with coarsening step: %d nodes, %d elements\n",
 560              myChunk,g.nnodes,g.nelems);
 561     if (1) { //Publish data to the net
 562       NetFEM n=NetFEM_Begin(myChunk,t,2,NetFEM_WRITE);
 563       int count=0;
 564       double *vcoord = new double[2*g.nnodes];
 565       double *vnodeid = new double[g.nnodes];
 566       int *maptovalid = new int[g.nnodes];
 567       for(int i=0;i<g.nnodes;i++){
 568         maptovalid[i] = -1;
 569         if(g.validNode[i]){
 570           vcoord[2*count] = ((double *)g.coord)[2*i];
 571           vcoord[2*count+1] = ((double *)g.coord)[2*i+1];
 572           maptovalid[i] = count;
 573           printf("node~~~~~~~ %d %d %.6lf %.6lf \n",count,i,vcoord[2*count],vcoord[2*count+1]);
 574           vnodeid[count] = i;
 575           count++;
 576         }
 577       }
 578       NetFEM_Nodes(n,count,(double *)vcoord,"Position (m)");
 579       NetFEM_Scalar(n,vnodeid,1,"Node ID");
 580       /*    NetFEM_Vector(n,(double *)g.d,"Displacement (m)");
 581             NetFEM_Vector(n,(double *)g.v,"Velocity (m/s)");*/
 582       count=0;
 583       int *vconn = new int[3*g.nelems];
 584       double *vid = new double[3*g.nelems];
 585       for(int i=0;i<g.nelems;i++){
 586         if(g.validElem[i]){
 587           vconn[3*count] = maptovalid[g.conn[3*i]];
 588           vconn[3*count+1] = maptovalid[g.conn[3*i+1]];
 589           vconn[3*count+2] = maptovalid[g.conn[3*i+2]];
 590           printf("~~~~~~~ %d %d < %d,%d %d,%d %d,%d >\n",count,i,vconn[3*count],g.conn[3*i],vconn[3*count+1],g.conn[3*i+1],vconn[3*count+2],g.conn[3*i+2]);
 591           vid[count]=count;
 592           count++;
 593         }
 594       }
 595       NetFEM_Elements(n,count,3,(int *)vconn,"Triangles");
 596       NetFEM_Scalar(n,vid,1,"Element ID");
 597       /*        NetFEM_Scalar(n,g.S22,1,"Y Stress (pure)");
 598                 NetFEM_Scalar(n,g.S12,1,"Shear Stress (pure)");*/
 599       NetFEM_End(n);
 600       delete [] vcoord;
 601       delete [] vconn;
 602       delete [] maptovalid;
 603       delete [] vid;
 604       delete [] vnodeid;
 605     }
 606   }
 607   if (CkMyPe()==0)
 608     CkPrintf("Driver finished\n");
 609 }
 610
 611