Geomesher and Pandamesh

[Huite]

Hi @cheginit,

I saw a while ago you forked pandamesh for geomesher -- happy to see it's been useful to you.

Recently, I've made a bunch of changes, see: https://deltares.github.io/pandamesh/api/changelog.html
Some might be of interest to you. The Gmsh fields stuff wasn't properly finished before and I've also increased the test coverage.

I think forking might not have been necessary: your changes / additions seem quite reasonable. For some cases, I've used different solutions e.g. the weighted area remapping logic is in https://github.com/Deltares/numba_celltree, exposed in a nicer way via https://deltares.github.io/xugrid/api/xugrid.OverlapRegridder.html#xugrid.OverlapRegridder

You might not be interested in all those dependencies; I'm fine with a shapely-backed function too (the reason I went through the numba_celltree effort was due to poor shapely intersection performance for large grids: see https://github.com/Huite/clipping-benchmarks)

I don't quite agree that triangle isn't maintained: it occasionally takes a while before the Python bindings are updated to newer Python versions, but it's used very widely and does very fast mesh generation (although with an unfortunate license situation). Currently I set gmsh as an optional (pip) dependency, but I could do the same for triangle and make both optional (though I generally distribute all dependencies via conda-forge so a user has access to everything; not sure if there's a good way to configure that).

Anyway, I'd be happy to discuss changes to suit your needs.

[cheginit]

Hi @Huite,

Thanks for reaching out!

The main reason that I forked it was that at the time I was planning to use it in part of my work and wanted to make further changes that would have meant making rewriting many parts of the code. However, later on I realized using triangulated mesh will not work out for my application and decided to work on curvilinear meshes instead. I published a package called CurviRiver. It hasn't been updated in a while and I have changed the methodology since then, but I will update it with the most recent methodology once I submit its accompanying paper. Parts of it are already public in CatSmoothing and PyWBT packages.

Regarding triangle, the main reason that I dropped it as a dep was that I looked at their repo and didn't see any recent activities, and I remember there were some dep conflicts with packages that I was using. Since gmsh was producing good results, I decided it's better to just keep that.

Currently I set gmsh as an optional (pip) dependency, but I could do the same for triangle and make both optional

I think it makes sense to have one default mesh generator as a hard dep and make the others optional. If you have some in-house mesh generator that does a decent job, using that as the hard dep and the rest as optional make sense too. But, I think, there should at least be one default mesh generator as a dep.

It's interesting that a couple of days ago, I stumbled upon numba_celltree. I was impressed by how effectively numba was used! There's this GitHub repo that profiles different regridding packages, I think it'd be interesting to see how your implementation in xugrid compares to those.

I looked at the clipping profiling that you linked, I wonder how shapely>=2 compares, since I think that's before its release. It's supposed to be much faster that shapely<2. So, it'd be interesting to see if that's really the case for this use case.

Anyway, I'd be happy to discuss changes to suit your needs.

Absolutely, if you have any ideas for curvilinear meshes that preserve anisotropy, I'd be interested in discussing it.

[Huite]

The main reason that I forked it was that at the time I was planning to use it in part of my work and wanted to make further changes that would have meant making rewriting many parts of the code. However, later on I realized using triangulated mesh will not work out for my application and decided to work on curvilinear meshes instead. I published a package called CurviRiver. It hasn't been updated in a while and I have changed the methodology since then, but I will update it with the most recent methodology once I submit its accompanying paper. Parts of it are already public in CatSmoothing and PyWBT packages.

Ah, right, that makes sense! Thanks for linking to Catsmoothing, it might solve a problem I'm dealing with!

It's interesting that a couple of days ago, I stumbled upon numba_celltree. I was impressed by how effectively numba was used! There's this GitHub repo that profiles different regridding packages, I think it'd be interesting to see how your implementation in xugrid compares to those.

Yeah, numba worked out pretty well there, I've had other attempts where I couldn't really fit it into a form that numba would deal with efficiently. But it's so easy to use from Python and distribute, compared to Cython, C++, etc.

Ah, interesting for the resample profiling repo. The focus seems on reprojection, but that can be done by reprojecting the UGRID topology to a common coordinate system, and then regridding. (It would also be worthwhile to compare regridding in the same coordinate system.)

I looked at the clipping profiling that you linked, I wonder how shapely>=2 compares, since I think that's before its release. It's supposed to be much faster that shapely<2. So, it'd be interesting to see if that's really the case for this use case.

As far as I know, the pygeos timings should be representative for shapely 2. I was quite surprised to see the performance difference, my guess is that the GEOS (C++) algorithms just have a lot more overhead (but also a lot more options and robustness) than a basic clipping algorithm. Back then, a colleague and I also did the benchmark in Julia, calling GEOS directly, and the timings were similar to the pygeos timings.

if you have any ideas for curvilinear meshes that preserve anisotropy

Gmsh has this concept of transfinite meshes, see this example. My feeling is that you still have to do most of the work though, you have to specify the opposite lines and the number of points. The primary usefulness is that it can be embedded in a larger triangular grid easily.

Some of my Deltares colleagues work on MeshKernel, with some Python bindings here. It requires clean input though.
(It's an abstract TO DO to make a nice Pythonic API for it, but I have many concrete TO DO's to tend to in the mean time...)

Not sure if I understand "preserve anisotropy" correctly, though!

[cheginit]

Ah, right, that makes sense! Thanks for linking to Catsmoothing, it might solve a problem I'm dealing with!

Awesome!

As far as I know, the pygeos timings should be representative for shapely 2.

I guess it shouldn't be far off, but they're also using numpy arrays for everything, so I am guessing that will also help.

Gmsh has this concept of transfinite meshes, see this example ... Some of my Deltares colleagues work on MeshKernel,

Ah, interesting. I think I worked with transfinite meshes while I was exploring Gmsh, but, as you said, the main issue is that you need to actually have the left and right banks, as well as the inflow and outflow edges, to being with. This is the main issue that makes it unsuitable for large scale applications and complex river geometries that can have holes too (braided rivers). The methodology that I developed is geared toward automating as much as possible with minimal user input and being efficient so it can be applied at large scales and works with braided systems.

Not sure if I understand "preserve anisotropy" correctly, though!

This is important for interpolating, for example, bathymetric survey data, onto river meshes. Preserving anisotropy ensures that changes in bathymetry in flow direction is preserved during interpolation. Essentially, you need to do the interpolation in the curvilinear mesh reference system (SN) instead of a Cartesian reference (XY). For example, you can check out this paper. Having accurate representation of a river centerline is the key to this regridding from XY to SN.