1 #include "pgm.h"
3 //Compute forces on constant-strain triangles:
4 void CST_NL(const vector2d *coor,const connRec *lm,vector2d *R_net,
5             const vector2d *d,const double *c,
6             int numnp,int numel,
7             double *S11o,double *S22o,double *S12o)
8 {
9   int n1,n2,n3,i;
10   double S11,S22,S12,u1,u2,u3,v1,v2,v3,x21,y21,x31,y31,x32,y32;
11   double aa,B1,B2,B3,B4,B5,B6,dudx,dvdy,dudy,dvdx;
12   double E11,E22,E12;
14   for (i=0;i<numel;i++) {
15     n1=lm[i][0];
16     n2=lm[i][1];
17     n3=lm[i][2];
18           u1 = d[n1].x;
19           u2 = d[n2].x;
20           u3 = d[n3].x;
21           v1 = d[n1].y;
22           v2 = d[n2].y;
23           v3 = d[n3].y;
25           x21 = coor[n2].x-coor[n1].x;
26           y21 = coor[n2].y-coor[n1].y;
27           x31 = coor[n3].x-coor[n1].x;
28           y31 = coor[n3].y-coor[n1].y;
29           x32 = coor[n3].x-coor[n2].x;
30           y32 = coor[n3].y-coor[n2].y;
32           aa = x21*y31-x31*y21;
33           B1 = -y32/aa;
34           B2 = x32/aa;
35           B3 = y31/aa;
36           B4 = -x31/aa;
37           B5 = -y21/aa;
38           B6 = x21/aa;
40           dudx = B1*u1 + B3*u2 + B5*u3;
41           dvdy = B2*v1 + B4*v2 + B6*v3;
42           dudy = B2*u1 + B4*u2 + B6*u3;
43           dvdx = B1*v1 + B3*v2 + B5*v3;
44           E11 = dudx + 0.5*(dudx*dudx + dvdx*dvdx);
45           E22 = dvdy + 0.5*(dvdy*dvdy + dudy*dudy);
46           E12 = dudy + dvdx + dudx*dudy + dvdy*dvdx;
48           // Calculate CST stresses
49           S11 = E11*c[0] + E22*c[1];
50           S22 = E11*c[1] + E22*c[2];
51           S12 = E12*c[3];
52           S11o[i]=S11;
53           S22o[i]=S22;
54           S12o[i]=S12;
56           // Calculate R-internal force vector
57           R_net[n1] -= aa*0.5*vector2d(
58                S11*B1*(1.0+dudx) +
59                S22*B2*dudy +
60                S12*(B2*(1.0+dudx) + B1*dudy)
61             ,
62                S11*B1*dvdx +
63                S22*B2*(1.0+dvdy) +
64                S12*(B1*(1.0+dvdy)+B2*dvdx)
65             );
66           R_net[n2] -= aa*0.5*vector2d(
67                S11*B3*(1.0+dudx) +
68                S22*B4*dudy +
69                S12*(B4*(1.0+dudx) + B3*dudy)
70             ,
71                S11*B3*dvdx +
72                S22*B4*(1.0+dvdy) +
73                S12*(B3*(1.0+dvdy)+B4*dvdx)
74             );
75           R_net[n3] -= aa*0.5*vector2d(
76                S11*B5*(1.0+dudx) +
77                S22*B6*dudy +
78                S12*(B6*(1.0+dudx) + B5*dudy)
79             ,
80                S11*B5*dvdx +
81                S22*B6*(1.0+dvdy) +
82                S12*(B5*(1.0+dvdy)+B6*dvdx)
83             );
84   }
85 }