MATLAB: Converting FORTRAN code for Finite Element Analysis – Dealing with a Go To Command

finite element analysisfinite element methodfortranfortran to matlabgo toMATLABtriangular elements

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

Best Answer

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
Related Question