//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.

// Clear any previous data in the memory of scilab console of your computer 
clc;
clear;

n=11; // No. of nodes
A=zeros(n,n); // Global coefficient matrix corresponding to Phi'' term (del^2)
D=zeros(n,n); // Global coefficient matrix corresponding to Phi term 
B=zeros(n,1); // Global coefficient matrix corresponding to column vector (for source/ boundary conditions)

// For simplicity Phi (unknown nodal potentials) variable is denoted with U

x=0:1/(n-1):1; // Discretizing into n-1 segments or elements, for example: if n =11, no. of segments: 1/(11-1) = 1/10 = 0.1

// Formation of connectivity matrix (t)
for i=1:n-1
    t(1,i)=i; // for example: element 1 is between nodes 1 and 2, element 2 between nodes 2 and 3, and so on
    t(2,i)=i+1;
end
//t:
// start node 1 2 3 5 ....
// end node   2 3 4 6 ....

for element = 1:n-1
    // Coordinates of the 2 nodes of each element 
    xn(1)=x(t(1,element));
    xn(2)=x(t(2,element));
    // for element 1: Xn(1) = 0, Xn(2) = 1/10, element 2: Xn(2) = 2/10, Xn(3) = 3/10, and so on 
    L=abs(xn(2)-xn(1)); // length of each element
    // formation of element level matrices
    for i=1:2 
        for j=1:2
            if(i==j) then // Diagonal element of 'a' matrix 
                a(element,i,j)=1/L;
            else // Off-diagonal element of 'a' matrix
                a(element,i,j)=-1/L;
            end
        end
    end
    d(element,:,:) = (L/6)*[2 1;1 2]; // Formation of element level matrix corresponding to Phi term 
    // Formation of element level matrix corresponding to source/boundary conditions
    b(element,1)=((xn(2)*((xn(2)^2)-(xn(1)^2))/2)-(((xn(2)^3)-(xn(1)^3))/3))/L; // for node 1 of each element
    b(element,2)=(-(xn(1)*((xn(2)^2)-(xn(1)^2))/2)+(((xn(2)^3)-(xn(1)^3))/3))/L; // for node 2 of each element
end

// Formation of global level matrices using the connectivity matrix
for element=1:n-1
    nodes=t(1:2,element); // global node numbers of each element 
    for i=1:2
        for j=1:2
        // 1, 2 are local node number of each element
            A(nodes(i),nodes(j))=a(element,i,j)+A(nodes(i),nodes(j)); 
            D(nodes(i),nodes(j))=d(element,i,j)+D(nodes(i),nodes(j)); 
            // 2 X 2 element matrices
        end
        B(nodes(i))=b(element,i)+B(nodes(i)); // 2 X 1 element matrices
    end
end

K=A-D; // A*U - D*U = B, (A-D)*U = B, K = A-D 

// Imposing boundary conditions
K(1,:)=zeros(1,n); // for the 1st node
K(1,1)=1;
B(1)=0;

K(n,:)=zeros(1,n); // for the end node
K(n,n)=1;
B(n)=0;

// for example
// 1 -1 0 0   U1 = B1
// -1 2 -1 0  U2 = B2
// 0 -1 2 -1  U3 = B3
// 0 0 -1 1   U4 = B4

// converts to 

// 1 0 0 0   U1  = 0   --> U1 = 0 
// -1 2 -1 0  U2 = B2
// 0 -1 2 -1  U3 = B3
// 0 0 0 1   U4  = B4  --> U4 = 0

U=K\B;

plot(x',U,'r') // here 'r' denotes the color of the curve is in red
// x' to make it column matrix since U is a column matrix 
//-----------------------------------------------------------------------
// Analytical solution or Exact solution
xan = 0:1/100:1; 
yan = (sin(xan)/sin(1))-xan;
plot(xan,yan,'b') // here 'b' denotes the color of the curve is in blue

y1=10*(xan.*(1-xan))/36;
plot(xan,y1,'k')

hl=legend(['1D FEM (10 elements)';'Exact solution';'2nd order approximation']);