Close 2D geometry curves by default

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

examples/preprocessing/complex_shape_2d.jl

@@ -7,24 +7,23 @@
 # 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
 using Plots
 
 particle_spacing = 0.05
 
-filename = "inverted_open_curve"
+filename = "star"
 file = joinpath("examples", "preprocessing", "data", filename * ".asc")
 
 geometry = load_geometry(file)
 
 trixi2vtk(geometry)
 
 point_in_geometry_algorithm = WindingNumberJacobson(; geometry,
-                                                    winding_number_factor=0.4,
                                                     hierarchical_winding=true)
 
 # Returns `InitialCondition`

src/preprocessing/geometries/geometries.jl

@@ -4,13 +4,71 @@ 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)
 
     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)},

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}

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.

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

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

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