finite element method with matrix error
Show older comments
I am trying to calculate energy using the folowing code , the matrices elements2s, nodes2s was the results of another script.
clc
clear all
% Initial data:
TubeLength = 3.38e-9; %m
TubeDiameter = 1.011e-9; %m
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 in MN 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 DOF):
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:
it=1;
while it<numIterations+1
load elements2s.txt
load nodes2s.txt
nNodes2s= 199;
nElements2s = 199;
nodesDef=zeros(1,nNodes2s);
k=zeros(nNodes2s*D,nNodes2s*D);
% Defects:
defectsPercentage=1; %per cent
numDef=round((nElements2s*defectsPercentage)/100);
for i=1:numDef
def=round((nElements2s-1)*rand)+1; %The addition of 1 and -1 is to avoid zero as the resulting random number. elements(def,:)=[];
nElements2s=nElements2s-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:nElements2s-1
n = elements2s(i,1);
m = elements2s(i,2);
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 = RR'*KBond*RR;
n = elements2s(i,1);
m = elements2s(i,2);
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(1:D, 1: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;
save k_real.mat k;
loads=zeros(1,D*nNodes);
for i=1:size(head(1,:),2)
for j=1:D
if head(j+1,i)~=0 % displacement prescribed
for m=1:D*nNodes
loads(m)=loads(m)-k(m,(head(1,i)-1)*D+j)*fload(j+1,i);
end
k((head(1,i)-1)*D+j,:)=zeros(1,D*nNodes);
k(:,(head(1,i)-1)*D+j)=zeros(D*nNodes,1);
k((head(1,i)-1)*D+j,(head(1,i)-1)*D+j)=1.0;
loads(1,(head(1,i)-1)*D+j)=fload(j+1,i);
if fload(j+1,i) ~= 0
kdof=j;
end
else % force prescribed
loads(1,(head(1,i)-1)*D+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;
load 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)*D+kdof);
avdisp=avdisp+disps((tip(i)-1)*D+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)*D+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;
end
matlab writng an error which is
{Index in position 1 is invalid. Array indices must be positive integers or logical values.
Error in Untitled (line 99)
k((m-1)*dof+1:m*dof,(m-1)*dof+1:m*dof)=k((m-1)*dof+1:m*dof, (m-1)*dof+1:m*dof)+KElement(1:dof, 1:dof);
>> }
can anyone help me please
many thanks in advance
15 Comments
Walter Roberson
on 26 Oct 2020
m = elements2s(i,2);
As outside observers, we have no reason to know that is going to be a positive integer.
(m-1)*dof+1:m*dof
If m is not an integer with minimum value 1, then (m-1)*dof+1 would not be an integer.
zina shadidi
on 26 Oct 2020
Walter Roberson
on 26 Oct 2020
That code would be a problem if n was ever 0 or negative.
Give the command
dbstop if error
and run the code. When the code stops, examine the values of n and D and any other variables mentioned on the line it is failing on.
zina shadidi
on 26 Oct 2020
Walter Roberson
on 26 Oct 2020
dbstop if error
and run the code. Then show us the output of
disp(n)
disp(D)
disp(size(k))
disp(m)
disp(dof)
zina shadidi
on 26 Oct 2020
zina shadidi
on 26 Oct 2020
Walter Roberson
on 27 Oct 2020
(7-1)*3+1:7*3 -> 19:21 . So by the time you reach n = 7, you would need k to be at least 21 x 21 in order to be able to index k(19:21, 19:21) on the right hand side of the assignment. But your k is only 18 x 18.
Your k needs to be at least max(n) * D on each side for the indexing to work.
If you want to automatically extend k with zeros if needed, then before that assignment line you could put in:
if size(k,1) < n*D | size(k,2) < n*D
k(n*D, n*D) = 0; %extends k with 0s
end
zina shadidi
on 28 Oct 2020
zina shadidi
on 28 Oct 2020
zina shadidi
on 28 Oct 2020
Walter Roberson
on 28 Oct 2020
if size(k,1) < n*D || size(k,2) < n*D
k(n*D, n*D) = 0; %extends k with 0s
end
zina shadidi
on 29 Oct 2020
Edited: zina shadidi
on 29 Oct 2020
Walter Roberson
on 29 Oct 2020
sz = size(k);
if sz(1) < n*d
k(sz(1)+1:n*d, :) = k(sz(1),:) + (1:n*d-sz(1));
end
if sz(2) < n*d
k(:, sz(2)+1:n*d) = (1:n*d-sz(2)).' + k(:, sz(2));
end
You only talked about extending according to the last value on each "line", but your matrices are short on both sides, so this code also extends by adding 1 along columns as needed.
I think you are talking the wrong approach: I think you should be initializing to the maximum size.
zina shadidi
on 29 Oct 2020
Answers (1)
Asad (Mehrzad) Khoddam
on 26 Oct 2020
0 votes
Can you show the values of the nodes and elements data?
Categories
Find more on Creating and Concatenating Matrices in Help Center and File Exchange
on 26 Oct 2020
on 29 Oct 2020
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
