River meshing

[SorooshMani-NOAA]

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.

[krober10nd]

Hi Soroosh, I didn't finish the CGAL call (in cpp) to support constraining edges in the triangulation. If you're interesting in learning more about that, I can point you in the right direction.

I don't agree that the meshers should just "mesh". For numerical simulation, the user needs to understand how best to construct the model. Besides, this requires a small number of parameter selections inspired by the downstream application.

Anyway, I believe that automation shouldn't result in a loss of knowledge.

[SorooshMani-NOAA]

Hi Keith, thank you for your response.

I started looking into using CGAL directly too. They have a nice SWIG binding package which covers most of the triangulation and meshing use cases (I still need to try it!).

In any case, I was wondering if oceanmesh already supports constrained triangulation-only use case? In my case there's an algorithm developed by VIMS team that provides optimal node location for river/watershed mesh. I just need a library to triangulate it. I failed to do so with jigsawpy through ocsmesh, so I tried oceanmesh where I was unsuccessful again. So I started looking into CGAL direct use.

If you believe this use case is not something relevant to oceanmesh I can close the ticket.

[krober10nd]

Soorosh,

oceanmesh uses CGAL to do the Delaunay triangulation already.

What I'm saying is that I did not finish adding the functionality to constrain edges in the triangulation, as noted.

In this folder you can see the existing interface with CGAL.
https://github.com/CHLNDDEV/oceanmesh/tree/master/oceanmesh/cpp

additional code would need to be written to constrain edges.

Anyway, I think meshing rivers would ideally be a supported use-case.

[SorooshMani-NOAA]

I think we're kind of saying the same thing, but there's a miscommunication:

Anyway, I think meshing rivers would ideally be a supported use-case.

The reason I thought you might not be considering this use case was what you said in the earlier comment that:

I don't agree that the meshers should just "mesh". For numerical simulation, the user needs to understand how best to construct the model. Besides, this requires a small number of parameter selections inspired by the downstream application.

Anyway, I believe that automation shouldn't result in a loss of knowledge.

oceanmesh uses CGAL to do the Delaunay triangulation already.

I know! And actually oceanmesh inspired me to look into using CGAL directly, and it looks like a great library with a lot of functionalities!

In this folder you can see the existing interface with CGAL.
https://github.com/CHLNDDEV/oceanmesh/tree/master/oceanmesh/cpp

additional code would need to be written to constrain edges.

I see that you're just inserting points in the triangulation, but not the constraint edges. So it looks relatively straightforward. However, is there any reason you wrote custom code to use CGAL? Was any functionality missing from the SWIG generated bindings?

[krober10nd]

However, is there any reason you wrote custom code to use CGAL? Was any functionality missing from the SWIG generated bindings?

I had a lot of difficulty implementing and deploying SWIG and had some prior experience with pybind11.