Converting FORTRAN code for Finite Element Analysis - Dealing with a Go To Command

9 ビュー (過去 30 日間)
Hello everyone,
I'm working on a problem and trying to come to grips with finite element analysis. I'm working from "A Simple Introduction to Finite Element Electromagnetic Problems", MATTHEW N. 0. SADIKU, IEEE TRANSACTIONS ON EDUCATION, VOL. 32, NO. 2, MAY 198. In it some sample FORTRAN code is provided for one of the example problems. My attempt is then to translate this program into MATLAB however I'm struggling to deal with the Go To statements.
The FORTRAN program is:
EDIT: Updated Attempt At Code:
%Nodes:
N(1,:) = [0,0];N(2,:) = [0.2,0];N(3,:) = [0.4,0];N(4,:) = [0.6,0];N(5,:) = [0.8,0];N(6,:) = [1.0,0];
N(7,:) = [0,0.2];N(8,:) = [0.2,0.2];N(9,:) = [0.4,0.2];N(10,:) = [0.6,0.2];N(11,:) = [0.8,0.2];
N(12,:) = [0,0.4];N(13,:) = [0.2,0.4];N(14,:) = [0.4,0.4];N(15,:) = [0.6,0.4];N(16,:) = [0,0.6];
N(17,:) = [0.2,0.6];N(18,:) = [0.4,0.6];N(19,:) = [0,0.8];N(20,:) = [0.2,0.8];N(21,:) = [0,1];
X = N(:,1); Y = N(:,2);
%Triangles:
TRI = [1,2,7;2,8,7;2,3,8;3,9,8;3,4,9;4,10,9;4,5,10;5,11,10;5,6,11;7,8,12;...
8,13,12;8,9,13;9,14,13;9,10,14;10,15,14;10,11,15;12,13,16;13,17,16;...
13,14,17;14,18,17;14,15,18;16,17,19;17,20,19;17,18,20;19,20,21];
%Fixed (perscribed) Potentials [Vf,Node Number]
Vp = [0,1;0,2;0,3;0,4;0,5;50,6;100,11;100,15;100,18;100,20;50,21;0,19;0,16;0,12;0,7];
ND = size(N,1); %Number of Nodes
NE = size(TRI,1); %Number of Elements
NP = size(Vp,1); %Number of Fixed Nodes
NF = ND-NP; %Number of Free Nodes
%Preallocations:
Cg = zeros(ND,ND);
B = zeros(ND,1); %-[Cfp][Vf] rearranged into a matrix with 1's at p nodes
for I = 1:NE %Over Each Element {I}
nodes = TRI(I,:);
Xe = X(nodes); Ye = Y(nodes);
%Local Coefficient Matrix:
P(1) = Ye(2) - Ye(3); Q(1) = Xe(3) - Xe(2);
P(2) = Ye(3) - Ye(1); Q(2) = Xe(1) - Xe(3);
P(3) = Ye(1) - Ye(2); Q(3) = Xe(2) - Xe(1);
A{I} = 0.5.*(P(2)*Q(3) - P(3)*Q(2)); %Saves Area of Specific Element
for i = 1:3
for j = 1:3
C(i,j) = (1./(4*A{I})).*(P(i).*P(j) +Q(i).*Q(j));
%Calculates Individual Element Coefficient Matrix
end
end
clearvars i j
%Tested & Working To Here
%Global Coefficient Matrix:
for j = 1:3 %Runs over local node numbers for rows:
for L = 1:3
IR = TRI(I,j); %Row Index
IC = TRI(I,L); %Column Index
for k = 1:NP
if IR == Vp(k,2) %If Row is at a fixed node:
Cg(IR,IR) = 1;
B(IR,1)=Vp(k,1);
end
end
for k = 1:NP
if IC == Vp(k,2) %If Column is at a fixed node:
B(IR,1) = B(IR,1) - C(j,L)*Vp(k,1);
%if IC & IR statements are not true do this
end
end
if IR ~= Vp(k,2) && IC ~= Vp(k,2)
Cg(IR,IC) = Cg(IR,IC) + C(j,L);
end
end
end
end
Vsol = Cg \ B;
The code above does not reproduce the expected potentials when the inverse is performed (V = Cg \ B). The resulting B matrix is however tested and correct. The issue comes therefore in the calculation of the global coefficient matrix, Cg.
For rows where IR == Vp(k,2) is true, the row
Cg(IR,:) = [0,...,Irth Position, 0,,,] e.g IR = 1, Cg(1,:) = [1,0,0,0] for a system with 4 nodes
And the program appears to set these 1's correctly but I'm also getting non-zero values where I don't expect them in the resulting matrix.
Thanks for the time
  2 件のコメント
Witold Waldman
Witold Waldman 2022 年 4 月 18 日
編集済み: Witold Waldman 2022 年 4 月 18 日
Is the completed MATLAB code for the sample FORTRAN program available?
Ben Barrowes
Ben Barrowes 2022 年 4 月 19 日
If you did have the compilable fortran code, I could use my f2matlab (fileexchange) then remgoto (homegrown tool) on your code (refactors all goto's).
Helping you is complicated by the fact that you don't use the same variable names and structure as the original code. But in general, "while", "break", and "continue" in matlab can reproduce the behavior of 100% of legal fortran goto's. PM me for more help. Otherwise, good luck!

サインインしてコメントする。

採用された回答

Walter Roberson
Walter Roberson 2019 年 9 月 23 日
The general pattern for needing to exit two layers of looping is:
need_to_exit_loop = false;
for first_variable = start_value : end_value
for second_variable = start_second : end_second
if condition_to_leave
need_to_exit_loop = true;
break;
end
end
if need_to_exit_loop
break;
end
end
  2 件のコメント
ADSW121365
ADSW121365 2019 年 9 月 23 日
編集済み: ADSW121365 2019 年 9 月 23 日
I'm not sure if this helps me or not. My understanding maybe the issue. Think the break command maybe the way forward I'm just unsure of the specific application.
for I = 1:NE %Over Each Element {I}
%Calculate Local Coefficient Matrix First C = size([ND x ND])
%Global Coefficient Matrix:
for j = 1:3 %Runs over local node numbers for rows:
IR = TRI(I,j);
for k = 1:NP
if IR == Vp(k,2)
%Do Calculations:
Cg(IR,IR) = 1;
B(IR,1) = Vp(k,1);
%At this point should break such j loop is closed I THINK
%k loop should run to completion such Cg, B are assigned for all IR
end
end
%if above isn't met for rows begin next loop for columns
% (i.e this only runs when IR =~ Vk(2,k))
for L = 1:3
IC = TRI(I,L);
for k = 1:NP
if IC == Vp(k,2)
B(IR,1) = B(IR,1) - C(j,L).*Vp(k,1);
%Again should break to close j loop
%k loop should run to completion such B is assigned for all IC
end
end
%if IC & IR statements are not true do this
Cg(IR,IC) = Cg(IR,IC) + C(j,L);
end
end
end
ADSW121365
ADSW121365 2019 年 9 月 23 日
The more I think about it, the more I think it wants me to skip to the second from final end ( the one closing the j loop) if either condition is met, after it does the calculation but I'm honestly not certain as FORTRAN is a language I have zero experience in.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeFortran with MATLAB についてさらに検索

製品


リリース

R2018b

Community Treasure Hunt

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

Start Hunting!

Translated by