Load_HCT_Sub_T3_SP.m

function [fe]=Load_HCT_Sub_T3_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;
% o: The centroid of triangle element

%                                          o-node (3)
%                                            o
%                                           / \

%								L2=0      /		\   L1=0

%										/ Sub-T_1 \

%                         1-node  (1) o-------------o 3-node (2)
%                                            L3=0 
% x0=(x1+x2+x3)/3;   y0=(y1+y2+y3)/3;
b20= y2-y0;                 b01=y0-y1;                  b12=y1-y2;
c20= x2-x0;                 c01=x0-x1;                  c12=x1-x2;
l20=sqrt(b20^2+c20^2);      l01=sqrt(b01^2+c01^2);      l12=sqrt(b12^2+c12^2);
% mu20=(l12^2-l01^2)/l20^2;   mu01=(l20^2-l12^2)/l01^2;   
mu12=(l01^2-l20^2)/l12^2;
A32=b20*(-c01)-(-c20)*b01;      % A22=2*area_Sub_triangle_3
A34=2*A32;
 
fe3 =[ (5*b01*c12 - 5*b12*c01 - 2*b01*c20 + 2*b20*c01 + 7*b12*c20 - 7*b20*c12 - b01*c12*mu12 + b12*c01*mu12 - b12*c20*mu12 + b20*c12*mu12)/(80*(b01*c12 - b12*c01)), ((A32*b01*l12^2)/120 - (7*A32*b12*l12^2)/240 + (A32*b12*l12^2*mu12)/240)/(l12^2*(b01*c12 - b12*c01)) + (A34*c12)/(120*l12^2), (A34*b12)/(120*l12^2) - ((A32*c01*l12^2)/120 - (7*A32*c12*l12^2)/240 + (A32*c12*l12^2*mu12)/240)/(l12^2*(b01*c12 - b12*c01)), (7*b01*c12 - 7*b12*c01 - 2*b01*c20 + 2*b20*c01 + 5*b12*c20 - 5*b20*c12 + b01*c12*mu12 - b12*c01*mu12 + b12*c20*mu12 - b20*c12*mu12)/(80*(b12*c20 - b20*c12)), ((7*A32*b12*l12^2)/240 - (A32*b20*l12^2)/120 + (A32*b12*l12^2*mu12)/240)/(l12^2*(b12*c20 - b20*c12)) + (A34*c12)/(120*l12^2), (A34*b12)/(120*l12^2) - ((7*A32*c12*l12^2)/240 - (A32*c20*l12^2)/120 + (A32*c12*l12^2*mu12)/240)/(l12^2*(b12*c20 - b20*c12)), -(b01*c12 - b12*c01 - b01*c20 + b20*c01 + b12*c20 - b20*c12)/(20*(b01*c20 - b20*c01)), (A32*(b01 - b20))/(60*(b01*c20 - b20*c01)), -(A32*(c01 - c20))/(60*(b01*c20 - b20*c01))];
 
fe=Load*A32*fe3';
end