finite element method with matrix error

2 views (last 30 days)
zina shadidi
zina shadidi on 26 Oct 2020
Commented: zina shadidi on 29 Oct 2020
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
Show 13 older comments
Walter Roberson
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
zina shadidi on 29 Oct 2020
Walter Roberson, I can understand you but I have Null matrices as a results of young modulus , and other variables , and this is a wrong results in physics . In fact I dont know what I will do , may be I should start with another code.
I do appreciate your help and thank you very much.

Sign in to comment.

Sign in to answer this question.

Answers (1)

Asad (Mehrzad) Khoddam
Asad (Mehrzad) Khoddam on 26 Oct 2020
Can you show the values of the nodes and elements data?

Categories

Find more on Logical in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!