Cleaned up code and use new interface FEM_Mesh_data().
[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   FEM_Get_node(&nnodes,&ignored);
144   g.coord=new vector2d[nnodes];
145   FEM_Get_node_data((double *)g.coord);  
146
147   FEM_Get_elem(0,&nelems,&ignored,&ignored);
148   g.nnodes=nnodes; g.nelems=nelems;
149   g.conn=new connRec[nelems];
150   g.S11=new double[nelems];
151   g.S22=new double[nelems];
152   g.S12=new double[nelems];
153   FEM_Get_elem_conn(0,(int *)g.conn);
154
155   //Initialize associated data
156   g.R_net=new vector2d[nnodes]; //Net force
157   g.d=new vector2d[nnodes];//Node displacement
158   g.v=new vector2d[nnodes];//Node velocity
159   g.a=new vector2d[nnodes];//Node acceleration
160   for (i=0;i<nnodes;i++)
161     g.R_net[i]=g.d[i]=g.v[i]=g.a[i]=vector2d(0.0);
162
163 //Apply a small initial perturbation to positions
164   for (i=0;i<nnodes;i++) {
165           const double max=1.0e-10/15.0; //Tiny perturbation
166           g.d[i].x+=max*(i&15);
167           g.d[i].y+=max*((i+5)&15);
168   }
169
170   int fid=FEM_Create_simple_field(FEM_DOUBLE,2);
171
172   //Timeloop
173   if (myId==0)
174     CkPrintf("Entering timeloop\n");
175   int tSteps=5000;
176   double startTime, totalStart;
177   startTime=totalStart=CkWallTimer();
178   for (int t=0;t<tSteps;t++) {
179     if (1) { //Structural mechanics
180
181     //Compute forces on nodes exerted by elements
182         CST_NL(g.coord,g.conn,g.R_net,g.d,matConst,nnodes,nelems,g.S11,g.S22,g.S12);
183
184     //Communicate net force on shared nodes
185         FEM_Update_field(fid,g.R_net);
186
187     //Advance node positions
188         advanceNodes(dt,nnodes,g.coord,g.R_net,g.a,g.v,g.d,0);
189     }
190
191     //Debugging/perf. output
192     double curTime=CkWallTimer();
193     double total=curTime-startTime;
194     startTime=curTime;
195     if (myId==0 && (t%64==0)) {
196             CkPrintf("%d %.6f sec for loop %d \n",CkNumPes(),total,t);
197             if (0) {
198               CkPrintf("    Triangle 0:\n");
199               for (int j=0;j<3;j++) {
200                     int n=g.conn[0][j];
201                     CkPrintf("    Node %d: coord=(%.4f,%.4f)  d=(%.4g,%.4g)\n",
202                              n,g.coord[n].x,g.coord[n].y,g.d[n].x,g.d[n].y);
203               }
204             }
205     }
206     /* perform migration-based load balancing */
207     if (0)
208       FEM_Migrate();
209     
210     if (t==4608) { //Publish data to the net
211             NetFEM n=NetFEM_Begin(FEM_My_partition(),t,2,NetFEM_POINTAT);
212             
213             NetFEM_Nodes(n,nnodes,(double *)g.coord,"Position (m)");
214             NetFEM_Vector(n,(double *)g.d,"Displacement (m)");
215             NetFEM_Vector(n,(double *)g.v,"Velocity (m/s)");
216             
217             NetFEM_Elements(n,nelems,3,(int *)g.conn,"Triangles");
218                 NetFEM_Scalar(n,g.S11,1,"X Stress (pure)");
219                 NetFEM_Scalar(n,g.S22,1,"Y Stress (pure)");
220                 NetFEM_Scalar(n,g.S12,1,"Shear Stress (pure)");
221             
222             NetFEM_End(n);
223     }
224   }
225
226   if (myId==0) {
227     double elapsed=CkWallTimer()-totalStart;
228     CkPrintf("Driver finished: average %.6f s/step\n",elapsed/tSteps);
229   }
230 }
231
232
233 // A PUP function to allow for migration and load balancing of mesh partitions.
234 // The PUP function is not needed if no migration or load balancing is desired.
235 void pup_myGlobals(pup_er p,myGlobals *g) 
236 {
237   FEM_Print("-------- called pup routine -------");
238   pup_int(p,&g->nnodes);
239   pup_int(p,&g->nelems);
240   int nnodes=g->nnodes, nelems=g->nelems;
241   if (pup_isUnpacking(p)) {
242     g->coord=new vector2d[nnodes];
243     g->conn=new connRec[nelems];
244     g->R_net=new vector2d[nnodes]; //Net force
245     g->d=new vector2d[nnodes];//Node displacement
246     g->v=new vector2d[nnodes];//Node velocity
247     g->a=new vector2d[nnodes];
248         g->S11=new double[nelems];
249         g->S22=new double[nelems];
250         g->S12=new double[nelems];
251   }
252   pup_doubles(p,(double *)g->coord,2*nnodes);
253   pup_ints(p,(int *)g->conn,3*nelems);
254   pup_doubles(p,(double *)g->R_net,2*nnodes);
255   pup_doubles(p,(double *)g->d,2*nnodes);
256   pup_doubles(p,(double *)g->v,2*nnodes);
257   pup_doubles(p,(double *)g->a,2*nnodes);
258   pup_doubles(p,(double *)g->S11,nelems);
259   pup_doubles(p,(double *)g->S22,nelems);
260   pup_doubles(p,(double *)g->S12,nelems);
261   if (pup_isDeleting(p)) {
262     delete[] g->coord;
263     delete[] g->conn;
264     delete[] g->R_net;
265     delete[] g->d;
266     delete[] g->v;
267     delete[] g->a;
268         delete[] g->S11;
269         delete[] g->S22;
270         delete[] g->S12;
271   }
272 }
273
274
275
276 //Update node position, velocity, acceleration based on net force.
277 void advanceNodes(const double dt,int nnodes,const vector2d *coord,
278                   vector2d *R_net,vector2d *a,vector2d *v,vector2d *d,bool dampen)
279 {
280   const double nodeMass=1.0e-6; //Node mass, kilograms (HACK: hardcoded)
281   const double xm=1.0/nodeMass; //Inverse of node mass
282   const vector2d z(0,0);
283
284   const double shearForce=1.0e-6/(dt*dt*xm);
285
286   bool someNaNs=false;
287   int i;
288   for (i=0;i<nnodes;i++) {
289     vector2d R_n=R_net[i];
290 #if NANCHECK
291     if (((R_n.x-R_n.x)!=0)) {
292             CkPrintf("%d (%.4f,%.4f)   ",i,coord[i].x,coord[i].y);
293             someNaNs=true;
294     }
295 #endif
296     R_net[i]=z;
297 //Apply boundary conditions (HACK: hardcoded!)
298     if (1) {
299        if (coord[i].x<0.00001)
300                continue; //Left edge will NOT move
301        if (coord[i].y>0.02-0.00001)
302                R_n.x+=shearForce; //Top edge pushed hard to right
303     }
304 //Update displacement and velocity
305     vector2d aNew=R_n*xm;
306     v[i]+=(dt*0.5)*(aNew+a[i]);
307     d[i]+=dt*v[i]+(dt*dt*0.5)*aNew;
308     a[i]=aNew;   
309     //if (coord[i].y>0.02-0.00001) d[i].y=0.0; //Top edge in horizontal slot
310   }
311   if (dampen)
312     for (i=0;i<nnodes;i++)
313           v[i]*=0.9; //Dampen velocity slightly (prevents eventual blowup)
314   if (someNaNs) {
315           CkPrintf("Nodes all NaN!\n");
316           CkAbort("Node forces NaN!");
317   }
318 }