Added points which are on the boundary are added to the constrained segments

[Kevin-Mattheus-Moerman]

@DanielVandH I am working on a minimal example (only have a clunky one currently), but I've noticed that when one has added interior points that happen to be exactly on a boundary edges, they are then included/dded to the constrained segments provided. This seems to me to be undesirable behaviour. It would be better to exclude both "in" and "on" points from the triangulation.

[DanielVandH]

I'm not sure why this is not ideal. The point being added is indeed on a segment and so the segment should get split into two. This is a very important feature in e.g. mesh refinement to avoid long edges. Similarly in e.g. cell simulations it's not uncommon to have vertices collinear with each other initially, where again this feature is important. Ultimately if you would like to avoid this feature, it's the users responsibility. I do not intend to change triangulate to ignore user input in this way.

…

On Wed, 24 Apr 2024, 2:59 pm Kevin Mattheus Moerman, < ***@***.***> wrote: @DanielVandH <https://github.com/DanielVandH> I am working on a minimal example (only have a clunky one currently), but I've noticed that when one has added interior points that happen to be exactly on a boundary edges, they are then included/dded to the constrained segments provided. This seems to me to be undesirable behaviour. It would be better to exclude both "in" and "on" points from the triangulation. — Reply to this email directly, view it on GitHub <#105>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AWZPH4HBT4BN4GZQGBXQHVTY663D3AVCNFSM6AAAAABGW7NZQ6VHI2DSMVQWIX3LMV43ASLTON2WKOZSGI3DCMZXG4ZDGMQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[Kevin-Mattheus-Moerman]

@DanielVandH thanks. I understand what you are saying. However, I think when one uses constraints, I interpret them as "edges" e.g. I want an edge from point n to point m. However if a new point p is added on the edge n<->m then instead we get the edge n<->p and p<->m. I think you what you've implemented is boundary_nodes which form the boundary, rather than boundary edges. So since n and m will be part of the output boundary as per the constrain, you'd say that is fine.
Is there a way to perhaps add an option though to say that these are the sole boundary points, i.e. to see them as contrained edges? Perhaps it is a simple optional parameter, e.g. constraint_edges = true? That way it would not alter the behaviour for the other applications you mentioned.
I've implemented a fix now (albeit inefficient), so this is not too urgent. But thanks anyway for this great package and your time looking into this.

[Kevin-Mattheus-Moerman]

Below are some test cases for what I built with this in case you are interested:

see also regiontrimesh here.

[DanielVandH]

Hm. I see your point, yes I agree that this could be nice. I think I have an efficient way to implement this. Let me try and get an example working soon and I'll ask for your input.

[DanielVandH]

Haven't been able to get much progress with this since there's a lot of technical details here to do this efficiently (and avoid the annoying amount of edge cases - this package could be so much smaller if not for edge cases :) ), but I did get some ideas down that maybe you are interested in going forward with. I don't really have much time to delve into it much more than this though, sorry.

The idea would be that you first triangulate, and then afterwards you just call enforce_constraints!(tri, segments) (or a keyword argument, as you suggest, would be used to just call enforce_constraints! at the end of triangulate I suppose).

using DelaunayTriangulation: integer_type,
    contains_unoriented_edge,
    is_boundary_node, 
    get_boundary_chain

"""
    fuse_edge!(tri::Triangulation, i, j, r)

Given two vertices `(i, j)` separated only by a single vertex `r`, deletes
that vertex and joins the two vertices with an edge. 
"""
function fuse_edge!(tri::Triangulation, i, j, r) # this is the inverse of split_edge!
    # Problem with this function is that it assumes that the vertex adjoining
    # the edges (i, r) and (r, j) are the same, and similarly for the edges 
    # (j, r) and (r, i). This is not always the case. We need to modify 
    # delete_point! so that it works for boundary points (and for 
    # segment vertices) so that we can just do 
    #   delete_point!(tri, r)
    #   add_segment!(tri, i, j)
    w₁ = get_adjacent(tri, i, r) # assumption: get_adjacent(tri, i, r) == get_adjacent(tri, r, j)
    w₂ = get_adjacent(tri, j, r) # assumption: get_adjacent(tri, j, r) == get_adjacent(tri, r, i)
    delete_triangle!(tri, i, r, w₁)
    delete_triangle!(tri, r, j, w₁)
    delete_triangle!(tri, j, r, w₂)
    delete_triangle!(tri, r, i, w₂)
    add_triangle!(tri, i, j, w₁)
    add_triangle!(tri, j, i, w₂) # need to make sure ghost triangles are handled correctly here 
    return tri
end

"""
    get_segment_chain(tri::Triangulation, i, j) -> Vertices

Given two vertices `i` and `j`, returns the chain of vertices that connects them.
"""
function get_segment_chain(tri::Triangulation, i, j)
    if is_boundary_node(tri, i)[1] && is_boundary_node(tri, j)[1] # that this implies (i, j) should be a boundary edge might not always be true, 
                                                                  # i.e. there might be a segment chain that includes corner points and interior points. 
                                                                  # For the application in this issue, this is fine I believe.
        g = is_boundary_node(tri, i)[2] # the ghost vertex associated with the boundary edge
        chain = get_boundary_chain(tri, i, j, g) # assumes that (i, j) follows the orientation of the boundary - see ?get_boundary_chain
    else 
        chain = get_interior_segment_chain(tri, i, j)
    end     
    return chain
end
function get_interior_segment_chain(tri::Triangulation, i, j)
    I = integer_type(tri)
    chain = [I(i)]
    w = last(chain)
    while w ≠ j # assumes that the chain ends in j. Not always true - e.g. 
                # i could adjoin many segments and the below code will fail
        for k in get_neighbours(tri, w)
            if k ≠ i 
                if contains_segment(tri, w, k) # For the example mentioned above that this could be problematic
                                               # if w adjoins many segments (in which case we may never find j), 
                                               # one way to resolve this for your application would be to check for collinearity, e.g. force 
                                               #     is_collinear(cert),
                                               # where 
                                               #     cert = point_position_relative_to_line(tri, w, k, j)
                                               # If it's possible that the collinearity isn't exact, though, this can 
                                               # still be problematic. You could just try and find all chains until all 
                                               # possible chains end (being careful to not loop the chains), and then 
                                               # take the one that ends in j. Maybe that's what you do in regiontrimesh? I didn't look too much into the details for yours.
                    push!(chain, k)
                    break 
                end
            end # If there are no other segments, note that this will just cause the loop to go on forever.
        end
        w = last(chain)
    end
    return chain
end

"""
    enforce_constraint!(tri::Triangulation, i, j)

Given a segment `(i, j)` represented in `tri` by a sequence of 
collinear segments, removes the vertices defined by the collinear segments 
so that only the segment `(i, j)` remains.
"""
function enforce_constraint!(tri::Triangulation, i, j)
    chain = get_segment_chain(tri, i, j)
    # Note that chain[1] = i and chain[end] = j
    n = length(chain)
    for k in 2:(n-1) # this loop does nothing if contains_unoriented_edge(tri, i, j) is already true, since n = 2
        fuse_edge!(tri, i, chain[k+1], chain[k])
        # Here, note that we originally start with 
        #   i -- chain[2] -- chain[3] -- ... -- chain[end-1] -- j
        # and we want to end with
        #   i -- j
        # So, we start with fusing the chain 
        #   i -- chain[2] -- chain[3]
        # which leaves 
        #   i -- chain[3]
        # or 
        #   i -- chain[3] -- ... -- chain[end-1] -- j
    end
    # Will need to also take some care here, 
    # making sure that all the other segments are removed from 
    #   get_all_segments(tri)
    #   get_interior_segments(tri) (if appropriate)
    # and that (i, j) is added to the appropriate lists.
    return tri 
end

"""
    enforce_constraints!(tri::Triangulation, segments)

Given a list of segments `segments`, enforces the constraints. See [`enforce_constraint`](@ref).
"""
function enforce_constraints!(tri::Triangulation, segments)
    for e in each_edge(segments)
        i, j = edge_vertices(e)
        enforce_constraint!(tri, i, j)
    end
    return tri
end

It would be nice if delete_point! worked for segment vertices for this, but it doesn't at the moment due to there being too many edge cases that I kept encountering I had to give up on it.

[DanielVandH]

Did any of the above work out for you @Kevin-Mattheus-Moerman?

[DanielVandH]

Going to close this due to inactivity and I'm not sure about this feature. The functions I give above maybe work but they won't in general - maybe for your usecase you can take those ideas though.

Feel free to reopen if needed.