clear; 
% DEFINE BOUNDARIES/PARAMETERS 
Lb=1;
D=5;
Q=0;
tm=3e7;bsb=0.167;
fc=20000;
ft=2000;
%平面應變問題
Dmat0=(tm/(1-2*bsb)/(1+bsb))*[(1-bsb) bsb 0;bsb (1-bsb) 0;0 0 (1-2*bsb)/2];




% SET UP NODAL COORDINATES
ndivl=8;
ndivw=20;
[x,conn,numcell,numnod] = mesh2(Lb,D,ndivl,ndivw);
numnod=length(x);

% SET UP QUADRATURE CELLS
ndivlq = 8;ndivwq = 20;
[xc,conn,numcell,numq] = mesh2(Lb,D,ndivlq,ndivwq);
% DETERMINE DOMAINS OF INFLUENCE - UNIFORM NODAL SPACING
dmax=3.5;
xspac = Lb/ndivl;
yspac = D/ndivw;
dm(1,1:numnod)=dmax*xspac*ones(1,numnod);
dm(2,1:numnod)=dmax*yspac*ones(1,numnod);

% SET UP QUADRATURE 求積分 CELLS
ndivlq = 8;ndivwq = 20;
[xc,conn,numcell,numq] = mesh2(Lb,D,ndivlq,ndivwq);
stress(1:3,:)=0;
stress1(1:3,:)=0;
for ix=1:10
    Q=Q+10;
% SET UP GAUSS POINTS, WEIGHTS, AND JACOBIAN FOR EACH CELL
quado = 4;
[gauss] = gauss2(quado);
numq2 = numcell*quado^2;
gs = zeros(4,numq2);
[gs] = egauss(xc,conn,gauss,numcell);
ngs=length(gs);
stress(1:3,:)=stress;
gs(5,:)=stress(1,:);
gs(6,:)=stress(2,:);
gs(7,:)=stress(3,:);

% LOOP OVER GAUSS POINTS TO ASSEMBLE DISCRETE EQUATIONS
k = sparse(numnod*2,numnod*2);
for gg=gs
    gpos=gg(1:2,:);
    weight=gg(3);
    jac=gg(4);
    stressx(1,1)=gs(5);
stressx(2,1)=gs(6);
stressx(3,1)=gs(7);
[tm,bsb]=fxx(stressx,tm,bsb,fc,ft);
% PLANE STRESS DMATRIX 
Dmat =(tm/(1-2*bsb)/(1+bsb))*[(1-bsb) bsb 0;bsb (1-bsb) 0;0 0 (1-2*bsb)/2];

% DETERMINE NODES IN NEIGHBORHOOD OF GAUSS POINT
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
Bmat=zeros(3,2*L);
for j=1:L
Bmat(1:3,(2*j-1):2*j) = [dphix(j) 0;0 dphiy(j);dphiy(j) dphix(j)];
end
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
end
k(en,en) = k(en,en)+sparse((weight*jac)*Bmat'*Dmat*Bmat);
end

f=21:21:21*9;
for i=1:length(f)
    qn(2*i-1)=2*f(i)-1;
    qn(2*i)=2*f(i);
    k(qn,qn)=k(qn,qn)*1e005;
end


% DETERMINE NODES ON BOUNDARY, SET UP BC'S
ind1 = 0;ind2 = 0;ind3=0;ind4=0;
for j=1:numnod
if(x(1,j)==-Lb/2)
ind1=ind1+1;
nnu(1,ind1) = x(1,j);
nnu(2,ind1) = x(2,j);
no1(ind1)=j;
end
if(x(1,j)==Lb/2)
ind2=ind2+1;
nt(1,ind2) = x(1,j);
nt(2,ind2) = x(2,j);
no2(ind2)=j;
end
if(x(2,j)==D/2)
ind3=ind3+1;
nt2(1,ind3) = x(1,j);
nt2(2,ind3) = x(2,j);
no3(ind3)=j;
end
if(x(2,j)==-D/2)
ind4=ind4+1;
nt1(1,ind4) = x(1,j);
nt1(2,ind4) = x(2,j);
no4(ind4)=j;
end
end
lthu = length(nnu);
ltht = length(nt);
ltht2=length(nt2);
ltht1=length(nt1);
f = zeros(1,numnod*2);


%計算zuo節點荷載
qn=Q/D*yspac;
qs=1/6*qn*yspac;
qx=2*qs;
qk=0.5*qn*yspac;

degao=zeros(3,lthu);
for i=1:lthu
    degao(1,i)=qs;
end
for i=1:lthu
    degao(2,i)=qx;
end
for i=1:lthu
    degao(3,i)=qk;
end



ind=1;
suan=0;
jisuan=zeros(4,lthu);
for i=2:lthu-1
    ind=ind+1;
    suan=suan+degao(3,i);
    jisuan(1,ind)=suan;
end
ind=1;
for i=2:lthu-1
    ind=ind+1;
    jisuan(2,ind)= jisuan(1,ind)+ jisuan(1,ind-1);
end
for i=1:lthu-1
    jisuan(3,i)=degao(1,i);
end
for i=2:lthu
    jisuan(4,i)=degao(2,i);
end

suanmain=zeros(2,lthu);
for i=1:lthu
    suanmain(1,i)=jisuan(2,i)+jisuan(3,i)+jisuan(4,i);
end

for i=1:lthu
    f(2*i-1)=suanmain(1,i);
    f(2*i)=suanmain(2,i);
end
f(41)=jisuan(1,20)+qx;
f

%求解節點位移
d=k\f';
u=d(1:2*numnod);
for i=1:numnod
    uf(1,i)=u(2*i-1);
    uf(2,i)=u(2*i);
end
u;
%求解高斯點的位移
% SOLVE FOR OUTPUT VARIABLES - DISPLACEMENTS
displ=zeros(1,2*numnod);
ind=0;
for gg=x
ind = ind+1;
gpos = gg(1:2);
v = domain(gpos,x,dm,numnod);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
displ(2*ind-1)= phi*uf(1,v)';
displ(2*ind)=phi*uf(2,v)';
end
disp1=displ';
for i=1:numnod
    disp2(1,i)=displ(2*i-1);
    disp2(2,i)=displ(2*i);
end

%求解高斯點的應力
ind = 0;
enorm=0;
L=length(gs);
stress(1:3,:)=stress;
gs(5,:)=stress(1,:);
gs(6,:)=stress(2,:);
gs(7,:)=stress(3,:);
for gg=gs
ind = ind+1;
gpos = gg(1:2);
weight = gg(3);
jac = gg(4);
stressx(1,1)=gs(5);
stressx(2,1)=gs(6);
stressx(3,1)=gs(7);
[tm,bsb]=fxx(stressx,tm,bsb,fc,ft);
% PLANE STRESS DMATRIX 
Dmat =(tm/(1-2*bsb)/(1+bsb))*[(1-bsb) bsb 0;bsb (1-bsb) 0;0 0 (1-2*bsb)/2];
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
Bmat=zeros(3,2*L);
for j=1:L
Bmat(1:3,(2*j-1):2*j) = [dphix(j) 0;0 dphiy(j);dphiy(j) dphix(j)];
end
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
end
stress(1:3,ind) = Dmat*Bmat*disp1(en);
end


%求解節點應力
ind=0;
enorm=0;
sx=x;
L1=length(x);
stress1(1:3,:)=stress1;
sx(3,:)=stress1(1,:);
sx(4,:)=stress1(2,:);
sx(5,:)=stress1(3,:);
for gg=sx
ind=ind+1;
 stressx(1,1)=sx(5);
stressx(2,1)=sx(6);
stressx(3,1)=sx(7);
[tm,bsb]=fxx(stressx,tm,bsb,fc,ft);
% PLANE STRESS DMATRIX 
Dmat =(tm/(1-2*bsb)/(1+bsb))*[(1-bsb) bsb 0;bsb (1-bsb) 0;0 0 (1-2*bsb)/2];
gpos=gg(1:2);
v = domain(gpos,x,dm,numnod);
L=length(v);
en=zeros(1,2);
[phi,dphix,dphiy]=shape(gpos,dmax,x,v,dm);
Bmat=zeros(3,2);
for j=1:L
Bmat(1:3,(2*j-1):2*j)=[dphix(j) 0;0 dphiy(j);dphiy(j) dphix(j)];
en(2*j-1)=2*v(j)-1;
en(2*j)=2*v(j);
stress1(1:3,ind)=Dmat*Bmat*u(en);
end
end
end
x;
stress1;
stress;
displ;



%繪制擋土墻體應力云圖
figure(1)
x1=-Lb/2:xspac:Lb/2;
x2=linspace(D/2,-D/2,ndivw+1);
z12=stress1(2,:);
z12=reshape(z12,ndivw+1,length(x1));
pcolor(x1,x2,z12);
shading interp;
xlabel('擋土墻體各點X坐標');
ylabel('擋土墻體各點Y坐標');
title('擋土墻體各點σy應力云圖');
axis equal;


figure(2)
x3=-Lb/2:xspac:Lb/2;
x4=linspace(D/2,-D/2,ndivw+1);
z34=stress1(1,:);
z34=reshape(z34,ndivw+1,length(x1));
pcolor(x3,x4,z34);
shading interp;
xlabel('擋土墻體各點X坐標');
ylabel('擋土墻體各點Y坐標');
title('擋土墻體各點σx應力云圖');
axis equal;


figure(3)
xb=[-0.5 0.5 0.5 -0.5 -0.5];
yb=[-2.5 -2.5 2.5 2.5 -2.5];
plot(xb,yb,'k','LineWidth',1.5);
hold on
xb1=x(:,no1)+500*uf(:,no1);
xb2=x(:,no2)+500*uf(:,no2);
xb3=x(:,no3)+500*uf(:,no3);
xb4=x(:,no4)+500*uf(:,no4);
plot(xb1(1,:),xb1(2,:),'r',xb2(1,:),xb2(2,:),'r',xb3(1,:),xb3(2,:),'r',xb3(1,:),xb3(2,:),'r','LineWidth',1.5);
axis equal;