finite element method with matrix error
이전 댓글 표시
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
Walter Roberson
2020년 10월 26일
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
2020년 10월 26일
Walter Roberson
2020년 10월 26일
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
2020년 10월 26일
Walter Roberson
2020년 10월 26일
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
2020년 10월 26일
zina shadidi
2020년 10월 26일
Walter Roberson
2020년 10월 27일
(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
2020년 10월 28일
zina shadidi
2020년 10월 28일
zina shadidi
2020년 10월 28일
Walter Roberson
2020년 10월 28일
if size(k,1) < n*D || size(k,2) < n*D
k(n*D, n*D) = 0; %extends k with 0s
end
zina shadidi
2020년 10월 29일
편집: zina shadidi
2020년 10월 29일
Walter Roberson
2020년 10월 29일
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
2020년 10월 29일
답변 (1개)
Asad (Mehrzad) Khoddam
2020년 10월 26일
0 개 추천
Can you show the values of the nodes and elements data?
카테고리
도움말 센터 및 File Exchange에서 Creating and Concatenating Matrices에 대해 자세히 알아보기
2020년 10월 26일
2020년 10월 29일
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
