Charm++: Charm Source Code Documentation

00001 /*
00002 Triangle Mesh Refinement & Coarsening (TMRC2D) demo code.
00003 
00004 Reads in a Triangle mesh.
00005 Passes it to the REFINE2D subsystem.
00006 Remeshes via a sequence of refinements and coarsenings.
00007 */
00008 #include <stdlib.h>
00009 #include <unistd.h>
00010 #include <stdio.h>
00011 #include <math.h>
00012 #include "charm++.h"
00013 #include "fem.h"
00014 /*My HEADERS*/
00015 #include "ckvector3d.h"
00016 #include "charm-api.h"
00017 #include "fem_mesh.h"
00018 #include "netfem.h"
00019 #include "refine.h"
00020 #include "femrefine.h"
00021 #include "pgm.h"
00022 
00023 //The material constants c, as computed by fortran mat_const
00024 // I think the units here are Pascals (N/m^2)
00025 const double matConst[4]={3.692e9,  1.292e9,  3.692e9,  1.200e9 };
00026 
00027 //Material density, Kg/m^3
00028 const double density=5.0*1000.0;
00029 //Plate thickness, meters
00030 const double thickness=0.0001;
00031 
00032 //The timestep, in seconds
00033 // This aught to be adjusted when the mesh changes!
00034 const double dt=1.0e-12;
00035 
00036 static void die(const char *str) {
00037   CkError("Fatal error: %s\n",str);
00038   CkExit();
00039 }
00040 
00042 void resize_nodes(void *data,int *len,int *max);
00043 void resize_elems(void *data,int *len,int *max);
00044 void resize_edges(void *data,int *len,int *max);
00045 
00046 #define NANCHECK 1 /*Check for NaNs at each timestep*/
00047 
00048 extern "C" void
00049 init(void)
00050 {
00051   CkPrintf("init started\n");
00052 
00053   const char *eleName="out.1024.ele";
00054   const char *nodeName="out.1024.node";
00055   const char *edgeName="out.1024.edge";
00056   int nPts=0; //Number of nodes
00057   vector2d *pts=0; //Node coordinates
00058 
00059   CkPrintf("Reading node coordinates from %s\n",nodeName);
00060   //Open and read the node coordinate file
00061   {
00062     char line[1024];
00063     FILE *f=fopen(nodeName,"r");
00064     if (f==NULL) die("Can't open node file!");
00065     fgets(line,1024,f);
00066     if (1!=sscanf(line,"%d",&nPts)) die("Can't read number of points!");
00067     pts=new vector2d[nPts];
00068     for (int i=0;i<nPts;i++) {
00069       int ptNo;
00070       if (NULL==fgets(line,1024,f)) die("Can't read node input line!");
00071       if (3!=sscanf(line,"%d%lf%lf",&ptNo,&pts[i].x,&pts[i].y)) 
00072     die("Can't parse node input line!");
00073       pts[i].y = 4.0*pts[i].y;
00074     }
00075     fclose(f);
00076   }
00077   CkPrintf("Passing node coords to framework\n");
00078   
00079   // Register the node entity and its data arrays that will be used later. This
00080   // needs to be done so that the width of the data segments are set correctly 
00081   // and can be used later
00082   FEM_Register_entity(FEM_Mesh_default_write(),FEM_NODE,NULL,nPts,nPts,resize_nodes);
00083   for(int k=0;k<=4;k++){
00084     if(k != 0){
00085       vector2d *t = new vector2d[nPts];
00086       FEM_Register_array(FEM_Mesh_default_write(),FEM_NODE,FEM_DATA+k,t,FEM_DOUBLE,2);
00087     }else{
00088       FEM_Register_array(FEM_Mesh_default_write(),FEM_NODE,FEM_DATA+k,pts,FEM_DOUBLE,2);
00089     }
00090   }
00091   double *td = new double[nPts];
00092   FEM_Register_array(FEM_Mesh_default_write(),FEM_NODE,FEM_DATA+5,td,FEM_DOUBLE,1);
00093   int *validNodes = new int[nPts];
00094   for(int ii=0;ii<nPts;ii++){
00095     validNodes[ii]=1;
00096   }
00097   FEM_Register_array(FEM_Mesh_default_write(),FEM_NODE,FEM_VALID,validNodes,FEM_INT,1);
00098   
00099   int nEle=0;
00100   int *ele=NULL;
00101   CkPrintf("Reading elements from %s\n",eleName);
00102   //Open and read the element connectivity file
00103   {
00104     char line[1024];
00105     FILE *f=fopen(eleName,"r");
00106     if (f==NULL) die("Can't open element file!");
00107     fgets(line,1024,f);
00108     if (1!=sscanf(line,"%d",&nEle)) die("Can't read number of elements!");
00109     ele=new int[3*nEle];
00110     for (int i=0;i<nEle;i++) {
00111       int elNo;
00112       if (NULL==fgets(line,1024,f)) die("Can't read element input line!");
00113       if (4!=sscanf(line,"%d%d%d%d",&elNo,&ele[3*i+0],&ele[3*i+1],&ele[3*i+2])) 
00114     die("Can't parse element input line!");
00115       ele[3*i+0]--; //Fortran to C indexing
00116       ele[3*i+1]--; //Fortran to C indexing
00117       ele[3*i+2]--; //Fortran to C indexing
00118     }
00119     fclose(f);
00120   }
00121   
00122   CkPrintf("Passing elements to framework\n");
00123   // Register the Element entity and its connectivity array. Register the
00124   // data arrays to set up the widths correctly at the beginning
00125   FEM_Register_entity(FEM_Mesh_default_write(),FEM_ELEM,NULL,nEle,nEle,resize_elems);
00126   FEM_Register_array(FEM_Mesh_default_write(),FEM_ELEM,FEM_CONN,ele,FEM_INDEX_0,3);
00127   
00128   for(int k=0;k<3;k++){
00129     void *t = new double[nEle];
00130     FEM_Register_array(FEM_Mesh_default_write(),FEM_ELEM,FEM_DATA+k,t,FEM_DOUBLE,1);
00131   }
00132   int *validElem = new int[nEle];
00133   for(int ii=0;ii<nEle;ii++){
00134     validElem[ii]=1;
00135   }
00136   FEM_Register_array(FEM_Mesh_default_write(),FEM_ELEM,FEM_VALID,validElem,FEM_INT,1);
00137   
00138   /*Build the ghost layer for refinement border*/
00139   FEM_Add_ghost_layer(2,0); /*2 nodes/tuple, do not add ghost nodes*/
00140   //FEM_Add_ghost_layer(2,1); /*2 nodes/tuple, do not add ghost nodes*/
00141   const static int tri2edge[6]={0,1, 1,2, 2,0};
00142   FEM_Add_ghost_elem(0,3,tri2edge);
00143 
00144   //open and read the .edge (edge connectivity file) needed for boundary values
00145   int nEdge;
00146   int *edgeConn;
00147   int *edgeBoundary;
00148     int *validEdge;
00149   {
00150     char line[1024];
00151     FILE *f=fopen(edgeName,"r");
00152     if (f==NULL) die("Can't open edge file!");
00153     fgets(line,1024,f);
00154     if (1!=sscanf(line,"%d",&nEdge)) die("Can't read number of elements!");
00155     edgeConn = new int[2*nEdge];
00156     edgeBoundary = new int[nEdge];
00157         validEdge = new int[nEdge];
00158     for(int i=0;i<nEdge;i++){
00159       int edgeNo;
00160       if (NULL==fgets(line,1024,f)) die("Can't read edge input line!");
00161       if (4 != sscanf(line,"%d%d%d%d",&edgeNo,&edgeConn[i*2+0],&edgeConn[i*2+1],&edgeBoundary[i])){
00162     die("Can't parse edge input line!");
00163       }
00164       edgeConn[i*2+0]--;
00165       edgeConn[i*2+1]--;        
00166             validEdge[i] = 1;
00167     }
00168     fclose(f);
00169   }
00170   printf("Number of edges %d \n",nEdge);
00171   FEM_Register_entity(FEM_Mesh_default_write(),FEM_SPARSE,NULL,nEdge,nEdge,resize_edges);
00172   FEM_Register_array(FEM_Mesh_default_write(),FEM_SPARSE,FEM_CONN,edgeConn,FEM_INDEX_0,2);
00173   FEM_Register_array(FEM_Mesh_default_write(),FEM_SPARSE,FEM_BOUNDARY,edgeBoundary,FEM_INT,1);
00174     FEM_Register_array(FEM_Mesh_default_write(),FEM_SPARSE,FEM_VALID,validEdge,FEM_INT,1);
00175   CkPrintf("Finished with init\n");
00176 }
00177 
00178 struct myGlobals {
00179   int nnodes,maxnodes;
00180   int nelems,maxelems;
00181   int nedges,maxedges;
00182   int *conn; //Element connectivity table
00183   vector2d *coord; //Undeformed coordinates of each node
00184   vector2d *R_net, *d, *v, *a; //Physical fields of each node
00185   double *m_i; //Inverse of mass at each node
00186   int m_i_fid; //Field ID for m_i
00187   int *validNode,*validElem;
00188   int *nodeBoundary; //node boundary
00189   double *S11, *S22, *S12; //Stresses for each element
00190   int *edgeConn;//edge Connectivity table
00191   int *edgeBoundary;//edge boundary value
00192   int *validEdge;
00193 };
00194 
00195 void pup_myGlobals(pup_er p,myGlobals *g) 
00196 {
00197   FEM_Print("-------- called pup routine -------");
00198   pup_int(p,&g->nnodes);
00199   pup_int(p,&g->nelems);
00200   pup_int(p,&g->nedges);
00201   pup_int(p,&g->maxelems);
00202   pup_int(p,&g->maxnodes);
00203   pup_int(p,&g->maxedges);
00204 
00205   int nnodes=g->nnodes, nelems=g->nelems;
00206   if (pup_isUnpacking(p)) {
00207     g->coord=new vector2d[g->maxnodes];
00208     g->conn=new int[3*g->maxelems];
00209     g->R_net=new vector2d[g->maxnodes]; //Net force
00210     g->d=new vector2d[g->maxnodes];//Node displacement
00211     g->v=new vector2d[g->maxnodes];//Node velocity
00212     g->a=new vector2d[g->maxnodes];
00213     g->m_i=new double[g->maxnodes];
00214     g->S11=new double[g->maxelems];
00215     g->S22=new double[g->maxelems];
00216     g->S12=new double[g->maxelems];
00217     g->validNode = new int[g->maxnodes];
00218     g->validElem = new int[g->maxelems];
00219     g->edgeConn = new int[2*g->maxedges];
00220     g->edgeBoundary = new int[g->maxedges];
00221         g->validEdge = new int[g->maxedges];
00222   }
00223   pup_doubles(p,(double *)g->coord,2*nnodes);
00224   pup_ints(p,(int *)g->conn,3*nelems);
00225   pup_doubles(p,(double *)g->R_net,2*nnodes);
00226   pup_doubles(p,(double *)g->d,2*nnodes);
00227   pup_doubles(p,(double *)g->v,2*nnodes);
00228   pup_doubles(p,(double *)g->a,2*nnodes);
00229   pup_doubles(p,(double *)g->m_i,nnodes);
00230   pup_doubles(p,(double *)g->S11,nelems);
00231   pup_doubles(p,(double *)g->S22,nelems);
00232   pup_doubles(p,(double *)g->S12,nelems);
00233   pup_ints(p,(int *)g->validNode,nnodes);
00234   pup_ints(p,(int *)g->validElem,nelems);
00235   pup_ints(p,(int *)g->edgeConn,2*g->nedges);
00236   pup_ints(p,(int *)g->edgeBoundary,g->nedges);
00237     pup_ints(p,(int *)g->validEdge,g->nedges);
00238 
00239   if (pup_isDeleting(p)) {
00240     delete[] g->coord;
00241     delete[] g->conn;
00242     delete[] g->R_net;
00243     delete[] g->d;
00244     delete[] g->v;
00245     delete[] g->a;
00246     delete[] g->m_i;
00247     delete[] g->S11;
00248     delete[] g->S22;
00249     delete[] g->S12;
00250     delete[] g->validNode;
00251     delete[] g->validElem;
00252     delete[] g->edgeConn;
00253     delete[] g->edgeBoundary;
00254         delete[] g->validEdge;
00255   }
00256 }
00257 
00258 //Return the signed area of triangle i
00259 double calcArea(myGlobals &g, int i)
00260 {
00261   int n1=g.conn[3*i+0];
00262   int n2=g.conn[3*i+1];
00263   int n3=g.conn[3*i+2];
00264   vector2d a=g.coord[n1];
00265   vector2d b=g.coord[n2];
00266   vector2d c=g.coord[n3];
00267   c-=a; b-=a;
00268   double area=0.5*(b.x*c.y-c.x*b.y);
00269   return area;
00270 }
00271 
00272 // Check the quality of triangle i
00273 void checkTriangle(myGlobals &g, int i)
00274 {
00275   double area=calcArea(g,i);
00276   if (area<0) {
00277     CkError("Triangle %d of chunk %d is inverted! (area=%g)\n",
00278         i,FEM_My_partition(),area);
00279     CkAbort("Inverted triangle");
00280   }
00281   if (area<1.0e-15) {
00282     CkError("Triangle %d of chunk %d is a sliver!\n",i,FEM_My_partition());
00283     CkAbort("Sliver triangle");
00284   }
00285 }
00286 
00287 //Compute forces on constant-strain triangles:
00288 void CST_NL(const vector2d *coor,const int *lm,vector2d *R_net,
00289         const vector2d *d,const double *c,
00290         int numnp,int numel,
00291         double *S11o,double *S22o,double *S12o);
00292 
00293 //Update node position, velocity, accelleration based on net force.
00294 void advanceNodes(const double dt,int nnodes,const vector2d *coord,
00295           vector2d *R_net,vector2d *a,vector2d *v,vector2d *d,
00296           const double *m_i,bool dampen)
00297 {
00298   const vector2d z(0,0);
00299   const double shearForce=1.0e-11/(dt*dt);
00300   bool someNaNs=false;
00301   int i;
00302   for (i=0;i<nnodes;i++) {
00303     vector2d R_n=R_net[i];
00304 #if NANCHECK
00305     if (((R_n.x-R_n.x)!=0)) {
00306       CkPrintf("R_net[%d]=NaN at (%.4f,%.4f)   ",i,coord[i].x,coord[i].y);
00307       CmiAbort("nan node");
00308       someNaNs=true;
00309     }
00310     if (fabs(d[i].x)>1.0) {
00311       CkPrintf("d[%d] %f large at (%.4f,%.4f)   ",i,d[i].x,coord[i].x,coord[i].y);
00312       someNaNs=true;
00313     }
00314 #endif
00315     R_net[i]=z;
00316     //Apply boundary conditions (HACK: hardcoded!)
00317     if (1) {
00318       if (coord[i].x<0.00001)
00319     R_n.y+=shearForce/m_i[i]; //Bottom edge pushed hard down
00320       if (coord[i].y>0.02-0.00001)
00321     R_n=z; //Top edge held in place
00322     }
00323     //Update displacement and velocity
00324     vector2d aNew=R_n*m_i[i];
00325     v[i]+=(dt*0.5)*(aNew+a[i]);
00326     d[i]+=dt*v[i]+(dt*dt*0.5)*aNew;
00327     a[i]=aNew;   
00328     //if (coord[i].y>0.02-0.00001) d[i].y=0.0; //Top edge in horizontal slot
00329   }
00330   if (dampen)
00331     for (i=0;i<nnodes;i++)
00332       v[i]*=0.9; //Dampen velocity slightly (prevents eventual blowup)
00333   
00334   if (someNaNs) {
00335     CkPrintf("Nodes all NaN!\n");
00336     CkAbort("Node forces NaN!");
00337   }
00338 }
00339 
00340 //Fill out the m_i array with the inverse of the node masses
00341 void calcMasses(myGlobals &g) {
00342   int i;
00343   double *m_i=g.m_i;
00344   //Zero out node masses
00345   for (i=0;i<g.nnodes;i++) m_i[i]=0.0;
00346   //Add mass from surrounding triangles:
00347   for (i=0;i<g.nelems;i++) {
00348     if(g.validElem[i]){
00349       int n1=g.conn[3*i+0];
00350       int n2=g.conn[3*i+1];
00351       int n3=g.conn[3*i+2];
00352       double area=calcArea(g,i);
00353       //        if (1 || i%100==0) CkPrintf("Triangle %d (%d %d %d) has area %.3g\n",i,n1,n2,n3,area);
00354       double mass=0.333*density*(thickness*area);
00355       m_i[n1]+=mass;
00356       m_i[n2]+=mass;
00357       m_i[n3]+=mass;
00358     }
00359   }
00360   //Include mass from other processors
00361   FEM_Update_field(g.m_i_fid,m_i);
00362   //Invert masses to get m_i
00363   for (i=0;i<g.nnodes;i++) {
00364     double mass=m_i[i];
00365     if (mass<1.0e-10) m_i[i]=1.0; //Disconnected node (!)
00366     else m_i[i]=1.0/mass;
00367   }
00368 }
00369 
00370 void init_myGlobal(myGlobals *g){
00371   g->coord = g->R_net = g->d = g->v = g->a = NULL;
00372   g->m_i = NULL;
00373   g->conn = NULL;
00374   g->S11 = g->S22 = g->S12 = NULL;
00375 }
00376 
00377 
00378 void resize_nodes(void *data,int *len,int *max){
00379   printf("[%d] resize nodes called len %d max %d\n",FEM_My_partition(),*len,*max);
00380   FEM_Register_entity(FEM_Mesh_default_read(),FEM_NODE,data,*len,*max,resize_nodes);
00381   myGlobals *g = (myGlobals *)data;
00382   vector2d *coord=g->coord,*R_net=g->R_net,*d=g->d,*v=g->v,*a=g->a;
00383   double *m_i=g->m_i;
00384   int *validNode = g->validNode;
00385   
00386   g->coord=new vector2d[*max];
00387   g->coord[0].x = 0.9;
00388   g->coord[0].y = 0.8;
00389   g->maxnodes = *max;
00390   g->R_net=new vector2d[g->maxnodes]; //Net force
00391   g->d=new vector2d[g->maxnodes];//Node displacement
00392   g->v=new vector2d[g->maxnodes];//Node velocity
00393   g->a=new vector2d[g->maxnodes];//Node accelleration
00394   g->m_i=new double[g->maxnodes];//Node mass
00395   g->validNode = new int[g->maxnodes]; //is the node valid
00396   
00397   if(coord != NULL){
00398     for(int k=0;k<*len;k++){
00399       printf("before resize node %d ( %.6f %.6f ) \n",k,coord[k].x,coord[k].y);
00400     }
00401   } 
00402   
00403   FEM_Register_array(FEM_Mesh_default_read(),FEM_NODE,FEM_DATA,(void *)g->coord,FEM_DOUBLE,2);
00404   FEM_Register_array(FEM_Mesh_default_read(),FEM_NODE,FEM_DATA+1,(void *)g->R_net,FEM_DOUBLE,2);
00405   FEM_Register_array(FEM_Mesh_default_read(),FEM_NODE,FEM_DATA+2,(void *)g->d,FEM_DOUBLE,2);
00406   FEM_Register_array(FEM_Mesh_default_read(),FEM_NODE,FEM_DATA+3,(void *)g->v,FEM_DOUBLE,2);
00407   FEM_Register_array(FEM_Mesh_default_read(),FEM_NODE,FEM_DATA+4,(void *)g->a,FEM_DOUBLE,2);
00408   FEM_Register_array_layout(FEM_Mesh_default_read(),FEM_NODE,FEM_DATA+5,(void *)g->m_i,g->m_i_fid);
00409   FEM_Register_array(FEM_Mesh_default_read(),FEM_NODE,FEM_VALID,(void *)g->validNode,FEM_INT,1);
00410   
00411   for(int k=0;k<*len;k++){
00412     printf("after resize node %d ( %.6f %.6f )\n",k,g->coord[k].x,g->coord[k].y);
00413   }
00414   
00415   if(coord != NULL){
00416     delete [] coord;
00417     delete [] R_net;
00418     delete [] d;
00419     delete [] v;
00420     delete [] a;
00421     delete [] m_i;
00422     delete [] validNode;
00423   }
00424 };
00425 
00426 void resize_elems(void *data,int *len,int *max){
00427   printf("[%d] resize elems called len %d max %d\n",FEM_My_partition(),*len,*max);
00428   FEM_Register_entity(FEM_Mesh_default_read(),FEM_ELEM,data,*len,*max,resize_elems);
00429   myGlobals *g = (myGlobals *)data;
00430   int *conn=g->conn;
00431   double *S11 = g->S11,*S22 = g->S22,*S12 = g->S12;
00432   int *validElem = g->validElem;
00433   
00434   g->conn = new int[3*(*max)];
00435   g->maxelems = *max;
00436   g->S11=new double[g->maxelems];
00437   g->S22=new double[g->maxelems];
00438   g->S12=new double[g->maxelems];
00439   g->validElem = new int[g->maxelems];
00440   
00441   FEM_Register_array(FEM_Mesh_default_read(),FEM_ELEM,FEM_CONN,(void *)g->conn,FEM_INDEX_0,3);  
00442   CkPrintf("Connectivity array starts at %p \n",g->conn);
00443   FEM_Register_array(FEM_Mesh_default_read(),FEM_ELEM,FEM_DATA,(void *)g->S11,FEM_DOUBLE,1);    
00444   FEM_Register_array(FEM_Mesh_default_read(),FEM_ELEM,FEM_DATA+1,(void *)g->S22,FEM_DOUBLE,1);  
00445   FEM_Register_array(FEM_Mesh_default_read(),FEM_ELEM,FEM_DATA+2,(void *)g->S12,FEM_DOUBLE,1);  
00446   FEM_Register_array(FEM_Mesh_default_read(),FEM_ELEM,FEM_VALID,(void *)g->validElem,FEM_INT,1);    
00447   
00448   if(conn != NULL){
00449     delete [] conn;
00450     delete [] S11;
00451     delete [] S22;
00452     delete [] S12;
00453     delete [] validElem;
00454   }
00455 };
00456 
00457 void resize_edges(void *data,int *len,int *max){
00458   printf("[%d] resize edges called len %d max %d\n",FEM_My_partition(),*len,*max);
00459   FEM_Register_entity(FEM_Mesh_default_read(),FEM_SPARSE,data,*len,*max,resize_edges);
00460   myGlobals *g = (myGlobals *)data;
00461   
00462   int *conn = g->edgeConn;
00463   int *bound = g->edgeBoundary;
00464     int *validEdge = g->validEdge;
00465   g->maxedges = *max;   
00466   g->edgeConn = new int[2*(*max)];
00467   g->edgeBoundary = new int[(*max)];
00468     g->validEdge = new int[(*max)];
00469   
00470   FEM_Register_array(FEM_Mesh_default_read(),FEM_SPARSE,FEM_CONN,(void *)g->edgeConn,FEM_INDEX_0,2);    
00471   FEM_Register_array(FEM_Mesh_default_read(),FEM_SPARSE,FEM_BOUNDARY,(void *)g->edgeBoundary,FEM_INT,1);
00472   FEM_Register_array(FEM_Mesh_default_read(),FEM_SPARSE,FEM_VALID,(void *)g->validEdge,FEM_INT,1);
00473   if(conn != NULL){
00474     delete [] conn;
00475     delete [] bound;    
00476         delete [] validEdge;
00477   }
00478 }
00479 
00480 void repeat_after_split(void *data){
00481   myGlobals *g = (myGlobals *)data;
00482   g->nelems = FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_ELEM);
00483   g->nnodes = FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_NODE);
00484   for(int k=0;k<g->nnodes;k++){
00485     if(g->validNode[k]){
00486       printf(" node %d ( %.6f %.6f )\n",k,g->coord[k].x,g->coord[k].y);
00487     }   
00488   }
00489   calcMasses(*g);
00490 };
00491 
00492 
00493 void publishMeshToNetFEM(myGlobals &g,int myChunk,int t);
00494 int countValidEntities(int *validData,int total);
00495 
00496 
00497 extern "C" void
00498 driver(void)
00499 {
00500   int ignored;
00501   int i, count;  
00502   int myChunk=FEM_My_partition();
00503   
00504   /*Add a refinement object to FEM array*/
00505   CkPrintf("[%d] begin init\n",myChunk);
00506   FEM_REFINE2D_Init();
00507   CkPrintf("[%d] end init\n",myChunk);
00508   
00509   myGlobals g;
00510   FEM_Register(&g,(FEM_PupFn)pup_myGlobals);
00511   init_myGlobal(&g);
00512   
00513   g.nnodes = FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_NODE);
00514   int maxNodes = g.nnodes;
00515   g.maxnodes=2*maxNodes;
00516   g.m_i_fid=FEM_Create_field(FEM_DOUBLE,1,0,sizeof(double));
00517   resize_nodes((void *)&g,&g.nnodes,&maxNodes);
00518   int nghost=0;
00519   g.nelems=FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_ELEM);
00520   g.maxelems=g.nelems;
00521   resize_elems((void *)&g,&g.nelems,&g.maxelems);
00522 
00523   FEM_REFINE2D_Newmesh(FEM_Mesh_default_read(),FEM_NODE,FEM_ELEM);
00524   
00525   //Initialize associated data
00526   for (i=0;i<g.maxnodes;i++) {
00527     g.R_net[i]=g.d[i]=g.v[i]=g.a[i]=vector2d(0.0);
00528   }
00529   
00530   //Apply a small initial perturbation to positions
00531   for (i=0;i<g.nnodes;i++) {
00532     const double max=1.0e-15/15.0; //Tiny perturbation
00533     g.d[i].x+=max*(i&15);
00534     g.d[i].y+=max*((i+5)&15);
00535   }
00536   
00537   int fid=FEM_Create_field(FEM_DOUBLE,2,0,sizeof(vector2d));
00538   
00539   for (i=0;i<g.nelems;i++){
00540     checkTriangle(g,i);
00541   } 
00542   sleep(5);
00543   //Timeloop
00544   if (CkMyPe()==0){
00545     CkPrintf("Entering timeloop\n");
00546   } 
00547   //  int tSteps=0x70FF00FF;
00548   int tSteps=4;
00549   int z=13;
00550   calcMasses(g);
00551   double startTime=CkWallTimer();
00552   double curArea=2.5e-5/1024;
00553   int t = 0;
00554 
00555   // THIS IS THE INITIAL MESH SENT TO NetFEM
00556   if (1) { //Publish data to the net
00557     publishMeshToNetFEM(g,myChunk,t);
00558   }
00559   double desiredArea;
00560   /* 
00561   //should not be necessary as it would have been set in the init
00562   for (i=0; i<g.nnodes; i++) {
00563     g.validNode[i] = 1;
00564   }
00565   for (i=0; i<g.nelems; i++) {
00566     g.validElem[i] = 1;
00567   }*/
00568   double avgArea = 0.0;
00569   for (i=0;i<g.nelems;i++) {
00570     avgArea += calcArea(g, i);
00571   }
00572   avgArea /= g.nelems;
00573   for (t=1;t<=tSteps;t++) {
00574     /*    if (1) { //Structural mechanics
00575     //Compute forces on nodes exerted by elements
00576     CST_NL(g.coord,g.conn,g.R_net,g.d,matConst,g.nnodes,g.nelems,g.S11,g.S22,g.S12);
00577     //Communicate net force on shared nodes
00578     FEM_Update_field(fid,g.R_net);
00579     //Advance node positions
00580     advanceNodes(dt,g.nnodes,g.coord,g.R_net,g.a,g.v,g.d,g.m_i,(t%4)==0);
00581     }*/
00582     
00583     //Debugging/perf. output
00584     double curTime=CkWallTimer();
00585     double total=curTime-startTime;
00586     startTime=curTime;
00587     vector2d *loc;
00588     double *areas;
00589     
00590     // prepare to coarsen
00591     loc=new vector2d[2*g.nnodes];
00592     for (i=0;i<g.nnodes;i++) {
00593       loc[i]=g.coord[i];//+g.d[i];
00594     }
00595     areas=new double[g.nelems];
00596     for (i=0;i<g.nelems;i++) {
00597       areas[i] = avgArea;
00598     }
00599     //coarsen one element at a time
00600     //int coarseIdx = (23  + 4*t)%g.nnodes;
00601     //areas[coarseIdx] = calcArea(g,coarseIdx)*2.5;
00602     
00603     CkPrintf("[%d] Starting coarsening step: %d nodes, %d elements\n", myChunk,countValidEntities(g.validNode,g.nnodes),countValidEntities(g.validElem,g.nelems));
00604     FEM_REFINE2D_Coarsen(FEM_Mesh_default_read(),FEM_NODE,(double *)g.coord,FEM_ELEM,areas,FEM_SPARSE);
00605     repeat_after_split((void *)&g);
00606     g.nelems = FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_ELEM);
00607     g.nnodes = FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_NODE);
00608     CkPrintf("[%d] Done with coarsening step: %d nodes, %d elements\n",
00609          myChunk,countValidEntities(g.validNode,g.nnodes),countValidEntities(g.validElem,g.nelems));
00610     delete [] loc;
00611     delete[] areas;
00612     // THIS IS THE COARSENED MESH SENT TO NetFEM
00613     if (1) { //Publish data to the net
00614       publishMeshToNetFEM(g,myChunk,2*t-1);
00615     }
00616 
00617     //prepare to refine
00618     loc=new vector2d[2*g.nnodes];
00619     for (i=0;i<g.nnodes;i++) {
00620       loc[i]=g.coord[i];//+g.d[i];
00621     }
00622     
00623     areas=new double[g.nelems];
00624     for (i=0;i<g.nelems;i++) {
00625       areas[i] = avgArea;
00626     }
00627     //refine one element at a time
00628     //int refIdx = (13  + 3*t)%g.nnodes;
00629     //areas[refIdx] = calcArea(g,refIdx)/1.5;
00630     
00631     CkPrintf("[%d] Starting refinement step: %d nodes, %d elements\n", myChunk,countValidEntities(g.validNode,g.nnodes),countValidEntities(g.validElem,g.nelems));
00632     FEM_REFINE2D_Split(FEM_Mesh_default_read(),FEM_NODE,(double *)loc,FEM_ELEM,areas,FEM_SPARSE);
00633     repeat_after_split((void *)&g);
00634     
00635     g.nelems = FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_ELEM);
00636     g.nnodes = FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_NODE);
00637     CkPrintf("[%d] Done with refinement step: %d nodes, %d elements\n",
00638          myChunk,countValidEntities(g.validNode,g.nnodes),countValidEntities(g.validElem,g.nelems));
00639     delete [] loc;
00640     delete[] areas;
00641     // THIS IS THE REFINED MESH SENT TO NetFEM
00642     if (1) { //Publish data to the net
00643       publishMeshToNetFEM(g,myChunk,2*t);
00644     }
00645   }
00646   if (CkMyPe()==0)
00647     CkPrintf("Driver finished\n");
00648 }
00649 
00650 int countValidEntities(int *validData,int total){
00651     int sum =0 ;
00652     for(int i=0;i<total;i++){
00653         sum += validData[i];
00654     }
00655     return sum;
00656 }
00657 
00658 
00659 void publishMeshToNetFEM(myGlobals &g,int myChunk,int t){
00660     NetFEM n=NetFEM_Begin(myChunk,t,2,NetFEM_WRITE);
00661     int count=0;
00662     double *vcoord = new double[2*g.nnodes];
00663     double *vnodeid = new double[g.nnodes];
00664     int *maptovalid = new int[g.nnodes];
00665     for(int i=0;i<g.nnodes;i++){
00666       if(g.validNode[i]){
00667     vcoord[2*count] = ((double *)g.coord)[2*i];
00668     vcoord[2*count+1] = ((double *)g.coord)[2*i+1];
00669     maptovalid[i] = count;
00670     printf("~~~~~~~ %d %d %.6lf %.6lf \n",count,i,vcoord[2*count],vcoord[2*count+1]);
00671     vnodeid[count] = i;
00672     count++;    
00673       }
00674     }
00675     NetFEM_Nodes(n,count,(double *)vcoord,"Position (m)");
00676     NetFEM_Scalar(n,vnodeid,1,"Node ID");
00677     /*    NetFEM_Vector(n,(double *)g.d,"Displacement (m)");
00678       NetFEM_Vector(n,(double *)g.v,"Velocity (m/s)");*/
00679     count=0;
00680     int *vconn = new int[3*g.nelems];
00681     double *vid = new double[3*g.nelems];
00682     for(int i=0;i<g.nelems;i++){
00683       if(g.validElem[i]){
00684     vconn[3*count] = maptovalid[g.conn[3*i]];
00685     vconn[3*count+1] = maptovalid[g.conn[3*i+1]];
00686     vconn[3*count+2] = maptovalid[g.conn[3*i+2]];
00687     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]);
00688     vid[count]=count;
00689     count++;    
00690       }
00691     }
00692     NetFEM_Elements(n,count,3,(int *)vconn,"Triangles");
00693     NetFEM_Scalar(n,vid,1,"Element ID");
00694     /*  NetFEM_Scalar(n,g.S22,1,"Y Stress (pure)");
00695     NetFEM_Scalar(n,g.S12,1,"Shear Stress (pure)");*/
00696     NetFEM_End(n);
00697     delete [] vcoord;
00698     delete [] vconn;
00699     delete [] maptovalid;
00700     delete [] vid;
00701     delete [] vnodeid;
00702 }
libs/ck-libs/TMRC2D/pgm-remesh.C
