clear; 
% DEFINE BOUNDARIES/PARAMETERS 
Lb=1;
D=5;
Q=0;
Q1=0;
Q2=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


% 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(numnod*2,1);

%SET UP GAUSS POINTS ALONG TRACTION BOUNDARY
xs=linspace(-Lb/2,Lb/2,9);
xz=linspace(D/2,-D/2,21);
xshang=[xs;(D/2)*ones(1,length(xs))];
xxia=[xs;(-D/2)*ones(1,length(xs))];
xzuo=[(-Lb/2)*ones(1,length(xz));xz];
xyou=[(Lb/2)*ones(1,length(xz));xz];
lxs=length(xshang);
lxx=length(xxia);
lxz=length(xzuo);
lxy=length(xyou);
ind=0;
gauss=gauss2(quado);
for i=1:(lxy-1)
ycen=(xyou(2,i+1)+xyou(2,i))/2;
jcob = abs((xyou(2,i+1)-xyou(2,i))/2);
for j=1:quado
mark(j) = ycen-gauss(1,j)*jcob;
ind = ind+1;
gst(1,ind)=xyou(1,i);
gst(2,ind)=mark(j);
gst(3,ind)=gauss(2,j);
gst(4,ind)=jcob;
end
end

%SET UP GAUSS POINTS ALONG DISPLACEMENT BOUNDARY
gsu=gst;
gsu(1,1:ind)=ones(1,ind)*(-Lb/2);

%INTEGRATE FORCES ALONG BOUNDARY
t1=1:4:4*(lxy-1);
t2=4:4:4*(lxy-1);
t3=2:4:4*(lxy-1);
t4=3:4:4*(lxy-1);
gst1=zeros(4,length(t1));
gst2=zeros(4,length(t2));
gst3=zeros(4,length(t3));
gst4=zeros(4,length(t4));
for i=1:length(t1)
       gst1(:,i)=gst(:,t1(i));
end
for i=1:length(t2)
       gst2(:,i)=gst(:,t2(i));
end
for i=1:length(t3)
       gst3(:,i)=gst(:,t3(i));
end
for i=1:length(t4)
       gst4(:,i)=gst(:,t4(i));
end
gst12=[gst1,gst2];
gst34=[gst3,gst4];
for gt=gst12
gpos = gt(1:2);
weight=gt(3);
jac = gt(4);
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
force=zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
tx=Q1*0.1997*abs(xyou(2,1)-xyou(2,2))/(weight*jcob);
ty=0;
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
force(2*i-1)=tx*phi(i);
force(2*i) = ty*phi(i);
end
f(en) = f(en) + jac*weight*force';
end
for gt=gst34
gpos = gt(1:2);
weight=gt(3);
jac = gt(4);
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
force=zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
tx=Q1*0.3003*abs(xyou(2,1)-xyou(2,2))/(weight*jcob);
ty=0;
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
force(2*i-1)=tx*phi(i);
force(2*i) = ty*phi(i);
end
f(en) = f(en) + jac*weight*force';
end
gsu1=zeros(4,length(t1));
gsu2=zeros(4,length(t2));
gsu3=zeros(4,length(t3));
gsu4=zeros(4,length(t4));
for i=1:length(t1)
       gsu1(:,i)=gsu(:,t1(i));
end
for i=1:length(t2)
       gsu2(:,i)=gsu(:,t2(i));
end
for i=1:length(t3)
       gsu3(:,i)=gsu(:,t3(i));
end
for i=1:length(t4)
       gsu4(:,i)=gsu(:,t4(i));
end
gsu12=[gsu1,gsu2];
gsu34=[gsu3,gsu4];
for gt=gsu12
gpos = gt(1:2);
weight=gt(3);
jac = gt(4);
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
force=zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
tx=-Q1*0.1997*abs(xzuo(2,1)-xzuo(2,2))/(weight*jcob);
ty=0;
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
force(2*i-1)=tx*phi(i);
force(2*i) = ty*phi(i);
end
f(en) = f(en) + jac*weight*force';
end
for gt=gsu34
gpos = gt(1:2);
weight=gt(3);
jac = gt(4);
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
force=zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
tx=-Q1*0.3003*abs(xzuo(2,1)-xzuo(2,2))/(weight*jcob);
ty=0;
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
force(2*i-1)=tx*phi(i);
force(2*i) = ty*phi(i);
end
f(en) = f(en) + jac*weight*force';
end


%shijia sanjiao hezai
gstsh=gsu;
ind=0;
for i=1:length(gstsh)
   ind=ind+1;
   gstshang(:,ind)=gstsh(:,ind);
end


%INTEGRATE FORCES ALONG BOUNDARY
t1=1:4:length(gstshang);
t2=4:4:length(gstshang);
t3=2:4:length(gstshang);
t4=3:4:length(gstshang);

degao=zeros(3,length(gstshang));
degao(1,t1)=0.1997*ones(1,length(t1))*abs(xyou(2,1)-xyou(2,2));
degao(1,t2)=0.1997*ones(1,length(t1))*abs(xyou(2,1)-xyou(2,2));
degao(1,t3)=0.3003*ones(1,length(t1))*abs(xyou(2,1)-xyou(2,2));
degao(1,t4)=0.3003*ones(1,length(t1))*abs(xyou(2,1)-xyou(2,2));

ind=0;
suan=0;
for i=1:length(gstshang)
 ind=ind+1;
 suan=suan+degao(1,i);
 jisuan(ind)=suan;
end

ind=0;
for i=1:length(gstshang)
 ind=ind+1;
 suan=Q*jisuan(i)/D;
 degao(2,ind)=suan;
end

ind=0;
suanmian(1)=0;
suanmian(2:length(gstshang)+1)=degao(2,:);
for i=1:length(gstshang)
 ind=ind+1;
 suan=(suanmian(i)+suanmian(i+1))*degao(1,i)/2;
 degao(3,ind)=suan;
end

ind=0;
for gt=gstshang
 ind=ind+1;
gpos = gt(1:2);
weight=gt(3);
jac = gt(4);
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
force=zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
tx=degao(3,ind)/(weight*jcob);
ty=0;
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
force(2*i-1)=tx*phi(i);
force(2*i) = ty*phi(i);
end
f(en) = f(en) + jac*weight*force';
end

%SET UP GAUSS POINTS ALONG TRACTION BOUNDARY
ind=0;
gauss=gauss2(quado);
gaussang1=gauss(1,:);
gaussang1=(-1)*gauss(1,:);
gaussang=[gaussang1;gauss(2,:)];
for i=1:(lxs-1)
ycen=(xshang(1,i+1)+xshang(1,i))/2;
jcob = abs((xshang(1,i+1)-xshang(1,i))/2);
for j=1:quado
mark(j) = ycen-gaussang(1,j)*jcob;
ind = ind+1;
gsth(1,ind)=mark(j);
gsth(2,ind)=xshang(2,i);
gsth(3,ind)=gaussang(2,j);
gsth(4,ind)=jcob;
end
end
gsuh=gsth;
gsuh(2,:)=-1*gsth(2,:);
%INTEGRATE FORCES ALONG BOUNDARY
t1=1:4:4*(lxs-1);
t2=4:4:4*(lxs-1);
t3=2:4:4*(lxs-1);
t4=3:4:4*(lxs-1);
gst1=zeros(4,length(t1));
gst2=zeros(4,length(t2));
gst3=zeros(4,length(t3));
gst4=zeros(4,length(t4));
for i=1:length(t1)
       gst1(:,i)=gsth(:,t1(i));
end
for i=1:length(t2)
       gst2(:,i)=gsth(:,t2(i));
end
for i=1:length(t3)
       gst3(:,i)=gsth(:,t3(i));
end
for i=1:length(t4)
       gst4(:,i)=gsth(:,t4(i));
end
gst12=[gst1,gst2];
gst34=[gst3,gst4];
for gt=gst12
gpos = gt(1:2);
weight=gt(3);
jac = gt(4);
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
force=zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
tx=0;
ty=-Q2*0.1997*abs(xshang(1,1)-xshang(1,2))/(weight*jcob);
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
force(2*i-1)=tx*phi(i);
force(2*i) = ty*phi(i);
end
f(en) = f(en) + jac*weight*force';
end
for gt=gst34
gpos = gt(1:2);
weight=gt(3);
jac = gt(4);
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
force=zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
tx=0;
ty=-Q2*0.3003*abs(xshang(1,1)-xshang(1,2))/(weight*jcob);
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
force(2*i-1)=tx*phi(i);
force(2*i) = ty*phi(i);
end
f(en) = f(en) + jac*weight*force';
end


% INTEGRATE G MATRIX AND Q VECTOR ALONG DISPLACEMENT BOUNDARY
qk = zeros(1,2*lxx);
GG = zeros(numnod*2,lxx*2);
ind1=0;ind2=0;
for i=1:(lxx-1)
ind1=ind1+1;
m1 = ind1; m2 = m1+1;
x1 = xxia(1,m1); x2 =xxia(1,m2);
len = x1-x2;
for j=1:quado               
ind2=ind2+1;
gpos = gsuh(1:2,ind2);
weight = gsuh(3,j);
jac = gsuh(4,j);
xp1 = gpos(1,1);
yp1 = gpos(2,1);
uyex1 = 0;
uxex1=0;
N1 = (gpos(1,1)-x2)/len; N2 = 1-N1;
qk(2*m1-1) = qk(2*m1-1)-weight*jac*N1*uxex1;
qk(2*m1) = qk(2*m1) - weight*jac*N1*uyex1;
qk(2*m2-1) = qk(2*m2-1) -weight*jac*N2*uxex1;
qk(2*m2) = qk(2*m2) - weight*jac*N2*uyex1;
v = domain(gpos,x,dm,numnod);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
L = length(v);
for n=1:L
G1 = -weight*jac*phi(n)*[N1 0;0 N1];
G2 = -weight*jac*phi(n)*[N2 0;0 N2];
c1=2*v(n)-1;c2=2*v(n);c3=2*m1-1;c4=2*m1;
c5=2*m2-1;c6=2*m2;
GG(c1:c2,c3:c4)=GG(c1:c2,c3:c4)+ G1;
GG(c1:c2,c5:c6)=GG(c1:c2,c5:c6)+G2;
end
end
end

% ENFORCE BC'S USING LAGRANGE MULTIPLIERS
f = [f;zeros(lxx*2,1)];
f(numnod*2+1:numnod*2+2*lxx,1) = qk';
m = sparse([k GG;GG' zeros(lxx*2)]);
d=m\f;
u=d(1:2*numnod);
for i=1:numnod
u2(1,i) = u(2*i-1);
u2(2,i) = u(2*i);
end



% 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*u2(1,v)';
displ(2*ind) = phi*u2(2,v)';
end
for i=1:numnod
    disp2(1,i)=displ(2*i-1);
    disp2(2,i)=displ(2*i);
end

% SOLVE FOR STRESSES AT GAUSS POINTS
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*u(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;
gpos=gg(1:2);
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];
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
disp2


%繪制擋土墻體應力云圖
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)+400*uf(:,no1);
xb2=x(:,no2)+400*uf(:,no2);
xb3=x(:,no3)+400*uf(:,no3);
xb4=x(:,no4)+400*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;