clear;
Lb=48;
D=12;
young=30e6;nu=0.3;
P=1000;
Dmat=(young/(1-nu^2))*[1 nu 0;nu 1 0;0 0 (1-nu)/2];
ndivl=10;ndivw=4;
[x,conn1,conn2,numnod]=mesh2(Lb,D,ndivl,ndivw);
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);
ndivlq=10;ndivwq=4;
[xc,conn,numcell,numq]=mesh2(Lb,D,ndivlq,ndivwq);
quado=4;
[gauss]=gauss2(quado);
numq2=numcell*quado^2;
gs=zeros(4,numq2);
[gs]=egauss(xc,conn,gauss,numcell);
k=sparse(numnod*2,numnod*2);
for gg=gs
gpos=gg(1:2);
weight=gg(3);
jac=gg(4);
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
ind1=0;ind2=0;
for j=1:numnod
if(x(1,j)==0.0)
ind1=ind1+1;
nnu(1,ind1)=x(1,j);
nnu(2,ind1)=x(2,j);
end
if(x(1,j)==Lb)
ind2=ind2+1;
nt(1,ind2)=x(1,j);
nt(2,ind2)=x(2,j);
end
end
lthu=length(nnu);
ltht=length(nt);
ubar=zeros(lthu*2,1);
f=zeros(numnod*2,1);
ind=0;
gauss=gauss2(quado);
for i=1:(ltht-1)
ycen=(nt(2,i+1)+nt(2,i))/2;
jcob=abs((nt(2,i+1)-nt(2,i))/2);
for j=1:quado
mark(j)=ycen-gauss(1,j)*jcob;
ind=ind+1;
gst(1,ind)=nt(1,i);
gst(2,ind)=mark(j);
gst(3,ind)=gauss(2,j);
gst(4,ind)=jcob;
end
end
gsu=gst;
gsu(1,1:ind)=zeros(1,ind);
qk=zeros(1,2*lthu);
Imo=(1/12)*D^3;
for gt=gst
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=-(P/(2*Imo))*((D^2)/4-gpos(2,1)^2);
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
GG=zeros(numnod*2,lthu*2);
ind1=0;ind2=0;
for i=1:(lthu-1)
ind1=ind1+1;
m1=ind1;m2=m1+1;
y1=nnu(2,m1);y2=nnu(2,m2);
len=y1-y2;
for j=1:quado
ind2=ind2+1;
gpos=gsu(1:2,ind2);
weight=gsu(3,j);
jac=gsu(4,j);
fac11=(-P*nnu(2,m1))/(6*young*Imo);
fac2=P/(6*young*Imo);
xp1=gpos(1,1);
yp1=gpos(2,1);
uxex1=(6*Lb-3*xp1)*xp1+(2+nu)*(yp1^2-(D/2)^2);
uxex1=uxex1*fac11;
uyex1=3*nu*yp1^2*(Lb-xp1)+0.25*(4+5*nu)*xp1*D^2+(3*Lb-xp1)*xp1^2;
uyex1=uyex1*fac2;
N1=(gpos(2,1)-y2)/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
f=[f;zeros(lthu*2,1)];
f(numnod*2+1:numnod*2+2*lthu,1)=-qk';
m=sparse([k GG;GG' zeros(lthu*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
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
ind=0;
enorm=0;
for gg=gs
ind=ind+1;
gpos=gg(1:2);
weight=gg(3);
jac=gg(4);
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);
stressex(1,ind)=(1/Imo)*P*(Lb-gpos(1,1))*gpos(2,1);
stressex(2,ind)=0;
stressex(3,ind)=-0.5*(P/Imo)*(0.25*D^2-gpos(2,1)^2);
err=stress(1:3,ind)-stressex(1:3,ind);
err2=weight*jac*(0.5*(inv(Dmat)*err)'*(err));
enorm=enorm+err2;
end
enorm=sqrt(enorm);
displx=displ(1:2:109);
disply=displ(2:2:110);
displtot=[displx;disply];
function[x,conn,numcell,numq]=mesh2(length,height,ndivl,ndivw)
numcell=ndivw*ndivl;
numq=(ndivl+1)*(ndivw+1);
for i=1:(ndivl+1)
for j=1:(ndivw+1)
x(1,((ndivw+1)*(i-1)+j))=(length/ndivl)*(i-1);
x(2,((ndivw+1)*(i-1)+j))=-(height/ndivw)*(j-1)+height/2;
end
end
for j=1:ndivl
for i=1:ndivw
elemn=(j-1)*ndivw+i;
nodet(elemn,1)=elemn+(j-1);
nodet(elemn,2)=nodet(elemn,1)+1;
nodet(elemn,3)=nodet(elemn,2)+ndivw+1;
nodet(elemn,4)=nodet(elemn,3)-1;
end
end
conn=nodet';
function [gs]=egauss(xc,conn,gauss,numcell)
index=0;
one=ones(1,4);
psiJ=[-1,+1,+1,-1];etaJ=[-1,-1,+1,+1];
l=size(gauss);
l=l(2);
for e=1:numcell
for j=1:4
je=conn(j,e);xe(j)=xc(1,je);ye(j)=xc(2,je);
end
for i=1:l
for j=1:l
index=index+1;
eta=gauss(1,i);psi=gauss(1,j);
N=.25*(one+psi*psiJ).*(one+eta*etaJ);
NJpsi=.25*psiJ.*(one+eta*etaJ);
NJeta=.25*etaJ.*(one+psi*psiJ);
xpsi=NJpsi*xe';ypsi=NJpsi*ye';xeta=NJeta*xe';yeta=NJeta*ye';
jcob=xpsi*yeta-xeta*ypsi;
xq=N*xe';yq=N*ye';
gs(1,index)=xq;
gs(2,index)=yq;
gs(3,index)=gauss(2,i)*gauss(2,j);
gs(4,index)=jcob;
end
end
end
function v=gauss2(k)
v(1,1)=-.861136311594052575224;
v(1,2)=-.339981043584856264803;
v(1,3)=-v(1,2);
v(1,4)=-v(1,1);
v(2,1)=.347854845137453857373;
v(2,2)=.652145154862546142627;
v(2,3)=v(2,2);
v(2,4)=v(2,1);
function v=domain(gpos,x,dm,nnodes)
dif=gpos*ones(1,nnodes)-x;
a=dm-abs(dif);
b=[a;zeros(1,nnodes)];
c=all(b>=-100*eps);
v=find(c);
function [phi,dphix,dphiy]=shape(gpos,dmax,x,v,dm)
L=length(v);
won=ones(1,L);
nv=x(1:2,v);
p=[won;nv];
dif=gpos*won-nv;
t=dm(1:2,v)/dmax;
[w,dwdx,dwdy]=cubwgt(dif,t,v,dmax,dm);
B=p.*[w;w;w];
pp=zeros(3);
aa=zeros(3);
daax=zeros(3);
daay=zeros(3);
for i=1:L
pp=p(1:3,i)*p(1:3,i)';
aa=aa+w(1,i)*pp;
daax=daax+dwdx(1,i)*pp;
daay=daay+dwdy(1,i)*pp;
end
pg=[1 gpos'];
[L,U,PERM]=lu(aa);
for i=1:3
if i==1
C=PERM*pg';
elseif i==2
C=PERM*([0 1 0]'-daax*gam(1:3,1));
elseif i==3
C=PERM*([0 0 1]'-daay*gam(1:3,1));
end
D1=C(1);
D2=(C(2)-L(2,1)*D1);
D3=(C(3)-L(3,1)*D1-L(3,2)*D2);
gam(3,i)=D3/U(3,3);
gam(2,i)=(D2-U(2,3)*gam(3,i))/(U(2,2));
gam(1,i)=(D1-U(1,2)*gam(2,i)-U(1,3)*gam(3,i))/U(1,1);
end
phi=gam(1:3,1)'*B;
dbx=p.*[dwdx;dwdx;dwdx];
dby=p.*[dwdy;dwdy;dwdy];
dphix=gam(1:3,2)'*B+gam(1:3,1)'*dbx;
dphiy=gam(1:3,3)'*B+gam(1:3,1)'*dby;
function [w,dwdx,dwdy]=cubwgt(dif,t,v,dmax,dm)
l=length(v);
for i=1:l
drdx=sign(dif(1,i))/dm(1,v(i));
drdy=sign(dif(2,i))/dm(2,v(i));
rx=abs(dif(1,i))/dm(1,v(i));
ry=abs(dif(2,i))/dm(2,v(i));
if rx>0.5
wx=(4/3)-4*rx+4*rx*rx-(4/3)*rx^3;
dwx=(-4+8*rx-4*rx^2)*drdx;
elseif rx<=0.5
wx=(2/3)-4*rx*rx+4*rx^3;
dwx=(-8*rx+12*rx^2)*drdx;
end
if ry>0.5
wy=(4/3)-4*ry+4*ry*ry-(4/3)*ry^3;
dwy=(-4+8*ry-4*ry^2)*drdy;
elseif ry<=0.5
wy=(2/3)-4*ry*ry+4*ry^3;
dwy=(-8*ry+12*ry^2)*drdy;
end
w(i)=wx*wy;
dwdx(i)=wy*dwx;
dwdy(i)=wx*dwy
end
