why i get this ERROR???

V

[ammar ahmed]

clc;
clear all;
E=200e9;
p=0.3
d=1;
t=1; 
L=5;
endload=-1000
nmeshx=32;
nmeshy=8;
%........%
Xspacing=(L/nmeshx)/2;
yspacing=(d/nmeshy)/2;
nel=nmeshx*nmeshy;
nnode=(nmeshx*2+1)*(nmeshy+1)+nmeshy*(nmeshx+1);
%.........coordinates 
nodecoord=zeros(nnode,2);
c=0;
for i=1:(nmeshy*2+1)
    if mod(i,2)==1
        for j=1:(2*nmeshx+1)
            c=c+1;
      xcoord=Xspacing*(j-1);
            ycoord=yspacing*(i-1);
            nodecoord(c,:) = [xcoord ycoord];
        end
    else
    for j=1:(nmeshx+1)
        c=c+1;
        xcoord=(2*Xspacing)*(j-1);
        ycoord=yspacing*(i-1);
        nodecoord(c,:)=[xcoord ycoord];
      end
    end
end
%..........conenectivy
conn=zeros(nel,4);
rdiff=nmeshx*2+nmeshx+2;
rdiffsmall=nmeshx*2+1;
c=0;
for j=0:nmeshy-1
    for i=0:nmeshx-1
        c=c+1;
        c1=(2*i+1)+rdiff*j;
        c2=c1+2;
        c4=c1+rdiff;
        c3=c4+2;
        conn(c,:)=[c1;c2;c3;c4];
    end
end
%..........load matrix
Load = zeros(nmeshy*2+1,3);
laod(1,:)=[rdiffsmall 0 endload/(nmeshy*2+1)];
for i=1:nmeshy*2
    Load(i+1,:)=[Load(i,1)+round((nmeshx+1)*(abs(cos((i+1)*pi()/2)))+rdiffsmall*(abs(cos(i*pi()/2 )))) 0 endload/(nmeshy*2+1)];
end
nLoad=size(Load,1);
%..................support
s=zeros(nmeshy*2+1,3);
s(1,:)=[1 1 1];
for i=1:nmeshy*2                 %s=[nodelid ux uy ] 1..fixed 0...free
    s(i+1,:)=[Load(i,1)+ 1 1 1];
end
numsupport=size(s,1);
c=0;
for i=1:nel
     nodecoordelem=zeros(4,2);
     for j=1:4
   nodecoordelem(j,:)=nodecoord(conn(i,j),:);
     end
     c=c+1;
     nodecoordelemlist(c).nodecoordelem=nodecoordelem;
end
EM=E/(1-p^2).*[1 p 0;p 1 0; 0 0 (1-p)/2];

%shape functions
syms ks ta
    n1 = 0.25*(1-ks)*(1-ta);
    n2 = 0.25*(1+ks)*(1-ta);
    n3 = 0.25*(1+ks)*(1+ta);
    n4 = 0.25*(1-ks)*(1+ta);
N8=[n1;n2;n3;n4];
dnks = [-0.25*(1-ta) 0.25*1*(1-ta)  0.25*1*(1+ta)  0.25*-1*(1+ta)];
dnta = [-0.25*(1-ks) 0.25*-1*(1+ks) 0.25*1*(1+ks) 0.25*(1-ks)];
D8=[dnks 
    dnta];
%element stiffness matrix
for g=1:nel
jac=D8*nodecoordelemlist(g).nodecoordelem; %jacobian matrix
capgamma = inv(jac);      %inverse jacobian matrix
B=sym(zeros(3,8)); 
for k=0:3
m=2*k+1;
n=m+1;
B(1,m) = capgamma(1,1)*dnks(k+1)+capgamma(1,2)*dnta(k+1);
B(2,n) = capgamma(2,1)*dnks(k+1)+capgamma(2,2)*dnta(k+1);
B(3,m) = B(2,n);
B(3,n) = B(1,m);
end
kelemlist(g).kelem = int(int(transpose(B)*EM*B*t*det(jac),ks,-1,1),ta,-1,1);
                                 
kelemlist(g).B = B;
end

%global stiffness matrix
kk = zeros(nnode*2,nnode*2);
for k=1:nel
 kelem=kelemlist(k).kelem;
for j=1:4
 for i=1:4
    n=conn(k,i);
    m=conn(k,j);
kk(2*n-1,2*m-1)=kk(2*n-1,2*m-1)+kelem(2*i-1,2*j-1);
kk(2*n-1,2*m)=kk(2*n-1,2*m)+kelem(2*i-1,2*j);
kk(2*n,2*m-1)=kk(2*n,2*m-1)+kelem(2*i,2*j-1);
kk(2*n,2*m)=kk(2*n,2*m)+kelem(2*i,2*j);
end
end
end
%apply the boundry condition 
kkWosup=kk;
for i=1:numsupport
    n=s(i,1);
    if(s(i,2)==1)
        kk(2*n-1,:)=0;
        kk(:,2*n-1)=0;
        kk(2*n-1,2*n-1)=1;
    end
    if(s(i,3)==1)
        kk(2*n,:)=0;
        kk(:,2*n)=0;
        kk(2*n,2*n)=1;
        end
end
    %load vector
    Loadvector = zeros(nnode*2,1);
for i=1:nLoad
n =Load(i,1);
Loadvector(2*n-1)=Load(i,2);
Loadvector(2*n) = Load(i,3);
end
%displacement computation
u = kk\Loadvector;
for i=1:nnode
unodelist(i).node = [u(2*i-1);u(2*i)];
end
for i=1:nel
temp = [];
for j=1:4
temp = [temp ; unodelist(conn(i,j)).node];
end
uelemlist(i).elem = temp;
end
beamtheory = (endload*L^3)/(3*E*(t*d^3/12)+(endload*L))/((E/2*(1-p))*(5/6)*d*t);

Stress=beamtheory/(2*Xspacing*yspacing*-1e-2);
%stress computation
scp = [-1 -1;1 -1;1 1;-1 1];
    for i=1:nel
uelem = uelemlist(i).elem';
temp = [];
felemlist(i).felem = double(kelemlist(i).kelem*uelem');
sigma = EM*kelemlist(i).B*uelem';
temp = [];
for j=1:4
temp = [temp double(subs(subs(sigma,ks,scp(j,1)),ta,scp(j,2)))];
end
sigmaelemlist(i).sigma = temp;
 end
fkk = kkWosup*u;

[ammar ahmed]

the code is above pleace help me as fast as you can i was like spending to much time trying to figur out the the problem

Answer 1

[Star Strider]

I suspect it is because ‘(2*n-1)’ is less than or equal to 0.

I have no idea what ‘n’ is, so it could also be the problem if it is not an integer.

Since apparently ‘Load’ is not throwing an error, I suspect it is a function you (or someone else) wrote. The MATLAB function has a lower-case ‘L’, and since MATLAB is case-sensitive, the difference is not the problem.

[Steven Lord]

The Workspace component in the picture of the MATLAB Desktop indicates n is 0. That means the code was trying to assign to element -1 of Loadvector. Arrays in MATLAB don't have a -1st element.

[ammar ahmed]

so what i have to do?? when i chanched n ny i it give me another kind of Error !!

[Walter Roberson]

It appears that you have some 0 in the first column of your array Load

[ammar ahmed]

[ammar ahmed]

i figerd out the firest problem but i have got another one

[Walter Roberson]

As outside observers, we have no reason to expect that the kk matrix you are creating will be non-singular.

Typically when a stiffness matrix is singular, the implication is that there are too many degrees of freedom, and that additional entries need to be added to pin down some part that is too free to move.

>> size(kk)
ans =
        1698        1698

>> rank(kk)
ans =
   591

Not even close :( You are missing pinning down a lot of values !

[Star Strider]

Steven — I didn’t see ‘n’ in the screencap. Thank you for discovering it. (My Answer remains valid.)

Walter — Thank you!