Rebuild geometry data after face deletion

docs/src/preprocessing/preprocessing.md

@@ -31,9 +31,9 @@ triangle = [125.0 375.0 250.0 125.0;
             175.0 175.0 350.0 175.0]
 
 # Delete all edges but one
-edge1 = deleteat!(TrixiParticles.Polygon(triangle), [2, 3])
-edge2 = deleteat!(TrixiParticles.Polygon(triangle), [1, 3])
-edge3 = deleteat!(TrixiParticles.Polygon(triangle), [1, 2])
+edge1 = delete_faces(TrixiParticles.Polygon(triangle), [2, 3])
+edge2 = delete_faces(TrixiParticles.Polygon(triangle), [1, 3])
+edge3 = delete_faces(TrixiParticles.Polygon(triangle), [1, 2])
 
 algorithm = WindingNumberJacobson()
 

src/TrixiParticles.jl

@@ -101,7 +101,7 @@ export trixi2vtk, vtk2trixi
 export RectangularTank, RectangularShape, SphereShape, ComplexShape
 export ParticlePackingSystem, SignedDistanceField
 export WindingNumberHormann, WindingNumberJacobson
-export VoxelSphere, RoundSphere, reset_wall!, extrude_geometry, load_geometry,
+export VoxelSphere, RoundSphere, reset_wall!, extrude_geometry, load_geometry, delete_faces,
        sample_boundary, planar_geometry_to_face
 export SourceTermDamping
 export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection,

src/preprocessing/geometries/polygon.jl

@@ -95,6 +95,40 @@ struct Polygon{NDIMS, ELTYPE}
         return new{NDIMS, ELTYPE}(vertices, edge_vertices, vertex_normals, edge_normals,
                                   edge_vertices_ids, min_corner, max_corner)
     end
+
+    function Polygon{NDIMS, ELTYPE}(vertices, edge_vertices, vertex_normals,
+                                    edge_normals, edge_vertices_ids,
+                                    min_corner, max_corner) where {NDIMS, ELTYPE}
+        return new{NDIMS, ELTYPE}(vertices, edge_vertices, vertex_normals, edge_normals,
+                                  edge_vertices_ids, min_corner, max_corner)
+    end
+end
+
+function rebuild_polygon_from_edges(edge_vertices, vertex_normals, edge_normals)
+    NDIMS = length(first(edge_normals))
+    ELTYPE = eltype(first(edge_normals))
+    vertices = SVector{NDIMS, ELTYPE}[]
+    vertex_ids = Dict{SVector{NDIMS, ELTYPE}, Int}()
+
+    edge_vertices_ids = map(edge_vertices) do edge
+        v1, v2 = edge
+        id1 = get!(vertex_ids, v1) do
+            push!(vertices, v1)
+            return length(vertices)
+        end
+        id2 = get!(vertex_ids, v2) do
+            push!(vertices, v2)
+            return length(vertices)
+        end
+
+        return (id1, id2)
+    end
+
+    min_corner = SVector([minimum(v[i] for v in vertices) for i in 1:NDIMS]...)
+    max_corner = SVector([maximum(v[i] for v in vertices) for i in 1:NDIMS]...)
+
+    return Polygon{NDIMS, ELTYPE}(vertices, edge_vertices, vertex_normals, edge_normals,
+                                  edge_vertices_ids, min_corner, max_corner)
 end
 
 function Base.show(io::IO, geometry::Polygon)
@@ -119,14 +153,25 @@ end
 
 @inline Base.eltype(::Polygon{NDIMS, ELTYPE}) where {NDIMS, ELTYPE} = ELTYPE
 
-@inline function Base.deleteat!(polygon::Polygon, indices)
-    (; edge_vertices, edge_normals, edge_vertices_ids) = polygon
+"""
+    delete_faces(geometry, indices)
+
+Return a geometry with the faces at `indices` removed and derived geometry data rebuilt.
+"""
+@inline function delete_faces(polygon::Polygon, indices)
+    edge_vertices = copy(polygon.edge_vertices)
+    vertex_normals = copy(polygon.vertex_normals)
+    edge_normals = copy(polygon.edge_normals)
 
     deleteat!(edge_vertices, indices)
-    deleteat!(edge_vertices_ids, indices)
+    deleteat!(vertex_normals, indices)
     deleteat!(edge_normals, indices)
 
-    return polygon
+    if isempty(edge_vertices)
+        throw(ArgumentError("cannot delete all polygon edges"))
+    end
+
+    return rebuild_polygon_from_edges(edge_vertices, vertex_normals, edge_normals)
 end
 
 @inline nfaces(mesh::Polygon) = length(mesh.edge_normals)

src/preprocessing/geometries/triangle_mesh.jl

@@ -171,15 +171,19 @@ end
 
 @inline face_normal(triangle, geometry::TriangleMesh) = geometry.face_normals[triangle]
 
-@inline function Base.deleteat!(mesh::TriangleMesh, indices)
-    (; face_vertices, face_vertices_ids, face_edges_ids, face_normals) = mesh
+@inline function delete_faces(mesh::TriangleMesh, indices)
+    face_vertices = copy(mesh.face_vertices)
+    face_normals = copy(mesh.face_normals)
 
     deleteat!(face_vertices, indices)
-    deleteat!(face_vertices_ids, indices)
-    deleteat!(face_edges_ids, indices)
     deleteat!(face_normals, indices)
 
-    return mesh
+    if isempty(face_vertices)
+        throw(ArgumentError("cannot delete all triangle mesh faces"))
+    end
+
+    vertices = collect(Iterators.flatten(face_vertices))
+    return TriangleMesh(face_vertices, face_normals, vertices)
 end
 
 @inline nfaces(mesh::TriangleMesh) = length(mesh.face_normals)

test/preprocessing/geometries/geometries.jl

@@ -54,6 +54,34 @@
         end
     end
 
+    @testset verbose=true "`delete_faces` Rebuilds Derived Data" begin
+        triangle = [0.0 1.0 0.5 0.0;
+                    0.0 0.0 0.7 0.0]
+
+        edge_only = TrixiParticles.delete_faces(TrixiParticles.Polygon(triangle), [1, 2])
+
+        @test TrixiParticles.nfaces(edge_only) == 1
+        @test length(edge_only.vertices) == 2
+        @test length(edge_only.vertex_normals) == 1
+        @test edge_only.min_corner == min.(edge_only.edge_vertices[1]...)
+        @test edge_only.max_corner == max.(edge_only.edge_vertices[1]...)
+
+        A = SVector(0.0, 0.0, 0.0)
+        B = SVector(1.0, 0.0, 0.0)
+        C = SVector(0.0, 1.0, 0.0)
+        D = SVector(1.0, 1.0, 0.0)
+        face_vertices = [(A, B, C), (B, D, C)]
+        face_normals = [SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, 1.0)]
+        mesh = TrixiParticles.TriangleMesh(face_vertices, face_normals, [A, B, C, D])
+
+        mesh = TrixiParticles.delete_faces(mesh, 1)
+
+        @test TrixiParticles.nfaces(mesh) == 1
+        @test length(mesh.vertices) == 3
+        @test length(mesh.edge_normals) == 3
+        @test mesh.face_vertices == [face_vertices[2]]
+    end
+
     @testset verbose=true "Real World Data" begin
         data_dir = pkgdir(TrixiParticles, "examples", "preprocessing", "data")
         validation_dir = pkgdir(TrixiParticles, "test", "preprocessing", "data")

test/preprocessing/packing/nhs_faces.jl

@@ -4,7 +4,7 @@
                     0.0 0.0 0.7 0.0]
 
         # Only use the third edge of the triangle, i.e. the edge from [0.1, 0.0] to [0.0, 0.0]
-        edge_aligned = deleteat!(TrixiParticles.Polygon(triangle), [1, 2])
+        edge_aligned = TrixiParticles.delete_faces(TrixiParticles.Polygon(triangle), [1, 2])
         edge_id = 1 # Only one edge in `Polygon`
 
         cell_sizes = [1.0 + sqrt(eps()), 0.1]
@@ -27,7 +27,8 @@
         end
 
         # Only use the first edge of the triangle, i.e. the edge from [0.0, 0.0] to [0.5, 0.7]
-        edge_arbitrary = deleteat!(TrixiParticles.Polygon(triangle), [2, 3])
+        edge_arbitrary = TrixiParticles.delete_faces(TrixiParticles.Polygon(triangle),
+                                                     [2, 3])
         edge_id = 1 # Only one edge in `Polygon`
 
         expected_ncells_bbox = [(1, 1), (6, 7)]