Fixed some more bugs. Now crack2D runs on multiple processors.
[charm.git] / examples / fem / crack2D / driver.C
1 #include <fstream.h>
2 #include <stddef.h>
3 #include "crack.h"
4
5 extern "C" void
6 init(void)
7 {
8 }
9
10 static void
11 initNodes(GlobalData *gd)
12 {
13   int idx;
14   for(idx=0; idx<gd->nn; idx++) {
15     Node *np = &(gd->nodes[idx]);
16     // Zero out node accumulators, update node positions
17     np->Rco.y = np->Rco.x = np->Rin.y =  np->Rin.x = 0.0;
18     np->disp.x += gd->delta2*np->accel.x + gd->delta*np->vel.x;
19     np->disp.y += gd->delta2*np->accel.y + gd->delta*np->vel.y;
20   }
21 }
22
23 static void 
24 updateNodes(GlobalData *gd, double prop, double slope)
25 {
26   for(int idx=0; idx<gd->nn; idx++) {
27     Node *np = &(gd->nodes[idx]);
28     if(!np->isbnd) {
29       double aX, aY;
30       aX = (np->Rco.x-np->Rin.x)*np->xM.x;
31       aY = (np->Rco.y-np->Rin.y)*np->xM.y;
32       np->vel.x += (gd->delta*(np->accel.x+aX)*(double) 0.5);
33       np->vel.y += (gd->delta*(np->accel.y+aY)*(double) 0.5);
34       np->accel.x = aX;
35       np->accel.y = aY;
36     } else {
37       double acc;
38       if (!(np->id1)) {
39         np->vel.x = (np->r.x)*prop;
40         np->accel.x = (np->r.x)*slope;
41       } else {
42         acc = (np->r.x*prop+ np->Rco.x - np->Rin.x)*np->xM.x;
43         np->vel.x += (gd->delta*(np->accel.x+acc)*0.5);
44         np->accel.x = acc;
45       }
46       if (!(np->id2)) {
47         np->vel.y = np->r.y*prop;
48         np->accel.y = np->r.y*slope;
49       } else {
50         acc = (np->r.y*prop+ np->Rco.y - np->Rin.y)*np->xM.y;
51         np->vel.y = (np->vel.y + gd->delta*(np->accel.y+acc)*0.5);
52         np->accel.y = acc;
53       }
54     }
55   }
56 }
57
58 extern "C" void
59 driver(int nn, int *nnums, int ne, int *enums, int npere, int *conn)
60 {
61   GlobalData *gd = new GlobalData;
62   Node *nodes = new Node[nn];
63   Element *elements = new Element[ne];
64   gd->myid = FEM_My_Partition();
65   gd->nn = nn;
66   gd->nnums = nnums;
67   gd->ne = ne;
68   gd->enums = enums;
69   gd->npere = npere;
70   gd->conn = conn;
71   gd->nodes = nodes;
72   gd->elements = elements;
73   readFile(gd);
74   int massfield = FEM_Create_Field(FEM_DOUBLE, 2, offsetof(Node, xM), 
75                                    sizeof(Node));
76   int rfield = FEM_Create_Field(FEM_DOUBLE, 4, offsetof(Node, Rin), 
77                                 sizeof(Node));
78   FEM_Update_Field(massfield, gd->nodes);
79   int i;
80   int kk = -1;
81   double prop, slope;
82   double stime, etime;
83   stime = CkTimer();
84   for(i=0;i<gd->nTime;i++)
85   {
86     if (gd->ts_proportion[kk+1] == i)
87     {
88       kk++;
89       prop = gd->proportion[kk];
90       slope = (gd->proportion[kk+1]-prop)/gd->delta;
91       slope /= (double) (gd->ts_proportion[kk+1]- gd->ts_proportion[kk]);
92     }
93     else
94     {
95       prop = (double)(i - gd->ts_proportion[kk])*
96                     slope*gd->delta+gd->proportion[kk];
97     }
98     initNodes(gd);
99     lst_NL(gd);
100     lst_coh2(gd);
101     updateNodes(gd, prop, slope);
102     FEM_Update_Field(rfield, gd->nodes);
103   }
104   etime = CkTimer();
105   if(gd->myid == 0)
106     CkPrintf("Time per iteration = %lf seconds\n", (etime-stime)/gd->nTime);
107 }
108
109 extern "C" void
110 finalize()
111 {
112 }