triangle_voltage_distribution_fem.m

Voltage Distribution of a Triangle (FEM & FDFD Method)

%%%%%%%%%%%%%%%%%%%%%%% In The Name of God %%%%%%%%%%%%%%%%%%%%%%
% 
% Numerical Method in Electromagnetic
% Dr. M.S Abrishamian
% 
% Mohsen Bahrami Panah        8600374                14/2/87
% 
% Problem 20 -- *** Voltage Distribution of a Triangle (FEM & FDFD) ***

clear all
close all
clc
clf
opengl neverselect;
f1=figure(1);
set(f1,'position',[1 29 1024 672],'color',[1 0.9 1]);
axis off
hold on

L = 1;
DL = 0.2;
NL = L/DL;
Nnod = (L/DL)+1;
x = 0:DL:L;
y = 0:DL:L;
Nn = 0;
for J = 1:Nnod
    for I = 1:Nnod-J+1
        Nn = Nn+1;        
        NodNum(Nnod-J+1,I) = Nn;
    end
end

F = 0;
for J = 1:Nnod
    for I = 1:Nnod-J+1
        F = F+1;                
        X(F) = x(I);
        Y(F) = y(J);
        plot(X(F),Y(F),'oc')
        text(X(F)+DL/10,Y(F),num2str(F),'color','r')
    end
end

Ne = 0;
N = Nnod;
for I = 0:Nnod-2
    N = N-1;
    S = 0;
    for J = 1:2*N-1
        Ne = Ne+1;            
        if rem(J,2)==1
            S = S+1;
            M(Ne,1) = NodNum(Nnod-I,S);
            M(Ne,2) = NodNum(Nnod-I,S+1);
            M(Ne,3) = NodNum(Nnod-I-1,S);
        else
            M(Ne,1) = NodNum(Nnod-I,S+1);
            M(Ne,2) = NodNum(Nnod-I-1,S+1);
            M(Ne,3) = NodNum(Nnod-I-1,S);
        end
    end
end
for I = 1:Ne
    Xn = [X(M(I,1)) X(M(I,2)) X(M(I,3))];
    Xm = sum(Xn(1:end))/3;
    XX(I,:) = Xn;
    Yn = [Y(M(I,1)) Y(M(I,2)) Y(M(I,3))];
    Ym = sum(Yn(1:end))/3;
    YY(I,:) = Yn;
    TRI = delaunay(Xn,Yn);
    figure(1)
    triplot(TRI,Xn,Yn,'m')
    text(Xm,Ym,num2str(I),'color','b')
    axis([0 L 0 L])
    title('\bf\it Triangle Mesh Generation','fontsize',10,'color','b',...
        'BackgroundColor','c')
    pause(0.01)
end
% line([X(1) Y(1)],[X(Nnod) Y(Nnod)],'linewidth',1,'color','m')
% print -djpeg -r500 Triangle_Mesh_Generation.jpg
% pause(3)
% close all


V=zeros(F,1);
% p = zeros(Ne,1);
for I=1:F
    if X(I)==0 || Y(I)==0
        V(I,1) = 0;
        P(I) = 1;
    elseif abs(L-Y(I)-X(I))<=DL/100
        V(I,1) = 100;
        P(I) = 1;
    end
    if I==Nnod || I==F
        V(I,1) = 50;
        P(I) = 1;
    end
end
C = zeros(F,F);
for I = 1:Ne
    p(1) = Y(M(I,2)) - Y(M(I,3));
    p(2) = Y(M(I,3)) - Y(M(I,1));
    p(3) = Y(M(I,1)) - Y(M(I,2));
    Q(1) = X(M(I,3)) - X(M(I,2));
    Q(2) = X(M(I,1)) - X(M(I,3));
    Q(3) = X(M(I,2)) - X(M(I,1));
    A = 0.5 * (p(2)*Q(3) - p(3)*Q(2));
    for J = 1:3
        for K = 1:3
            Ce(J,K) = 1/(4*A)*(p(J)*p(K)+Q(J)*Q(K));
        end
    end
    for i = 1:3
        B = M(I,i);
        for j = 1:3
            G = M(I,j);
            C(B,G) = C(B,G) + Ce(i,j);
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%<<  Starting Iteration Method  >>%%%%%%%%%%%%%%%%%%%%%
Ni = 25;
for i=1:Ni
    for j=1:F
        if P(j)==0
            V1 = 0;
            for k=1:F
                if j~=k
                V1 = V1+V(k,1)*C(j,k);
                end
            end
            if j~=Ne
                V(j,1) = -V1/C(j,j);
            end
        end
    end
end
F = 0;
for j=NL+1:-1:1
    for i=1:j
        F = F+1;
        V_FEM(j,i) = V(F);
    end
end
f2=figure(2);
set(f2,'position',[1 29 1024 672],'color',[1 0.9 1]);
surf(0:DL:L,0:DL:L,V_FEM)
view([1,1,1])
colormap(hsv(500))
title('\bf\it Voltage Distribution of a Triangle (Finite Element Method)',...
    'fontsize',10,'color','b','BackgroundColor','c')
% print -djpeg -r300 Voltage_Distribution_Triangle.jpg
% pause(3)
% close all



V_FDFD = FDFD(Nnod);

fprintf('\n')
fprintf('Node      X         Y       Potential(FEM)    Potential(FDFD)\n')
fprintf('______________________________________________________________')
fprintf('\n')
fprintf('\n')
for I=1:NodNum
    if I>=1 & I<10
        fprintf('0')
    end
    fprintf('%1.0f) ',I)
    fprintf('     %3.2f ',X(I))
    fprintf('    %3.2f ', Y(I))
    fprintf('     V=%8.4f ', V(I))
    fprintf('       V=%8.4f ', V_FDFD(I))
    fprintf('\n')
end


fprintf('\n')
fprintf('\n')
Difference = V' - V_FDFD