% Mesh Generation Control File Gold Coast/ Northern NSW
clearvars; clc;
addpath(genpath('utilities/'))
addpath(genpath('datasets/'))
addpath(genpath('m_map/'))
%% STEP 1: Set mesh extents and set parameters for mesh.
min_el    = 100;  	% Minimum resolution in meters.
max_el    = 7500; 	% Maximum resolution in meters. 
max_el_ns = 400;     % Maximum resolution nearshore in meters.
dt        = 2;          % Encourage mesh to be stable at a 1 s timestep
grade     = 0.25; 	% Mesh grade in decimal percent.
R         = 2;    	% Number of elements to resolve feature width.    

min_el_in    = 25;  	% Minimum resolution in meters.
max_el_in    = 4500; 	% Maximum resolution in meters. 
max_el_ns_in = 50;        % Maximum resolution nearshore in meters.
dt_in        = 0;          % Encourage mesh to be stable at a 1 s timestep
grade_in     = 0.25; 	% Mesh grade in decimal percent.
R_in         = 3;    	% Number of elements to resolve feature width.
slp_in       = 5; 
%% STEP 2: Specify geographical datasets and process the geographical data 
%% to be used later with other OceanMesh classes.
coastline = 'coast_contour_dissolve_v3';
dem       = 'bathy_final_v2.nc';
bbox = [153.35, 154.0;
         -28.90, -27.1]
     
bbox_inner = [153.35, 153.8;
         -28.6, -27.3]     
     

%%
gdat = geodata('shp',coastline,...
               'dem',dem,...
               'bbox',bbox,...
               'h0',min_el); 
           
%%
plot(gdat, 'dem')
%% STEP 3: Create an edge function class.
fh = edgefx('geodata',gdat,...
            'fs',R,'max_el_ns',max_el_ns,...
            'max_el',max_el,'dt',dt,'g',grade);
%%
gdat_inner = geodata('shp',coastline,...
               'dem',dem,...
               'bbox',bbox_inner,...
               'h0',min_el_in); 
%% STEP 3: Create an edge function class.
fh_inner = edgefx('geodata',gdat_inner,...
            'fs',R,'max_el_ns',max_el_ns_in,...
            'max_el',max_el_in,'dt',dt_in,'g',grade_in, 'slp',slp_in);
        
%% STEP 4: Pass your edgefx class object along with some meshing options and
% build the mesh...
mshopts = meshgen('ef',{fh,fh_inner},'bou',{gdat,gdat_inner},'proj','utm','plot_on',1 );
mshopts = mshopts.build; 
%% STEP 5: Plot it and write a triangulation fort.14 compliant file to disk.
% Get out the msh class and put on bathy and nodestrings
m = mshopts.grd;

m = interp(m,gdat); % interpolate bathy to the mesh with minimum depth of 1 m
m = makens(m,'auto',gdat);  % make the nodestring boundary conditions 
                                 % with min depth of 5 m on open boundary
plot(m,'bd'); plot(m,'blog');    % plot triangulation, and bathy on log scale
write(m,'Gold_Coast_v6');                % write to ADCIRC compliant ascii file