Finished changes to use newer API with its FEM_Mesh_data() calls instead of those...
[charm.git] / examples / fem / simple2D / pgm.C
1 /*
2  Charm++ Finite-Element Framework Program:
3    Performs simple 2D structural simulation on Triangle-style inputs.
4    
5  init: 
6     read node input file
7     pass nodes into FEM framework
8     read element input file
9     pass elements into FEM framework
10  
11  driver:
12     extract mesh chunk from framework
13     calculate masses of nodes
14     timeloop
15       compute forces within my chunk
16       communicate
17       apply boundary conditions and integrate
18       pass data to NetFEM
19  
20  Among the hideous things about this program are:
21    -Hardcoded material properties, timestep, and boundary conditions
22    -Limited to 2D
23  
24  Converted from f90 Structural Materials program by 
25         Orion Sky Lawlor, 2001, olawlor@acm.org
26
27  Updated to new FEM interface by
28     Isaac Dooley, 2005
29  */
30
31 #include "pgm.h"
32
33
34 extern "C" void
35 init(void)
36 {
37   CkPrintf("init started\n");
38   double startTime=CmiWallTimer();
39   const char *eleName="xxx.1.ele";
40   const char *nodeName="xxx.1.node";
41   int nPts=0; //Number of nodes
42   vector2d *pts=0; //Node coordinates
43
44   CkPrintf("Reading node coordinates from %s\n",nodeName);
45   //Open and read the node coordinate file
46   {
47     char line[1024];
48     FILE *f=fopen(nodeName,"r");
49     if (f==NULL) die("Can't open node file!");
50     fgets(line,1024,f);
51     if (1!=sscanf(line,"%d",&nPts)) die("Can't read number of points!");
52     pts=new vector2d[nPts];
53     for (int i=0;i<nPts;i++) {
54       int ptNo;
55       if (NULL==fgets(line,1024,f)) die("Can't read node input line!");
56       if (3!=sscanf(line,"%d%lf%lf",&ptNo,&pts[i].x,&pts[i].y)) 
57         die("Can't parse node input line!");
58     }
59     fclose(f);
60   }
61  
62   int nEle=0;
63   connRec *ele=NULL;
64   CkPrintf("Reading elements from %s\n",eleName);
65   //Open and read the element connectivity file
66   {
67     char line[1024];
68     FILE *f=fopen(eleName,"r");
69     if (f==NULL) die("Can't open element file!");
70     fgets(line,1024,f);
71     if (1!=sscanf(line,"%d",&nEle)) die("Can't read number of elements!");
72     ele=new connRec[nEle];
73     for (int i=0;i<nEle;i++) {
74       int elNo;
75       if (NULL==fgets(line,1024,f)) die("Can't read element input line!");
76       if (4!=sscanf(line,"%d%d%d%d",&elNo,&ele[i][0],&ele[i][1],&ele[i][2])) 
77         die("Can't parse element input line!");  
78       ele[i][0]--; //Fortran to C indexing
79       ele[i][1]--; //Fortran to C indexing
80       ele[i][2]--; //Fortran to C indexing
81     }
82     fclose(f);
83   }
84  
85
86
87   int fem_mesh=FEM_Mesh_default_write(); // Tell framework we are writing to the mesh
88
89   CkPrintf("Passing node coords to framework\n");
90
91   /*   Old versions used FEM_Set_node() and FEM_Set_node_data()
92    *   New versions use the more flexible FEM_Set_Data()
93    */
94
95   FEM_Mesh_data(fem_mesh,        // Add nodes to the current mesh
96                 FEM_NODE,        // We are registering nodes
97                 FEM_DATA+0,      // Register the point locations which are normally 
98                                  // the first data elements for an FEM_NODE
99                 (double *)pts,   // The array of point locations
100                 0,               // 0 based indexing
101                 nPts,            // The number of points
102                 FEM_DOUBLE,      // Coordinates are doubles
103                 2);              // Points have dimension 2 (x,y)
104  
105
106   CkPrintf("Passing elements to framework\n");
107
108   /*   Old versions used FEM_Set_elem() and FEM_Set_elem_conn() 
109    *   New versions use the more flexible FEM_Set_Data()
110    */
111
112   FEM_Mesh_data(fem_mesh,      // Add nodes to the current mesh
113                 FEM_ELEM+0,      // We are registering elements with type 0
114                                  // The next type of element could be registered with FEM_ELEM+1
115                 FEM_CONN,        // Register the connectivity table for this
116                                  // data elements for this type of FEM entity
117                 (int *)ele,      // The array of point locations
118                 0,               // 0 based indexing
119                 nEle,            // The number of elements
120                 FEM_INDEX_0,     // We use zero based node numbering
121                 3);              // Elements have degree 3, since triangles are defined 
122                                  // by three nodes
123  
124    
125   delete[] ele;
126   delete[] pts;
127
128   CkPrintf("Finished with init (Reading took %.f s)\n",CmiWallTimer()-startTime);
129
130 }
131
132
133 // A driver() function 
134 // driver() is required in all FEM programs
135 extern "C" void
136 driver(void)
137 {
138   int nnodes,nelems,ignored;
139   int i, myId=FEM_My_partition();
140   myGlobals g;
141   FEM_Register(&g,(FEM_PupFn)pup_myGlobals);
142   
143   
144   int mesh=FEM_Mesh_default_read(); // Tell framework we are reading data from the mesh
145   
146   // Get node data
147   nnodes=FEM_Mesh_get_length(mesh,FEM_NODE); // Get number of nodes
148   g.coord=new vector2d[nnodes];
149   // Get node positions
150   FEM_Mesh_data(mesh, FEM_NODE, FEM_DATA+0, (double*)g.coord, 0, nnodes, FEM_DOUBLE, 2);  
151
152
153   // Get element data
154   nelems=FEM_Mesh_get_length(mesh,FEM_ELEM+0); // Get number of elements
155   g.nnodes=nnodes; g.nelems=nelems;
156   g.conn=new connRec[nelems];
157   g.S11=new double[nelems];
158   g.S22=new double[nelems];
159   g.S12=new double[nelems];
160   // Get connectivity for elements
161   FEM_Mesh_data(mesh, FEM_ELEM+0, FEM_CONN, (int *)g.conn, 0, nelems, FEM_INDEX_0, 3);  
162
163
164   //Initialize associated data
165   g.R_net=new vector2d[nnodes]; //Net force
166   g.d=new vector2d[nnodes];//Node displacement
167   g.v=new vector2d[nnodes];//Node velocity
168   g.a=new vector2d[nnodes];//Node acceleration
169   for (i=0;i<nnodes;i++)
170     g.R_net[i]=g.d[i]=g.v[i]=g.a[i]=vector2d(0.0);
171
172 //Apply a small initial perturbation to positions
173   for (i=0;i<nnodes;i++) {
174           const double max=1.0e-10/15.0; //Tiny perturbation
175           g.d[i].x+=max*(i&15);
176           g.d[i].y+=max*((i+5)&15);
177   }
178
179   int fid=FEM_Create_simple_field(FEM_DOUBLE,2);
180
181   //Timeloop
182   if (myId==0)
183     CkPrintf("Entering timeloop\n");
184   int tSteps=5000;
185   double startTime, totalStart;
186   startTime=totalStart=CkWallTimer();
187   for (int t=0;t<tSteps;t++) {
188     if (1) { //Structural mechanics
189
190     //Compute forces on nodes exerted by elements
191         CST_NL(g.coord,g.conn,g.R_net,g.d,matConst,nnodes,nelems,g.S11,g.S22,g.S12);
192
193     //Communicate net force on shared nodes
194         FEM_Update_field(fid,g.R_net);
195
196     //Advance node positions
197         advanceNodes(dt,nnodes,g.coord,g.R_net,g.a,g.v,g.d,0);
198     }
199
200     //Debugging/perf. output
201     double curTime=CkWallTimer();
202     double total=curTime-startTime;
203     startTime=curTime;
204     if (myId==0 && (t%64==0)) {
205             CkPrintf("%d %.6f sec for loop %d \n",CkNumPes(),total,t);
206             if (0) {
207               CkPrintf("    Triangle 0:\n");
208               for (int j=0;j<3;j++) {
209                     int n=g.conn[0][j];
210                     CkPrintf("    Node %d: coord=(%.4f,%.4f)  d=(%.4g,%.4g)\n",
211                              n,g.coord[n].x,g.coord[n].y,g.d[n].x,g.d[n].y);
212               }
213             }
214     }
215     /* perform migration-based load balancing */
216     if (t%1024==0)
217       FEM_Migrate();
218     
219     if (t%1024==0) { //Publish data to the net
220             NetFEM n=NetFEM_Begin(FEM_My_partition(),t,2,NetFEM_POINTAT);
221             
222             NetFEM_Nodes(n,nnodes,(double *)g.coord,"Position (m)");
223             NetFEM_Vector(n,(double *)g.d,"Displacement (m)");
224             NetFEM_Vector(n,(double *)g.v,"Velocity (m/s)");
225             
226             NetFEM_Elements(n,nelems,3,(int *)g.conn,"Triangles");
227                 NetFEM_Scalar(n,g.S11,1,"X Stress (pure)");
228                 NetFEM_Scalar(n,g.S22,1,"Y Stress (pure)");
229                 NetFEM_Scalar(n,g.S12,1,"Shear Stress (pure)");
230             
231             NetFEM_End(n);
232     }
233   }
234
235   if (myId==0) {
236     double elapsed=CkWallTimer()-totalStart;
237     CkPrintf("Driver finished: average %.6f s/step\n",elapsed/tSteps);
238   }
239 }
240
241
242 // A PUP function to allow for migration and load balancing of mesh partitions.
243 // The PUP function is not needed if no migration or load balancing is desired.
244 void pup_myGlobals(pup_er p,myGlobals *g) 
245 {
246   FEM_Print("-------- called pup routine -------");
247   pup_int(p,&g->nnodes);
248   pup_int(p,&g->nelems);
249   int nnodes=g->nnodes, nelems=g->nelems;
250   if (pup_isUnpacking(p)) {
251     g->coord=new vector2d[nnodes];
252     g->conn=new connRec[nelems];
253     g->R_net=new vector2d[nnodes]; //Net force
254     g->d=new vector2d[nnodes];//Node displacement
255     g->v=new vector2d[nnodes];//Node velocity
256     g->a=new vector2d[nnodes];
257         g->S11=new double[nelems];
258         g->S22=new double[nelems];
259         g->S12=new double[nelems];
260   }
261   pup_doubles(p,(double *)g->coord,2*nnodes);
262   pup_ints(p,(int *)g->conn,3*nelems);
263   pup_doubles(p,(double *)g->R_net,2*nnodes);
264   pup_doubles(p,(double *)g->d,2*nnodes);
265   pup_doubles(p,(double *)g->v,2*nnodes);
266   pup_doubles(p,(double *)g->a,2*nnodes);
267   pup_doubles(p,(double *)g->S11,nelems);
268   pup_doubles(p,(double *)g->S22,nelems);
269   pup_doubles(p,(double *)g->S12,nelems);
270   if (pup_isDeleting(p)) {
271     delete[] g->coord;
272     delete[] g->conn;
273     delete[] g->R_net;
274     delete[] g->d;
275     delete[] g->v;
276     delete[] g->a;
277         delete[] g->S11;
278         delete[] g->S22;
279         delete[] g->S12;
280   }
281 }
282
283
284
285 //Update node position, velocity, acceleration based on net force.
286 void advanceNodes(const double dt,int nnodes,const vector2d *coord,
287                   vector2d *R_net,vector2d *a,vector2d *v,vector2d *d,bool dampen)
288 {
289   const double nodeMass=1.0e-6; //Node mass, kilograms (HACK: hardcoded)
290   const double xm=1.0/nodeMass; //Inverse of node mass
291   const vector2d z(0,0);
292
293   const double shearForce=1.0e-6/(dt*dt*xm);
294
295   bool someNaNs=false;
296   int i;
297   for (i=0;i<nnodes;i++) {
298     vector2d R_n=R_net[i];
299 #if NANCHECK
300     if (((R_n.x-R_n.x)!=0)) {
301             CkPrintf("%d (%.4f,%.4f)   ",i,coord[i].x,coord[i].y);
302             someNaNs=true;
303     }
304 #endif
305     R_net[i]=z;
306 //Apply boundary conditions (HACK: hardcoded!)
307     if (1) {
308        if (coord[i].x<0.00001)
309                continue; //Left edge will NOT move
310        if (coord[i].y>0.02-0.00001)
311                R_n.x+=shearForce; //Top edge pushed hard to right
312     }
313 //Update displacement and velocity
314     vector2d aNew=R_n*xm;
315     v[i]+=(dt*0.5)*(aNew+a[i]);
316     d[i]+=dt*v[i]+(dt*dt*0.5)*aNew;
317     a[i]=aNew;   
318     //if (coord[i].y>0.02-0.00001) d[i].y=0.0; //Top edge in horizontal slot
319   }
320   if (dampen)
321     for (i=0;i<nnodes;i++)
322           v[i]*=0.9; //Dampen velocity slightly (prevents eventual blowup)
323   if (someNaNs) {
324           CkPrintf("Nodes all NaN!\n");
325           CkAbort("Node forces NaN!");
326   }
327 }