%% Get mesh info
% Example is taken from the quadrilateral duct homework

% 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);
in = zeros(4,ne);
ibc1 = zeros(nbc1);
ibc2 = zeros(3,nbc2);
bc1 = zeros(nbc1);
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);
    in(1,j) = fscanf(fid, '%d',1);
    in(2,j) = fscanf(fid, '%d',1);
    in(3,j) = fscanf(fid, '%d',1);
    in(4,j) = 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(1,j) = fscanf(fid,'%d',1);
    ibc2(2,j) = fscanf(fid,'%d',1);
    ibc2(3,j) = fscanf(fid,'%d',1);
    bc2(j) = fscanf(fid,'%f',1);
end
status = fclose(fid);


%%  Build the vectors for sparse matrix storage

NinL = 4;
neighbor = 12;

Dof = zeros(nn,1) + 3; % For example u, v, w displacements
DofSummed = zeros(nn,1);
DofSummed(1) = Dof(1);
for i = 2:nn
    DofSummed(i) = DofSummed(i-1) + Dof(i-1);
end

% Single DOF example consistent with Conjugate Gradient ICCGLU solver
ia = zeros(nn,1);
ja = zeros(neighbor*nn,1);

ia(1) = 1;

for i = 1:nn
    start = ia(i);
    ja(start) = i;
    extent = start + 1;
    for L = 1:ne
        val_1 = ismember(i, in(:,L));
        if val_1 == 1
            for j = 1:NinL
                jj = in(j,L);
                val_2 = ismember(jj,ja(start:extent));
                if val_2 == 0
                    ja(extent) = jj;
                    extent = extent + 1;
                end
            end
        end
    end
    ia(i+1) = extent;
end

%% Write out values

[fname,pname] = uiputfile('*.mat','Please save the Output Vector file');

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

save(filename,'ia','ja');