545b4595739a1234c1af65ab66e336a8a0024b88
[charm.git] / examples / fem / crack2D / lst_NL.C
1 /*Subroutine: 
2   lst_NL:  converted from Fortran on 9/7/99 by Orion Sky Lawlor
3            Computes each node's Rin.
4 */
5
6 #include "crack.h"
7
8 //BtoR: sum the "B" spatial derivative array into the
9 // "R" partial internal force array.
10
11 static void BtoR(const Coord *B/*6-array*/,
12     Coord *R/*6-array*/,Node **n/*6-array*/,
13     Element *v,const VolMaterial *vm,
14     int pkCount,double iaa)
15
16 {
17   double dudx,dvdy,dudy,dvdx;      //    ! partial derivatives of disp
18   double E11,E22,E12;              //            ! strains
19   int k;                         //Loop index
20
21 //C-----Calculate displacement gradient (H)
22   dudx=dvdy=dudy=dvdx=0.0;
23   for (k=0;k<6;k++) {
24     dudx+=B[k].x*n[k]->disp.x *iaa;
25     dvdy+=B[k].y*n[k]->disp.y *iaa;
26     dudy+=B[k].y*n[k]->disp.x *iaa;
27     dvdx+=B[k].x*n[k]->disp.y *iaa;
28   }
29
30 //C-----Calculate Lagrangian strain (E)
31   
32   E11 = dudx + 0.5*(dudx*dudx + dvdx*dvdx);
33   E22 = dvdy + 0.5*(dvdy*dvdy + dudy*dudy);
34   E12 = dudy + dvdx + dudx*dudy + dvdy*dvdx;
35   
36 //C-----Calculate the pkCount'th Piola-Kirchhoff stress (S)
37   
38   double s11c,s12c,s22c;
39   s11c=v->v.s11l[pkCount] = E11*vm->c[0] + E22*vm->c[1];
40   s22c=v->v.s22l[pkCount] = E11*vm->c[1] + E22*vm->c[2];
41   s12c=v->v.s12l[pkCount] = E12*vm->c[3];
42   
43 //Update R
44   for (k=0;k<6;k++) {
45     R[k].x+=s11c*B[k].x*(1.0+dudx)  +  s22c*B[k].y* dudy+
46       s12c*(B[k].y*(1.0+dudx) + B[k].x*dudy);
47     R[k].y+=s11c*B[k].x*dvdx        +  s22c*B[k].y*(1.0+dvdy)+
48       s12c*(B[k].y*     dvdx  + B[k].x*(1.0+dvdy));
49   }
50
51 }
52
53
54 void
55 lst_NL(GlobalData *gd)
56 {
57   int idx;
58   for(idx=gd->svol; idx<gd->evol; idx++) {
59     Element *v = &(gd->elements[idx]);
60     Node *n[6];
61     int k;
62     for(k=0;k<6;k++)
63       n[k] = &(gd->nodes[gd->conn[idx*6+k]]);
64
65     const double c5v3 = 1.66666666666667;
66     const double c1v3 = 0.33333333333333;
67     const double c2v3 = 0.66666666666667;
68     const double c1v6 = 0.166666666666667;
69   
70     double x21,y21,x31,y31,x32,y32;         //! coor(1,n2)-coor(1,n1)
71     double iaa;                       //! inverse of twice the element area
72       
73     Coord B[6];                           //! spacial derivatives
74     Coord R[6];                           //! partial internal force
75       
76     //Nodes 0, 1, and 2 are the outer corners;
77     //Nodes 4, 5, and 6 are on the midpoints of the sides.
78       
79     VolMaterial *vm = &(gd->volm[v->material]); //Current material
80       
81     for (k=0;k<6;k++) {
82         R[k].x=R[k].y=0.0;
83     }
84       
85     x21 = n[1]->pos.x-n[0]->pos.x;
86     y21 = n[1]->pos.y-n[0]->pos.y;
87     x31 = n[2]->pos.x-n[0]->pos.x;
88     y31 = n[2]->pos.y-n[0]->pos.y;
89     x32 = n[2]->pos.x-n[1]->pos.x;
90     y32 = n[2]->pos.y-n[1]->pos.y;
91     iaa = 1.0/(x21*y31-x31*y21);
92       
93     //Perform the three stages of Piola-Kirchoff stress
94     //First stage:
95     B[0].x = -y32*c5v3;             //                !   B1 = aa*dN1/dx
96     B[0].y = x32*c5v3;              //                !   B2 = aa*dN1/dy
97     B[1].x = -y31*c1v3;             //                !   B3 = aa*dN2/dx
98     B[1].y = x31*c1v3;              //                !   B4 = aa*dN2/dy
99     B[2].x = y21*c1v3;              //                !   B5 = aa*dN3/dx
100     B[2].y = -x21*c1v3;             //                !   B6 = aa*dN3/dy
101     B[3].x = (4.0*y31 - y32)*c2v3;  //                !   B7 = aa*dN4/dx
102     B[3].y = (x32 - 4.0*x31)*c2v3;  //                !   B8 = aa*dN4/dy
103     B[4].x = (y31 - y21)*c2v3;      //                !   B9 = aa*dN5/dx
104     B[4].y = (x21 - x31)*c2v3;      //                !   B10 = aa*dN5/dy
105     B[5].x = -(y32 + 4.0*y21)*c2v3; //                !   B11 = aa*dN6/dx
106     B[5].y = (x32 + 4.0*x21)*c2v3;  //                !   B12 = aa*dN6/dy
107           
108     BtoR(B,R,n,v,vm,0,iaa);
109           
110     //Second stage:
111     B[0].x = y32*c1v3;               //   !   B1 = aa*dN1/dx
112     B[0].y = -x32*c1v3;              //   !   B2 = aa*dN1/dy
113     B[1].x = y31*c5v3;               //   !   B3 = aa*dN2/dx
114     B[1].y = -x31*c5v3;              //   !   B4 = aa*dN2/dy
115     B[2].x = y21*c1v3;               //   !   B5 = aa*dN3/dx
116     B[2].y = -x21*c1v3;              //   !   B6 = aa*dN3/dy
117     B[3].x = (y31 - 4.0*y32)*c2v3;   //   !   B7 = aa*dN4/dx
118     B[3].y = (4.0*x32 - x31)*c2v3;   //   !   B8 = aa*dN4/dy
119     B[4].x = (y31 - 4.0*y21)*c2v3;   //   !   B9 = aa*dN5/dx
120     B[4].y = (4.0*x21 - x31)*c2v3;   //   !   B10 = aa*dN5/dy
121     B[5].x = -(y32 + y21)*c2v3;      //   !   B11 = aa*dN6/dx
122     B[5].y = (x32 + x21)*c2v3;       //   !   B12 = aa*dN6/dy
123     
124     BtoR(B,R,n,v,vm,1,iaa);
125     
126     //Third stage:
127     B[0].x = y32*c1v3;               //  !   B1 = aa*dN1/dx
128     B[0].y = -x32*c1v3;              //  !   B2 = aa*dN1/dy
129     B[1].x = -y31*c1v3;              //  !   B3 = aa*dN2/dx
130     B[1].y = x31*c1v3;               //  !   B4 = aa*dN2/dy
131     B[2].x = -y21*c5v3;              //  !   B5 = aa*dN3/dx
132     B[2].y = x21*c5v3;               //  !   B6 = aa*dN3/dy
133     B[3].x = (y31 - y32)*c2v3;       //  !   B7 = aa*dN4/dx
134     B[3].y = (x32 - x31)*c2v3;       //  !   B8 = aa*dN4/dy
135     B[4].x = (4.0*y31 - y21)*c2v3;   //  !   B9 = aa*dN5/dx
136     B[4].y = (x21 - 4.0*x31)*c2v3;   //  !   B10 = aa*dN5/dy
137     B[5].x = -(4.0*y32 + y21)*c2v3;  //  !   B11 = aa*dN6/dx
138     B[5].y = (4.0*x32 + x21)*c2v3;   //  !   B12 = aa*dN6/dy
139     
140     BtoR(B,R,n,v,vm,2,iaa);
141     
142     //Now we've calculated R.  Update each node`s Rin with the new values.
143     for (k=0;k<6;k++) {
144       n[k]->Rin.x+=c1v6*R[k].x;
145       n[k]->Rin.y+=c1v6*R[k].y;
146     }
147   }
148 }