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