[Proof of Concept] NHS only include intersecting faces

[LasNikas]

As discussed with @efaulhaber, there's no real improvement in checking whether a face intersects a cell in the face-embedding cell grid. However, this PR is used to store the methods.

To show that it doesn't really improve performance, here's a quick showcase. In both cases, 2D and 3D, all the faces are almost the same size.

3D (geometry with ca. 25k faces)

julia> analyze_face_surplus_in_nhs(geometry, ratios)
10.51 % when face is ca. 4.0 times the search radius
mean_face_size ≈ 0.211, search_radius ≈ 0.05

3.14 % when face is ca. 2.0 times the search radius
mean_face_size ≈ 0.211, search_radius ≈ 0.11

0.73 % when face is ca. 1.0 times the search radius
mean_face_size ≈ 0.211, search_radius ≈ 0.21

0.28 % when face is ca. 0.67 times the search radius
mean_face_size ≈ 0.211, search_radius ≈ 0.32

0.12 % when face is ca. 0.5 times the search radius
mean_face_size ≈ 0.211, search_radius ≈ 0.42

0.05 % when face is ca. 0.33 times the search radius
mean_face_size ≈ 0.211, search_radius ≈ 0.63

2D (geometry with ca. 63 edges)

julia> analyze_face_surplus_in_nhs(geometry, ratios)
17.7 % when face is ca. 4.0 times the search radius
mean_face_size ≈ 0.1, search_radius ≈ 0.02

7.67 % when face is ca. 2.0 times the search radius
mean_face_size ≈ 0.1, search_radius ≈ 0.05

2.96 % when face is ca. 1.0 times the search radius
mean_face_size ≈ 0.1, search_radius ≈ 0.1

1.09 % when face is ca. 0.67 times the search radius
mean_face_size ≈ 0.1, search_radius ≈ 0.15

1.27 % when face is ca. 0.5 times the search radius
mean_face_size ≈ 0.1, search_radius ≈ 0.2

0.0 % when face is ca. 0.33 times the search radius
mean_face_size ≈ 0.1, search_radius ≈ 0.3

julia> particle_spacing = 0.01
0.01

julia> @benchmark SignedDistanceField($geometry, $particle_spacing; always_true=$true)
BenchmarkTools.Trial: 4602 samples with 1 evaluation.
 Range (min … max):  790.334 μs …   7.876 ms  ┊ GC (min … max):  0.00% … 86.66%
 Time  (median):     979.104 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):     1.086 ms ± 639.109 μs  ┊ GC (mean ± σ):  10.07% ± 13.43%

   ▄██▂                                                         ▁
  █████▆▆▃▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▁▃▃▄▃▃▃▃▄▄▅▃▆▆▅▄▅▄▆▅▅▁▃▆▄▅▅▆▃▅▅ █
  790 μs        Histogram: log(frequency) by time       4.77 ms <

 Memory estimate: 2.68 MiB, allocs estimate: 2892.

julia> @benchmark SignedDistanceField($geometry, $particle_spacing; always_true=$false)
BenchmarkTools.Trial: 5184 samples with 1 evaluation.
 Range (min … max):  702.583 μs …   9.633 ms  ┊ GC (min … max):  0.00% … 89.30%
 Time  (median):     863.395 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   963.647 μs ± 615.840 μs  ┊ GC (mean ± σ):  10.48% ± 13.39%

  ▁▄█▆▁                                                         ▁
  █████▆▆▅▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▁▁▁▃▅▁▄▄▅▄▄▄▃▄▅▄▅▅▅▅▄▄▅▅▄▅▄▅▄▄▁▄▅▅ █
  703 μs        Histogram: log(frequency) by time       4.43 ms <

 Memory estimate: 2.60 MiB, allocs estimate: 2457.

function bbox_size(face, geometry::TrixiParticles.TriangleMesh)
    min_corner = minimum(stack(geometry.face_vertices[face]), dims=2)
    max_corner = maximum(stack(geometry.face_vertices[face]), dims=2)

    return max_corner .- min_corner
end

function bbox_size(edge, geometry::TrixiParticles.Polygon)
    min_corner = minimum(stack(geometry.edge_vertices[edge]), dims=2)
    max_corner = maximum(stack(geometry.edge_vertices[edge]), dims=2)

    return max_corner .- min_corner
end


function analyze_face_surplus_in_nhs(geometry, ratios)
    mean_face_size = 0.0
    for face in TrixiParticles.eachface(geometry)
        face_size = bbox_size(face, geometry)
        mean_face_size += TrixiParticles.norm(face_size)
    end

    mean_face_size /= TrixiParticles.nfaces(geometry)

    for search_radius_to_face_size in ratios

        particle_spacing = mean_face_size * search_radius_to_face_size
        search_radius = mean_face_size * search_radius_to_face_size

        nhs = TrixiParticles.FaceNeighborhoodSearch{ndims(geometry)}(; search_radius)
        TrixiParticles.initialize!(nhs, geometry; always_true=true)

        nfaces_always_true = sum(keys(nhs.neighbors.hashtable)) do key
            length(nhs.neighbors[key])
        end


        nhs = TrixiParticles.FaceNeighborhoodSearch{ndims(geometry)}(; search_radius)
        TrixiParticles.initialize!(nhs, geometry; always_true=false)

        nfaces_not_always_true = sum(keys(nhs.neighbors.hashtable)) do key
            length(nhs.neighbors[key])
        end

        ratio = (nfaces_always_true - nfaces_not_always_true) / nfaces_always_true * 100

        println("$(round(ratio, digits=2)) % when face is ca. $(round(1/search_radius_to_face_size, digits=2)) times the search radius")
        println("mean_face_size ≈ $(round(mean_face_size, digits=3)), search_radius ≈ $(round(search_radius, digits=2))")
        println("")

    end
end

[efaulhaber]

In our discussion, we also noted that plane_intersects_cell only checks if the plane of the triangle intersects the cell, not if the cell lies inside the triangle or outside, so there are always some unnecessary extra cells that are not filtered out correctly.