diff --git a/docs/src/preprocessing/preprocessing.md b/docs/src/preprocessing/preprocessing.md index 0518574b73..e3d0a3f3d8 100644 --- a/docs/src/preprocessing/preprocessing.md +++ b/docs/src/preprocessing/preprocessing.md @@ -269,6 +269,15 @@ For example: 0.0 1.0 ``` It is the user’s responsibility to ensure the points are ordered correctly. +For 2D `.asc` and `.dxf` files, `load_geometry` appends the first point by default +when it is not already repeated. This is only a convenience for complete, ordered +boundaries that omit the final duplicate point; it does not repair missing +segments, gaps, self-intersections, or incorrectly ordered points. Use +`load_geometry(file; close_curve=false)` for intentional open curves. Operations +that sample or classify a region, such as [`ComplexShape`](@ref), [`intersect`](@ref), +and [`setdiff`](@ref), require closed geometries. Boundary packing with +[`SignedDistanceField`](@ref) also requires a closed geometry, since it needs a +well-defined outside region. This format is easy to generate and inspect manually. ## DXF Format (.dxf) – recommended diff --git a/examples/preprocessing/complex_shape_2d.jl b/examples/preprocessing/complex_shape_2d.jl index fa2b762988..62de84ad49 100644 --- a/examples/preprocessing/complex_shape_2d.jl +++ b/examples/preprocessing/complex_shape_2d.jl @@ -7,8 +7,8 @@ # 3. Utilize the Winding Number algorithm to determine if points are inside or outside. # 4. Visualize the sampled particles and the winding number field. # -# The example uses an "inverted_open_curve" geometry, where standard inside/outside -# definitions might be ambiguous without a robust point-in-polygon test like winding numbers. +# The example uses a polygonal star geometry, where standard inside/outside +# definitions benefit from a robust point-in-polygon test like winding numbers. # ========================================================================================== using TrixiParticles @@ -16,7 +16,7 @@ using Plots particle_spacing = 0.05 -filename = "inverted_open_curve" +filename = "star" file = joinpath("examples", "preprocessing", "data", filename * ".asc") geometry = load_geometry(file) @@ -24,7 +24,6 @@ geometry = load_geometry(file) trixi2vtk(geometry) point_in_geometry_algorithm = WindingNumberJacobson(; geometry, - winding_number_factor=0.4, hierarchical_winding=true) # Returns `InitialCondition` diff --git a/src/preprocessing/geometries/geometries.jl b/src/preprocessing/geometries/geometries.jl index 228b21cc22..4fc9fd145e 100644 --- a/src/preprocessing/geometries/geometries.jl +++ b/src/preprocessing/geometries/geometries.jl @@ -4,6 +4,63 @@ include("io.jl") @inline eachface(mesh) = Base.OneTo(nfaces(mesh)) +""" + is_closed_geometry(geometry) + +Return `true` if a polygon or triangle mesh forms a closed region or surface. +""" +function is_closed_geometry(polygon::Polygon) + vertices_close = isapprox(first(polygon.vertices), last(polygon.vertices)) + vertices_close || return false + + expected_edges = count(1:(length(polygon.vertices) - 1)) do i + !isapprox(polygon.vertices[i], polygon.vertices[i + 1]) + end + + return nfaces(polygon) == expected_edges +end + +function is_closed_geometry(mesh::TriangleMesh) + return all(==(2), edge_face_counts(mesh)) +end + +function require_closed_geometry(geometry, operation) + is_closed_geometry(geometry) && return nothing + + msg = "`$operation` requires a closed geometry. " * + closure_error_detail(geometry) + + throw(ArgumentError(msg)) +end + +function closure_error_detail(polygon::Polygon) + if !isapprox(first(polygon.vertices), last(polygon.vertices)) + return "The first and last polygon vertices are different. " * + "If the vertices already trace a complete 2D boundary, construct or load " * + "the geometry with `close_curve=true`; otherwise provide a closed boundary." + end + + return "The polygon edge list does not form a complete closed curve." +end + +function closure_error_detail(mesh::TriangleMesh) + invalid_edges = count(!=(2), edge_face_counts(mesh)) + + return "Found $invalid_edges mesh edges with an incident-face count different from 2." +end + +function edge_face_counts(mesh::TriangleMesh) + edge_face_counts = zeros(Int, length(mesh.edge_vertices_ids)) + + for face_edges in mesh.face_edges_ids + edge_face_counts[face_edges[1]] += 1 + edge_face_counts[face_edges[2]] += 1 + edge_face_counts[face_edges[3]] += 1 + end + + return edge_face_counts +end + function Base.setdiff(initial_condition::InitialCondition, geometries::Union{Polygon, TriangleMesh}...) geometry = first(geometries) @@ -11,6 +68,7 @@ function Base.setdiff(initial_condition::InitialCondition, if ndims(geometry) != ndims(initial_condition) throw(ArgumentError("all passed geometries must have the same dimensionality as the initial condition")) end + require_closed_geometry(geometry, "setdiff") coords = reinterpret(reshape, SVector{ndims(geometry), eltype(initial_condition.coordinates)}, @@ -38,6 +96,7 @@ function Base.intersect(initial_condition::InitialCondition, if ndims(geometry) != ndims(initial_condition) throw(ArgumentError("all passed geometries must have the same dimensionality as the initial condition")) end + require_closed_geometry(geometry, "intersect") coords = reinterpret(reshape, SVector{ndims(geometry), eltype(initial_condition.coordinates)}, diff --git a/src/preprocessing/geometries/io.jl b/src/preprocessing/geometries/io.jl index bbccc276ec..a7c6106827 100644 --- a/src/preprocessing/geometries/io.jl +++ b/src/preprocessing/geometries/io.jl @@ -1,5 +1,5 @@ """ - load_geometry(filename; element_type=Float64) + load_geometry(filename; element_type=Float64, close_curve=true) Load file and return corresponding type for [`ComplexShape`](@ref). Supported file formats are `.stl`, `.asc` and `dxf`. @@ -18,16 +18,20 @@ For comprehensive information about the supported file formats, refer to the doc # Keywords - `element_type`: Element type (default is `Float64`) +- `close_curve`: Close 2D `.asc` and `.dxf` curves by appending the first point + when it is not already repeated. This assumes the vertices already + trace a complete, ordered boundary. Set this to `false` for intentional + open curves. Region sampling and classification reject open geometries. """ -function load_geometry(filename; element_type=Float64) +function load_geometry(filename; element_type=Float64, close_curve=true) ELTYPE = element_type file_extension = splitext(filename)[end] if file_extension == ".asc" - geometry = load_ascii(filename; ELTYPE, skipstart=1) + geometry = load_ascii(filename; ELTYPE, skipstart=1, close_curve) elseif file_extension == ".dxf" - geometry = load_dxf(filename; ELTYPE) + geometry = load_dxf(filename; ELTYPE, close_curve) elseif file_extension == ".stl" geometry = load(FileIO.query(filename); ELTYPE) else @@ -37,21 +41,21 @@ function load_geometry(filename; element_type=Float64) return geometry end -function load_ascii(filename; ELTYPE=Float64, skipstart=1) +function load_ascii(filename; ELTYPE=Float64, skipstart=1, close_curve=true) # Read the data from the ASCII file in as a matrix of coordinates. # Ignore the first `skipstart` lines of the file (e.g. headers). points = DelimitedFiles.readdlm(filename, ' ', ELTYPE, '\n'; skipstart)[:, 1:2] - return Polygon(copy(points')) + return Polygon(copy(points'); close_curve) end -function load_dxf(filename; ELTYPE=Float64) +function load_dxf(filename; ELTYPE=Float64, close_curve=true) points = Tuple{ELTYPE, ELTYPE}[] load_dxf!(points, filename) - return Polygon(stack(points)) + return Polygon(stack(points); close_curve) end function load_dxf!(points::Vector{Tuple{T, T}}, filename) where {T} diff --git a/src/preprocessing/geometries/polygon.jl b/src/preprocessing/geometries/polygon.jl index c56315b729..6215deef2e 100644 --- a/src/preprocessing/geometries/polygon.jl +++ b/src/preprocessing/geometries/polygon.jl @@ -8,21 +8,30 @@ struct Polygon{NDIMS, ELTYPE} min_corner :: SVector{NDIMS, ELTYPE} max_corner :: SVector{NDIMS, ELTYPE} - function Polygon(vertices) + function Polygon(vertices; close_curve=true) NDIMS = size(vertices, 1) - return Polygon{NDIMS}(vertices) + return Polygon{NDIMS}(vertices; close_curve) end # Function barrier to make `NDIMS` static and therefore `SVector`s type-stable - function Polygon{NDIMS}(vertices_) where {NDIMS} - n_vertices = size(vertices_, 2) + function Polygon{NDIMS}(vertices_; close_curve=true) where {NDIMS} ELTYPE = eltype(vertices_) - min_corner = SVector{NDIMS}(minimum(vertices_, dims=2)) - max_corner = SVector{NDIMS}(maximum(vertices_, dims=2)) + vertices = collect(reinterpret(reshape, SVector{NDIMS, ELTYPE}, vertices_)) - vertices = reinterpret(reshape, SVector{NDIMS, ELTYPE}, vertices_) + if length(vertices) < 3 + throw(ArgumentError("polygon requires at least three vertices")) + end + + if close_curve && !isapprox(first(vertices), last(vertices)) + push!(vertices, first(vertices)) + end + + n_vertices = length(vertices) + + 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]...) # Sum over all the edges and determine if the vertices are in clockwise order # to make sure that all normals pointing outwards. @@ -63,6 +72,10 @@ struct Polygon{NDIMS, ELTYPE} push!(edge_normals, edge_normal) end + if length(edge_vertices) < 3 + throw(ArgumentError("polygon requires at least three non-degenerate edges")) + end + vertex_normals = Vector{NTuple{2, SVector{NDIMS, ELTYPE}}}() # Calculate vertex pseudo-normals. diff --git a/src/preprocessing/particle_packing/signed_distance.jl b/src/preprocessing/particle_packing/signed_distance.jl index 01e862365f..79326de2a5 100644 --- a/src/preprocessing/particle_packing/signed_distance.jl +++ b/src/preprocessing/particle_packing/signed_distance.jl @@ -21,6 +21,8 @@ to this surface. - `use_for_boundary_packing`: Set to `true` if [`SignedDistanceField`] is used to pack a boundary [`ParticlePackingSystem`](@ref). Use the default of `false` when packing without a boundary. + This requires a closed geometry, since boundary packing + needs a well-defined outside region. """ struct SignedDistanceField{ELTYPE, P, N, D} positions :: P @@ -38,6 +40,11 @@ function SignedDistanceField(geometry, particle_spacing; NDIMS = ndims(geometry) ELTYPE = eltype(particle_spacing) + if use_for_boundary_packing + require_closed_geometry(geometry, + "SignedDistanceField with `use_for_boundary_packing=true`") + end + sdf_factor = use_for_boundary_packing ? 2 : 1 search_radius = sdf_factor * max_signed_distance diff --git a/src/setups/complex_shape.jl b/src/setups/complex_shape.jl index 6a78b412e2..0f44b9df8a 100644 --- a/src/setups/complex_shape.jl +++ b/src/setups/complex_shape.jl @@ -52,6 +52,8 @@ function ComplexShape(geometry; particle_spacing, density, throw(ArgumentError("`WindingNumberHormann` only supports 2D geometries")) end + require_closed_geometry(geometry, "ComplexShape") + if grid_offset < 0.0 throw(ArgumentError("only a positive `grid_offset` is supported")) end diff --git a/test/preprocessing/geometries/geometries.jl b/test/preprocessing/geometries/geometries.jl index 4f428d266c..203f6e77e6 100644 --- a/test/preprocessing/geometries/geometries.jl +++ b/test/preprocessing/geometries/geometries.jl @@ -54,12 +54,66 @@ end end + @testset verbose=true "Open Polygon Closure" begin + open_square = [1.0 2.0 2.0 1.0; + 1.0 1.0 2.0 2.0] + + geometry = TrixiParticles.Polygon(open_square) + + @test TrixiParticles.nfaces(geometry) == 4 + @test first(geometry.vertices) == last(geometry.vertices) + @test TrixiParticles.volume(geometry) ≈ 1.0 + + mktempdir() do dir + filename = joinpath(dir, "open_square.asc") + open(filename, "w") do io + println(io, "# ASCII") + for vertex in eachcol(open_square) + println(io, vertex[1], " ", vertex[2]) + end + end + + geometry_from_file = load_geometry(filename) + + @test TrixiParticles.nfaces(geometry_from_file) == 4 + @test first(geometry_from_file.vertices) == last(geometry_from_file.vertices) + @test TrixiParticles.volume(geometry_from_file) ≈ 1.0 + end + end + + @testset verbose=true "Closed Geometry Detection" begin + open_square = [1.0 2.0 2.0 1.0; + 1.0 1.0 2.0 2.0] + + closed_polygon = TrixiParticles.Polygon(open_square) + open_polygon = TrixiParticles.Polygon(open_square; close_curve=false) + partial_polygon = deleteat!(TrixiParticles.Polygon(open_square), 2) + + @test TrixiParticles.is_closed_geometry(closed_polygon) + @test !TrixiParticles.is_closed_geometry(open_polygon) + @test !TrixiParticles.is_closed_geometry(partial_polygon) + + shape = RectangularShape(0.5, (2, 2), (1.0, 1.0), density=1.0) + @test_throws ArgumentError intersect(shape, open_polygon) + @test_throws ArgumentError setdiff(shape, open_polygon) + + file = pkgdir(TrixiParticles, "test", "preprocessing", "data") + planar_geometry = load_geometry(joinpath(file, "inflow_geometry.stl")) + closed_mesh = extrude_geometry(planar_geometry, 0.8) + open_mesh = extrude_geometry(planar_geometry, 0.8; omit_top_face=true) + + @test !TrixiParticles.is_closed_geometry(planar_geometry) + @test TrixiParticles.is_closed_geometry(closed_mesh) + @test !TrixiParticles.is_closed_geometry(open_mesh) + end + @testset verbose=true "Real World Data" begin data_dir = pkgdir(TrixiParticles, "examples", "preprocessing", "data") validation_dir = pkgdir(TrixiParticles, "test", "preprocessing", "data") @testset verbose=true "2D" begin files = ["hexagon", "circle", "inverted_open_curve"] + close_curves = [true, true, false] n_edges = [6, 63, 240] volumes = [2.5980750000000006, 3.1363805763454, 2.6153740535469048] @@ -74,7 +128,8 @@ points = vcat((data.var"Points:0")', (data.var"Points:1")') - geometry = load_geometry(joinpath(data_dir, files[i] * ".asc")) + geometry = load_geometry(joinpath(data_dir, files[i] * ".asc"); + close_curve=close_curves[i]) @test TrixiParticles.nfaces(geometry) == n_edges[i] @@ -202,47 +257,22 @@ omit_bottom_face=true) winding_number_factor = 0.2 - @testset verbose=true "Omit Top Face" begin - expected_min_corner = [-0.036399998962879196; 0.24624998748302457; -0.5233639197487431;;] - expected_max_corner = [0.38360000103712083; 1.1462499874830245; -0.07336391974874301;;] - - ic_1 = ComplexShape(geometry_extruded_1; particle_spacing=0.03, density=1.0, - point_in_geometry_algorithm=WindingNumberJacobson(; - geometry=geometry_extruded_1, - winding_number_factor)) - - @test nparticles(ic_1) == 2994 - @test isapprox(maximum(ic_1.coordinates, dims=2), expected_max_corner) - @test isapprox(minimum(ic_1.coordinates, dims=2), expected_min_corner) - end - - @testset verbose=true "Omit Bottom Face" begin - expected_min_corner = [-0.0663999989628792; 0.1562499874830246; -0.49336391974874305;;] - expected_max_corner = [0.38360000103712083; 1.0562499874830245; -0.07336391974874301;;] - ic_2 = ComplexShape(geometry_extruded_2; particle_spacing=0.03, density=1.0, - point_in_geometry_algorithm=WindingNumberJacobson(; - geometry=geometry_extruded_2, - winding_number_factor)) - - @test nparticles(ic_2) == 2988 - @test isapprox(maximum(ic_2.coordinates, dims=2), expected_max_corner) - @test isapprox(minimum(ic_2.coordinates, dims=2), expected_min_corner) - end - - @testset verbose=true "Omit Both" begin - expected_min_corner = [-0.0663999989628792; 0.1562499874830246; -0.5233639197487431;;] - expected_max_corner = [0.38360000103712083; 1.1462499874830245; -0.07336391974874301;;] - - ic_3 = ComplexShape(geometry_extruded_3; particle_spacing=0.03, density=1.0, - point_in_geometry_algorithm=WindingNumberJacobson(; - geometry=geometry_extruded_3, - winding_number_factor)) - - @test nparticles(ic_3) == 3258 - @test isapprox(maximum(ic_3.coordinates, dims=2), expected_max_corner) - @test isapprox(minimum(ic_3.coordinates, dims=2), expected_min_corner) - end + @test_throws ArgumentError ComplexShape(geometry_extruded_1; + particle_spacing=0.03, density=1.0, + point_in_geometry_algorithm=WindingNumberJacobson(; + geometry=geometry_extruded_1, + winding_number_factor)) + @test_throws ArgumentError ComplexShape(geometry_extruded_2; + particle_spacing=0.03, density=1.0, + point_in_geometry_algorithm=WindingNumberJacobson(; + geometry=geometry_extruded_2, + winding_number_factor)) + @test_throws ArgumentError ComplexShape(geometry_extruded_3; + particle_spacing=0.03, density=1.0, + point_in_geometry_algorithm=WindingNumberJacobson(; + geometry=geometry_extruded_3, + winding_number_factor)) end end diff --git a/test/preprocessing/packing/signed_distance.jl b/test/preprocessing/packing/signed_distance.jl index cb02753b26..969d54292f 100644 --- a/test/preprocessing/packing/signed_distance.jl +++ b/test/preprocessing/packing/signed_distance.jl @@ -44,6 +44,16 @@ @test repr("text/plain", signed_distance_field) == show_box end + @testset verbose=true "Open Geometry Validation" begin + open_square = [0.0 1.0 1.0 0.0; + 0.0 0.0 1.0 1.0] + geometry = TrixiParticles.Polygon(open_square; close_curve=false) + + @test SignedDistanceField(geometry, 0.1) isa SignedDistanceField + @test_throws ArgumentError SignedDistanceField(geometry, 0.1; + use_for_boundary_packing=true) + 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/point_in_poly/winding_number_jacobson.jl b/test/preprocessing/point_in_poly/winding_number_jacobson.jl index ee2119a62d..a6136cff37 100644 --- a/test/preprocessing/point_in_poly/winding_number_jacobson.jl +++ b/test/preprocessing/point_in_poly/winding_number_jacobson.jl @@ -30,4 +30,17 @@ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", winding) == show_box end + + @testset verbose=true "Open Geometry Validation" begin + open_square = [0.0 1.0 1.0 0.0; + 0.0 0.0 1.0 1.0] + geometry = TrixiParticles.Polygon(open_square; close_curve=false) + points = [SVector(0.5, 0.5)] + + jacobson = WindingNumberJacobson(; hierarchical_winding=false) + hormann = WindingNumberHormann() + + @test jacobson(geometry, points)[1] isa Vector{Bool} + @test hormann(geometry, points)[1] isa Vector{Bool} + end end diff --git a/test/setups/complex_shape.jl b/test/setups/complex_shape.jl index 2917661225..1235e95889 100644 --- a/test/setups/complex_shape.jl +++ b/test/setups/complex_shape.jl @@ -41,7 +41,7 @@ end @testset verbose=true "Real World Data" begin - files = ["hexagon", "circle", "inverted_open_curve"] + files = ["hexagon", "circle"] algorithms = [ WindingNumberHormann(), WindingNumberJacobson(; hierarchical_winding=false) @@ -72,7 +72,8 @@ # See https://docs.julialang.org/en/v1/base/base/#var%22name%22 coords = vcat((data.var"Points:0")', (data.var"Points:1")') - geometry = load_geometry(joinpath(data_dir, files[j] * ".asc")) + geometry = load_geometry(joinpath(data_dir, files[j] * ".asc"); + close_curve=true) shape_sampled = ComplexShape(geometry; particle_spacing=0.05, density=1.0, point_in_geometry_algorithm) @@ -82,6 +83,15 @@ end end + @testset verbose=true "Open Geometry Validation" begin + open_square = [0.0 1.0 1.0 0.0; + 0.0 0.0 1.0 1.0] + geometry = TrixiParticles.Polygon(open_square; close_curve=false) + + @test_throws ArgumentError ComplexShape(geometry; particle_spacing=0.1, + density=1.0) + end + @testset verbose=true "Intersect of Overlapping Shapes and Geometries" begin shape = RectangularShape(0.1, (10, 10), (0.0, 0.0), density=1.0) geometry = load_geometry(joinpath(data_dir, "circle.asc"))