Issues with refining mesh in channels

[soloyant]

Hello,

I am trying to build a mesh that needs to be refined at a h0=10m in resolution in a network of channels while the main domain's resolution should be bounded to h0=30m. The domain being large, I can't afford to lower h0 to 10m everywhere, the number of elements gets too high. I am a bit stuck and would appreciate any advice or suggestions you may have to try to achieve this goal.

So far, I have tried 3 approaches:

1. Defining sub regions with h0=10m, but there are 2 problems with this approach: a) Each subregion is likely to require (or not) to be flipped using the inpoly_flip switch. This seems to be happening at random, and requires manual checking which limits the number of subregions (ideally, this mesh could use hundreds of such polygons, my tries were limited to using the 3 largest ones). I am not sure to understand how this in/out poly works in the program, and wasn't able to fix the issue ahead in the process, I tried using clockwise and counterclockwise arrays of polygon coordinates but it didn't seem have any effect. b) The build algorithm seems to get confused and ignores parts of the domain in the process (see animation here bellow). The final result after cleaning is a somewhat destroyed mesh with lots of holes, including pfix points.

2. Using only one main domain and lowering the h0 argument down to 10m everywhere in the domain while using the slope as a refinement guide with a low slp value of 3 to 5 (which is far from ideal according to the manual), in order to limit the number of elements. In this case the build method iterates correctly but the mesh gets destroyed again during the cleaning process.

3. Using the polyline function fed with the centerlines of the network of channels. My understanding is that this method defines channel width as a function of depth, but this hypothesis is not valid for the region of interest, so it is probably not ideal. In any case, successful runs of the build method showed that the mesh size along centerlines seems to be bounded by the main domain h0 = 30 despite using the argument min_el_ch = 10m, therefore making it impossible to refine channels more than their surrounding environment. At this point, approach 3 seems to be the most promising one, and I was hoping to find a way to make the min_el_ch argument taking over h0 within channels, but I'm not sure how to yet, nor if this is doable.

Script of approach 3:

% Load domain data
dem = 'dem.nc';
domain = readtable('domain.csv');
domain.Properties.VariableNames = {'lat', 'long'};
channels = shaperead('channels_centerlines.shp');
pts = cellmat();
for i = 1:size(channels,1)
    pts{i} = [channels(i).X; channels(i).Y]';
end

% Mesh settings - main domain:
bbox = [domain.long, domain.lat]; % polygon boubox
min_el    = 30;                   % minimum resolution in meters.
max_el    = inf;                  % maximum resolution in meters. 
max_el_ns = 300;                  % maximum nearshore resolution in meters. 
grade     = 0.2;                  % mesh grade in decimal percent. 
R         = 3;                    % Number of elements to resolve feature.
dt        = 0.5;                  % Automatically set timestep based on nearshore res
wl        = 10;                   % Minimum number of triangles to resolve tidal wavelengths;
slp       = 10;
pfix = [domain.long, domain.lat];

% Mesh settings - channels:
angleOfReslope = 157 ;            % Control width of channel by changing angle of reslope.
ch = 0.4 ;                        % Scale resolution propotional to depth nearby thalweg.
min_el_ch = 10;

gdat{1} = geodata('dem', dem, ...
                  'h0', min_el, ...
                  'bbox', bbox);
fh{1} = edgefx('geodata', gdat{1}, ...
               'fs', R, ...
               'max_el', max_el, ...
               'g', grade, ...
               'dt', dt, ...
               'wl', wl, ...
               'slp', slp,...
               'ch',ch,...
               'AngOfRe',angleOfReslope,...
               'Channels',pts,...
               'min_el_ch',min_el_ch);

if gdat{1}.inpoly_flip == 1
    gdat{1}.inpoly_flip = 0;
else
    gdat{1}.inpoly_flip = 1;
end

mshopts = meshgen('ef', fh, ...
                  'bou', gdat, ...
                  'plot_on', 1, ...
                  'proj', 'Mercator', ...
                  'itmax', 9999, ...
                  'memory_gb', 60, ...
                  'pfix', pfix, ...
                  'nscreen', 1, ...
                  'dj_cutoff', 0);
mshopts = mshopts.build;

Thanks!

https://github.com/CHLNDDEV/OceanMesh2D/assets/36890211/b8bf8372-d49e-4aee-a3b4-4c2e764314ed

[krober10nd]

To save in memory I would divide the domain into four vertical equal sized slices and mesh each one independently. The slices should be slightly overlapping.

For refinement, I would make a shapefile filled with those green lines and l use this only in an edgefx definition with a distance based mesh function. Not the geodata definition. These things can be separated.

Then I would merge the slices using msh.plus.

[krober10nd]

The pull request up for review can address this provided of course the vector does not contain features smaller than the 2x or 3x the minimum element size. In practice, you build the mesh once and then write to a .2dm file (m.write('my_mesh.2dm','2dm') and in QGIS with the mesh and vector file overlaid, edit the vector file based on how it meshed around the structures.

https://github.com/CHLNDDEV/OceanMesh2D/pull/264

[GHGouthamHari]

i am new here in oceanmesh 2d user. i have created mesh for oceanpart so far, but to do same for some extend in land boundary, i am not able to do so, any suggesition????

[krober10nd]

You can make the domain any polygon this polygon could represent a floodplain extent instead of a shoreline.

[GHGouthamHari]

sir i am trying examples of oceanmesh 2d, but in Example_6b_GBAY_w_floodplain.m i am facing problem about dem file, if i use srtm data instead of that file mentioned there, will there be an issue??

thank you in advance.

[WPringle]

sir i am trying examples of oceanmesh 2d, but in Example_6b_GBAY_w_floodplain.m i am facing problem about dem file, if i use srtm data instead of that file mentioned there, will there be an issue??

thank you in advance.

Of course you can use any DEM file you choose just need to keep in mind the limitations of that particular DEM with regards to resolution and accuracy

[GHGouthamHari]

I am trying to create a mesh using the GEBCO-2023 gridded data which has both bathymetry and topography information. I have used code "Example_7_Global.m" for the mesh generation. However, even though the ocean part is coming properly, the bathymetry in the mesh is shown till 1m depth. There is no topography information in the output mesh. All land points and bathymetry between 0-1m depth are replaced by constant 1m depth.

Can we generate mesh using GEBCO gridded data, for both ocean and land part using this example code?

[krober10nd]

This issue is not the appropriate place to post this question. Please create a new issue.

[GHGouthamHari]

Thank you, I will do that.

On Wed, Jul 31, 2024 at 5:20 PM Keith Roberts @.***> wrote:

This issue is not the appropriate place to post this question. Please create a new issue.

— Reply to this email directly, view it on GitHub https://github.com/CHLNDDEV/OceanMesh2D/issues/291#issuecomment-2260338014, or unsubscribe https://github.com/notifications/unsubscribe-auth/BKBXGSPIDX6T5OPJTLYDUX3ZPDFRDAVCNFSM6AAAAAAZ7JUIICVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDENRQGMZTQMBRGQ . You are receiving this because you commented.Message ID: @.***>