Enforcing nodes along the line?

[ivicajan]

Dear all,
is it possible to add (or use) feature similar to
hfun.add_feature which would enforce that nodes are exactly on the "line"?
In other words, I am trying to create mesh which will have nodes (hence sides of the triangles) along the coastline. It is important as I have river and would like to capture exactly river bank as well as some man-made structures. Seems to me it is softly imposed and not in a hard way. Would then another approach of generating the grid using "roi" exactly following the coastline as "wet" grid and then another grid from the coastline up to dry region (i.e. isoline). At the end option to merge those 2 grids might result in the perfectly aligned nodes (sides) along the wet/dry "coastline"?

Thanks in advance
Ivica

[SorooshMani-NOAA]

@ivicajan, Jigsaw has the capability to specify initial nodes that should not be touch. It doesn't prevent it from adding more nodes, but maybe that's good enough for your usecase. OCSMesh doesn't directly support it, but if in the last step instead of calling JigsawDriver from OCSMesh, you directly call JigsawPy, you can achieve this:

geom = ...  # Like before
hfun = ...  # Like before

geom_jig = geom.msh_t()
hfun_jig = hfun.msh_t()

# TODO: Make sure geom_jig.crs and hfun_jig.crs are the same

opts = jigsawpy.jigsaw_jig_t()
opts.hfun_scal = "absolute"
opts.hfun_hmin = np.min(hfun_jig.value)
opts.hfun_hmax = np.max(hfun_jig.value)
opts.mesh_dims = +2

# Initial mesh (e.g. points you'd like to embed in final mesh)
embedding_points_xy = ...  # List of (x,y) values for the points of interest, e.g. boundary nodes.
init_jig = jigsawpy.jigsaw_msh_t()
init_jig.mshID = 'euclidean-mesh'
init_jig.ndims = 2
init_jig.crs = init_jig.crs
init_jig.vert2 = np.array(
    [(coord, -1) for id, coord in embedding_points_xy],
    dtype=jigsawpy.jigsaw_msh_t.VERT2_t
)

# Empty mesh object
mesh_jig= jigsawpy.jigsaw_msh_t()
mesh_jig.mshID = 'euclidean-mesh'
mesh_jig.ndims = 2
mesh_jig.crs = init_jig.crs

jigsawpy.lib.jigsaw(
    opts, geom_jig, mesh_jig,
    init=init_jig,
    hfun=hfun_jig)

ocsmesh.utils.finalize_mesh(mesh_jig, sieve=None)

mesh = ocsmesh.Mesh(mesh_jig)
mesh.write(...)

I usually used this approach with a whole init_jig mesh instead of only passing the points of interest. But I think passing the points only should work too. Please let me know how it goes. Note that here the -1 in definition of points is crucial. It specifies to Jigsaw that these points are final and should not be touched (or at least that's how I understand it!)

[SorooshMani-NOAA]

You can specify fixed edges with the same -1 tag, but with Jigsaw it's not guaranteed to be fixed (at least in my experience).

The script above is basically the simplest for of example I think, I actually copy pasted the relevant parts from the OCSMesh's CLI scripts where this approach is used, e.g. https://github.com/noaa-ocs-modeling/OCSMesh/blob/main/ocsmesh/cli/remesh.py.

This is an example I worked on a while ago: local_remesh_example.zip. I don't have a self contained example, so you need to replace the DEMs and meshes with your own.

[ivicajan]

Thanks for reply

I tried your approach but for some strange reason I can't get it right. Normal meshing is fine,
however, jig mesh with the additional points is empty.
Not sure if I put points correctly in .vert2

Please, if you could have a look and enlighten me what I am doing wrong.
I put Jupyter notebook and all needed files at the same place, in the:
https://github.com/ivicajan/OCSmesh_sandpit/blob/main/sandpit_v1.ipynb

Thanks a lot!

[SorooshMani-NOAA]

@ivicajan after running the notebook, it seems to me that the initial points are not at the same location as the domain (maybe wrong crs?):

Mesh (without init)	Init points

[ivicajan]

Thanks a lot, you are the legend!
The issue was within rasters, and for some strange reason I was not able to use standard GDA EPSG:3112 projection but needed all into ZONE 50S EPSG:32750.
However, after doing that another transformation they are all sitting on the same page and are used inside the meshing.

Now I am trying to impose edge as well.
If the coastline is simple segment type then I should connect them in a consequence ordered way [1,2], [2,3] etc?

I updated notebook and rasters for your reference as well.

Cheers,
Ivica

[SorooshMani-NOAA]

No problem, I'm glad that it helped! Yes for segments of the boundary you just need to set sequence of nodes as you said to the edge2 of the msh_t object.