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