zabaqer asked . 2020-11-16

how to add matrices to global matrix

I have a problem with matrices in this code clc 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; 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 ; 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 ; 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 ]; 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; 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 ]*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 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; else % force prescribed load(1,(head(1,i)-1)*dof+j)=fload(j+1,i); if fload(j+1,i) ~= 0 kdof=j; end 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; the error was Index in position 2 is invalid. Array indices must be positive integers or logical values. Error in Untitled1111111 (line 129) load(m)=load(m)-k(m,(((head(1,i)-1)*dof+j)*fload(j+1,i)));

numerical analysis