Deltares / xugrid

Xarray and unstructured grids
https://deltares.github.io/xugrid/
MIT License
64 stars 8 forks source link

Add utility to remove small faces #281

Closed Huite closed 2 months ago

Huite commented 2 months ago

I've already written a snap_nodes function. My use case was to pre-process some data in some way, but it can be directly used to get rid of tiny faces.

  1. snap the nodes
  2. replace the node indices in the face node connectivity
  3. remove resulting repeated vertices
  4. remove resulting 0-area'd faces
Huite commented 2 months ago

Not quite as simple, of course:

Current snapping has a "communicative" property, see #282, fixed by #283...:

image

Also highlights that the result is not necessarily convex:

image

Huite commented 2 months ago

A triangular mesh might be guaranteed convex, but might still have sliver faces.

Huite commented 2 months ago

For what it's worth, a basic implementation:

node_x = grid.node_x
node_y = grid.node_y
face_nodes = grid.face_node_connectivity
fill_value = grid.fill_value
inverse, new_x, new_y = xu.ugrid.snapping.snap_nodes(node_x, node_y, max_distance=3.0)
simplified = inverse[face_nodes]
simplified[face_nodes == fill_value] = fill_value
closed, isfill = xu.ugrid.connectivity.close_polygons(simplified, fill_value)

# Remove repeated nodes
valid = (closed[:, :-1] != closed[:, 1:]) & (simplified != fill_value)
m_per_row = (valid).sum(axis=1)

n, m = simplified.shape
index = xu.ugrid.connectivity.ragged_index(n, m, m_per_row)
new_face_nodes = np.full((n, m), fill_value)
new_face_nodes[index] = simplified[valid]

# Cases:
# m == 0 or 1: face has been contracted to a single point, safe to remove.
# m == 2: face has been contracted to a single edge, safe to remove.
# m >= 3: face remains valid.
out = new_face_nodes[m_per_row >= 3]

new_grid = xu.Ugrid2d(
    node_x=new_x,
    node_y=new_y,
    fill_value=fill_value,
    face_node_connectivity=out,
)

This'll simplify to some satisfaction in some cases:

image

But will generate slivers, possibly zero area triangles if the vertices are colinear:

image

Anything more sided than triangles is hopeless regardless.

A better solution is to take a stepwise approach to grid generation instead. Generate the grid, then collect the nodes, snap them to get rid of the smallest distances, then regenerate the grid with the vertices and the bounding polygon.

Thanks to Teun:

Before:

image

After:

image