chlnddev/oceanmesh2d

A two-dimensional triangular mesh generator with pre- and post-processing utilities written in pure MATLAB (no toolboxes required, some support for Octave) designed specifically to build models that solve shallow-water equations or wave equations in a coastal environment (ADCIRC, FVCOM, WaveWatch3, SWAN, SCHISM, Telemac, etc.).

Precise distance-based two-dimensional automated mesh generation toolbox intended for coastal ocean/shallow water flow models

Note

This is the default and recommended PROJECTION branch. Please use it unless you otherwise require legacy (MASTER branch) or the absolute newest features (DEV branch).

Warning

The boundary of the meshing domain must be a polygon (first point equals the last and non-self intersecting) but it does not need to be simplified. Read the user guide for more information about the inputs.

Table of contents

• Getting help
• Contributing
• Code framework
• Starting out
• Testing
• References
• Gallery
• Changelog

OceanMesh2D is a set of user-friendly MATLAB functions to generate two-dimensional (2D) unstructured meshes for coastal ocean circulation problems. These meshes are based on a variety of feature driven geometric and bathymetric mesh size functions, which are generated according to user-defined parameters. Mesh generation is achieved through a force-balance algorithm combined with a number of topological improvement strategies aimed at improving the worst case triangle quality. The software embeds the mesh generation process into an object-orientated framework that contains pre- and post-processing workflows, which makes mesh generation flexible, reproducible, and script-able.

Getting help

Tip

PLEASE READ THE USER GUIDE! A recent pdf of the user guide is located in this branch. For a continually updated version click here (wait for compilation and then click download PDF)

Besides posting issues with the code on Github, you can also ask questions via our Slack channel here.

Note: If the slack link invite isn't working, please send either one of us and an email and we'll fix it. By default, the invite link expires every 30 days.

Otherwise please reach out to either Dr. William Pringle ([email protected]) or Dr. Keith Roberts ([email protected]) with questions or concerns or feel free to start an Issue in the issues tab above.

Contributing

All contributions are welcome!

To contribute to the software:

1. Fork the repository.
2. Clone the forked repository, add your contributions and push the changes to your fork.
3. Create a Pull request

Before creating the pull request, make sure that the examples pass...

Some things that will increase the chance that your pull request is accepted:

• Write minimal working examples that demonstrate the functionality.
• Write good commit and pull request messages.

Code framework

OceanMesh2D consists of four standalone classes that are called in sequence. It requires no paid toolboxes to build meshes and has been tested to work with a trial version of MATLAB.

OceanMesh2D::
├── geodata -- process geospatial data.
├── edgefx  -- build mesh size functions.
├── meshgen -- generate mesh based on mesh size functions and boundaries.
└── msh     -- store, write, read, inspect, and visualize meshes and their
               axuillary components for numerical simulation.

Starting Out

1. Clone or download and unzip the current repository.

2. Run the setup.sh bash script or setup.bat batch script to download the required m_map package and the base datasets:

   • GSHHG global shoreline
   • SRTM15_PLUS global topobathy DEM

3. Create or modify startup.m file located under your home folder by executing command in MATLAB:

   edit(fullfile(userpath,'startup.m'))
   

   And add the line:

   run(<PATH_TO_OCEANMESH2D>/setup_oceanmesh2d.m)
   

   Replace <PATH_TO_OCEANMESH2D> with the folder name where you cloned/unzipped the repository.

4. Restart MATLAB to get paths added.

Additional data required for some of the following examples must be downloaded manually from here. Specifically, Examples 2, 3, 4, 5 and 5b require additional datasets from the google drive folder while base datasets are sufficient for the other examples.

Featured in  ┌── Examples/Example_1_NZ.m   -- A simple mesh around South Island
             |                                New Zealand that uses GSHHS shoreline.
user guide   ├── Examples/Example_2_NY.m   -- A high-resolution mesh around
             |                                the New York/Manhattan area that uses a DEM
             |                                created from LiDAR data.
             └── Examples/Example_3_ECGC.m -- Builds a mesh for the western North Atlantic
                                              with a local high-resolution nest around New York

Featured in         ┌── Examples/Example_4_PRVI.m -- Builds a mesh for the western North Atlantic
                    |                                with three high-resolution nests around
Geoscientific Model |                                Puerto Rico and US Virgin Islands
                    ├── Examples/Example_5_JBAY.m -- An extremely high-fidelity (15-m) mesh from
Development paper[1]|                                LiDAR data around Jamaica Bay with CFL-limiting.
                    └── Examples/Example_6_GBAY.m -- An example of the polyline/thalweg mesh size
                                                     function along the Houston Ship Channel.

See Testing to test OceanMesh2D on your system.

Testing

To ensure the software is fully functional on your system before building some crazy meshes it is strongly recommended to run the tests (RunTests.m) in the Tests/ directory.

We test all pull requests using this test suite on a local host before accepting. For substantial pull requests we will also test the Examples from the Examples/ directory.

References!

If you make use of OceanMesh2D please include a reference to [1], and to any of [2]-[5] if pertinent (latex .bib file). We would also appreciate using our logo in a presentation featuring OceanMesh2D.

[1] - Roberts, K. J., Pringle, W. J., and Westerink, J. J., 2019.
      OceanMesh2D 1.0: MATLAB-based software for two-dimensional unstructured mesh generation in coastal ocean modeling,
      Geosci. Model Dev. (GMD), https://doi.org/10.5194/gmd-12-1847-2019.
[2] - Roberts, K. J., Pringle, W. J, 2018.
      OceanMesh2D: User guide - Precise distance-based two-dimensional automated mesh generation toolbox intended for coastal
      ocean/shallow water. https://doi.org/10.13140/RG.2.2.21840.61446/2.
[3] - Roberts, Keith J. Unstructured Mesh Generation and Dynamic Load Balancing for Coastal Ocean Hydrodynamic Simulation, 2019.
      PhD Thesis, University of Notre Dame. https://curate.nd.edu/show/4q77fr0022c.
[4] - Roberts, Keith J., Pringle W.J., Westerink J. J. Contreras, M.T., Wirasaet, D., 2019.
      On the automatic and a priori design of unstructured mesh resolution for coastal ocean circulation models,
      Ocean Modelling, 144, 101509. https://doi.org/10.1016/j.ocemod.2019.101509.
[5] - Pringle, W. J., Wirasaet, D., Roberts, K. J., and Westerink, J. J., 2021.
      Global Storm Tide Modeling with ADCIRC v55: Unstructured Mesh Design and Performance,
      Geoscientific Model Development, 14(2), 1125-1145. https://doi.org/10.5194/gmd-14-1125-2021.

In addition, best practice when using software in a scientific publication is to cite the permanent doi corresponding to the version used (e.g., for reproducibility). All our releases are archived at the following Zenodo repository doi link.

Authors (202X). CHLNDDEV/OceanMesh2D: OceanMesh2D VX.X. Zenodo. https://doi.org/10.5281/zenodo.1341384

Please fill in the version (VX.X), author list and year corresponding to the version used.

We would also like to acknowledge various scripts and algorithms from mesh2d included in OceanMesh2D that have been developed by @dengwirda. Please also see JIGSAW-GEO:

[i] - Engwirda, D., 2017.
      JIGSAW-GEO (1.0): Locally orthogonal staggered unstructured grid generation for general circulation modelling on the sphere.
      Geoscientific Model Development, 10(6), 2117–2140. https://doi.org/10.5194/gmd-10-2117-2017.

GALLERY:

• These images can be made by running the examples

                                               

• The following images are provided from happy users. Please send us your meshes.

Jiangchao Qiu from his publication.

       

Changelog

The format is based on Keep a Changelog

Unreleased (on current HEAD of the Projection branch)

[6.0.0] - 2024-02-28

Added

• Added ability to generate constrainsts for shorelines automatically using the new high-fidelity method. #264
• auxiliary file reader function that can be used standalone e.g., m = m.read({'fort.13','fort.15'})), as well as called from within msh() function. #287
• swanoutput namelist for SWAN model outputs. #287
• Added new function in msh called remesh_patch to remesh within specified polygonal domains and insert back into parent mesh. #301
• Read and write to 2dm format. #300
• namelist and RSTIMNC input arguments for Make_f15.m fort.15 generator. updated the help message for all input argumebts to Make_f15. #283
• New stability namelist options to Make_f15.m fort.15 generator. #283
• Optionality for mesh2dgen to choose the method of mesh generation kind, and the maximum iteration count iter. #272
• Option improve_with_reduced_quality to meshgen for allowing for mesh improvements even when quality is decreasing or large number of nodes eliminated, which sometimes is necessary to force the advancement in mesh quality.
• Option delaunay_elim_on_exit to meshgen to skip the last call to delaunay_elim to potentially avoid deleting boundary elements.
• Geoid offset nodal attribute in Calc_f13 subroutine. #251
• Support for writing Self Attraction and Loading (SAL) files in NetCDF for the ADCIRC model. #231
• Added 'high-fidelity' option for automatically forming and constraining edges into the mesh. #264

Changed

• treating logicals in namelists as a logicals type within MATLAB but ouput as a string. #287
• Calc_Sponge.m to allow for sponge in polygon region and coefficents based on spatially varying depth. #287
• removed namelists from default setup in Make_f15.m fort.15 generator to be invoked by user as an input argument instead. #283
• Use implicit smoother (ds=1) in msh.clean when fix points are present. #283
• Default filename for the dynamicWaterLevelCorrection is now null so that it is not evoked by default. #272
• Default mesh improvement strategy is ds 2.
• Retrieve boundary indices in msh.get_boundary_of_mesh method. #259
• msh.offset63 struct and associated write/make routines for dynamicwaterlevel offset functionality. #259
• dynamicWaterLevelCorrection to fort.15 namelist, and PRBCKGRND option to met fort.15 namelist. #261
• Improved code readability by auto-indenting core functions and reworking 'setup' and 'Examples' functions, including adding necessary files for separation from core functions. #308
• Renamed functions for clarity and fixed multiple small bugs; changes are compatible with OpenEarthTools but distinct. #308

Fixed

• catch for offset_nodes of int type in writeoffset63.m. #287
• writing out the wave coupling timestep RSTIMNC on the WTIMNC line of fort.15 when NWS > 300. #283
• make_bc call with empty varargin, e.g., when calling inner. #283
• Make_offset63 call with constant offset (length 1 time_vector). #283
• ourKNNsearch call in nanfill option of Griddata (msh.interp). #283
• m_map link in setup.sh Issue #277. #283
• Updated Calc_f13.m to avoid an "Unrecognized variable" error by ensuring "broken" is always defined. #282
• Fixed test for likely geographic coordinates in Make_f15.m. #282
• updated Gridded_to_Mesh_SeaBed_DepthAveraged.m to fix the infinite loop in using Cal_IT_Fric.m by filling in the NaNs at greater depths with values from above. #280
• Recursive cleaning issues: infinite loop and preservation of fixed points.
• msh.interp method for K argument of length 1, and for the test to determine whether the bathymetry grid is irregular. #259
• Printing of namelist character strings or numbers. #261
• Make_offset63.m time interval computation. #261 and #272
• Removed dependency on statistics toolbox when using the 'nanfill' option in msh.interp. #269
• Missing routines for reading in elvstaname and velstaname in readfort15.m by adding readlinevecname() method. #281
• Incorrect reference to ibtype was changed to ibtypee in map_mesh_properties function. #298

[5.0.0] - 2021-07-29

Added

• meshgen.build() now will rewind the iteration set in the case mesh improvement cannot improve the qualities. #234
• msh.plot() now has a cmap in which the user can specify their any cmocean colormap
• radius_separated_points function that trims the points in the mesh to have a specified resolution that can be used before m_quiver so that vectors are evenly plotted. #225
• Deleting boundary conditions by specifyng their indices in msh.object.bd field. See #205
• Ability for user to set their own axis limits when plotting with msh.plot(). #224

Fixed

• Minor fix to msh.make_bc using the auto method. #237
• correction in setting stereographic projection bounds in setProj to make sure points are not pushed outside and become NaNs (was limited to radius of 178 deg but made sure can go up to full 180 deg). #225
• Correctly deleting weirs from boundary object through make_bc delete method. See #205
• Array format fix for reading in ibtype and nvell from fort.14 file and when executing carry_over_weirs. See #206
• Fix for irregular grid spacings in DEMs. See #204
• tidal constituents for Make_f15 can now contain "major8" in addition to other constituents in the string/cell array #221
• Correctly collect NDBC and NOS stations in mesh when creating fort15 file using Make_f15 for meteorological, velocity and elevation records #242

Changed

• msh.plot() using type bd option now creates a legend for the different boundary condition types. #247
• forcing facecolor to white in m_trimesh so that it does not intefere with background color option. #245
• made topographic elevation bound option for max_ele, wl, slp, and g edgefx kwargs consistent and added explanation of this option is included in edgefx help. #230
• m_plot() function calls m_grid() with background color input kwarg (if backcolor option used) instead of manual application. #225
• tidal_data_to_ob function called from Make_f15 populates boundary condition tidal constituents that do not exist in the tidal database with zero values so that user can add user-defined values later (previously did not populate). #225
• Improved cmocean for pivot handling with discrete colormap. #225
• Renamed Calc_NLCD_Mannings to Calc_Mannings_Landcover and making option for 'ccap' landcover type in addition to 'nlcd' (default) and added the ability to using user specified inteprolation (e.g., nearest, linear, cell-averaging, etc.) of the landcover data to the mesh vertices. #221

[4.0.0] - 2021-03-14

Added

• mesh2d interface improvements to filter small polygons.
• Support for creation of fort.20 files for forcing rivers by @Jiangchao3
• Mesh "cleaning" modes moderate and aggressive transfer nodal attributes via improvements to msh.map_mesh_properties
• msh.remove_attribute() method to remove f13 attribute(s)
• new msh.interp() slope_calc kwarg option to set the method of computing the topographic gradients (slopes), either rms [default] or abs
• new extract_subdomain() keep_numbering kwarg option to keep the full mesh triangulation numbering on the subdomain [off by default].

Changed

• msh.plot() overhaul. All options specified via kwarg.
• msh.plot() option subset option is now called subdomain
• msh.plot() arbitary f13 option now utilizes colormap kwarg
• utilities/extract_subdomain now is called with kwargs.
• Cleaning up msh.bound_courant_number() to use msh.get_boundary_of_mesh() for Delaunay-triangulation and allowing msh.clean() to do the transfer of attributes automatically.
• msh.plus(obj1,obj2) can now carry over obj2 f13 attributes if also exist in obj1
• msh() more efficient storing of boundary conditions read in from fort.xx files, and msh.write() can write out arbitrary vertex indices (instead of just 1 to NP).

Fixed

• Boundary labeling fix
• Prompt when labeling bcs using outer kwarg in make_bc
• fix for boundary condition mapping in msh.map_mesh_properties() especially for weirs/barriers
• fix for barrier mapping in msh.plus() routine
• fix for msh.make_bc(m,auto,gdat) where gdat is empty. In this case it uses the depths on the mesh to determine the open boundaries.
• check for poly2ccw mapping toolbox function in kml2struct
• fix for msh.plot() on log colormap when plotting f13 attributes

Deleted

• Deprecating msh.CheckTimestep() for msh.bound_courant_number. Added error message and instruction in the CheckTimestep help

[3.3.0] - 2020-12-21

Fixed

• Users without mapping toolbox could not read in shapefiles because of a bug that made them required to have a 3d shapefiles.
• plotting gdat with no shoreline.
• plotting a mesh's bathymetry with a non-zero datum using cmocean.
• cell-averaging interpolation method in msh.interp fixed for unequal lon-lat DEM grid spacings

Added

• Mesh patch smoother
• Ability to remesh abritary patches of elements within the domain while respecting user-defined mesh sizes and the patches boundaries.
• Ability to use the TPXO9 Atlas for the tidal bcs and sponge (inside tidal_data_to_ob.m and Calc_Sponge.m) by using '**' wildcards in place of the constituent name within the tidal atlas filename (the atlas has an individual file for each constituent).
• Introducing 'auto_outer' option for the make_bc msh method which populates the bc for the outermost mesh boundary polygon (ignores islands)
• Changelog to README
• "mapMeshProperties" msh method ports over mesh properties for a mesh subset
• 'invert' option in the msh.interp method to turn off the DEM value inversion typically performed

Changed

• for the make_bc msh method 'auto'/'auto_outer' options, allowing for the 'depth' method of classification to use the interpolated depths on the mesh if gdat is empty.
• improving help for make_bc msh method, Make_f15.m and Calc_Sponge.m
• renamed "ExtractSubDomain.m" to "extract_subdomain.m"
• improving "extract_subdomain.m" help and facilitating NaN-delimited polygons
• ability to return boundary as a cell in "get_boundary_of_mesh" msh method
• "Example_1_NZ.m" includes example of plotting bcs of a msh subset
• using "mapMeshproperties" method in "fixmeshandcarry"
• using "fixmeshandcarry" in the "cat" msh method
• improving warning and error messages for the "interp" msh method
• adding geofactor into "writefort15" for the GAHM vortex model

oceanmesh2d's Issues

speed up?

The test to determine whether a given point is in/out of the meshing domain does not require the densified PSLG. Because the PSLG is by definition linear, the densification of points to h0/2 along the PSLG can only create new straight lines between two of the original points on the pslg.

With this said, the densified PSLG is still necessary for the nearest distance calculation. However, a large savings may be obtained through the suggested optimization by reducing the inpoly execution time. This could especially benefit the multiscale method.

Errors in user guide

In \subsection{Around the world: \textit{Example\_1\_NZ.m}} code example, you call the variable min_res wich does not exist.

In \subsubsection{Methods}, the third element in list has a broken reference Section \ref(sec:).

Add comments to 'auto' method of makens

It needs some comments on the logic applied as it's nearly impossible to tell how things are working after taking a break and coming back to it.

User-specified edge function

Hi @krober10nd and @WPringle,

I was wondering if you guys would be interested in adding the ability to input user-defined edge functions. I wrote a bit of crude code to do this in OceanMesh2D and am happy to share it back to the repo. If you're interested, I can open a PR. If not, no big deal at all!

Improve boundary conformity

• On termination, it's possible to repeatedly project only the boundary vertices back onto the 0-level set for better boundary conformity.

• See here what I did in Python.

https://github.com/krober10nd/SeismicMesh/blob/fa4cbc97a987d3468e6362fce6c76efe05f47bdf/SeismicMesh/generation/mesh_generator.py#L436

• It removes some of the "dents" that happen when meshing with an implicit domain.

Here's the function call:
https://github.com/krober10nd/SeismicMesh/blob/fa4cbc97a987d3468e6362fce6c76efe05f47bdf/SeismicMesh/generation/mesh_generator.py#L629

Mesh for SWAN with refinement nearshore

Hi,

I know this software is mainly designed for circulation model, but I want use it for SWAN that mainly refines the mesh nearshore. Hence, I only turn on distance function and closed options of wavelength, slope and polyline.

After that, the software mainly refines mesh with 20 km offshore (from red line to blue line), how can I extend such refinement to 50 km offshore (from red line to green line)? I know it can be refined in a small rectangular region, but not like offshore distance-based region in the figure. I also try to revise d-value in distance function, but it doesn't work. Could you help me to fix it? Thanks for your help.

Kind regards,
Jin

plot(m, 'b') fails due to missing projtype

(latest checkout)

>> plot(m, 'b')
Cell contents reference from a non-cell array object.

Error in setProj (line 52)
            projtype = projtype{1};

Error in msh/plot (line 274)
            del = setProj(obj,proj,projtype) ;

@msh/msh.m line 245 should catch this, but doesn't.
if proj && isempty(projtype)

In my test m.coord=[] and proj=1.
m was generated with m = msh('mesh14file');

thanks,
Hamish

Change to fixed points and edges in msh class

Currently the msh class contains fixed points (m.pfix) and fixed edges (m.egfix) variables.
m.pfix variables are just coordinates that presumably also exist in the main mesh point space (m.p). And m.edgfix refer to indices of m.pfix.

This is a somewhat clumsy way of treating fixed points and edges which makes it a bit annoying to manage them, and may require us to do a knnsearch to find the indices that match between m.p and m.pfix. For the meshgen it makes sense that we pass fixed points coordinates, but as an output into the msh class, the fixed points should just refer to indices in m.p.

I propose that we eliminate m.pfix and only keep m.egfix which should just be indices of m.p. That way we can also use the m.egfix in the constraints option for the DelaunayTriangulation class directly as is. m.pfix can be trivially recovered by m.p(unique(m.egfix(:)),:). This change will require us to treat the fixed points differently in the cleaning and merging routines etc, but it will be much neater imo.

multi-pass mesh generation using a fort.64 seed

Hi, an idea for your consideration--

I'm wondering if it could be useful to do multi-pass mesh generation by:

• using the existing M2 wavelength mesh size function with your input DEM in the usual manner,
• run ADCIRC with the resulting mesh and output a fort.64 file (Depth-averaged Velocity Time Series at All Nodes in the Model Grid),
• then read that fort.64.nc file back into OM2D,
• and finally run a new mesh size function based on observed velocities.

The idea being to better focus the mesh effort into areas of known high flow that the M2 shallow water wave eq. might not specifically capture, e.g. tidal straits or the inlet of an estuary which includes a large river flow as a boundary condition.

i.e. use a 'real' value for |u| in the CFL calculation instead of just estimating it.

At minimum it could be used to simply plot a map view highlighting areas where your mesh is either too coarse or needlessly fine, based on initial fort.64 results.

Distorted elements for constant resolution mesh

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:

Fiji:

Alaska:

South Africa:

Triangle mesh for WAVEWATCH III

Hi my friend,

I was wondering if OceanMsh2D can generate mesh file that can be read by WAVEWATCH III directly? Or I need to revise the mesh format (fort.14) to fit WAVEWATCH III mesh format (.msh)? Thanks.

Kind regards,
Jin

Plotting small meshes `plot`

• Often times small meshes are not immediately visible when plotting (i.e., plot).
• The user has to scroll around randomly to find it on the screen.
• Is this related to the axis sizing in simpplot?

Overloaded plot hold optional

Would be nice to have plot have one more option to ‘hold on’. By default all plot options should create a new figure but sometimes you want to plot on top of different plots.

bug

hello,my friend:
I downloaded the compression package and unzipped it under my MATLAB/Toolbox folder. I opened MATLAB and added the path. The software showed an error: the folder starting with @ could not be added to the path.
So I changed all the folder names with @ (such as @ann to Ann), but could not add them, showing that MATLAB could not save the changes to the modified path,pathdef.m may be a read-only file.I added each folder separately and found that it was the Ann folder that made the error.
I tried many methods but could not solve it. I installed other MATLAB tools normally, and the MATLAB I used was 2017A version. Hope to get your answer or email, thank you very much.

fort.15 cosmetics

Hi,

two small cosmetic issues with the generation of the fort.15 files, the first is mostly cosmetic, the second has functionality implications.

1. @msh/private/writefort15.m

diff --git a/@msh/private/writefort15.m b/@msh/private/writefort15.m
index 0d536df..6e108a0 100644
--- a/@msh/private/writefort15.m
+++ b/@msh/private/writefort15.m
@@ -148,7 +148,7 @@ end
 % Open boundary harmonic forcing  
 for k = 1: f15dat.nbfr
     fprintf(fid, '%s \n', f15dat.opealpha(k).name  ) ; 
-    fprintf(fid, '%16.9e %16.9e \n', f15dat.opealpha(k).val' ) ; 
+    fprintf(fid, '%16.9e %16.10g \n', f15dat.opealpha(k).val' ) ; 
 end
 
 %  ANGINN

The phase angle is given in engineering style, which makes it annoying to parse as regex in Perl. As engineering units are nice for values which vary over many orders of magnitude, but don't really add much for values which are evenly distributed over the range of 0-360, it's nice for humans reading it to see the regular floating point representation and easier to parse as well.
%16.10g preserves all the precision of %16.9e.
I don't know if anyone else beyond me has to deal with the Perl regex matching, but the human readability aspect might make the change of general interest.

2. utilities/Make_f15.m

diff --git a/utilities/Make_f15.m b/utilities/Make_f15.m
index 1893f37..aad52be 100644
--- a/utilities/Make_f15.m
+++ b/utilities/Make_f15.m
@@ -192,7 +192,7 @@ obj.f15.dtdp = dt ;
 % Runday
 obj.f15.rndy = datenum(te)-datenum(ts);
 
-obj.f15.extraline(10).msg = ts;
+obj.f15.extraline(10).msg = datestr(ts, 'yyyy-mm-dd HH:MM:SS UTC');
 
 % Checking f11
 if ~isempty(obj.f11)

The ts and te time variables are suggested to be in the format of '01-Jan-2018 03:00' in the Make_f15.m header-help and LaTeX docs* and that gets written to the NetCDF footer in fort.15. The parameter notes at adcirc.org state:

NCDATE = The format of NCDATE must be as follows so that ADCIRC can create netcdf files that comply with the CF standard: yyyy-MM-dd hh:mm:ss tz. For example, if the cold start date/time of the run is midnight UTC on 1 May 2010, the NCDATE parameter should be set to 2010-05-01 00:00:00 UTC.

and besides that it means you can't use GNU date to do timestep math on the result, e.g.:
$ date -u -d '2018-01-01 00:00:00 UTC + 3600 sec'
works but
$ date -u -d '01-Jan-2018 03:00 + 3600 sec'
doesn't.
(the above trick with date is a very handy for making time stamps for animations of fort.63 frames, given the start time from the .nc file and 1680000 or whatever seconds of the output timestep)

Matlab's datestr() does not like time zone info like "UTC" or "+0000", so that has to be added by hand. (and assumed to be true!)

[*] n.b. the underscores of Make_f15 links in the LaTeX doc have been UTF-8'd and so become hard to ctrl-f searched for

thanks,
Hamish

error in ParseShoreline

I am testing out the mesh generation software and I am getting the following error when running Example_1_NZ using m_map 1.4 and Matlab 2017a on a windows laptop. Any help is much appreciated.

Reading shapefile with m_shaperead
179837 Records in file, by 500:

.......................................................................................................................................................................................................................................................................................................................................................................Partitioning the boundary into islands, mainland, ocean
Index exceeds matrix dimensions.
Error in Read_shapefile (line 149)
tmpCC{nn,2} = tmpC{ii,2};
Error in geodata/ParseShoreline (line 276)
polygon_struct = Read_shapefile( obj.contourfile, [], ...
Error in geodata (line 240)
obj = ParseShoreline(obj) ;
Error in Example_1_NZ (line 17)
gdat = geodata('shp',coastline,'bbox',bbox,'h0',min_el);

Handle DEM as scatter points..

G'day @krober10nd and @WPringle,

I have bathymetry (*.mat) data in scatter points (xyz, three columns) , now I want to use it in geodata instead of using "dem" as *.nc file.

Is there any way to handle it?

Thank you very much!
Huy

Issue running examples (edgefx/featfx error).

Hi,

I have been following the guide and I am trying to run through the examples.
I have installed m_map in the root folder and put the 'GSHHS_f_L1' shapefile in the datasets/ folder.
However I am running into an issue, the following error occurs when running Example_1_NZ.m:

Reading shapefile with m_shaperead

179837 Records in file, by 500:

.......................................................................................................................................................................................................................................................................................................................................................................Partitioning the boundary into islands, mainland, ocean
Smoothing coastline with 5 point window
Reading shapefile with m_shaperead
5706 Records in file, by 500:

...........Partitioning the boundary into islands, mainland, ocean
'Read in meshing boundary: ' 'GSHHS_f_L1'

Building feature size function...
Warning: No mainland or inner

In edgefx/featfx (line 262)
In edgefx (line 212)
In Example_1_NZ (line 20)
Error using &
Matrix dimensions must agree.

Error in edgefx/finalize (line 775)
hh_m(nearshore & hh_m > obj.max_el_ns) = obj.max_el_ns;

Error in edgefx (line 246)
obj = finalize(obj,feat);

Error in Example_1_NZ (line 20)
fh =

edgefx('geodata',gdat,...

Height in shapefiles

Some shapefiles for example those from grass can be made 3D. In this case, they contain a 3rd dimension height. However, triggers the usage of toolbox functionality in geodata ('polybool') which is incredibly slow. Can we activate this through some other means?

Vert2Vert is wrong

It should be something like this:

% Generate triangulation
dtr = delaunayTriangulation(x,y);

% 1. Get all the triangles attached to a particular vertex in the
% triangulation.  
attachedTriangles = vertexAttachments(dtr);


for i = 1: size(x,1)
    % 2. Use the connectivity list to get the vertex indices of all these
    % triangles
    verticesOfTI         = dtr.ConnectivityList(attachedTriangles{i},:);
    % 3. Find all the unique vertices and remove the current vertex
    neighbours{i} = setdiff(unique(verticesOfTI), i);
end

Output msh as vtk/vtu files

A new file format for write should be added called vtu and vtk. This enables users to visualize their meshes inside Paraview and for solvers like FireDrake. The former visualization software is quite efficient compared to MATLAB's visualization capabilities.

For example this MathWork's submission is particularly good. https://www.mathworks.com/matlabcentral/fileexchange/47814-vtkwrite-exports-various-2d-3d-data-to-paraview-in-vtk-file-format

Search for 'longitude' and 'latitude' variables when reading .nc DEM files

(Sorry for being so annoying - I only want to improve this project)

When the program reads the netCDF DEM file, only checks for the existence of 'lon', 'lat', 'x' or 'y' variables. In some datasets (including the SRTM 15 Plus), those variables are called 'longitude' and 'latitude'.

geodata.m at line 161 must try to find those nomenclatures before throwing an error.

Mainland/sea boundary detection

Is there any way to force the program to mesh a selected side of the shapefile?

In some cases, I got meshed the land surface instead of the water one.

Smooth(er) outer buffer?

OceanMesh2D/@meshgen/private/smooth_outer.m

Line 10 in f0d0f6e

inx = x >= min(efs{bn}.F.GridVectors{1}) - 1*dx & ...

Increasing the multiplier on dx and dy makes resolution transitions smoother between nests appeasing some our our more critical users.

string quoting bugs in msh.m

Hi, on the commit 3 weeks ago " is used for quoting error message strings not ', and matlab is unhappy with that.

Error using msh
Error: File: /.../@msh/msh.m Line: 1176 Column: 23
Creating a string using double quotes is not supported.  Use the string function.

the following patch fixes it:

diff --git a/@msh/msh.m b/@msh/msh.m
index 71b3c97..673490f 100644
--- a/@msh/msh.m
+++ b/@msh/msh.m
@@ -1173,7 +1173,7 @@ classdef msh
             % varargins: 
             % none
             if nargin < 2
-                error("Needs type: one of 'auto', 'outer', 'inner', 'delete', or 'weirs'")
+                error('Needs type: one of "auto", "outer", "inner", "delete", or "weirs"')
             end
             if iscell(varargin{1})
                 varargin = varargin{1};
@@ -1183,11 +1183,11 @@ classdef msh
                 case('auto')
                     % automatically applies elevation and no-flux bcs to mesh based on geodata class 
                     if isempty(varargin)
-                        error("must supply a geodata class for 'auto' as third input")
+                        error('must supply a geodata class for "auto" as third input')
                     end
                     gdat = varargin{1};
                     if ~isa(gdat,'geodata')
-                        error("third input must be a geodata class for 'auto'")
+                        error('third input must be a geodata class for "auto"')
                     end
                     depth_lim = -10;
                     dist_lim = 10*gdat.gridspace;
@@ -1202,14 +1202,14 @@ classdef msh
                             end
                         case('depth')
                             if isempty(gdat.Fb)
-                                error("No DEM info to use 'depth' classification criteria")
+                                error('No DEM info to use "depth" classification criteria')
                             end
                             if length(varargin) > 2 && ~isempty(varargin{3})
                                 depth_lim = -varargin{3};
                             end
                         case('both')
                             if isempty(gdat.Fb)
-                                error("No DEM info to use 'both' classification criteria")
+                                error('No DEM info to use "both" classification criteria')
                             end
                             if length(varargin) > 2 && ~isempty(varargin{3})
                                 dist_lim = varargin{3};
@@ -1224,7 +1224,7 @@ classdef msh
                     Le = find(obj.p(:,1) < -179, 1);
                     Ri = find(obj.p(:,1) > 179, 1);
                     if ~isempty(Le) && ~isempty(Ri)
-                        error("Detected global mesh, cannot apply 'auto' makens. Try 'inner' option")
+                        error('Detected global mesh, cannot apply "auto" makens. Try "inner" option')
                     end
                     
                     % Get the boundary edges and mid-points
@@ -1391,7 +1391,7 @@ classdef msh
                     Le = find(obj.p(:,1) < -179, 1);
                     Ri = find(obj.p(:,1) > 179, 1);
                     if ~isempty(Le) && ~isempty(Ri)
-                        error("Detected global mesh, cannot apply 'auto' makens. Try 'inner' option")
+                        error('Detected global mesh, cannot apply "auto" makens. Try "inner" option')
                     end
                     
                     % Get the boundaries

gdat plot

• plotting the DEM with no shoreline results in a blank figure.

dem = geodata('dem',DEMNAME,'h0',min_el); 
plot(gdat,'dem');

error in Calc_tau0

Hi @krober10nd and @WPringle,

I am following the typical workflow provided in the userguide to create input files (fort.15) for ADCIRC simulation. When I type the code m = Calc_tau0(m), an error occurs and prompts "No bathymetry data in the grid to calculate the tau0 coefficients". However, I have complete the interpolation by carrying out the previous step m = interp(m, DEMFILE, 'type', 'depth'). In the msh class, the attribute of b is not empty and the value of b is not always 0. I don't know why this happened. The following is my code,

% Example: Mesh the PRD
clearvars; clc;
addpath(genpath('utilities/'))
addpath(genpath('datasets/'))
addpath(genpath('m_map/'))
%% STEP 1: set mesh extents and set parameters for mesh.
bbox = [111.5 116.5; % lon_min lon_max
18 24]; % lat_min lat_max
min_el = 150; % minimum resolution in meters.
max_el = 5e4; % maximum resolution in meters.
max_el_ns = 200; % maximum resolution nearshore in meters.
grade = 0.3; % mesh grade in decimal percent.
R = 3; % number of elements to resolve feature width.
%% STEP 2: specify geographical datasets and process the geographical data
%% to be used later with other OceanMesh classes...
dem = 'topo15_compressed.nc'; coastline = 'GSHHS_f_L1';
gdat = geodata('shp',coastline,'dem',dem,'h0',min_el,...
'bbox',bbox);

% NOTE: You can plot the shapefile with bounding box by using the
% overloaded plot command:
% plot(gdat,'shp'); plot(gdat,'dem');
%% STEP 3: create an edge function class
fh = edgefx('geodata',gdat,...
'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',fh,'bou',gdat,'plot_on',1,'nscreen',1);
mshopts = mshopts.build;
m = mshopts.grd;
% plot(m,'tri');

DT = 2 ; % Stable timestep
TS = '01-Aug-2012 00:00' ; % start time of simulation
TE = '31-Nov-2012 00:00' ; % end time of simulation
DEMFILE = 'topo15_compressed.nc';
BUOYFILE = 'Gridded_N_values.mat';
TPXO9 = 'h_tpxo9.v1.nc';
CONST = 'major8';

m = interp(m, DEMFILE, 'type','depth');
plot(m,'b')

m = CheckTimestep(m,DT);

m = renum(m);

m = make_bc(m,'auto',gdat); % make the nodestring boundary conditions
% plot(m,'bd',0);
% makens(m,'outer',1)

m = Calc_tau0(m);

Looking forward to your reply.

Best wishes.
Jiangchao

SRTM15+ change version

SRTM15+ changed version, from v2, to v2.1. You may correct the link in your README.md to:

ftp://topex.ucsd.edu/pub/srtm15_plus/SRTM15+V2.1.nc

Add checker to deal with cases where mesh is unable to become traversable

I noticed that errors with traversability even after passing Make_Mesh_Traversable can occur, but this typically only happens after complicated processes such as mesh merging and never from a mesh that has been just generated by meshgen.build.

I found that this traversability problem seems to becoming from the fact that a triangle is specified as being made up only two unique nodes (i.e., one of the nodes is entered twice in a triangle).

As Keith has mentioned previously, I agree that we need to introduce some sort of parser/listener to automatically deal with and/or notify about such strange triangulations or NaNs being introduced into the mesh.

Mesh Corner Elements coordinate error

Hi,

Apologies for raising another issue, so far I have been pretty impressed with the adaptability, I have read through the manual and cannot seem to find anything that I think might help other than perhaps using pfix in meshgen at the corners.

So I started by using example_2_NY as a starting point and have been progressively worked from there add elements I want. However all the nested meshes I produce either one or both of the eastern boundaries seem to end up with some sort of error, whereby the coordinates for those elements seem to be 0,0.

Attached is my matlab file:
mesh_build_gc_nest.txt

The reason I think that it is a coordinate error is when I import the file to Blue Kenue it looks like the below.

Thanks,

Drak

Mesh boundaries doesn't assign correctly using interpFP

Hi,
I have been trying to create the mesh boundaries using "makens" at the end of mesh workflow and noticed too issues/glitches , may be you guys can take a look at it.

1. For the entire North Atlantic mesh we noticed that it create 1200+ mainland nodestrings using the "auto" option and the only possible way around we figured out was to use "islands" option. That way we will load the mesh in SMS and delete the most outer island to recreate mainland and ocean boundary only.

2. Using the interpFP interpolation scheme and later assigning the mesh boundaries using "makens" for the small rivers like Potomac, it will not recognize the ocean and mainland boundary. And it assigns only the island boundaries by using "auto" or "islands".
   meshFull_FP = interpFP(meshFull,gdatovFull,muwFull,gdatuwFull,0.05,0);
   It will be great if you kindly take a look at it.
   Thanks

user guide needs updating

Need to add BoundCr method details and updated details about MakeBC...

Msh.cleaning destroys previous bathy

If you execute msh.clean, the routine will not correctly map you bathymetry back onto the msh vertices correctly.

Slack invite needs refresh

Apparently it dies every 30 days unless we pay?!

Question about automatically determining the time step

hi @krober10nd

I am confused about the method of automatically determining the time step.

In Example_3_ECGC, the time step was set to 0.

When I followed this setting, I found that the time step dt calculated based on the nearshore resolution is very small.

Therefore, the mesh resolution calculated according to the original code of Example 3 is very small, and it runs very slowly.

When I set dt to 2, as in Example_2_NY (Encourage mesh to be stable at a 2 s time step), It turned out to look good.

If using the automatic step size method, is such a small value of dt reasonable to use in the simulation? Or should I specify it myself?

Best wishes
Jiangchao

Cfl warning when writing

If a mesh grossly violates the cfl condition, let's output a warning when writing.

Inverse CFL handling

Mainly for functionality for SCHISM, it seems that it would be good to include the inverse of the current 'dt' function in edgefx and CheckTimestep functionality that will make the mesh finer where the CFL is too small (They suggest something like CFL > 0.5 everywhere for SCHISM so that the Lagrangian advection scheme stays accurate).

Ideally we can should be able to set a band so that CFL should be at least 0.5 and less than say 6 or whatever.

For example in CheckTimestep we can add points somewhere close to where the CFL is too small in addition to deleting points where CFL is too high like we do now.

Add a tool to remove pits from the mesh bathymetry

Pits in the bathymetry can make the wetting/drying portion go unstable.
Here's a little code contribution to replace pits with the average bathymetry
of the neighboring nodes, if it may be useful.

% find pits in OceanMesh2D mesh bathymetry
%   Hamish Bowman 10 Sept. 2020

m = msh('fort14_with_pits');
m_orig.b = m.b;

% find node points likely to see some wetting and drying
% may want to get rid of upper bound (here 5m above MSL)
intertidal = find(m.b < 3 & m.b > -5);

sz = size(m.t);
is_pit = zeros(size(m.b));  % init array
mean_neighbors = nan(size(m.b));  % init array

% brute force
% note: this modifies along the way, not atomically at the end
%    so pits may come and go as the neighboring pit nodes get filled
tic
for n = [intertidal']
   idx = find(m.t == n);
   [rows cols] = ind2sub(sz, idx);

   hits = m.t(rows, :);
   neighbors = unique(hits(find(hits ~= n)));

   % ignore spur nodes along the boundary (maybe skip all bdy nodes?)
   %   or run QA code to nuke spur elements first?
   if(length(neighbors) > 2 & m.b(n) > max(m.b(neighbors)))
      disp(['node ' num2str(n) ' (' num2str(m.b(n), '%.3f') ...
         ') is an intertidal pit with ' ...
         num2str(length(neighbors)) ' neighbors who average ' ...
	 num2str(mean(m.b(neighbors)), '%.3f') ' meters depth'])
      is_pit(n) = 1;
      mean_neighbors(n) = mean(m.b(neighbors));

      % replace if pit is more than 10 cm deep (does size of pit matter?)
      if(m.b(n) - max(m.b(neighbors)) > 0.10)
         disp(['  * diff: ' num2str(m.b(n) - max(m.b(neighbors)))])
         m.b(n) = mean(m.b(neighbors));
      end
   end
end
toc

disp([num2str(length(find(is_pit))) ' pits found in mesh bathymetry'])


% or to fill in pits after the scan,  (e.g. if scan is multi-threaded)
%for n = [find(is_pit)']
%   % replace if pit is more than 10 cm deep (does size of pit matter?)
%   if(m.b(n) - mean_neighbors(n) > 0.10)
%      m.b(n) = mean_neighbors(n);
%   end
%end


%%% you may want to iterate through this 2-10 times.

Meaning of the colors

First of all, thank you very much for this tool you developed! It is great!

I am trying to generate my first grid to later be used in SWAN. Everything is working great so far. However, I am a bit confused about the meaning of the colors when the program shows the mesh. For example the difference between the red and blue here:

I tried to find this information in the manual without success. Hopefully, I am not missing something obvious.

Thanks again.

Almir

CheckTimestep without pfix

Hi guys,

I was attempting to run CheckTimestep on a mesh without any fixed points. This gives me an error when it tries to clean up the mesh:
Error in msh/CheckTimestep (line 1962)
obj = clean(obj,'passive','proj',0,'pfix',pf,'con',con);

I fixed this on my own by just adding a

if ~isempty(obj.pfix)
            [pf(:,1),pf(:,2)] = m_ll2xy(obj.pfix(:,1),obj.pfix(:,2));

else
pf = [];
end

around line 1936 as a quick fix, but I'm not sure if this is going to mess anything up within the code.

Using mesh in SWAN with vertices numbering problem

Hi Keith, William, Joannnes,

Thanks for the smart mesh generator in coastal region. I was wondering if there is any program to sort the vertices ? Since when run swan using the generated mesh, it shows "numbering of vertices is not sequential in Triangle file". Vertices and triangles file in SWAN have to meet its numbering requirement. Thanks a lot for your patient explanation.

Kind regards,
Jin

Smooth outer does nothing

• Minimal working example

bbox = [0, 1; 0,1]; 
boubox = [0,0; 1,0; 1,1; 0,1; 0,0; NaN NaN ];
min_el = 1e3; 

gdat1 = geodata('pslg',boubox,'bbox',bbox,'h0',min_el);
fh1 = edgefx('geodata',gdat1,'g',0.15);

bbox2 = [-1, 2; -1,2]; 
boubox2 = [-1,-1; 2,-1; 2,2; -1,2; -1,-1; NaN NaN ];
min_el2 = min_el*5; 

gdat2 = geodata('pslg',boubox2,'bbox',bbox2,'h0',min_el2);
fh2 = edgefx('geodata',gdat2,'g',0.10);

mshopts = meshgen('ef',{fh2, fh1},'bou',{gdat2,gdat1},'plot_on',1,'itmax',5,'cleanup',0);
mshopts = mshopts.build;

m = mshopts.grd; 
plot(m)

Whether or not I comment out the call to smooth_outer in meshgen produces the same results. There should be a smooth transition out of the nest.

Excessive coastline smoothing

[dev BRANCH]

Using a shapefile containing a single poligon, I get some coastline areas better smoothed than others. A picture is worth a thousand words:

Is there any way to fix it?

Mixed triangular-quad elements & curved open boundary

Hi @krober10nd and @WPringle,

I am wondering if oceanmesh2D can handle mixed (tri and quad) grids and curved lines for open boundaries (instead of a rectangular box) as done in SMS? or is there any idea to make it possible?

Thanks a lot :)

Cheers,
Huy

user guide corrections for Calc_IT_Fric() and write() usage

Hi,

In the pdf user guide's Typical Workflow Example (sec. 6.9), line 30 of the listing says:

m = Calc_IT_Fric(m,'Nfname',BUOYFILE) ;

but the current version of the function does not have 'Nfname' as a varargin pair, it's mandatory. So if you try to run the code as presented it looks for a file called 'Nfname.mat', unsuccessfully.

Also TPXO is now up to version 10, you may want to update that link in the text of the section.

(by the way, the final version of WOA 2018 is now available for download since Aug 29th, I don't think it is announced anywhere but the netCDF internal metadata will show it)

---

In the 7.1.2 Methods help for write() it says:

however, the user can optionally decide
to write only a select few of these files by passing a comma-delimited list of the post-fix name
types (e.g., 14,15 to denote fort.14 or fort.15).

it looks for a string not an int, so: e.g '14','15' to denote fort.14 or fort.15 helps the reader avoid a few moments of confusion.

thanks,
Hamish

Block DEM reading

As in meshgen with the formation of points, it’d be useful to read in the DEM in blocks in geodata based on say an assumption the user has 8 gab of RAM and downsample to avoid overloading the RAM each block as they are read in

Rename makens to make_boundaries

I propose to rename the msh method makens to make_boundaries for easier understanding.

shoreline constraints

Hi Keith, William,

I am making a mesh over the wharf like the land, but I want to keep its shape (e.g. straight line, not weirs) An example is attached. Is OceanMesh2D capable of doing that?

Thanks, mates
Huy

Weir struct file missing

OceanMesh2D/Examples/Example_5b_JBAY_w_weirs.m

Line 23 in 2047bed

load weirs_struct

Add example with mesh addition/plus

...to highlight the patching idea rather than using the multiscale approach.

In this example, we should highlight the lock_dis parameter to show that the outside region is unchanged.