clear; 
% DEFINE BOUNDARIES/PARAMETERS 
Lb=1;
D=5;
Q=100;
as=1.8e-5;
h1=0.25;
tm=0;bsb=0;
tm0=3e7;bsb0=0.167;
stm=2.1e8; sbsb=0.3;
tm1=(tm0*(h1-as)+stm*as)/h1;
bsb1=(bsb0*(h1-as)+sbsb*as)/h1;
tm2=(tm0+tm1)/2;
bsb2=(bsb0+bsb1)/2;


%平面應(yīng)變問題
Dmat=(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);

%ind=0;
%for i=1:length(x)
%if (x(1,i)<-0.375)
 %   tm=tm1;
%    bsb=bsb1;
%end
%if(x(1,i)=-0.375)
 %tm=tm2;
  %  bsb=bsb2;
  %end   
%if(x(1,i)> -0.375)
 %tm=tm0;
 %  bsb=bsb0;
 %end   




% 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 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);



% 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);
% DETERMINE NODES IN NEIGHBORHOOD OF GAUSS POINT
v = domain(gpos,x,dm,numnod);
for i=1:length(v)
if (x(1,v(i))<-0.375)
    tm=tm1;
    bsb=bsb1;
end
if(x(1,v(i))==-0.375)
 tm=tm2;
    bsb=bsb2;
end   
if(x(1,v(i))> -0.375)
 tm=tm0;
    bsb=bsb0;
end 
end
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;
    a=qn(2*i-1);
   k(a,a)=k(a,a)*100000;
end
for i=1:length(f)
    qn(2*i)=2*f(i);
    a=qn(2*i);
   k(a,a)=k(a,a)*100000;
end
 for i=1:length(f)
    qn(2*i)=2*f(i);
     qn(2*i-1)=2*f(i)-1;
    a=qn(2*i);
    b=qn(2*i-1);
   k(a,b)=k(a,b)*100000;
 end
for i=1:length(f)
    qn(2*i)=2*f(i);
     qn(2*i-1)=2*f(i)-1;
    a=qn(2*i);
    b=qn(2*i-1);
   k(b,a)=k(b,a)*100000;
 end

k;

%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);


%計(jì)算zuo節(jié)點(diǎn)荷載
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

%求解節(jié)點(diǎn)位移
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
uf


%求解節(jié)點(diǎn)應(yīng)力
ind=0;
enorm=0;
for gg=x
ind=ind+1;
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
x;
stress1



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


figure(2)
x3=-Lb/2:xspac:Lb/2;
x4=linspace(D/2,-D/2,ndivw+1);
z34=stress1(2,:);
z34=reshape(z34,ndivw+1,length(x1));
pcolor(x3,x4,z34);
shading interp;
xlabel('擋土墻體各點(diǎn)X坐標(biāo)');
ylabel('擋土墻體各點(diǎn)Y坐標(biāo)');
title('擋土墻體各點(diǎn)σx應(yīng)力云圖');
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)+1000*uf(:,no1);
xb2=x(:,no2)+1000*uf(:,no2);
xb3=x(:,no3)+1000*uf(:,no3);
xb4=x(:,no4)+1000*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);