Load_HCT_Sub_T2_SP.m

function [fe]=Load_HCT_Sub_T2_SP(Load,x,y)
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;

 
fe2 =[ -(2*b03*c10 - 2*b10*c03 - 5*b03*c31 + 5*b31*c03 + 7*b10*c31 - 7*b31*c10 + b03*c31*mu31 - b31*c03*mu31 - b10*c31*mu31 + b31*c10*mu31)/(80*(b03*c31 - b31*c03)), ((A22*b03*l31^2)/120 - (7*A22*b31*l31^2)/240 + (A22*b31*l31^2*mu31)/240)/(l31^2*(b03*c31 - b31*c03)) + (A24*c31)/(120*l31^2), (A24*b31)/(120*l31^2) - ((A22*c03*l31^2)/120 - (7*A22*c31*l31^2)/240 + (A22*c31*l31^2*mu31)/240)/(l31^2*(b03*c31 - b31*c03)), (2*b03*c10 - 2*b10*c03 - 7*b03*c31 + 7*b31*c03 + 5*b10*c31 - 5*b31*c10 - b03*c31*mu31 + b31*c03*mu31 + b10*c31*mu31 - b31*c10*mu31)/(80*(b10*c31 - b31*c10)), (A24*c31)/(120*l31^2) - ((7*A22*b31*l31^2)/240 - (A22*b10*l31^2)/120 + (A22*b31*l31^2*mu31)/240)/(l31^2*(b10*c31 - b31*c10)), ((7*A22*c31*l31^2)/240 - (A22*c10*l31^2)/120 + (A22*c31*l31^2*mu31)/240)/(l31^2*(b10*c31 - b31*c10)) + (A24*b31)/(120*l31^2), (b03*c10 - b10*c03 - b03*c31 + b31*c03 + b10*c31 - b31*c10)/(20*(b03*c10 - b10*c03)), (A22*(b03 - b10))/(60*(b03*c10 - b10*c03)), -(A22*(c03 - c10))/(60*(b03*c10 - b10*c03))];

fe=Load*A22*fe2';
end