#include <unistd.h> for sleep
[charm.git] / examples / fem / matrix / pgm.C
1 /*
2 Test out iterative (matrix-based) FEM interface.
3
4 Orion Sky Lawlor, olawlor@acm.org, 1/27/2003 
5 */
6 #include <stdlib.h>
7 #include <stdio.h>
8 #include <math.h>
9 #include <unistd.h>
10 #include "charm++.h"
11 #include "tcharmc.h"
12 #include "fem.h"
13 #include "netfem.h"
14 #include "ifemc.h"
15
16 #if CMK_CC_PGCC
17 #undef offsetof
18 #define offsetof(TYPE, MEMBER) ((size_t) &((TYPE *)0)->MEMBER)
19 #endif
20
21 //Number of time steps to simulate
22 int tsteps=10;
23 int dim=10;//Elements per side of the FEM mesh
24 const int np=4; //Nodes per element for a quad
25 const int width=3; //Number of dimensions to solve for (1, 2, or 3)
26
27 //Index of a node at x,y
28 int nodeDex(int x,int y) {return x+y*(dim+1);}
29 //Index of element connecting x,y to x+1,y+1
30 int elemDex(int x,int y) {return x+y*dim;}
31
32 // An imposed displacment boundary condition
33 struct boundaryCondition {
34         int node; //node this BC applies to
35         int dof; //Coordinate axis this BC applies to
36         double disp; //Imposed displacement value
37 };
38 int dof(const boundaryCondition &b) {
39   return b.node*width+b.dof;
40 }
41
42 // Save or restore this list of boundary conditions to/from the FEM framework
43 void fem_boundaryConditions(int mesh,int entity,int n,boundaryCondition *bc) 
44 {
45         int distance=sizeof(boundaryCondition);
46         FEM_Mesh_data_offset(mesh,entity,FEM_CONN, bc, 0,n, 
47                 IDXL_INDEX_0, 1,offsetof(boundaryCondition,node),distance,0);
48         
49         FEM_Mesh_data_offset(mesh,entity,FEM_DATA+0, bc, 0,n, 
50                 IDXL_INT, 1,offsetof(boundaryCondition,dof),distance,0);
51         
52         FEM_Mesh_data_offset(mesh,entity,FEM_DATA+1, bc, 0,n, 
53                 IDXL_DOUBLE, 1, offsetof(boundaryCondition,disp),distance,0);
54 }
55
56 extern "C" void
57 init(void)
58 {
59   CkPrintf("init called\n");
60   int argc=CkGetArgc();
61   char **argv=CkGetArgv();
62   if (argc>1) dim=atoi(argv[1]);
63   int nelems=dim*dim, nnodes=(dim+1)*(dim+1);
64   CkPrintf("Generating %d elements, %d node serial mesh\n",nelems,nnodes);
65   
66   //Describe the nodes and elements
67   int mesh=FEM_Mesh_default_write();
68   
69   int x,y,e;
70   
71 //Create node coordinates and initial conditions
72   double *nodes=new double[3*nnodes];
73   double domainX=2.0; //Size of domain, meters
74   double domainY=2.0;
75   for(y=0;y<dim+1;y++) 
76   for (x=0;x<dim+1;x++) {
77     nodes[nodeDex(x,y)*width+0]=domainX/float(dim)*x;
78     if (width>1) 
79       nodes[nodeDex(x,y)*width+1]=domainY/float(dim)*y;
80     if (width>2)
81       nodes[nodeDex(x,y)*width+2]=0.1;
82   }
83   FEM_Mesh_data(mesh,FEM_NODE,FEM_DATA, nodes,
84         0,nnodes, FEM_DOUBLE,width);
85
86   int i,j,n,c;
87   int length=width*nnodes;
88   
89   // Prepare net forces (boundary conditions):
90   double *netforce=new double[length];
91   double domainDensity=0.7*1.0e3; //Domain's density, Kg/m^3
92   double domainThickness=0.02; // Domain thickness, meters
93   double domainVolume=domainThickness*domainX*domainY;
94   double domainMass=domainDensity*domainVolume; //Entire domain's mass, Kg
95   double nodeMass=domainMass/((dim)*(dim)); //Node's mass, Kg
96   for (n=0;n<nnodes;n++) {
97     for (c=0;c<width;c++) {
98       double f=0;
99       if (c==width-1) f=9.8*nodeMass; //Weight of node, in Newtons
100       netforce[n*width+c]=f;
101     }
102   }
103   FEM_Mesh_data(mesh,FEM_NODE,FEM_DATA+1, netforce,
104         0,nnodes, FEM_DOUBLE,width);
105   
106   // Prepare imposed displacements (boundary conditions)
107   int nBC=2*width;
108   boundaryCondition *bc=new boundaryCondition[nBC];
109   for (c=0;c<width;c++) {
110     bc[0*width+c].node=0; //Left end: all displacments==0.1
111     bc[0*width+c].dof=c;
112     bc[0*width+c].disp=0.1;
113     bc[1*width+c].node=nnodes-1; //Right end: all == 0.2
114     bc[1*width+c].dof=c;
115     bc[1*width+c].disp=0.2;
116   }
117   fem_boundaryConditions(mesh,FEM_SPARSE+0,nBC,bc);
118
119 //Create the connectivity array
120   int *conn=new int[nelems*np];
121   for (y=0;y<dim;y++) 
122   for (x=0;x<dim;x++) {
123            e=elemDex(x,y);
124            conn[e*np+0]=nodeDex(x  ,y  );
125            conn[e*np+1]=nodeDex(x+1,y  );
126            conn[e*np+2]=nodeDex(x+1,y+1);
127            conn[e*np+3]=nodeDex(x  ,y+1);
128   }
129   FEM_Mesh_conn(mesh,FEM_ELEM+0, conn,
130         0,nelems, np);
131   delete[] conn;
132
133 }
134
135 extern "C"
136 int compareInts(const void *a,const void *b) {
137   return *(int *)a - *(int *)b;
138 }
139
140 /// User-defined class: contains my part of the FEM problem
141 class myMesh {
142   int myId;
143   int mesh;
144
145   int nnodes;
146   double *coord; //Node coordinates
147   double *netforce; //(knowns) Node net forces (plus forces from fixed displacements)
148   double *disp; //(unknowns) Node displacements 
149   
150   int nBC; //Number of fixed displacement boundary conditions
151   boundaryCondition *bc;
152   
153   int nelems;
154   int *conn; //Element->node mapping
155   
156   int mvps; //counter for matrix-vector-products
157   
158   // Return the undisplaced distance between these two nodes
159   double dist(int n1,int n2) {
160     double sum=0;
161     for (int i=0;i<width;i++) {
162       double del=coord[n1*width+i]-coord[n2*width+i];
163       sum+=del*del;
164     }
165     return sqrt(sum);
166   }
167 public:
168   myMesh(int mesh_)
169         :mesh(mesh_)
170   {
171     myId = FEM_My_partition();
172     nnodes=FEM_Mesh_get_length(mesh,FEM_NODE);
173     coord=new double[width*nnodes];
174     FEM_Mesh_data(mesh,FEM_NODE,FEM_DATA, coord,
175         0,nnodes, FEM_DOUBLE,width);
176     
177     netforce=new double[width*nnodes];
178     FEM_Mesh_data(mesh,FEM_NODE,FEM_DATA+1, netforce,
179         0,nnodes, FEM_DOUBLE,width);
180     disp=new double[width*nnodes];
181     
182     nelems=FEM_Mesh_get_length(mesh,FEM_ELEM);
183     conn=new int[np*nelems];
184     FEM_Mesh_conn(mesh,FEM_ELEM+0, conn,
185         0,nelems, np);
186
187     nBC=FEM_Mesh_get_length(mesh,FEM_SPARSE+0);
188     bc=new boundaryCondition[nBC];
189     fem_boundaryConditions(mesh,FEM_SPARSE+0,nBC,bc);
190     
191     mvps=0;
192   }
193   ~myMesh() {
194     delete[] coord;
195     delete[] netforce;
196     delete[] disp;
197     delete[] bc;
198     delete[] conn;
199   }
200   
201   // Solve our system using this solver
202   void solve(ILSI_Solver s);
203   
204   // Compute the forces we exert under these displacements.
205   //  If skipBoundaries is true, it's like zeroing out the rows and 
206   //  columns corresponding to boundary nodes.
207   void applyElementForces(const double *disp, double *force) 
208   {
209      if ((myId==0) && (mvps++)%32==0) //FEM_My_partition()==0) 
210         CkPrintf("MatrixVectorProduct (d=%g,%g,%g)\n",disp[0],disp[1],disp[2]);
211      
212      //Zero out node forces
213      int i,length=width*nnodes; //i loops over entries of our vectors
214      for (i=0;i<length;i++) force[i]=0;
215      
216      /**
217       * Add in (stupid) forces from local elements:
218       *    n1 - A - n2
219       *    |        |
220       *    D        B
221       *    |        |
222       *    n4 - C - n3
223       */
224      for (int e=0;e<nelems;e++) {
225        int n1=conn[np*e+0];
226        int n2=conn[np*e+1];
227        int n3=conn[np*e+2];
228        int n4=conn[np*e+3];
229        const double k=1.0e2; //Spring stiffness, N/m^2 (?)
230        double kA=k*(1.0/dist(n1,n2)); //Spring constant, N/m
231        double kB=k*(1.0/dist(n2,n3)); //Spring constant, N/m
232        double kC=k*(1.0/dist(n3,n4)); //Spring constant, N/m
233        double kD=k*(1.0/dist(n4,n1)); //Spring constant, N/m
234        for (int c=0;c<width;c++) {
235          double f;
236          f=-kA*(disp[n1*width+c]-disp[n2*width+c]);
237          force[n1*width+c]+=f; force[n2*width+c]-=f;
238          
239          f=-kB*(disp[n2*width+c]-disp[n3*width+c]);
240          force[n2*width+c]+=f; force[n3*width+c]-=f;
241          
242          f=-kC*(disp[n3*width+c]-disp[n4*width+c]);
243          force[n3*width+c]+=f; force[n4*width+c]-=f;
244          
245          f=-kD*(disp[n4*width+c]-disp[n1*width+c]);
246          force[n4*width+c]+=f; force[n1*width+c]-=f;
247        }
248      }
249      
250      // Communicate forces from remote elements
251      IDXL_Layout_t layout=IDXL_Layout_create(IDXL_DOUBLE,width);
252      FEM_Update_field(layout,force);
253      IDXL_Layout_destroy(layout);
254   }
255   
256   void netfem(int ts) {
257     NetFEM f=NetFEM_Begin(FEM_My_partition(),ts,width,NetFEM_POINTAT);
258       NetFEM_Nodes(f,nnodes,coord,"Node Locs");
259         NetFEM_Vector(f,disp,"Displacement");
260         NetFEM_Vector(f,netforce,"Net force");
261       NetFEM_Elements(f,nelems,np,conn,"Elements");
262     NetFEM_End(f);
263   }
264 };
265
266 extern "C" 
267 void mesh_matrix_product(void *ptr,
268         int length,int width,const double *src, double *dest)
269 {
270         myMesh *m=(myMesh *)ptr;
271         m->applyElementForces(src,dest);
272 }
273
274
275 // Solve our system using this solver
276 void myMesh::solve(ILSI_Solver s) {
277 //Split up boundary conditions into DOF and value arrays:
278   int *bcDOF=new int[nBC];
279   double *bcValue=new double[nBC];
280   int i;
281   for (i=0;i<nBC;i++) {
282     bcDOF[i]=dof(bc[i]);
283     bcValue[i]=bc[i].disp;
284   }
285   
286 //Pick a reasonable initial guess
287   for (i=0;i<width*nnodes;i++) disp[i]=0;
288
289 //Solve the system
290   ILSI_Param param;
291   ILSI_Param_new(&param);
292   param.maxResidual=1.0e-6;
293   // param.maxIterations=10000;
294   if (myId==0) CkPrintf("Solving...\n");
295   IFEM_Solve_shared_bc(ILSI_CG_Solver,&param,
296         mesh, FEM_NODE, nnodes,width,
297         nBC,bcDOF,bcValue,
298         mesh_matrix_product,this,netforce,disp);
299   
300   delete[] bcDOF; delete[] bcValue;
301   
302 //Compute the actual forces applied to each node:
303   applyElementForces(disp,netforce);
304   
305   if (myId==0) {
306     CkPrintf("Solved-- %d iterations, %g residual\n",
307         (int)param.iterations,param.residual);
308     if (nnodes<10)
309       for (int n=0;n<nnodes;n++) {
310         int c;
311         printf("  node %d  disp= ",n);
312         for (c=0;c<width;c++) printf(" %g",disp[width*n+c]);
313         printf("  ");
314         printf("  force= ",n);
315         for (c=0;c<width;c++) printf(" %g",netforce[width*n+c]);
316         printf("\n");
317       }
318   }
319 }
320
321 extern "C" void
322 driver(void)
323 {
324   FEM_Print("Starting driver...");
325   int myId = FEM_My_partition();
326
327   int fem_mesh=FEM_Mesh_default_read();
328   myMesh mesh(fem_mesh);
329   
330   // Solve a little problem on the mesh:
331   mesh.solve(ILSI_CG_Solver);
332
333   if (0) { //Loop for netfem access
334     if (myId==0)
335       FEM_Print("Waiting for NetFEM client to connect (hit ctrl-c to exit)");
336     int ts=0;
337     while(1) {
338       sleep(1);
339       mesh.netfem(ts);
340       ts++;
341       FEM_Barrier();
342     }
343   }
344 }
345
346 extern "C" void
347 mesh_updated(int param)
348 {
349   CkPrintf("mesh_updated(%d) called.\n",param);
350 }