extending matrix with numbers (not zero)

1 view (last 30 days)
zina shadidi
zina shadidi on 17 Nov 2020
Edited: KSSV on 17 Nov 2020
In the folowing code; I have error
Index in position 1 exceeds array bounds (must not exceed 4).
Error in Untitled1111111 (line 137)
load(1,(head(1,i)-1)*(dof+j))=fload(j+1,i);
the matrices dimention is:
load(1x300), head (4x30), fload(4x30)
and the cod is:
close all;
clear all
dbstop if error
% Initial data:
TubeLength =2.38e-9;
%%TubeLength = input ('TubeLength'); %m
TubeDiameter=1.62e-9;
%%%TubeDiameter = input ('TubeDiameter'); %m
r=TubeDiameter/2;
A=pi*r^2;
D = 3; %degrees of freedom at each node
presc_disp=1e-9; %m
L=0.141e-9; %m Length of a C-C bond
numIterations=1; %number of iterations
numNaN=0;
% %%%%% the ref .
%%Force constants from Tserpes et al. (2005):
% Kr is EA=6.52e-7 %N·nm^-1
% ktheta is EI=8.76e-10 %N·nm·rad^-2
% ktao is GJ=2.78e-10 %N·nm·rad^-2
EA=L*652 %N
EI=L*8.76e-19 %N·m^2·rad^-2
GJ=L*2.78e-19; %N·m^2·rad^-2
% Transformation to nondimensional space:
% [F]=EI/L^2; [L]=L. Then, multiply the resulting forces with
% EI/L^2 and displacements with L to get forces and
% displacements in nm.
Lreal=L;
EIreal=EI;
presc_disp_real=presc_disp;
EA=EA/EI*L^2;
GJ=GJ/EI;
presc_disp=presc_disp/L;
L=1;
EI=1;
% Stiffness matrix of a C-C bond (adapted to 3 degree of freedom):
KBond=[EA/L 0 0 -EA/L 0 0
0 12*EI/L^3 0 0 -12*EI/L^3 0 ;
0 0 12*EI/L^3 0 0 -12*EI/L^3;
-EA/L 0 0 EA/L 0 0 ;
0 -12*EI/L^3 0 0 12*EI/L^3 0 ;
0 0 -12*EI/L^3 0 0 12*EI/L^3]
% Calculations:
[elementss1]=xlsread('elementss1.xlsx');
[nodess1]=xlsread('nodess1.xlsx');
%%nnodess = input ('number of nodes');
%%nElementss = input ('number of elements');
nnodess1= 200; %%% number of nodes
nElementss2 = 200; %%%% number of elements
nodesDef=zeros(1,nnodess1);
k=zeros(nnodess1*D,nnodess1*D);
% Defects:
defectsPercentage=1; %per cent
numDef=round((nnodess1*defectsPercentage)/100);
for i=1:numDef
def=round((nodess1-1)*rand)+1; %The addition of 1 and -1 is to avoid zero as the resulting random number. elements(D,:)=[];
nnodess1=nnodess1-1;
end
zeroMatrix = [0 0 0; 0 0 0; 0 0 0];
vi = [1 0 0];
vj = [0 1 0];
vk = [0 0 1];
for i=1:nnodess1-1
n = nodess1(2,i);
m = nodess1(3,i)
localZ = [1-2, 3-2, 3-4];
localZ = localZ/sqrt(localZ*localZ');
% Find local y
if not(localZ(3)==0)
y3 = (-localZ(1)-2*localZ(2))/localZ(3);
localY = [1 2 y3];
elseif not(localZ(2)==0)
y2 = (-localZ(1)-2*localZ(3))/localZ(2);
localY = [1 y2 2];
elseif not(localZ(1)==0)
y1 = (-localZ(2)-2*localZ(3))/localZ(1);
localY = [y1 1 2];
end
localY = localY/sqrt(localY*localY');
% Find local z
localX = cross(localY, localZ);
localX = localX/sqrt(localX*localX');
% Values for rotation matrix
cosZI = localZ*vi';
cosZJ = localZ*vj';
cosZK = localZ*vk';
cosYI = localY*vi';
cosYJ = localY*vj';
cosYK = localY*vk';
cosXI = localX*vi';
cosXJ = localX*vj';
cosXK = localX*vk';
rotationSmall = [
cosZI cosZJ cosZK;
cosYI cosYJ cosYK;
cosXI cosXJ cosXK];
RR = [rotationSmall zeroMatrix;zeroMatrix rotationSmall];
KElement = abs(RR'*KBond*RR)
n = nodess1(2,i);
m = nodess1(3,i);
k((n-1)*D+1:n*D,(n-1)*D+1:n*D)=k((n-1)*D+1:n*D, (n-1)*D+1:n*D)+KElement(1:D, 1:D);
k((m-1)*D+1:m*D,(m-1)*D+1:m*D)=k((m-1)*D+1:m*D, (m-1)*D+1:m*D)+KElement(D+1:2*D, D+1:2*D);
k((n-1)*D+1:n*D,(m-1)*D+1:m*D)=k((n-1)*D+1:n*D, (m-1)*D+1:m*D)+KElement(1:D, D+1:D*2);
k((m-1)*D+1:m*D,(n-1)*D+1:n*D)=k((m-1)*D+1:m*D, (n-1)*D+1:n*D)+KElement(D+1:D*2, 1:D);
end
head = [1 4 5 16 17 34 35 51 52 72 73 86 87 476 477 478 479 482 483 488 489 490 491 492 493 494 500 501 503 504;
1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
];
fload = [1 4 5 16 17 34 35 51 52 72 73 86 87 476 477 478 479 482 483 488 489 490 491 492 493 494 500 501 503 504;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
]*presc_disp;
dof=3;nNodes=100;
save k_real.mat k;
load=zeros(1,dof*nNodes);
for i=1:size(head(1,:),2)
for j=1:dof;
if head(j+1,i)~=0 % displacement prescribed
for m=1:dof*nNodes;
load(m)=load(m)-k(m,(head(1,i)-1)*dof+j)*fload(j+1,i);
end
i=i+1;
j=j+1;
k((head(1,i))*(dof+j),:)=zeros(1,dof*2*nNodes);
k(:,(head(1,i)-1)*dof+j)=zeros(dof*2*nNodes,1);
k((head(1,i)-1)*dof+j,(head(1,i)-1)*dof+j)=1.0;
load(1,(head(1,i)-1)*(dof+j))=fload(j+1,i);
if fload(j+1,i) ~= 0
kdof=j;
end
end
end
end
clear RR KEelement;
clear elements;
disps=loads/k;
clear k;
loads k_real.mat
forces=k*disps'
clear k
tip=[476 477 478 479 482 483 488 489 490 491 492 493 494];
% 'tip'contains the node numbers at one tip
avdisp=0;
for i=1:size(tip,2)
tip_disps(i)=disps((tip(i)-1)*dof+kdof);
avdisp=avdisp+disps((tip(i)-1)*dof+kdof);
end
fixed=[1 4 5 16 17 34 35 51 52 72 73 86 87];
% 'fixed'contains the node numbers at the other tip
totforce=0;
for i=1:size(tip,2)
totforce=totforce+forces((tip(i)-1)*dof+kdof);
end
avdisp=avdisp/size(tip,2);
avdisp=avdisp*Lreal; % in meters
tip_disps=tip_disps*Lreal; % in meters
totforce = totforce * EIreal/Lreal^2; % in Newtons
YY=totforce/(TubeDiameter*avdisp/TubeLength*pi) % N/m=Pa*m
Young(it)=YY;
it=it+1;

Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!