elemK2D8.m

function [K,elemB,D] = elemK2D8(EX,mu,h,x,y,isStress,reduce)
% 计算单元刚度矩阵
% 8节点四边形单元
% 全积分高斯点位置+0.774596669241483,-0.774596669241483,0
 %权系数为0.55555555555556,0.55555555555556,0.888888888888888889
% 减缩积分高斯点位置为+0.577350269189626,-0.577350269189626，权系数为1
% mnode == 8
% 4--7--3
% |     |
% 8     6
% |     |
% 1--5--2
K = zeros(16,16);
if reduce ==1
    elemB = zeros(12,16);
else
    elemB = 0;
end

if isStress == 1  % 平面应力
    D = EX/(1-mu^2)*[1,mu,0;mu,1,0;0,0,(1-mu)/2];
else              % 平面应变
    D = EX*(1-mu)/((1+mu)*(1-2*mu))*[1,mu/(1-mu),0;mu/(1-mu),1,0;0,0,(1-2*mu)/(2*(1-mu))];
end

if reduce == 0 % 采用完全积分
    nip = 3; % 单维方向的积分点数
    ks = [-0.774596669241483,0,0.774596669241483]; % 积分点
    w = [0.55555555555556,0.888888888888888889,0.55555555555556]; % 权系数
else
    nip = 2;
    ks = [-0.577350269189626,0.577350269189626]; % 积分点
    w = [1,1]; % 权系数
end

for m = 1:nip
    for n = 1:nip
        [J,B] = elemB2D8(x,y,ks(n),ks(m));
        if reduce == 1
            elemB(6*m+3*n-8:6*m+3*n-6,:) = B;
        end
        K = K + w(m)*w(n)*B'*D*B*det(J)*h;
    end
end