Poor quality triangles at shoreline

[crestocean]

Hi there,

After using OceanMesh2D for a few weeks now, I've run into a problem that is hard to resolve. I'm getting poor quality triangles at the shoreline that are preventing me from using the mesh in SWAN.

There seems to be two types of scenarios where these triangles are popping up:

Along an open portion of coast:

At the back of a tributary in an estuary:

My script is below: rng(0) % set the random seed clearvars; clc; close all PREFIX = 'NZ';

%% STEP 1: set mesh extents and set parameters for mesh. bbox1 = [ 162 182 % lon_min lon_max -50 -31 % lat_min lat_max ];

bnd_dxdy = 0.5; [X,Y] = meshgrid(bbox1(1):bnd_dxdy:bbox1(3),bbox1(2):bnd_dxdy:bbox1(4)); k = boundary(X(:),Y(:)); boubox = [X(k),Y(k)]; % half degree boundary nodes

min_el = 500; % minimum resolution in meters. max_el = 28000; % maximum resolution in meters. max_el_ns = 500; % maximum resolution nearshore in meters. grade = 0.15; % mesh grade in decimal percent. slp = 1; % slope

%% STEP 2: specify geographical datasets and process the geographical data % to be used later with other OceanMesh classes... coastline1 = 'new_zealand_OSM_ALL_simple1.shp'; coastline2 = 'new_zealand_OSM_ALL.shp'; dem1 = 'NZ_coarse_gauss2.nc';

gdat1 = geodata('shp',coastline2,'dem',dem1,'bbox',bbox1,'h0',min_el,'boubox',boubox);

%% STEP 3: create an edge function class fh1 = edgefx('geodata',gdat1,'dis',grade,'max_el_ns',max_el_ns,'max_el',max_el,'g',grade,'slp',slp);

%% Repeat STEPS 1-3 for a higher resolution domain: 2 bbox2 = [174.3048516456 -35.3124074509 174.3187290543 -35.3187969099 174.3375567742 -35.3079450223 174.3637227578 -35.2944717390 174.3820823208 -35.2735015111 174.3926375652 -35.2427895558 174.3967704258 -35.2064540920 174.3973056312 -35.1700760583 174.3965488751 -35.1361321469 174.3941375576 -35.1191872100 174.3897990085 -35.1053951482 174.3840131484 -35.0983026491 174.3719624723 -35.0876642518 174.3473880492 -35.0750550149 174.3001969791 -35.0608578554 174.2333169884 -35.0474094562 174.1813786870 -35.0402455026 174.1342545847 -35.0362148980 174.0746108698 -35.0360649098 174.0442844617 -35.0403078639 174.0211515661 -35.0485061241 173.9999198601 -35.0602519850 173.9844460499 -35.0743789825 173.9727946865 -35.0932454961 173.9577253524 -35.1187949662 173.9565010140 -35.1664552444 173.9528809465 -35.2116113066 173.9531107459 -35.2557352412 174.0131040540 -35.3221200441 174.0951310279 -35.3407955428 174.1860181153 -35.3404633683 174.2622981938 -35.3405650531 174.2874217242 -35.3232641623 174.3048516456 -35.3124074509 ];

min_el = 100; % minimum resolution in meters. max_el = 4000; % maximum resolution in meters. max_el_ns = 200; % maximum resolution nearshore in meters. grade = 0.15; % mesh grade in decimal percent. R = 6; % number of elements to resolve feature width.

dem2 = '2_main.nc';

gdat2 = geodata('shp',coastline2,'dem',dem2,'h0',min_el,'bbox',bbox2); fh2 = edgefx('geodata',gdat2,'fs',R,'max_el_ns',max_el_ns,'max_el',max_el,'g',grade);

%% STEP 4: Pass your edgefx class object along with some meshing options and % build the mesh... mshopts = meshgen('ef',{fh1,fh2},... 'bou',{gdat1,gdat2},... 'pfix',boubox,... 'plot_on',1,... 'nscreen',5); 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 nodestrings m = mshopts.grd; m = make_bc(m,'auto',gdat1,'distance'); % make the boundary conditions plot(m,'type','bd');

% save the grid as a fig file savefig(gcf,sprintf('%s_mesh.fig',PREFIX))

%% Interp dem to grid m = interp(m,{dem1,dem2},'type','all'); m = lim_bathy_slope(m,0.05,0); m.b = m.b * -1; plot(m,'type','b');

%% Save mesh files % Save as a msh class save(sprintf('%s_msh.mat',PREFIX),'m'); % Write an ADCIRC fort.14 compliant file to disk. write(m,sprintf('%s_mesh',PREFIX)) % rename the bathy file to fort.14 copyfile(sprintf('%s_mesh.14',PREFIX), 'fort.14')

Any help will be greatly appreciated.

Thanks, Sam

[krober10nd]

Use the high fidelity option in geodata to improve the element quality along the boundary.

https://github.com/CHLNDDEV/OceanMesh2D/blob/Projection/Examples/Example_13_New_Zealand_High_Fidelity.m

Make sure you execute msh.clean to remove singly connected elements. See the docstring for it.

[WPringle]

Yes, using the msh.clean with options with higher quality tolerance for "bad" boundary elements and removing singly connected elements would fix both these issues.

[crestocean]

Thank you for these suggestions. The mesh clean looks to have solved the issue, however it seems to have introduced a new problem which I think is related to the singly connected elements flag in the cleaning step. My call looks like this:

m = mshopts.grd; m = m.clean('db',0.35,'pfix',boubox,'sc_maxit',50); m = make_bc(m,'auto',gdat1,'distance'); % make the boundary conditions

SWAN is throwing the error:

Terminating error: SwanBpntlist: list of boundary vertices could not be completed

Is 'sc_maxit' what you mean by 'removing singly connected elements'? I don't get this error when I use the same boundary files on a simpler mesh that doesn't contain a higher resolution domain, but is otherwise built with the same parameters (and doesn't include the cleaning step).

When I omit the 'sc_maxit' flag SWAN throws a severe error:

forrtl: severe (157): Program Exception - access violation

I know that SWAN errors should probably be discussed on other forums but I'm hoping this is an issue that you still might be able to help me with as I think it is related to the arguments passed to the clean function.

Thank you.

[krober10nd]

sc_maxit = inf