//The scilab code is developed by 
//B. Sai Ram, Research Scholar, Dept. of Electrical Engineering, IIT Bombay
//for the MOOC course: Electrical Equipment and Machines: Finite Element Analysis
//by Prof. S. V. Kulkarni, IIT Bombay.

clc;
clear;

n=51; // number of nodes along x-axis
m=51; // number of nodes along y-axis
n_nodes=(n*m); // total number of nodes in the domain 


disp('intializing global coefficent matrix');
C=zeros(n_nodes,n_nodes);
BJ=zeros(n_nodes,1);

x=0:0.1/(n-1):0.1; // nodal points on x-axis 
y=0:0.1/(m-1):0.1; // nodal points on y- axis
// this forms a square with 0.1 X 0.1 dimensions in MKS system

disp('forming p matrix');
// formation of p matrix
// Given domain is discritized into squares and then each square is
// divided into 2 trainguler elements using diagonals of square
h=1; // h indicates node number (column index of p matrix)
for i=1:m // index for y
    for j=1:n // index for x
        // for each y 51 x coordinates are formed
        p(1,h)=x(j);
        p(2,h)=y(i);
        h=h+1;
    end
end

disp('----------------------------------------------------------------');
disp('forming t matrix');
// formation of t matrix, 3 gobal node numbers for each element and sub-domain number
h=1;// here, h indicates element number (column index of t matrix)
// forming t matix, of an element with an edge parallel to x-axis, 
for i=1:m-1
    for j=1:n-1
        t(2,h)=j+((i-1)*n); // for i = 1, j =1 , h =1, n=51, t(1,1) = 1
        t(3,h)=j+1+((i-1)*n);// t(2,1) = 2
        t(4,h)=j+1+((i)*n);// t(3,1) = 53
        // x coordinates of each node of the element
        Xc(1) = p(1,j+((i-1)*n));
        Xc(2) = p(1,j+1+((i-1)*n));
        Xc(3) = p(1,j+1+((i)*n));
        // y coordinates of each node of the element
        Yc(1) = p(2,j+((i-1)*n));
        Yc(2) = p(2,j+1+((i-1)*n));
        Yc(3) = p(2,j+1+((i)*n));
        // calculating the centroid of the element
        cX = sum(Xc)/3;
        cY = sum(Yc)/3;
        if (cX<=0.07 && cX>=0.03) && (cY<=0.07 && cY>=0.03) then
            t(1,h)=2; // same material everywhere
        else
            t(1,h)=1; // same material everywhere
        end
        h=h+1;// incrementing the element number
    end
end

for i=1:m-1
    for j=1:n-1

        t(2,h)=j+((i-1)*n);// for i = 1, j =1 , h =1, n=51, t(1,1) = 1
        t(3,h)=j+1+((i)*n);// t(2,1) = 53
        t(4,h)=j+((i)*n);// t(3,1) = 52
        // x coordinates of each node of the element
        Xc(1) = p(1,j+((i-1)*n));
        Xc(2) = p(1,j+1+((i)*n));
        Xc(3) = p(1,j+((i)*n));
        // y coordinates of each node of the element
        Yc(1) = p(1,j+((i-1)*n));
        Yc(2) = p(1,j+1+((i)*n));
        Yc(3) = p(1,j+((i)*n));
        // calculating the centroid of the element
        cX = sum(Xc)/3;
        cY = sum(Yc)/3;
        if (cX<=0.07 && cX>=0.03) && (cY<=0.07 && cY>=0.03) then
            t(1,h)=2; // same material everywhere
        else
            t(1,h)=1; // same material everywhere
        end
        h=h+1;// incrementing the element number
    end
end

disp('----------------------------------------------------------');
[A3 n_elements]=size(t);// A3=4, 
//size of t matrix is 4 X no. of elements and 4 gets assigned to A3
[A1 n_nodes]=size(p);// A1=2
//size of t matrix is 2 X no. of nodes and 2 gets assigned to A1
disp('number of elements');
disp(n_elements);
//for n = 51, m = 51, n_elements = 5000
disp('number of nodes');
disp(n_nodes);
//for n = 51, m = 51, n_nodes = 2601

//Forming element coeffcient matrices
disp('forming element coefficient matrix');
for element=1:n_elements
    nodes=t(2:4,element);// 1st three rows --> 3 global node numbers of each element
    //'nodes' is a 3 X 1 matrix containing global
    // node numbers of element under consideration  
    Xc=p(1,nodes');  //x-coordinates of nodes, Xc is a 3 X 1 matrix
    Yc=p(2,nodes');  //y-coordinates of nodes, Yc is a 3 X 1 matrix
    P=zeros(3,1);
    Q=zeros(3,1);
    P(1)=Yc(2)-Yc(3);
    P(2)=Yc(3)-Yc(1);
    P(3)=Yc(1)-Yc(2);
    Q(1)=Xc(3)-Xc(2);
    Q(2)=Xc(1)-Xc(3);
    Q(3)=Xc(2)-Xc(1);
    //finding x and y coordinates 
    //of the centroid 
    //see if the element lies within the 
    //rectangular conductor 
    cX = sum(Xc)/3;
    cY = sum(Yc)/3;
    if t(1,element) == 2 then
        Mu(element) = 4*%pi*1e-7;
        bj = 1e3; // J--> current density in the conductor
    else
        Mu(element) = 4*%pi*1e-7;
        bj = 0;
    end
    //eq. (15.21c) of 
    //M. N. O. Sadiku and S. V. Kulkarni, Principles of Electromagnetics,
    //Sixth Edition, Oxford University Press, India, 2015
    // area of triangular element
    delta= 0.5*abs((P(2)*Q(3))-(P(3)*Q(2)));//Absolute value is taken since 
    //the three nodes may not have been numbered in the anticlockwise 
     for i=1:3
        for j=1:3
            // 3 X 3 element coefficient matrix 
            c(element,i,j)=((P(i)*P(j))+(Q(i)*Q(j)))/(4*delta*Mu(element));
        end
     end
     // 3 X 1 element level source matrix
     be(element,:) = bj*(delta/3)*[1;1;1];
end

disp('forming global coefficient matrix');
for element=1:n_elements
    nodes=t(2:4,element);
    // Formation of the global coefficient matrix 
    for i=1:3
        for j=1:3            
            C(nodes(i),nodes(j))=C(nodes(i),nodes(j))+c(element,i,j);
        end
        BJ(nodes(i),1) = be(element,i)+ BJ(nodes(i),1);
    end
end

disp('applying boundary conditions');
// e matrix is not formed hence boundary conditions are applied manually
// 1st edge, bottom edge
for i=1:n
    C(i,:)=zeros(1,n_nodes);
    C(i,i)=1;
    BJ(i,1)=0;
end
// 2nd edge, vertical left side edge
for i=n+1:n:(n*m)-2*n+1
    C(i,:)=zeros(1,n_nodes);
    C(i,i)=1;
    BJ(i,1)=0;
end
// 3rd edge, vertical right side edge
for i=2*n:n:(n*m)
    C(i,:)=zeros(1,n_nodes);
    C(i,i)=1;
    BJ(i,1)=0;
end
// 4th edge, top edge
for i=(n*(m-1))+1:(n*m)-1
    C(i,:)=zeros(1,n_nodes);
    C(i,i)=1;
    BJ(i,1)=0;
end

disp('computing nodal potentials');
A=C\BJ;

//Here, A is a column vector - it is converted into n X m matrix,
//similarly, two rows of p matrix (x and y coordinates) are converted
//into n X m matrix
Vn=matrix(A(1:n*m),n,m);//potentials at each node of n X m grid
XX=matrix(p(1,1:n*m),n,m);// x coordinates of n X m grid
YY=matrix(p(2,1:n*m),n,m);// y coordinates of n X m grid

figure(1)
surf(XX,YY,Vn')

figure(2)
contour(x,y,Vn,25)