%用四節點平面板單元求解平板中心受集中載荷問題
%------------------------------------
%  輸入控制參數
%------------------------------------

clear                    %清除workplace殘留數據
nel=4;                   % 單元數
nnel=4;                  % 每個單元的節點數
ndof=3;                  % 每個節點的自由度數
nnode=9;                 % 系統總節點數
sdof=nnode*ndof;         % 系統總自由度數  
edof=nnel*ndof;          % 每個單元的自由度數
emodule=30e6;            % 楊氏彈性模量
poisson=0.3;             % 泊松比
t=0.1;                   % 薄板厚度
nglxb=2; nglyb=2;        % 彎曲對應的2x2高斯拉格朗日積分 
nglb=nglxb*nglyb;        % 彎曲對應的每個單元的高斯積分點數
nglxs=1; nglys=1;        % 剪切對應的1x1高斯拉格朗日積分 
ngls=nglxs*nglys;        % 剪切對應的每個單元的高斯積分點數

%---------------------------------------------
%  輸入節點坐標值
%  gcoord(i,j) i：節點號　j：x,y值
%---------------------------------------------

gcoord=[0.0  0.0; 2.5  0.0; 5.0  0.0; 
        0.0  2.5; 2.5  2.5; 5.0  2.5;
        0.0  5.0; 2.5  5.0; 5.0  5.0];

%---------------------------------------------------------
%  每個單元對應的節點號（逆時針排列）
%  nodes(i,j) i:節點號 j:對應的單元號
%---------------------------------------------------------

nodes=[1 2 5 4; 2 3 6 5; 4 5 8 7; 5 6 9 8];

%-------------------------------------
% 輸入邊界條件
%-------------------------------------

bcdof=[1 2 3 4 6 7 9 11 12 16 20 21 23 25 26];  % 約束的自由度
bcval=zeros(1,15);        % 對應的值

%----------------------------------------------
% 初始化矩陣和矢量
%----------------------------------------------

ff=zeros(sdof,1);       % 載荷矢量
kk=zeros(sdof,sdof);    % 系統剛度矩陣
disp=zeros(sdof,1);     % 系統位移矢量
index=zeros(edof,1);    % 每個單元所包含的自由度
kinmtpb=zeros(3,edof);   % 彎曲幾何函數矩陣
matmtpb=zeros(3,3);      % 彎曲材料系數矩陣
kinmtps=zeros(2,edof);   % 剪切幾何函數矩陣
matmtps=zeros(2,2);      % 剪切材料系數矩陣

%----------------------------
% 載荷矢量
%----------------------------

ff(27)=10;              % 結點9所受的集中載荷

%-----------------------------------------------------------------
%  單元剛度矩陣計算及其組合
%-----------------------------------------------------------------
%
%  彎曲相關計算
%
[pointb,weightb]=swp2(nglxb,nglyb);     % 積分點和權系數
matmtpb=sbm(emodule,poisson)*t^3/12;  %彎曲材料系數
%
%  剪切相關計算
%
[points,weights]=swp2(nglxs,nglys);     % 積分點和權系數
shearm=0.5*emodule/(1.0+poisson);          % 剪切模量
shcof=5/6;                                 % 剪切修正因數 
matmtps=shearm*shcof*t*[1 0; 0 1];         % 剪切材料系數矩陣

for iel=1:nel           % 對所有單元數的循環

for i=1:nnel
nd(i)=nodes(iel,i);         % 當前單元對應的節點
xcoord(i)=gcoord(nd(i),1);  % 節點對應的x坐標值
ycoord(i)=gcoord(nd(i),2);  % 節點對應的y坐標值
end

k=zeros(edof,edof);         % 初始化單元剛度矩陣
kb=zeros(edof,edof);        % 初始化彎曲剛度矩陣
ks=zeros(edof,edof);        % 初始化剪切剛度矩陣

%------------------------------------------------------
%  彎曲相關計算
%------------------------------------------------------

for intx=1:nglxb
x=pointb(intx,1);                  % x軸積分點坐標
wtx=weightb(intx,1);               % 權系數
for inty=1:nglyb
y=pointb(inty,2);                  % y軸積分點坐標
wty=weightb(inty,2) ;              % 權系數

[shape,dhdr,dhds]=ssf(x,y);     % 計算形函數和對其相應的求導

jacob2=sjc(nnel,dhdr,dhds,xcoord,ycoord);  % 計算雅可比行列式

detjacob=det(jacob2);                 % 計算雅可比行列式的值
invjacob=inv(jacob2);                 % 求雅可比行列式的逆

[dhdx,dhdy]=sxy(nnel,dhdr,dhds,invjacob); % 計算ddhdr,dhds在迪卡爾坐標下的值

kinmtpb=sbB(nnel,dhdx,dhdy);          % 計算彎曲幾何函數矩陣

%--------------------------------------------
%  計算彎曲剛度矩陣
%--------------------------------------------

kb=kb+kinmtpb'*matmtpb*kinmtpb*wtx*wty*detjacob;

end
end                      % 結束彎曲剛度矩陣的計算

%------------------------------------------------------
% 剪切相關計算
%------------------------------------------------------

for intx=1:nglxs
x=points(intx,1);                  % x軸積分點坐標
wtx=weights(intx,1);               % 權系數
for inty=1:nglys
y=points(inty,2);                  % y軸積分點坐標
wty=weights(inty,2) ;              % 權系數

[shape,dhdr,dhds]=ssf(x,y);     % 計算形函數和對其相應的教學求導

jacob2=sjc(nnel,dhdr,dhds,xcoord,ycoord);  % 計算雅可比行列式

detjacob=det(jacob2);                 % 計算雅可比行列式的值
invjacob=inv(jacob2);                 % 求雅可比行列式的逆

[dhdx,dhdy]=sxy(nnel,dhdr,dhds,invjacob); % 計算dhdr,dhds在迪卡爾坐標下的值

kinmtps=ssB(nnel,dhdx,dhdy,shape);        % 計算剪切幾何函數矩陣

%----------------------------------------
%  計算剪切剛度矩陣
%----------------------------------------

ks=ks+kinmtps'*matmtps*kinmtps*wtx*wty*detjacob;   

end
end                      % 結束剪切剛度矩陣的計算

%--------------------------------
% 計算單元  剛度矩陣
%--------------------------------

k=kb+ks;

index=etsd(nd,nnel,ndof);% 單元對應的系統自由度號

kk=ask(kk,k,index);  % 合成系統剛度矩陣

end

%-----------------------------
%   加載邊界條件
%-----------------------------

[kk,ff]=dbc(kk,ff,bcdof,bcval);

%----------------------------
% 求解
%----------------------------

disp=kk\ff;   

num=1:1:sdof;
nodedisp=[num' disp]                  % 輸出節點位移
%----------------------------
% 后處理
%----------------------------
result=zeros(25,3);
displace(75)=0;
a=re(1,0,disp);
a(75)=0;
displace=displace+a;
a=re(10,6,disp);
a(75)=0;
displace=displace+a;
a=re(19,12,disp);
a(75)=0;
displace=displace+a;
for i=1:15;
    displace(45+i)=displace(15+i);
    displace(60+i)=displace(i);
end
[result]=dtxy(displace);
for i=1:25
    displacex(i,:)=result(:,1);
    displacex(i,:)=result(:,1);
    displacey(i,:)=result(:,2);
    displacez(i,:)=result(:,3);
end
[X,Y] = meshgrid(0:10/24:10);
Z=sqrt(displacex.^2+displacey.^2+displacez.^2);
surf(Z);