gen14.m

%gen fort.14 
%
%help: http://adcirc.org/home/documentation/users-manual-v52/input-file-descriptions/adcirc-grid-and-boundary-information-file-fort-14/
%
dofs = @(p) (p+2)*(p+1)/2;
% script takes in reference compute weights for p = 1
ref.hght = 190;
ref.wdth = 190;
ref.comput = ref.hght*ref.wdth*dofs(1);
% and attempts to compute a rectangular inlet with an aspect ratio of 3
% to match the computational cost at a lower polynomial degree
p = 5;
% The equivalent mesh size would be
side = sqrt(ref.comput / dofs(p));
side = round(side); %to the nearest integer

fprintf('The equivalent mesh is a %4.d x %4.d mesh \n',side,side);
fprintf('Writing the mesh to file...\n');

x = linspace(0,3000,side + 1);
y = linspace(-1500,1500,side + 1);

bath = @(x,y) 18 - 8/3000*x;

fid = fopen('fort.14','w');
fprintf(fid,'Performance mesh for p = %2d\n',p);
%compute num nodes
nnodes = length(x) * length(y);
nele = (length(x) - 1)*(length(y) - 1)*2;

node_array = zeros(2,nnodes);

fprintf(fid,'%8d    %8d\n',nele, nnodes);

%print node locations
it = 1;
for ii = 1:length(x)
    for jj = 1: length(y)
        fprintf(fid,'%8d    %.15e    %.15e    %15e\n',...
            it,x(ii),y(jj),bath(x(ii),y(jj)));
        node_array(1,it) = x(ii);
        node_array(2,it) = y(jj);
        it = it + 1;
    end
end

%figure; hold on;
%scatter(node_array(1,:),node_array(2,:));

fprintf('Writing connectivity table...\n');
%print element connectivities
it = 1;
for ii = 1: length(x)-1
    for jj =  1:length(y) -1
        n2 = length(x)*(ii -1) + jj;
        n1 = n2 + 1;
        n3 = length(x)*ii + jj;
        [n1,n2,n3] = order_nodes(n1,n2,n3);
        fprintf(fid,'%8d  3 %8d %8d %8d\n',it,n1,n2,n3);
        %plot(node_array(1,[n1,n2,n3,n1]),node_array(2,[n1,n2,n3,n1]));
        it = it + 1;
        n2 = length(x)*(ii -1) + jj+1;
        n1 = n2 + length(x);
        n3 = n1 -1;
        [n1,n2,n3] = order_nodes(n1,n2,n3);
        fprintf(fid,'%8d  3 %8d %8d %8d\n',it,n1,n2,n3);
        %plot(node_array(1,[n1,n2,n3,n1]),node_array(2,[n1,n2,n3,n1]));
        it = it + 1;
    end
end
fprintf('Creating %8d edges\n\n',it -1);
fprintf('Writing boundary data...\n');
% print boundary data
fprintf(fid,'0 = Number of open boundaries\n');
fprintf(fid,'0 = Total number of open boundary nodes\n');
fprintf(fid,'1 = Number of land boundaries\n');
fprintf(fid,'%8d =  Total number of land boundary nodes\n',...
    2*length(x) + 2* length(y) - 4);
fprintf(fid,'%8d  10 = Number of nodes for land boundary 1\n',...
    2*length(x) + length(y) - 4);

for ii = 1:length(x)
    fprintf(fid,'%8d\n',(ii -1)*length(x) +1);
end

for ii = 2:length(y)-1
    fprintf(fid,'%8d\n',length(x)*(length(x) - 1) +ii);
end

for ii = 1:length(x)
    jj = length(x) - ii+1;
    fprintf(fid,'%8d\n',(jj -1)*length(x) + length(y));
end

%fprintf(fid,'%8d  12 = Number of nodes for land boundary 2\n',length(x));
for ii = 1:length(x)
    fprintf(fid,'%8d\n',length(x) + 1 - ii);
end

%done
fclose(fid);