Charm++: Charm Source Code Documentation

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