-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathStiff_HCT_Sub_T2_SP.m
58 lines (45 loc) · 5.02 KB
/
Stiff_HCT_Sub_T2_SP.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
function [ ke ] = Stiff_HCT_Sub_T2_SP( edof,H,x,y,gauss,we,ke )
L1=gauss(1); L2=gauss(2); L3=gauss(3);
x1=x(1); x2=x(2); x3=x(3); y1=y(1); y2=y(2); y3=y(3);
x0=(x1+x2+x3)/3; y0=(y1+y2+y3)/3;
% A22=2*area_Sub_triangle_2
% o: The centroid of triangle element
% 3-node (1)
% o
% / \
% L3=0 / \ L2=0
% / Sub-T_2 \
% 1-node (2) o-------------o 0-node (3)
% x0=(x1+x2+x3)/3; y0=(y1+y2+y3)/3;
b10= y1-y0; b03=y0-y3; b31=y3-y1;
c10= x1-x0; c03=x0-x3; c31=x3-x1;
l10=sqrt(b10^2+c10^2); l03=sqrt(b03^2+c03^2); l31=sqrt(b31^2+c31^2);
% mu10=(l31^2-l03^2)/l10^2; mu03=(l10^2-l31^2)/l03^2;
mu31=(l03^2-l10^2)/l31^2;
A22=b10*(-c03)-(-c10)*b03; % A22=2*area_Sub_triangle_2
A24=2*A22;
Q=[ b10 b03 b31 ;...
-c10 -c03 -c31];
% Qui uoc , w_x=-dw/dy, w_y=dw/dx
D2L11 =[ 6*L1 - (6*L3*(b03*c10 - b10*c03))/(b03*c31 - b31*c03) - (6*L2*(b10*c31 - b31*c10))/(b03*c31 - b31*c03), (2*A22*L3*b03)/(b03*c31 - b31*c03) - (2*A22*L2*b31)/(b03*c31 - b31*c03), (2*A22*L2*c31)/(b03*c31 - b31*c03) - (2*A22*L3*c03)/(b03*c31 - b31*c03), 0, 0, 0, 0, 0, 0];
D2L22 =[ 0, 0, 0, 6*L2 + (6*L3*(b03*c10 - b10*c03))/(b10*c31 - b31*c10) - (6*L1*(b03*c31 - b31*c03))/(b10*c31 - b31*c10), (2*A22*L3*b10)/(b10*c31 - b31*c10) - (2*A22*L1*b31)/(b10*c31 - b31*c10), (2*A22*L1*c31)/(b10*c31 - b31*c10) - (2*A22*L3*c10)/(b10*c31 - b31*c10), 0, 0, 0];
D2L33 =[ 0, 0, 0, 0, 0, 0, 6*L3 - (6*L1*(b03*c31 - b31*c03))/(b03*c10 - b10*c03) + (6*L2*(b10*c31 - b31*c10))/(b03*c10 - b10*c03), (2*A22*L1*b03)/(b03*c10 - b10*c03) - (2*A22*L2*b10)/(b03*c10 - b10*c03), (2*A22*L2*c10)/(b03*c10 - b10*c03) - (2*A22*L1*c03)/(b03*c10 - b10*c03)];
D2L12 =[ (3*L3*(2*b03*c10 - 2*b10*c03 + b03*c31 - b31*c03 - 3*b10*c31 + 3*b31*c10 - b03*c31*mu31 + b31*c03*mu31 + b10*c31*mu31 - b31*c10*mu31))/(2*(b03*c31 - b31*c03)) - (6*L1*(b10*c31 - b31*c10))/(b03*c31 - b31*c03), (A24*L3*c31)/l31^2 - (A22*L3*(2*b03 + 3*b31 - b31*mu31))/(2*(b03*c31 - b31*c03)) - (2*A22*L1*b31)/(b03*c31 - b31*c03), (A22*L3*(2*c03 + 3*c31 - c31*mu31))/(2*(b03*c31 - b31*c03)) + (A24*L3*b31)/l31^2 + (2*A22*L1*c31)/(b03*c31 - b31*c03), - (3*L3*(2*b03*c10 - 2*b10*c03 + 3*b03*c31 - 3*b31*c03 - b10*c31 + b31*c10 + b03*c31*mu31 - b31*c03*mu31 - b10*c31*mu31 + b31*c10*mu31))/(2*(b10*c31 - b31*c10)) - (6*L2*(b03*c31 - b31*c03))/(b10*c31 - b31*c10), (A24*L3*c31)/l31^2 - (A22*L3*(2*b10 + 3*b31 + b31*mu31))/(2*(b10*c31 - b31*c10)) - (2*A22*L2*b31)/(b10*c31 - b31*c10), (A22*L3*(2*c10 + 3*c31 + c31*mu31))/(2*(b10*c31 - b31*c10)) + (A24*L3*b31)/l31^2 + (2*A22*L2*c31)/(b10*c31 - b31*c10), 0, 0, 0];
D2L13 =[ (3*L2*(2*b03*c10 - 2*b10*c03 + b03*c31 - b31*c03 - 3*b10*c31 + 3*b31*c10 - b03*c31*mu31 + b31*c03*mu31 + b10*c31*mu31 - b31*c10*mu31))/(2*(b03*c31 - b31*c03)) - (6*L1*(b03*c10 - b10*c03))/(b03*c31 - b31*c03), (A24*L2*c31)/l31^2 - (A22*L2*(2*b03 + 3*b31 - b31*mu31))/(2*(b03*c31 - b31*c03)) + (2*A22*L1*b03)/(b03*c31 - b31*c03), (A22*L2*(2*c03 + 3*c31 - c31*mu31))/(2*(b03*c31 - b31*c03)) + (A24*L2*b31)/l31^2 - (2*A22*L1*c03)/(b03*c31 - b31*c03), -(3*L2*(2*b03*c10 - 2*b10*c03 + 3*b03*c31 - 3*b31*c03 - b10*c31 + b31*c10 + b03*c31*mu31 - b31*c03*mu31 - b10*c31*mu31 + b31*c10*mu31))/(2*(b10*c31 - b31*c10)), (A24*L2*c31)/l31^2 - (A22*L2*(2*b10 + 3*b31 + b31*mu31))/(2*(b10*c31 - b31*c10)), (A22*L2*(2*c10 + 3*c31 + c31*mu31))/(2*(b10*c31 - b31*c10)) + (A24*L2*b31)/l31^2, -(6*L3*(b03*c31 - b31*c03))/(b03*c10 - b10*c03), (2*A22*L3*b03)/(b03*c10 - b10*c03), -(2*A22*L3*c03)/(b03*c10 - b10*c03)];
D2L23 =[ (3*L1*(2*b03*c10 - 2*b10*c03 + b03*c31 - b31*c03 - 3*b10*c31 + 3*b31*c10 - b03*c31*mu31 + b31*c03*mu31 + b10*c31*mu31 - b31*c10*mu31))/(2*(b03*c31 - b31*c03)), (A24*L1*c31)/l31^2 - (A22*L1*(2*b03 + 3*b31 - b31*mu31))/(2*(b03*c31 - b31*c03)), (A22*L1*(2*c03 + 3*c31 - c31*mu31))/(2*(b03*c31 - b31*c03)) + (A24*L1*b31)/l31^2, (6*L2*(b03*c10 - b10*c03))/(b10*c31 - b31*c10) - (3*L1*(2*b03*c10 - 2*b10*c03 + 3*b03*c31 - 3*b31*c03 - b10*c31 + b31*c10 + b03*c31*mu31 - b31*c03*mu31 - b10*c31*mu31 + b31*c10*mu31))/(2*(b10*c31 - b31*c10)), (A24*L1*c31)/l31^2 - (A22*L1*(2*b10 + 3*b31 + b31*mu31))/(2*(b10*c31 - b31*c10)) + (2*A22*L2*b10)/(b10*c31 - b31*c10), (A22*L1*(2*c10 + 3*c31 + c31*mu31))/(2*(b10*c31 - b31*c10)) + (A24*L1*b31)/l31^2 - (2*A22*L2*c10)/(b10*c31 - b31*c10), (6*L3*(b10*c31 - b31*c10))/(b03*c10 - b10*c03), -(2*A22*L3*b10)/(b03*c10 - b10*c03), (2*A22*L3*c10)/(b03*c10 - b10*c03)];
% Extract from formulation (4.58), FEM - Zienkiewicz vol.2
for i=1:edof
DNL=zeros(3,3);
DNL=[D2L11(i) D2L12(i) D2L13(i);...
D2L12(i) D2L22(i) D2L23(i);...
D2L13(i) D2L23(i) D2L33(i)];
% Second derivatives of shape function coresponding to x,y in (4.58), pp.132
DNxy=Q*DNL*Q';
B(1,i)=DNxy(1,1);
B(2,i)=DNxy(2,2);
B(3,i)=2*DNxy(1,2);
end
B=(1/A22^2)*B;
% Formula (27.32) in IFEM,chapter 27.
ke=ke+0.5*A22*we*B'*H*B;
end