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