clear
elementX=input('Enter no of elements in X direction:');
elementY=input('Enter no of elements in Y direction:');
Xnode=elementX+1;
Ynode=elementY+1;
Xzero=0.0;
Yzero=0.0;
Xlength=input('Enter length in X direction:');
Ylength=input('Enter length in Y direction:');
eps=input('Enter eccentricity:');
Xnel=Xnode-1;
Ynel=Ynode-1;
nel=Xnel*Ynel;
nnode=Xnode*Ynode;
delX=Xlength/Xnel;
delY=Ylength/Ynel;
nodes=zeros(nel,4);
gcoord=zeros(nnode,2);
grhs=zeros(55,1);
gm=zeros(55,55);
gm1=zeros(55,55);
sol=zeros(55,1);
Q1=zeros(11);
E=0;
gcoord
Xcoord=gcoord(:,1);
Ycoord=gcoord(:,2);
for i=1:nel
for j=1:4
X(j)=Xcoord(nodes(i,j));
Y(j)=Ycoord(nodes(i,j));
end
Xvec=[X(1),X(2),X(3),X(4)];
Yvec=[Y(1),Y(2),Y(3),Y(4)];
plot(Xvec,Yvec);
hold on;
elX=mean(Xvec(1:4));
elY=mean(Yvec(1:4));
text(elX,elY,num2str(i));
end
xlabel('X axis');
ylabel('Y axis');
title(num2str(i));
for jj=1:nnode
text(gcoord(jj,1),gcoord(jj,2),num2str(jj));
end
zeta(1)=-1/1.732;
eta(1)=-1/1.732;
zeta(2)=1/1.732;
eta(2)=-1/1.732;
zeta(3)=1/1.732;
eta(3)=1/1.732;
zeta(4)=-1/1.732;
eta(4)=1/1.732;
for i1=1:nel;
for i=1:4;
Jacb=[dxzeta,dyzeta;dxeta,dyeta];
detJacob=det(Jacb);
detJacob
Jacobinv=inv(Jacb);
inv(Jacb)
for i=1:4;
dnx(i)=Jacobinv(1,1)*dzeta(i)+Jacobinv(1,2)*deta(i);
dny(i)=Jacobinv(2,1)*dzeta(i)+Jacobinv(2,2)*deta(i);
end
h=0;
dh=0;
%eps=0.6;
for k=1:4;
h=h+Shape(k)*(1+eps*cos(Xcoord(nodes(i1,k))));
dh=dh+dnx(k)*(1+eps*cos(Xcoord(nodes(i1,k))));
end
h1=h.^3/12;
for i=1:4;
rhs(i,1)=rhs(i,1)-Shape(i)*1/2*dh*detJacob;
for j=1:4;
em(i,j)=em(i,j)+h1*(dnx(i)*dnx(j)+dny(i)*dny(j))*detJacob;
end
end
end
em
for ij=1:4;
ii=nodes(i1,ij);
grhs(ii,1)=grhs(ii,1)+rhs(ij,1);
for ji=1:4;
jj=nodes(i1,ji);
gm(ii,jj)=gm(ii,jj)+em(ij,ji);
end
end
end
for m=1:55
for n=1:55
gm1(m,n)=gm(m,n);
end
end
gm1
gm;
grhs;
for i=1:11
for j=1:55
gm(i,j)=0;
end
gm(i,i)=1;
grhs(i,1)=0;
end
for i=45:55
for j=1:55
gm(i,j)=0; gm(i1,j)=0;
end
gm(i1,i1)=1;
grhs(i1,1)=0;
end
for i=1:55
disp(gm(i,i))
end
gm;
grhs;
sol=gm\grhs;
sol
xx=zeros(11);
yy=zeros(11);
for i=12:22
yy(i)=sol(i);
xx(i)=gcoord(i,1);
end
for i=23:33
yy(i)=sol(i);
xx(i)=gcoord(i,1);
end
for i=34:44
yy(i)=sol(i);
xx(i)=gcoord(i,1);
end
plot(xx,yy)
wzeta=0;
weta=0;
for ie=1:40
for ig=1:4;
%int=int+1;
Shape(1)=(1-zeta(ig))*(1-eta(ig))/4;
Shape(2)=(1+zeta(ig))*(1-eta(ig))/4;
Shape(3)=(1+zeta(ig))*(1+eta(ig))/4;
Shape(4)=(1-zeta(ig))*(1+eta(ig))/4;
dzeta(1)=-(1-eta(ig))/4;
dzeta(2)=(1-eta(ig))/4;
dzeta(3)=(1+eta(ig))/4;
dzeta(4)=-(1+eta(ig))/4;
deta(1)=-(1-zeta(ig))/4;
deta(2)=-(1+zeta(ig))/4;
deta(3)=(1+zeta(ig))/4;
deta(4)=(1-zeta(ig))/4;
dxzeta=0;dyzeta=0;dxeta=0;dyeta=0;
wzeta=wzeta+Shape(i)*sol(ii,1)*cos(gcoord(ii,1))*detJacob;
weta=weta+Shape(i)*sol(ii,1)*sin(gcoord(ii,1))*detJacob;
end
end
end
wload=sqrt(wzeta.^2+weta.^2)
wload
array indices must be positve
this error is coming please help
FIRST SET OF INPUT DATA:
Number of elements in X direction = 10
Number of elements in Y direction = 4
The length in X direction = 3.14
The length in Y direction = 2
Eccentricity ratio = 0.15
