% John Sullivan
% ME 515

clear all;


% Reading Input Mesh Information File
%   BE CAREFUL - ELEMENTS HAVE ONLY 3 NODES
%  and the code would continue...

%  MAKE TRIANGLES INTO QUADS.

for L = 1:ne;
    EL(L,4) = EL(L,3);
end

% Calculates the half-bandwidth (hb), diagonal (Diag), 
%         and bandwidth (Band) for any banded matrix.



% Once band issues determined, create matrices and vectors
B = zeros(nn,Band);
Rhs = zeros(nn,1);
Soln = zeros(nn,1);

% Declares local x,y arrays (length=4) and assign gauss values.
XL = zeros(4,1);
YL = zeros(4,1);
GPWeight = [1 1];
GPValue = [-sqrt(1.0/3.0) sqrt(1./3.0)];

% Five nested loops over all elements that calculates every coefficient
% of the banded matrix and the right hand side vector and places them in
% their correctly mapped locations.

for L = 1:ne
    XL(:) = x(EL(L,:));  % For each element put coordinates into
    YL(:) = y(EL(L,:));  %   the local (XL or YL) arrays
    for GP1 = 1:2
        W1 = GPWeight(GP1); % For one direction in perfect space
        XI = GPValue(GP1);  %   get a gauss point location
        for GP2 = 1:2
            W2 = GPWeight(GP2); % And do the same for the other space
            ETA = GPValue(GP2); %  Remember they are Nested!
            % This routine RETURNS the Ns (Basis), the dN/dx (dNx),
            %  the dN/dy (dNy), and the determinant of the Jacobian (DJ)
            %  at the SPECIFIC LOCATION Xi, Eta
            [Basis, dNx, dNy, DJ] = gp2dlf(XL, YL, XI, ETA);
            for i = 1:4
                IRow = EL(L,i);
                for j = 1:4
                    % This maps the columns into a banded matrix
                    JCol = Diag + (EL(L,j) - IRow);
                    % Enter your FEM formulation of the governing eqn
                    %   THIS LINE WOULD HOLD THE SS TUMOR EQUATION
                end
                % Rhs(IRow) = Rhs(IRow) + whatever belongs on right
                %        hand side from each element L
            end
        end
    end
end
% Applying Type 2 Boundary Conditions 
%  In this example they are all zero, so there is no action
for i = 1:nbc2

end

% Applying Type 1 Boundary Conditions
for i = 1:nbc1
    row = ibc1(i);
    B(row,:) = 0;
    B(row, Diag) = 1;
    Rhs(row) = bc1(i);
end

[B Soln] = slvdij(3, B, Rhs, nn, hb);

% Outputs the global, banded matrix and global right-hand side vector.
%  GIVE IT A NAME FOR THE TUMOR OUTPUT INSTEAD!
fid = fopen('TumorSoln.dat','wt');
fprintf(fid, '%12.8f\n', Soln);
fclose(fid);

% %  Create a Vertex Loop
V = zeros(nn,3);
for i = 1:nn
    V(i,1) = x(i);
    V(i,2) = y(i);
    V(i,3) = Soln(i);
end

patch('Faces',EL,'Vertices',V,'FaceVertexCData',Soln,'FaceColor','interp')

% Stream = zeros(nn,1);
% grid = 0.:.01:1.2;
% [X Y] = meshgrid(grid,grid);
% Stream = griddata(x, y, Soln, X, Y);
% contour(X, Y, Stream)