73f55ea3190beb90a11da81c5ecde41ffab9372a
[charm.git] / examples / fem / crack2D / lst_coh2.C
1 #include "crack.h"
2 #include <math.h>
3
4
5 //Update cohesive element traction itNo given Dn, Dt
6 static void 
7 updateTraction(double Dn, double Dt,Element *coh,CohMaterial *m,int itNo)
8 {
9   double Tn, Tt;
10   double c=coh->c.sidel[1],s=coh->c.sidel[2];
11   double x = 1.0 - sqrt(Dn*Dn+Dt*Dt);
12   
13   if (x <= 0.0) {
14     coh->c.Sthresh[itNo] = 0.0;
15   } else if (x <= coh->c.Sthresh[itNo])
16     coh->c.Sthresh[itNo] = x;
17   
18   if (Dn > 0.0)
19     Tn = coh->c.Sthresh[itNo]/(1.0-coh->c.Sthresh[itNo])*m->sigmax*Dn;
20   else
21     Tn = m->Sinit/(1.0-m->Sinit)*m->sigmax*Dn;
22   
23   Tt = coh->c.Sthresh[itNo]/(1.0-coh->c.Sthresh[itNo])*m->taumax*Dt;
24   coh->c.T[itNo].x = c*Tt - s*Tn;
25   coh->c.T[itNo].y = s*Tt + c*Tn;
26 }       
27
28 void
29 lst_coh2(GlobalData *gd) 
30 {
31   int idx;
32   for(idx=0; idx<gd->ne && gd->enums[idx]<gd->numCLST; idx++) {
33     Element *coh = &(gd->elements[idx]);
34     Node *n[6];
35     int k;
36     for(k=0;k<6;k++)
37       n[k] = &(gd->nodes[gd->conn[idx*6+k]]);
38
39     // local variables:
40     double deltn, deltt;  // the char length of current element
41     double length, c, s;
42     double Dx1, Dx2, Dx3;
43     double Dy1, Dy2, Dy3;
44     double Dn1, Dn2, Dn3;
45     double Dt1, Dt2, Dt3;
46     double Dn, Dt;
47     double x;
48     double Rx1,Rx2,Rx3;  // cohesive forces at nodes 1,2,5
49     double Ry1,Ry2,Ry3;  // cohesive forces at nodes 1,2,5
50   
51     // g1,g3: gauss quadrature points
52     // w1,w2,w3: weights; w1=w2, w3=w4
53     
54     CohMaterial *m = &(gd->cohm[coh->material]);
55       
56     deltn = (double)1.0 / m->deltan;
57     deltt = (double)1.0 / m->deltat;
58     length = coh->c.sidel[0] * (double)0.5;
59     c = coh->c.sidel[1];
60     s = coh->c.sidel[2];
61     
62     Dx1 = n[3]->disp.x - n[0]->disp.x;
63     Dy1 = n[3]->disp.y - n[0]->disp.y;
64     Dt1 =  (c*Dx1 + s*Dy1)*deltt;
65     Dn1 = (-s*Dx1 + c*Dy1)*deltn;
66     Dx2 = n[2]->disp.x - n[1]->disp.x;
67     Dy2 = n[2]->disp.y - n[1]->disp.y;
68     Dt2 =  (c*Dx2 + s*Dy2)*deltt;
69     Dn2 = (-s*Dx2 + c*Dy2)*deltn;
70     Dx3 = n[5]->disp.x - n[4]->disp.x;
71     Dy3 = n[5]->disp.y - n[4]->disp.y;
72     Dt3 =  (c*Dx3 + s*Dy3)*deltt;
73     Dn3 = (-s*Dx3 + c*Dy3)*deltn;
74       
75     // gauss point 1
76     Dt = g1*g1*0.5*(Dt1+Dt2-2.0*Dt3)+g1*0.5*(Dt2-Dt1)+Dt3;
77     Dn = g1*g1*0.5*(Dn1+Dn2-2.0*Dn3)+g1*0.5*(Dn2-Dn1)+Dn3;
78     updateTraction(Dn,Dt,coh,m,0);
79       
80     // gauss point 2
81     updateTraction(Dn3,Dt3,coh,m,1);
82   
83     // gauss point 3
84     Dt = g3*g3*0.5*(Dt1+Dt2-2.0*Dt3)+g3*0.5*(Dt2-Dt1)+Dt3;
85     Dn = g3*g3*0.5*(Dn1+Dn2-2.0*Dn3)+g3*0.5*(Dn2-Dn1)+Dn3;
86     updateTraction(Dn,Dt,coh,m,2);
87   
88     // cohesive forces
89     x = length*w1*g1*0.5;
90     Rx1 = (coh->c.T[0].x*(g1-1.0)+coh->c.T[2].x*(g1+1.0))*x;
91     Ry1 = (coh->c.T[0].y*(g1-1.0)+coh->c.T[2].y*(g1+1.0))*x;
92     Rx2 = (coh->c.T[0].x*(g1+1.0)+coh->c.T[2].x*(g1-1.0))*x;
93     Ry2 = (coh->c.T[0].y*(g1+1.0)+coh->c.T[2].y*(g1-1.0))*x;
94     Rx3 = ((coh->c.T[0].x+coh->c.T[2].x)*w1*(1.0-g1*g1)+coh->c.T[1].x*w2)*length;
95     Ry3 = ((coh->c.T[0].y+coh->c.T[2].y)*w1*(1.0-g1*g1)+coh->c.T[1].y*w2)*length;
96     n[0]->Rco.x += Rx1;
97     n[0]->Rco.y += Ry1;
98     n[3]->Rco.x -= Rx1;
99     n[3]->Rco.y -= Ry1;
100     n[1]->Rco.x += Rx2;
101     n[1]->Rco.y += Ry2;
102     n[2]->Rco.x -= Rx2;
103     n[2]->Rco.y -= Ry2;
104     n[4]->Rco.x += Rx3;
105     n[4]->Rco.y += Ry3;
106     n[5]->Rco.x -= Rx3;
107     n[5]->Rco.y -= Ry3;
108   }
109 }