Taking mesh right up to bbox

[edeyejedi]

Hi all, First post, and just looking for a pointer. I guess it may help someone else at some point as well. I cant seem to get the mesh so push right out to the specified bounding box - well the corners particularly. Understand this could be a limitation, but perhaps there are some settings/approaches that help with this.
Thanks for any input

ps. been experimenting with OceanMesh2D for the last 72 hours and find it brilliant and fascinating. Thanks all that have contributed.

[edeyejedi]

Example of the corner missing But also an empty value in the bathy output :( Not sure why that is happening - changes with resolutions, with high res resulting in more "holes". Bunch of missing vlaues on the southern boundary

[krober10nd]

Hey @edeyejedi thanks for the post! Could you please post the code here in a code formatted block so I can help you figure this out? Thank you!

[edeyejedi]

Hi Keith.

Thanks for your reply. I will look at putting something together for you. Just quickly, I am searching through the code as best I can trying to find the point at which triangles with angles < 5 or > 175 degrees are removed every iteration – SWAN has a couple of limitations on unstructured grids:

• The number of triangles that meet at each vertex inside the mesh should not be smaller than 4 or larger than 10.
• Angles inside each triangle should not be higher than 143 degrees. I think OM2D has the first point covered. Just need to change the upper angle limit. As per my first message, this package is awesome. Thanks for sharing it with everyone.

Regards Ed

From: Keith Roberts @.> Sent: Monday, 22 March 2021 10:25 pm To: CHLNDDEV/OceanMesh2D @.> Cc: Ed Atkin @.>; Mention @.> Subject: Re: [CHLNDDEV/OceanMesh2D] Taking mesh right up to bbox (#201)

Hey @edeyejedihttps://github.com/edeyejedi thanks for the post! Could you please post the code here in a code formatted block so I can help you figure this out? Thank you!

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/CHLNDDEV/OceanMesh2D/issues/201#issuecomment-803906406, or unsubscribehttps://github.com/notifications/unsubscribe-auth/APPMEIIDB4LEIAY26EX3USLTE4EFPANCNFSM4ZSHOAXA.

[krober10nd]

Thanks for that! There’s many optimization routines at the end and the way you build the mesh has a significant effect on the final quality (number of iterations, quality tolerance, mesh improvement rate, gradation bound)

In your case I would try lowering the gradation bound and/or enforcing a maximum element size.

Tô ensure the corner point is enforced, we can pass fixed points but this should happen automatically.

[krober10nd]

Also I would take care to resolve that large shelf break with the slp mesh heuristic. This will alter the wave transformation from deep to shallow water and especially in swan.

[krober10nd]

If you want to chat more about related topics we have a growing slack community

Join me on Slack -- it’s a faster, simpler way to work. Sign up here, from any device: https://join.slack.com/t/oceanmesh2d/shared_invite/zt-mxlhj1cl-i2TBm2SkBcCnDCTHJ9BgTw

[edeyejedi]

Thanks Keith. I have had a massive play around and managing to get in those corners most of the time. I have joined the slack thing as well. On the requirements for SWAN ( The number of triangles that meet at each vertex inside the mesh should not be smaller than 4 or larger than 10. Angles inside each triangle should not be higher than 143 degrees.) I have run a bunch of different tests, all different resolutions and settings but I cant seem to please SWAN. I started changing code to see if I can get the mesh to meet requirements, no luck. I could have it completely wrong though but here is some of the things I tried:

In meshgen, Line 942/943: % bad_ele = any(tq_n.vang < 1pi/180 | ... % tq_n.vang > 179pi/180,2); bad_ele = any(tq_n.vang < 1pi/180 | ... tq_n.vang > 143pi/180,2); Thing about this is that it makes no sense as the original code does not restrict angles to 175 degrees as per Roberts et al. 2019; so I could be off the mark here but I cant see where else angles are dealt with.

In get_small_connectivity in meshgen: % I = find(enum <= 4); I = find(enum <= 4 & enum > 9); Added the "greater than 9" to see if this would help but I don't think it is that end that is the problem.

I check the angles before writing the fort.14, sometimes this is good, and swan is happy with the angles, my changes to the meshgen code don't seem to make a difference though! The error in SWAN I cant shift is the smaller than 4 or larger than 10.
Any help would be appreciated.

Cheers

[krober10nd]

Hey Ed,

What's your full script that you're using to call the generator?

I can suggest how you can improve the final mesh if I can see this.

We have a module called msh.clean that is designed to clean up low quality elements and improve qualities.

For instance try

help msh.clean

you can try

m = clean(m,'ds',2,'conn',7); 

[edeyejedi]

Hi Keith,

Code below. Note that I have made this quite coarse (big value for min mesh res) so I can crank through testing.

%% STEP 1: set mesh extents and set parameters for mesh.
bbox1 = [166.1 169.5 -47.75 -44.1];
min_el=2000;  % Minimum mesh resolution in meters.
max_el=30e3;
max_el_ns=5000; 
R=5;  % Number of elements to resolve feature.
g=0.1; % Grade
%% STEP 2: specify geographical datasets and process the geographical
dem       = 'L2_gf1_posup.nc'; % local bathy
coastline = 'nz-coastlines-and-islands-polygons-topo-150k.shp'; % local coastline
gdat1 = geodata('shp',coastline,'dem',dem,'h0',min_el,'bbox',bbox1);
%% STEP 3: create an edge function class
fh1=edgefx('geodata',gdat1,'max_el',max_el,'fs',R,'max_el_ns',max_el_ns,'g',g);
mshopts1=meshgen('ef',fh1,'bou',gdat1,'plot_on',1,'proj','lam');
mshopts1 =mshopts1.build;
m1 = mshopts1.grd; 
m1=clean(m1,'ds',2,'con',7); % new recommendation from krober
m1 = interp(m1,gdat1,'mindepth',0); % interpolate bathy 
write(mshopts1.grd,'fort');

[krober10nd]

Thanks, so I would try increasing the itmax parameter to 200 or 300 and lower the termination quality tolerance firstly

mshopts1=meshgen('ef',fh1,'bou',gdat1,'plot_on',1,'proj','lam','itmax',200,'qual_tol',1e-3);

I'm in the process of depcretating interpolation from gdat as the DEM may be downsampled for runtime requirements. I'd rather you interpolate directly from the DEM using

m1 = interp(m1,dem,'mindepth',0); 

Also, for SWAN to work, I believe all boundaries need to be labeled but I may be mistaken. Check out help make_bc to do the boundary labeling.

[krober10nd]

The more iterations it does, the higher the final quality metrics should be (in general). You generally want to stop it on an evenly divisible number of 10 as the mesh improvement occurs every 30 iterations or so. You do not want to stop it after a mesh improvement cycle inside the meshgen imo.

[krober10nd]

Also, note that the result you get with a super coarse mesh will likely be worse than a much finer one for the same domain. The nature of the meshing algorithm tends to be more stable and work better for slowly graded and fine resolution meshes.

[krober10nd]

hey @edeyejedi can I close this issue? I think you said you ended up getting some favorable results after a lot of trial and error.

[edeyejedi]

Hi @krober10nd not sure. I am struggling again with the edges right now, not just he corners but along the edge of the bbox. most notably when I use combined edgefx and geodata obs (e.g. meshgen('ef',{fh1 fh2 fh3},'bou',{gdat1 gdat2 gdat3})), seems related to resolution (high res causing more problems). Have cranked up the itmax (500) and lowered the qual_tol (1e-6). Have not played with the gradation boundary (this is 'g' right?), and have it set the same in the edgefx for all 3 nests.

Perhaps we can close. There is probably more than enough here for others to read and figure out they need to do some tweaking of input parameters, and what those parameters might be.

[krober10nd]

I would try to comment this if clause out starting here

https://github.com/CHLNDDEV/OceanMesh2D/blob/f54a6d97e940734df94cfe76d02808c61959c783/%40meshgen/meshgen.m#L568

[krober10nd]

And re-run your case to see if it improves. I suspect it may.

[edeyejedi]

Copy. Will report back

[edeyejedi]

Distinct improvement by reducing 'g' from 0.1 to 0.05. Corner vertex still missing though (despite a pfix) Quick question, can you interp on a pair of DEM's? that is I have 2 DEMs, one high res and once low res,, would the following be effective: m0=interp(m0.grd,{dem1 dem2},'mindepth',0);
Can dive in a have a look, but just wondering if this has been accounted for (like in mshgen using multiple edgefx or geodata obs). Might also want to start new thread for that one!

[krober10nd]

Quick question, can you interp on a pair of DEM's? that is I have 2 DEMs, one high res and once low res,, would the following be effective: m0=interp(m0.grd,{dem1 dem2},'mindepth',0);

Yep.

[krober10nd]

I would just unpack the msh object from meshgen and do this:

m = m0.grd; 
m = interp(m,{dem1,dem2},'mindepth',0); 

see help msh.interp for more functionality.

and note #204 for DEMs with irregular grid spacing along an axis using cell-averaging interpolation.

[edeyejedi]

Yes re. the unpacking. Me not unpacking was a typo

[edeyejedi]

Still missing corners but not the end of the world for my application

[krober10nd]

Yea these corners have historically been an aesthetic issue but often I try to avoid corners in modeling domains and prefer arc'ed domains.. It's strange to me that pfix isn't keeping them put likely it's because one of the mesh improvement strategies may be touching it.

[krober10nd]

I'm going to close this one. Please start a new issue as issues arise!