% John Sullivan
% ME 515

clear all;


% Reading Input Mesh Information File
% John Sullivan
% ME 515
% This example opens, reads, and then closes 
% 03/21/10
% This example prompts the user for an input file
%     opens, reads, and then closes it
% The file expects to be formated as in "duct_quad.dat" 
% it reads the first few variables as integers ("%d")
% then within the various for loops it reads in the coordinates
% element connectivities and boundary conditions
% the "j = fscanf( etc." simply is used to match the format 
% of the input file.  Having the node numbers within the input 
% file helps when looking at just the input file listing for problems.
% Reads in space delimited data files and converts them to a matrix.

% Reading Input Mesh Information File
[fname,pname] = uigetfile('*.dat','Please select the Input Mesh file');

if (fname ~= 0)
    filename = sprintf('%s%s',pname,fname);
else
    disp('No file selected...Program will quit');
    return;
end

fid = fopen(filename);
nn = fscanf(fid, '%d', 1);
ne = fscanf(fid, '%d', 1);
nbc1 = fscanf(fid, '%d', 1);
nbc2 = fscanf(fid, '%d', 1);
x = zeros(nn,1);
y = zeros(nn,1);
EL = zeros(ne,4);
ibc1 = zeros(nbc1); bc1 = zeros(nbc1);
ibc2 = zeros(nbc2); bc2 = zeros(nbc2);
for i = 1:nn
    j = fscanf(fid, '%d',1);
    x(j) = fscanf(fid, '%f',1);
    y(j) = fscanf(fid, '%f',1);
    z = fscanf(fid, '%f',1);
end
for i = 1:ne
    j = fscanf(fid, '%d',1);
    EL(j,1) = fscanf(fid, '%d',1);
    EL(j,2) = fscanf(fid, '%d',1);
    EL(j,3) = fscanf(fid, '%d',1);
    EL(j,4) = fscanf(fid, '%d',1);
    mat = fscanf(fid, '%d',1);
end
for i = 1:nbc1
    j = fscanf(fid, '%d',1);
    ibc1(j) = fscanf(fid,'%d',1);
    bc1(j) = fscanf(fid,'%f',1);
end
for i = 1:nbc2
    j = fscanf(fid, '%d',1);
    ibc2(j) = fscanf(fid,'%d',1);
    bc2(j) = fscanf(fid,'%f',1);
end
status = fclose(fid);

%  and the code would continue...


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

fid = fopen('DuctQuadSoln.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)