chlnddev / oceanmesh Goto Github PK

View Code? Open in Web Editor NEW

46.0 16.0 14.0 16.23 MB

Automatic coastal ocean mesh generation in Python and C++. **under development**

License: GNU General Public License v3.0

Makefile 0.13% Python 91.97% CMake 0.38% C++ 6.34% Batchfile 1.18%

mesh-generation coastal-modelling python

oceanmesh's People

Contributors

Stargazers

Watchers

Forkers

mounirjaouhari stefanzieger ml14je dpzimmer alcoat tomas19 tomsail graig-sutherland pmav99 dong-dongd roy-shuvashish sgriffithjones jcharris

oceanmesh's Issues

Msh class

a msh class to provide basic meta data and storage capabilities
a port of the topological cleaning routine in the original version.

Support for reading in xy files

instead of shapefiles.

remove low valency elements

add to clean
https://github.com/CHLNDDEV/OceanMesh2D/blob/07f98a9ed3ee60d226f2aa099317c1fb9005a484/%40meshgen/meshgen.m#L969

Fix CI build system

ubuntu docker used to jump start the CI system no longer exists.

Installation error due to pybind11

When installing oceanmesh personalized python library with pip install -U -e . in the oceanmeshmater directorym I get a long big thread of errors due to pybind11.
I already installed the requirements sudo apt install libcgal-dev python3-pybind11.

The log with the errors from the installation of oceanmesh with pip is bellow. Did you get the same kind of errors? How could I solve this?

Installing collected packages: oceanmesh
  Running setup.py develop for oceanmesh
    ERROR: Command errored out with exit status 1:
     command: /home/nuno/resources/miniconda3/bin/python -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'/home/nuno/universidade/tese/unst_grid_generation/oceanmesh-master/setup.py'"'"'; __file__='"'"'/home/nuno/universidade/tese/unst_grid_generation/oceanmesh-master/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' develop --no-deps
         cwd: /home/nuno/universidade/tese/unst_grid_generation/oceanmesh-master/
    Complete output (107 lines):
    running develop
    running egg_info
    writing oceanmesh.egg-info/PKG-INFO
    writing dependency_links to oceanmesh.egg-info/dependency_links.txt
    writing requirements to oceanmesh.egg-info/requires.txt
    writing top-level names to oceanmesh.egg-info/top_level.txt
    reading manifest file 'oceanmesh.egg-info/SOURCES.txt'
    reading manifest template 'MANIFEST.in'
    writing manifest file 'oceanmesh.egg-info/SOURCES.txt'
    running build_ext
    -- OS detected: Linux
    -- CXX Compiler detected: GNU
    -- CMake additional search path for libraries:
    -- Using header-only CGAL
    -- Targetting Unix Makefiles
    -- Using /usr/bin/c++ compiler.
    -- Boost version: 1.65.1
    -- Boost include dirs: /usr/include
    -- Boost libraries:
    -- Using gcc version 4 or later. Adding -frounding-math
    -- CGAL version: 5.00.0.100
    -- Build type: Release
    -- USING CXXFLAGS = '-DVERSION_INFO=\"0.0.1\" -O3 -DNDEBUG'
    -- USING EXEFLAGS = ' '
    -- Requested component: MPFR
    -- Requested component: GMP
    -- Found pybind11 v2.0.1: /usr/include;/home/nuno/resources/miniconda3/include/python3.8
    -- Configuring done
    -- Generating done
    -- Build files have been written to: /home/nuno/universidade/tese/unst_grid_generation/oceanmesh-master/build/temp.linux-x86_64-3.8
    [ 25%] Building CXX object CMakeFiles/fast_geometry.dir/oceanmesh/cpp/fast_geometry.cpp.o
    [ 50%] Building CXX object CMakeFiles/delaunay_class.dir/oceanmesh/cpp/delaunay_class.cpp.o
    /home/nuno/universidade/tese/unst_grid_generation/oceanmesh-master/oceanmesh/cpp/fast_geometry.cpp: In function ‘pybind11::array unique_edges(pybind11::array_t<int, 17>)’:
    /home/nuno/universidade/tese/unst_grid_generation/oceanmesh-master/oceanmesh/cpp/fast_geometry.cpp:76:23: error: no matching function for call to ‘pybind11::buffer_info::buffer_info(int*, long unsigned int, std::__cxx11::string, int, std::vector<long int>&, std::vector<long int>&)’
                           ));
                           ^
    In file included from /usr/include/pybind11/pytypes.h:12:0,
                     from /usr/include/pybind11/cast.h:13,
                     from /usr/include/pybind11/attr.h:13,
                     from /usr/include/pybind11/pybind11.h:36,
                     from /usr/include/pybind11/complex.h:12,
                     from /home/nuno/universidade/tese/unst_grid_generation/oceanmesh-master/oceanmesh/cpp/fast_geometry.cpp:2:
    /usr/include/pybind11/common.h:287:5: note: candidate: pybind11::buffer_info::buffer_info(pybind11::buffer_info&&)
         buffer_info(buffer_info &&other) {
         ^~~~~~~~~~~
    /usr/include/pybind11/common.h:287:5: note:   candidate expects 1 argument, 6 provided
    /usr/include/pybind11/common.h:274:14: note: candidate: pybind11::buffer_info::buffer_info(Py_buffer*, bool)
         explicit buffer_info(Py_buffer *view, bool ownview = true)
                  ^~~~~~~~~~~
    /usr/include/pybind11/common.h:274:14: note:   candidate expects 2 arguments, 6 provided
    /usr/include/pybind11/common.h:270:5: note: candidate: pybind11::buffer_info::buffer_info(void*, pybind11::size_t, const string&, pybind11::size_t)
         buffer_info(void *ptr, size_t itemsize, const std::string &format, size_t size)
         ^~~~~~~~~~~
    /usr/include/pybind11/common.h:270:5: note:   candidate expects 4 arguments, 6 provided
    /usr/include/pybind11/common.h:262:5: note: candidate: pybind11::buffer_info::buffer_info(void*, pybind11::size_t, const string&, pybind11::size_t, const std::vector<long unsigned int>&, const std::vector<long unsigned int>&)
         buffer_info(void *ptr, size_t itemsize, const std::string &format, size_t ndim,
         ^~~~~~~~~~~
    /usr/include/pybind11/common.h:262:5: note:   no known conversion for argument 5 from ‘std::vector<long int>’ to ‘const std::vector<long unsigned int>&’
    /usr/include/pybind11/common.h:260:5: note: candidate: pybind11::buffer_info::buffer_info()
         buffer_info() { }
         ^~~~~~~~~~~
    /usr/include/pybind11/common.h:260:5: note:   candidate expects 0 arguments, 6 provided
    /home/nuno/universidade/tese/unst_grid_generation/oceanmesh-master/oceanmesh/cpp/fast_geometry.cpp: At global scope:
    /home/nuno/universidade/tese/unst_grid_generation/oceanmesh-master/oceanmesh/cpp/fast_geometry.cpp:79:16: error: expected constructor, destructor, or type conversion before ‘(’ token
     PYBIND11_MODULE(fast_geometry, m) {
                    ^
    CMakeFiles/fast_geometry.dir/build.make:62: recipe for target 'CMakeFiles/fast_geometry.dir/oceanmesh/cpp/fast_geometry.cpp.o' failed
    make[2]: *** [CMakeFiles/fast_geometry.dir/oceanmesh/cpp/fast_geometry.cpp.o] Error 1
    CMakeFiles/Makefile2:67: recipe for target 'CMakeFiles/fast_geometry.dir/all' failed
    make[1]: *** [CMakeFiles/fast_geometry.dir/all] Error 2
    make[1]: *** Waiting for unfinished jobs....
    /home/nuno/universidade/tese/unst_grid_generation/oceanmesh-master/oceanmesh/cpp/delaunay_class.cpp:66:16: error: expected constructor, destructor, or type conversion before ‘(’ token
     PYBIND11_MODULE(delaunay_class, m) {
                    ^
    CMakeFiles/delaunay_class.dir/build.make:62: recipe for target 'CMakeFiles/delaunay_class.dir/oceanmesh/cpp/delaunay_class.cpp.o' failed
    make[2]: *** [CMakeFiles/delaunay_class.dir/oceanmesh/cpp/delaunay_class.cpp.o] Error 1
    CMakeFiles/Makefile2:104: recipe for target 'CMakeFiles/delaunay_class.dir/all' failed
    make[1]: *** [CMakeFiles/delaunay_class.dir/all] Error 2
    Makefile:83: recipe for target 'all' failed
    make: *** [all] Error 2
    Traceback (most recent call last):
      File "<string>", line 1, in <module>
      File "/home/nuno/universidade/tese/unst_grid_generation/oceanmesh-master/setup.py", line 81, in <module>
        setup(
      File "/home/nuno/resources/miniconda3/lib/python3.8/site-packages/setuptools/__init__.py", line 153, in setup
        return distutils.core.setup(**attrs)
      File "/home/nuno/resources/miniconda3/lib/python3.8/distutils/core.py", line 148, in setup
        dist.run_commands()
      File "/home/nuno/resources/miniconda3/lib/python3.8/distutils/dist.py", line 966, in run_commands
        self.run_command(cmd)
      File "/home/nuno/resources/miniconda3/lib/python3.8/distutils/dist.py", line 985, in run_command
        cmd_obj.run()
      File "/home/nuno/resources/miniconda3/lib/python3.8/site-packages/setuptools/command/develop.py", line 34, in run
        self.install_for_development()
      File "/home/nuno/resources/miniconda3/lib/python3.8/site-packages/setuptools/command/develop.py", line 136, in install_for_development
        self.run_command('build_ext')
      File "/home/nuno/resources/miniconda3/lib/python3.8/distutils/cmd.py", line 313, in run_command
        self.distribution.run_command(command)
      File "/home/nuno/resources/miniconda3/lib/python3.8/distutils/dist.py", line 985, in run_command
        cmd_obj.run()
      File "/home/nuno/universidade/tese/unst_grid_generation/oceanmesh-master/setup.py", line 44, in run
        self.build_extension(ext)
      File "/home/nuno/universidade/tese/unst_grid_generation/oceanmesh-master/setup.py", line 76, in build_extension
        subprocess.check_call(
      File "/home/nuno/resources/miniconda3/lib/python3.8/subprocess.py", line 364, in check_call
        raise CalledProcessError(retcode, cmd)
    subprocess.CalledProcessError: Command '['cmake', '--build', '.', '--config', 'Release', '--', '-j2']' returned non-zero exit status 2.
    ----------------------------------------
ERROR: Command errored out with exit status 1: /home/nuno/resources/miniconda3/bin/python -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'/home/nuno/universidade/tese/unst_grid_generation/oceanmesh-master/setup.py'"'"'; __file__='"'"'/home/nuno/universidade/tese/unst_grid_generation/oceanmesh-master/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' develop --no-deps Check the logs for full command output.

ESPG=4326 is not recognised

This is one of the reasons that the CI actions don't pass:

Some of the examples in the README use EPSG = 4326 as an integer for the coordinate system.
However in edgefx.py, only strings are tested i.e. :

if crs == "EPSG:4326":

like in bathymetric_gradient_sizing_function() and wavelength_sizing_function()

if crs == "EPSG:4326" or crs == 4326:
seems to be a good fix

stereographic distortion close to the poles

This issue is related to what was mentioned in the PR for global meshes about the stereographic distortion.

Close to the poles, the triangles get distorted and don't have the right equilateralness.
As mentioned,

it has to do with the forces being too powerful

To solve this issue, I have implemented a parametric law in the mesh_generator():

def _parametric(lat):
    ones = np.ones(lat.shape)
    res = ((90 - lat) * 2 + 18) / 180 * np.pi
    return np.minimum(res, ones)

and illustrated it in the following python notebook

I will propose a PR with this parametric law. But I am still not convinced this is the best way to go..

The meshes look almost perfect closer to the poles now:

Keen to get another opinion on this.

Wavelength_sizing_function Error

The grid.values are given in SI units of metres in wavelength_sizing_function of oceanmesh/edgefx.py. These need to be translated to the grid coordinate system. Should be an easy fix?

elements to elements

This connectivity is required for the mesh cleaning step at the end. However, if there are disjoint/disconnected portions of the mesh my approach doesn't work...Any suggestions would be appreciated.

oceanmesh/oceanmesh/clean.py

Line 22 in c59f039

def _cell_to_cell(t):

def _cell_to_cell(t):
    """Cell to cell connectivity table.
        Cell `i` is connected to cells `ctoc[ix[i]:ix[i+1]]`
        By connected, I mean shares a mutual edge.
    Parameters
    ----------
    t: array-like
        Mesh connectivity table.
    Returns
    -------
    ctoc: array-like
        Cell numbers connected to cells.
    ix: array-like
        indices into `ctoc`
    """
    nt = len(t)
    t = np.sort(t, axis=1)
    # NB: we use order="F" to reshape because the `np.tile` command below
    e = t[:, [[0, 1], [0, 2], [1, 2]]].reshape((nt * 3, 2), order="F")
    trinum = np.tile(np.arange(nt), 3)
    j = _arg_sortrows(e)
    e = e[j, :]
    trinum = trinum[j]
    k = np.argwhere(~np.diff(e, axis=0).any(axis=1))
    ctoc = np.concatenate((trinum[k], trinum[k + 1]), axis=1)
    ctoc = np.append(ctoc, np.fliplr(ctoc), axis=0)
    ctoc = ctoc[np.argsort(ctoc[:, 0]), :]
    idx = np.argwhere(np.diff(ctoc[:, 0])) + 1
    idx = np.insert(idx, 0, 0)
    idx = np.append(idx, len(ctoc))
    return ctoc, idx

depth-based bound enforcement

Give the user the ability to bound mesh resolution as a function of 1 or more depth ranges.

Mesh resolution expanding from a point

and growing outward in the radial direction.

More edgefx

feature size
topographic length scale.
mesh size gradation.
cfl limiting

Implement CFL condition based on group velocity

for applications with phase resolving wave models.

Question: How the meshgen works

@WPringle, do any points get generated in the mesh gen code other than the ones initially created?

oceanmesh/oceanmesh/mesh_generator.py

Line 146 in b763b1a

def generate_mesh(domain, edge_length, **kwargs):

Refactor geodata

the class currently called Grid should be renamed Raster
a DEM is then a subclass of a raster.
Grid.project() should be named Raster.interpolate()
~~Geodata class should be called Vector and then is the parent class of Shoreline~~

black version update makes CI fail

@jcharris @krober10nd @sgriffithjones I see you had the same linting issue in the CI lint with black :

would reformat /home/runner/work/oceanmesh/oceanmesh/oceanmesh/idw.py
would reformat /home/runner/work/oceanmesh/oceanmesh/oceanmesh/geodata.py

Oh no! 💥 💔 💥
2 files would be reformatted, 29 files would be left unchanged.

I see 2 options here:

correct the scripts that fail
Downgrade black version for python 3.10 here

I guess 1 is the best one in the long run. I think we should also update black locally for python 3.11

Problem with distance_sizing_function

Hi guy, thanks for your work, please help me with this issue.
I try to run the examples you give in "README.md" or "test cases". Every time it runs into the code line with "distance_sizing_function" it has an error.

shoreline = om.Shoreline(fname, extent.bbox, min_edge_length)
edge_length = om.distance_sizing_function(shoreline,  #rate=0.15)

Error:

    raise ValueError("The points in dimension %d must be strictly "
ValueError: The points in dimension 1 must be strictly ascending

Slack link no longer active

Hi there. Looking to download and use oceanmesh for generating a triangular mesh for my DG-FEM code for modelling flow in and around submarine canyons, given by some fluid-depth function varying in space. I am however having issues with getting oceanmesh to install and work on my windows machine, and could benefit from help. Unfortunately, the slack channel link you posted seems to no longer be active.

Off-by-one error reading Etopo1 data

Hi,

I found an off-by-one error while trying to read Etopo1 data with DEM class. Well, first I need to convert Etopo data from 0-360 to the 180-180 referential, then I get this error upon reading:

bbox, min_edge_length = (-75.000, -70.001, 40.0001, 41.9000), 1e3  # from the example
dem = DEM("etopo1.nc", bbox)

Trace-back:

ValueError                                Traceback (most recent call last)
~/Github/ww3-brazil-wc/create_mesh.py in <module>
     22 
     23     fname = "/home/stringari/GIS/etopo1fix.nc"
---> 24     dem = DEM(fname, bbox)
     25 
     26     edge_length = distance_sizing_function(shore, max_edge_length=5e3)

~/Github/oceanmesh/oceanmesh/geodata.py in __init__(self, dem, bbox, verbose)
    501         dy = numpy.abs(lats[1] - lats[0])
    502         dx = numpy.abs(lons[1] - lons[0])
--> 503         super().__init__(
    504             bbox=bbox,
    505             dx=dx,

~/Github/oceanmesh/oceanmesh/grid.py in __init__(self, bbox, dx, dy, hmin, values, fill)
     66         self.nx = int(floor((self.bbox[1] - self.bbox[0]) / self.dx))
     67         self.ny = int(floor((self.bbox[3] - self.bbox[2]) / self.dy))
---> 68         self.values = values
     69         self.eval = None
     70         self.fill = fill

~/Github/oceanmesh/oceanmesh/grid.py in values(self, data)
    118             pass
    119         elif data.shape != (self.nx, self.ny):
--> 120             raise ValueError(
    121                 f"Shape of values {data.shape} does not match grid size {self.nx, self.ny}"
    122             )

ValueError: Shape of values (300, 114) does not match grid size (299, 113)

If I modify lines 66 and 67 of grid.py to

self.nx = int((self.bbox[1] - self.bbox[0]) // self.dx)
self.ny = int((self.bbox[3] - self.bbox[2]) // self.dy)

I can read both my file and the tests pass. Should I open a PR?

Caio

Handle non-square meshing domains

in which bbox is a 2-tuple.

River meshing

In order to get a better representation of rivers and watersheds in the coastal ocean model, I'm trying to generate a mesh for extracted river polygons:

These polygons are extracted from DEM and are not necessarily connected. Apart from that I'd like to use all the vertices and edges that the extraction algorithm has outputed to guide the meshing. I can use pfix for nodes, but is there anything similar for edges?

I know that probably I cannot rely on Shoreline and mesh sizing functions provided for "normal" coastal meshing, so I tried to hack it by defining mine (without worrying about performance for now). This is how I define the signed distance function

def get_dist(x):
    pts = gpd.points_from_xy(x[:, 0], x[:, 0])
    sign = pts.within(shp_geom) * -2 + 1
    dist = sign * pts.distance(shp_geom.boundary)
    return dist

sdf = om.Domain(extent.bbox, get_dist)

And this is how I define the edge length.

 edge_length = om.Grid(
    bbox=extent.bbox,
    dx=hmin,
    hmin=hmin,
    extrapolate=True,
    values=hmax,
    crs=gdf_river_polys.crs,
 )
edge_length.build_interpolant()

points, cells = om.generate_mesh(
    sdf,
    edge_length,
    pfix=df_river_coord[['lon', 'lat']].values,
    max_iter=20
)

I take the hmin to be the smallest of the initial edges and hmin to be the size of the largest edge. I'm still playing with what to specify at as edge_length function, but my ideal solution is to actually not have any size specified and the mesher to just triangulate the nodes specified as initial points pfix.

Issue with DEM()

Hi!

I am trying to figure out the different uses of DEM() since I have a netcdf that is not working. So I tried both the string and array options with an example that I know it works.

I have a netcdf, that if I load and plot with xarray:

plt.figure()
plt.contourf(ds.lon.data, ds.lat.data, ds.z.data)
plt.colorbar()

If I then use DEM() it all works correctly:

dem_etopo_built_nc = f"/home/filippog/Mesh/dem_etopo_built_{name}.nc"
dem_etopo_fine = om.DEM(dem_etopo_built_nc)
dem_etopo_fine.plot()

I wanted to get the same result giving an array as input.
1st ISSUE: you need to explicitly give bbox to DEM() since it fails to create it inside the function if not given (bbox=None) and then it crashes in lines 677-679 of geodata.py

bbox = (8, 10, 44, 45) # Liguria
dem_etopo_fine = om.DEM(ds.z.values, bbox=(8, 10, 44, 45))
dem_etopo_fine.plot()

I get a wrong dem:

After playing with the array I figured out that you need to give it as:
dem_etopo_fine = om.DEM(ds.z.values[::-1, :], bbox=(8, 10, 44, 45))
But cannot figure out why I need to invert the lon?

Any thoughts or ideas?

Thanks a lot!!

Meshgen

use serial and parallel mesh generator developed in SeismicMesh

rewrite dfs to avoid O(N^2)

so not quadratic complexity. Use NetworkX. https://networkx.org

def dfs(graph):
    seen = set()
    lifo = []
    start = random_node(graph)
    lifo.append(start)
    while lifo:
        n = lifo.pop()
        seen.add(n)
        for nei in neighbours(n, graph):
            if nei in seen:
                continue
            else:
                lifo.append(n)
    nflag = numpy.zeros(nnodes)
    nflag[list(seen)] = 1

multiscale meshing

We need Grid addition and min operations overloaded.
Mesh size grading

delete single connected elements

port a function that iteratively removes what we term "single connected elements" and friends
original implementation here: https://github.com/CHLNDDEV/OceanMesh2D/blob/Projection/utilities/Fix_single_connec_edge_elements.m

Update CI system

to use github workflows.

offset in sizing functions using dem (WL and bathy gradient)

The @setter for grid.values is misleading during the calculation of the sizing functions that require a dem.

The values grid.dx and grid.dy change here for the bathy gradient and here for the wavelength

Here an example for the world model where both the wavelength and bathy gradient have an offset:

The quick fix (without touching the Grid class) is to force in both functions the dx and dy in the sizing functions.

grid.dx = dem.dx
grid.dy = dem.dy

Here is the result after using the fix:

implement geolayer

User can specify the projection of the shapefile and DEM, and internally om will re-project this data into Mercator to perform the meshing.

Fix Mesh: Too many values to unpack

In the README.md, it is suggested that to ensure that the vertices of each element are arranged in a counter-clockwise fashion, one simply uses the command points, cells = fix_mesh(points, cells). I had the error of two many values to unpack. Upon investigation, I could see that fix_mesh returns three outputs. To fix this in my own code, I simply changed it to

points, cells, jx = fix_mesh(points, cells)

Sorry, just thought others would find that helpful. (I have changed that in mine, and will submit a pull request once I am satisfied with other edits I have made.