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