Automatically subdivide line elements based on distance to wells (or other vertices)

[Huite]

A relatively easy way to subdivide long lines based on distance to wells or other element vertices is by:

1. Projecting all (well) vertices onto a line segment and computing the (perpendicular) distance
2. Find where the projected point falls. Divide the segment in half if its length is larger than the distance to the vertex.
3. Repeat until all projections fall in a segment that is smaller than the distance to the vertex.

For example:

• V1 and V2 are the vertices of the line segment.
• P1 and P2 are the vertices to two wells
• The orange and green lines indicate the projection
• All other markers denote the created subdivision vertices

This has some nice properties: we don't have to introduce any logic when projections are close to each other, or when a vertex is far away:

def _halve_segments(vertices, segments, index):
    """Halve segments at indices provided by `index`."""
    halved = segments[index] * 0.5
    return np.insert(vertices, index, vertices[index] - halved)


def bisect(vertices, projections, distances):
    """
    Iteratively half a line segment, defined by vertices, until the line
    segments are smaller than the provided distances. The result is such that
    every projected point falls within a segment that is shorter than the
    distance to the point.
    
    Parameters
    ----------
    vertices: np.array of floats
        Location of the vertices, expressed as length along the line.
    projections: np.array of floats
        Location of the projected points, expressed as length along the line.
    distances: np.array of floats
        Distance of the points to the line.
    
    Returns
    -------
    refined_vertices: np.array of floats
        Location of new vertices along line.
    """
    def find_split(vertices, projections, distances):
        index = np.searchsorted(vertices, projections)
        segments = np.diff(vertices, prepend=0.0)
        to_split = segments[index] > distances
        return np.unique(index[to_split]), segments, to_split.any()
    
    index, segments, to_split = find_split(vertices, projections, distances)
    while to_split:
        vertices = _halve_segments(vertices, segments, index)
        index, segments, to_split = find_split(vertices, projections, distances)
    return vertices

However, just like in mesh generation with quadtrees or octrees, the line segments are not always "balanced" (a 2:1 ratio at max between neigbhors), see the first segment with the second:

Fortunately, this is not difficult to implement in 1D. I'm not sure if it's worthwhile...

def balance(vertices):
    """
    Iteratively half line segments until each segment is at most twice as big
    as a neighboring segment.

    Parameters
    ----------
    vertices: np.array of floats
        Location of the vertices, expressed as length along the line.

    Returns
    -------
    refined_vertices: np.array of floats
        Location of new vertices along line.
    """
    def find_split(vertices):
        segments = np.diff(vertices)
        left = segments[:-1]
        right = segments[1:]
        to_split = np.zeros(segments, dtype=bool)
        # 0.49 should be sufficient to avoid floating point roundoff
        to_split[:-1] = (0.49 * left) > right
        to_split[1:] = (0.49 * right) > left
        index = np.argwhere(to_split) + 1  # TODO: check
        return index, segments, to_split.any()
 
    index, segments, to_split = find_split(vertices)
    while to_split:
        vertices = _halve_segments(vertices, segments, index)
        index, segments, to_split = find_split(vertices)
    return vertices

Some other thoughts:

• I think it makes sense to define a "refine" method on the model object, because it should only be done after all elements are added to the model. Then:

• A collection of all vertices can be made.

• Project these on every line segment and compute distance.

• Split every segment with the logic above.

• Replace the element in the model by one with the additional vertices inserted (?)

• It's also possible to refine and "balance" a linestring in its entirety by expressing the vertices as distance measured along the linestring.

• Maybe (only) somewhat useful when two segments have very obtuse angles with each other?

• There might be an indexing trick to get this to work when the geometry forms a ring (e.g. numpy take has a "wrap" mode).

[dbrakenhoff]

Added link to this issue here: timflow-org/timflow#42