CHLNDDEV / OceanMesh2D

A two-dimensional triangular mesh generator with pre- and post-processing utilities written in pure MATLAB (no toolboxes required) designed specifically to build models that solve shallow-water equations or wave equations in a coastal environment (ADCIRC, FVCOM, WaveWatch3, SWAN, SCHISM, Telemac, etc.).
https://github.com/sponsors/krober10nd
GNU General Public License v3.0
188 stars 68 forks source link

Distorted elements for constant resolution mesh #82

Closed sbrus89 closed 4 years ago

sbrus89 commented 4 years ago

Hi @WPringle and @krober10nd,

Thanks for continuing to develop OceanMesh2D! I just started using it recently and noticed some distorted elements in a simple mesh I was making. I thought I'd ask for your advice. Maybe you have experienced something similar.

I modified Example 7 to create a constant resolution 50km global mesh and noticed that there some areas where elements get distorted. See below for examples. Other than these regions, the mesh looks great. I'm using ee76beadd.

New Zealand: new_zeland_issue

Fiji: fiji_issue

Alaska: alaska_issue

South Africa: south_africa_issue

sbrus89 commented 4 years ago

Here is my driver script:


clearvars; clc;
addpath(genpath('utilities/'));
addpath(genpath('datasets/'));
addpath(genpath('m_map/'));

%% STEP 1: set mesh extents and set parameters for mesh. 
%% The greater US East Coast and Gulf of Mexico region
min_el    = 50e3;                    % minimum resolution in meters.
bbox      = [-180 180; -88  85]; % lon min lon max; lat min lat max
max_el    = 50e3;                        % maximum resolution in meters. 
wl        = 30;                  % 30 elements resolve M2 wavelength.
dt        = 0;                   % Only reduces res away from coast
grade     = 0.25;                % mesh grade in decimal percent. 
R         = 3;                           % Number of elements to resolve feature.
slp       = 10;                  % slope of 10

outname = 'Global_50km'

%% STEP 2: specify geographical datasets and process the geographical data
%% to be used later with other OceanMesh classes...
dem       = 'bathy_SSG_1_120.nc';
coastline = {'GSHHS_i_L1','GSHHS_i_L6'};
gdat1 = geodata('shp',coastline,'dem',dem,...
                'bbox',bbox,'h0',min_el);

%% STEP 3: create an edge function class
fh1 = edgefx('geodata',gdat1,...
             'max_el',max_el);

%% STEP 4: Pass your edgefx class object along with some meshing options 
%% and build the mesh...
mshopts = meshgen('ef',fh1,'bou',gdat1,'nscreen',10,'plot_on',1,...
                  'proj','stereo');
mshopts = mshopts.build;

%% STEP 5: Match points and edges across boundary
m = mshopts.grd; % get out the msh object 

% Plotting the triangulation on Robinson
plot(m,'tri',1,'Robinson');
save([outname '.mat'],'m');
write(m,outname);```
krober10nd commented 4 years ago

Hey Steven! Interesting. I've seen this before. In the past, this distortion likely originates from the mesh topology cleaning step at the end of the generation stage. Specifically, what we call implicit smoothing. This particular step can be deactivated by deferring "mesh cleaning" to outside of the generation class.

  1. Mesh cleaning in the generation step can be deactivated by passing the name value argument "cleanup" 0 to the meshgen class constructor.
mshopts = meshgen('ef',fh1,'bou',gdat1,'nscreen',10,'plot_on',1,...
                  'proj','stereo','cleanup',0);
  1. After calling meshgen.build() and unpacking the mesh, write
m = mshopts.grd; 
m = clean(m,'ds',0); 

Let us know how that works. Could be something else, but easy enough to try I suppose.

Hope all is well!

sbrus89 commented 4 years ago

Thanks for the suggestion, Keith! I'll give it a try and let you know how it goes.

Hope all is well with you too!

sbrus89 commented 4 years ago

@krober10nd, that seems to have worked! Thanks for your help.

Here's (approximately) the same locations after re-running using your suggestion:

New Zealand: new_zealand

Fiji: fiji_fixed

Alaska: alaska_fixed

South Africa: south_africa_fixed

krober10nd commented 4 years ago

Cool glad it helped! I’ll take a closer look at the cause of the distorted elements this weekend.

On Wed, 13 May 2020 at 8:13 PM Steven Brus notifications@github.com wrote:

@krober10nd https://github.com/krober10nd, that seems to have worked! Thanks for your help.

Here's (approximately) the same locations after re-running using your suggestion:

New Zealand: [image: new_zealand] https://user-images.githubusercontent.com/3514023/81875077-b24a1f80-953c-11ea-94fe-c9320ced26e1.png

Fiji: [image: fiji_fixed] https://user-images.githubusercontent.com/3514023/81875092-bb3af100-953c-11ea-9c30-862e36e543a3.png

Alaska: [image: alaska_fixed] https://user-images.githubusercontent.com/3514023/81875112-c5f58600-953c-11ea-9150-3a9403860f17.png

South Africa: [image: south_africa_fixed] https://user-images.githubusercontent.com/3514023/81875131-d148b180-953c-11ea-8add-b78c9b0b14d1.png

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/CHLNDDEV/OceanMesh2D/issues/82#issuecomment-628293119, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEOBZ7B4KOVCKBQEVQUY6DDRRMSP5ANCNFSM4NAFAPZQ .

--

Keith Roberts Post-doctoral researcher University of São Paulo krober@usp.br krober10@nd.edu

krober10nd commented 4 years ago

looks like these distorted elements happen in the edge case when using uniform resolution, which we have not frequently tested in past development.

We'll make the mesh improvement hill-climbing in the sense that it will only accept the improvements to the connectivity if they do indeed improve the mesh OR the minimum quality is sufficiently low to warrant entering into the step in the first place. @WPringle

krober10nd commented 4 years ago

In #87, an option to use a different hill-climbing optimizer was added that should avoid potentially inverted elements and altering the underlying mesh resolution distribution by passing the name-value pair "ds",2 to the msh.clean routine.