setup.m

% Example mesh generation
clc; clearvars
close all;

addpath(genpath('/home/cliu/OneDrive/mesh/OceanMesh2D'));

%% STEP 1: set mesh extents and set parameters for mesh.
S_bnd = shaperead('bnd_latlon.shp');
bbox = cat(2,S_bnd.X', S_bnd.Y');
%bbox=load('../Fairfield_mesh/bnd_latlon.txt');


min_el    = 5.0;  		    % Minimum resolution in meters.
max_el    = 1e3;
dt        = 0.5;              % Ensure mesh is stable at a 2 s timestep
grade     = 0.2;           % Mesh grade in decimal percent.
R         = 3;              % Number of elements to resolve feature width.

%% STEP 2: specify geographical datasets and process the geographical data
%% to be used later with other OceanMesh classes...
coastline = 'Branford_0m_v2';
%dem       = '~/Documents/DATA/New_England_Topobathy_DEM_v1_78_merged.nc';
dem       = '/home/cliu/Documents/DATA/Branford_CoNED_latlon.nc';
% %% The weirs input is a structure with the following options
% %load weirs_struct.mat
% % Weirs is an array of
% %  structs each with fields:
% %              X: [N�1 double] % x georgraphic coordinates of crestline
% %              Y: [N�1 double] % y geographic coordinates of crestline
% %          width: 5 % seperation of front and back face in meters
% %        min_ele: 20 % minimum resolution along faces of weir in meters
% %    crestheight: 5 % in meters
% % weirs(1).X = [
% %   -73.9396
% %   -73.9304
% %   -73.9289];
% % weirs(1).Y = [
% %    40.5756
% %    40.5548
% %    40.5507];
% % weirs(1).width = 10;
% % weirs(1).min_ele = 100;
% % weirs(1).crestheight = 5;
% 
% T=xlsread('../CSTWP001_with gis contours-full plan profile/Plan and Profile Combined-images/digitization.xlsx','coordinates');
% weirs(1).X = T(:,2);
% weirs(1).Y = T(:,3);
% weirs(1).width = 10;
% weirs(1).min_ele = 5;
% T=xlsread('../CSTWP001_with gis contours-full plan profile/Plan and Profile Combined-images/digitization.xlsx','m');
% weirs(1).crestheight = T(:,6)+0.165;
% NOTE: weirs are modeled as "inner" geometry or artifical islands. A thin
% knife edge is added to both tips of the weirs to avoid bad element transistions.
% to nearby non-weir elements.
gdatuw = geodata('shp',coastline,...
               'dem',dem,...
               'bbox',bbox,...
               'h0',min_el);
           
S = shaperead('features.shp');
for i=1:length(S)
    pts{i} = cat(2,S(i).X', S(i).Y');
end


%% STEP 3: create an edge function class
fh = edgefx('geodata',gdatuw,...
            'fs',R,...
            'dt',dt,...
            'max_el',max_el,...
            'Channels',pts(3),...
            'ch',5,...
            'min_el_ch', 2,...
            'g',grade);
%% STEP 4: Pass your edgefx class object along with some meshing options and
% build the mesh...
mshopts = meshgen('ef',fh,'bou',gdatuw,'plot_on',1,'proj','trans');
% now build the mesh with your options and the edge function.
mshopts = mshopts.build;

%% STEP 5: Get fixed constraints and update gdat with overland meshing domain.
muw = mshopts.grd ;
muw = makens(muw,'auto',gdatuw) ; % apply so that extractFixedConstraints only grabs the shoreline constraints.

[pfix,egfix] = extractFixedConstraints(muw) ;

% 10-m contour extracted from the Coastal Relief Model.
coastline = 'Branford_5m';

gdat = geodata('shp',coastline,...
               'dem',dem,...
               'bbox',bbox,...
               'h0',min_el);
gdat.inpoly_flip =  mod(1,gdat.inpoly_flip) ; % if the meshing domain is inverted, you can always flip it .
           
% Here we pass our constraints to the mesh generator (pfix and egfix). 
mshopts = meshgen('ef',fh,'bou',gdat,'plot_on',1,'proj','lambert',...
                  'pfix',pfix,'egfix',egfix);

% now build the mesh with your options and the edge function.
mshopts = mshopts.build;

m = mshopts.grd ;



%% STEP 5: Manually specify open boundaries and weir crest heights
m = make_bc(m,'auto',gdat,'both',[],0.5); % auto make_bc with depth_lim = 5
%m.bd = [];                              % delete land bcs and replace with weirs below

%% STEP 6: interpolate bathy and plot and save the mesh
% gdat_ctdem = geodata('shp',coastline, 'dem','/home/cliu/Documents/DATA/DEM_NewHavenArea_4326.nc',...
%                     'h0',min_el, 'bbox',bbox);
%            
m = interp(m,gdat,'nan','fill'); % interpolate bathy to the
%                 % mesh with fill nan option to make sure corners get values
% m1 = interp(m,gdat_ctdem,'nan','fill');
% m1.b = m1.b * 0.3048;
% dep = m.b;
% dep(dep<=0 & m1.b>0) = m1.b(dep<=0 & m1.b>0);
% dep(m1.b<0) = m1.b(m1.b<0);
% m.b = dep;
% %m = interpFP(m,gdat,muw,gdatuw,-0.1,1);


% fix bathymetry under bridge

innodes = inpolygon(m.p(:,1), m.p(:,2), S(1).X', S(1).Y');
m.b(innodes) = max(0.5, m.b(innodes));
innodes = inpolygon(m.p(:,1), m.p(:,2), S(2).X', S(2).Y');
m.b(innodes) = min(-1.9, m.b(innodes));


% % raise bathymetry on bridge bank
% S = shaperead('../Branford_mesh/beach.shp');
% for i = 1:length(S)
%     innodes = inpolygon(m.p(:,1), m.p(:,2), S(i).X', S(i).Y');
%     m.b(innodes) = min(-3.0, m.b(innodes));
% end

% RENUMBER OB NODES
m = m.renum();
m = renumber_ob(m);

plot(m,'type','bmesh') % plot the bathy on the mesh
plot(m,'type','bd'); % visualize your boundaries
% plot(m,'bmesh'); % plot triangulation and bathy
caxis([-10 0]) ;
write(m,'Branford','14');
save('Branford.mat','m')

disp(['# nodes   : ',num2str(length(m.b))])
disp(['# elements: ',num2str(length(m.t))])
disp(['# ob nodes: ',num2str(m.op.neta)])

addpath(genpath('/home/cliu/Documents/fvcom-toolbox/'))
[x,y] = my_project(m.p(:,1), m.p(:,2), 'forward');
% dep=m.b;
% dep(isnan(dep))=-999;
write_SMS_2dm('Branford.2dm', m.t, m.p(:,1), m.p(:,2),m.b)
write_SMS_2dm('Branford_xy.2dm', m.t, x, y, m.b)


figure()
patch('Vertices',[m.p(:,1),m.p(:,2)], 'Faces',m.t,'Cdata',m.b,'edgecolor','interp','facecolor','interp')
caxis([-3 0]) ;
hold on
triplot(m.t, m.p(:,1),m.p(:,2), 'linewidth',0.1, 'color','r')