diff --git a/docs/src/preprocessing/preprocessing.md b/docs/src/preprocessing/preprocessing.md index 0518574b73..570c400279 100644 --- a/docs/src/preprocessing/preprocessing.md +++ b/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() diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 1829cf6977..919abe77d7 100644 --- a/src/TrixiParticles.jl +++ b/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, diff --git a/src/preprocessing/geometries/polygon.jl b/src/preprocessing/geometries/polygon.jl index c56315b729..6c0407a739 100644 --- a/src/preprocessing/geometries/polygon.jl +++ b/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) diff --git a/src/preprocessing/geometries/triangle_mesh.jl b/src/preprocessing/geometries/triangle_mesh.jl index ede02d4682..13512786b6 100644 --- a/src/preprocessing/geometries/triangle_mesh.jl +++ b/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) diff --git a/test/preprocessing/geometries/geometries.jl b/test/preprocessing/geometries/geometries.jl index 4f428d266c..77cc5b4003 100644 --- a/test/preprocessing/geometries/geometries.jl +++ b/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") diff --git a/test/preprocessing/packing/nhs_faces.jl b/test/preprocessing/packing/nhs_faces.jl index 292522a91b..a9b038f2ed 100644 --- a/test/preprocessing/packing/nhs_faces.jl +++ b/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)]