Rework particle packing

examples/preprocessing/complex_shape_2d.jl

@@ -21,9 +21,10 @@ shape_sampled = ComplexShape(geometry; particle_spacing, density=1.0,
 
 trixi2vtk(shape_sampled.initial_condition)
 
+coordinates = stack(shape_sampled.grid)
 # trixi2vtk(shape_sampled.signed_distance_field)
-# trixi2vtk(shape_sampled.grid, w=shape_sampled.winding_numbers)
+# trixi2vtk(coordinates, w=shape_sampled.winding_numbers)
 
 # Plot the winding number field
-plot(InitialCondition(; coordinates=shape_sampled.grid, density=1.0, particle_spacing),
+plot(InitialCondition(; coordinates, density=1.0, particle_spacing),
      zcolor=shape_sampled.winding_numbers)

examples/preprocessing/packing_2d.jl

@@ -1,25 +1,25 @@
 using TrixiParticles
-using OrdinaryDiffEq
+using OrdinaryDiffEq, Plots
 
 filename = "circle"
 file = pkgdir(TrixiParticles, "examples", "preprocessing", "data", filename * ".asc")
 
 # ==========================================================================================
 # ==== Packing parameters
-tlsph = true
+tlsph = false
 
 # ==========================================================================================
 # ==== Resolution
-particle_spacing = 0.03
+particle_spacing = 0.1
 
 # The following depends on the sampling of the particles. In this case `boundary_thickness`
 # means literally the thickness of the boundary packed with boundary particles and *not*
 # how many rows of boundary particles will be sampled.
-boundary_thickness = 8 * particle_spacing
+boundary_thickness = 5 * particle_spacing
 
 # ==========================================================================================
 # ==== Load complex geometry
-density = 1000.0
+density = 1.0
 
 geometry = load_geometry(file)
 
@@ -44,16 +44,18 @@ trixi2vtk(boundary_sampled, filename="boundary")
 # ==========================================================================================
 # ==== Packing
 
-# Large `background_pressure` can cause high accelerations. That is, the adaptive
-# time-stepsize will be adjusted properly. We found that the following order of
-# `background_pressure` result in appropriate stepsizes.
-background_pressure = 1e6 * particle_spacing^ndims(geometry)
+# A larger `background_pressure` makes the packing happen faster in physical time,
+# which results in a correspondingly smaller time step.
+# Essentially, the `background_pressure` just scales the physical time,
+# and can therefore arbitrarily be set to 1.
+background_pressure = 1.0
 
-packing_system = ParticlePackingSystem(shape_sampled;
+smoothing_length = 0.8 * particle_spacing
+packing_system = ParticlePackingSystem(shape_sampled; smoothing_length=smoothing_length,
                                        signed_distance_field, tlsph=tlsph,
                                        background_pressure)
 
-boundary_system = ParticlePackingSystem(boundary_sampled;
+boundary_system = ParticlePackingSystem(boundary_sampled; smoothing_length=smoothing_length,
                                         is_boundary=true, signed_distance_field,
                                         tlsph=tlsph, boundary_compress_factor=0.8,
                                         background_pressure)
@@ -67,8 +69,8 @@ tspan = (0, 10.0)
 ode = semidiscretize(semi, tspan)
 
 # Use this callback to stop the simulation when it is sufficiently close to a steady state
-steady_state = SteadyStateReachedCallback(; interval=1, interval_size=10,
-                                          abstol=1.0e-5, reltol=1.0e-3)
+steady_state = SteadyStateReachedCallback(; interval=10, interval_size=200,
+                                          abstol=1.0e-7, reltol=1.0e-6)
 
 info_callback = InfoCallback(interval=50)
 
@@ -77,7 +79,11 @@ saving_callback = save_intervals ?
                   SolutionSavingCallback(interval=10, prefix="", ekin=kinetic_energy) :
                   nothing
 
-callbacks = CallbackSet(UpdateCallback(), saving_callback, info_callback, steady_state)
+pp_cb_ekin = PostprocessCallback(; ekin=kinetic_energy, interval=1,
+                                 filename="kinetic_energy", write_file_interval=50)
+
+callbacks = CallbackSet(UpdateCallback(), saving_callback, info_callback, steady_state,
+                        pp_cb_ekin)
 
 sol = solve(ode, RDPK3SpFSAL35();
             save_everystep=false, maxiters=1000, callback=callbacks, dtmax=1e-2)
@@ -87,3 +93,9 @@ packed_boundary_ic = InitialCondition(sol, boundary_system, semi)
 
 trixi2vtk(packed_ic, filename="initial_condition_packed")
 trixi2vtk(packed_boundary_ic, filename="initial_condition_boundary_packed")
+
+shape = Plots.Shape(stack(geometry.vertices)[1, :], stack(geometry.vertices)[2, :])
+p1 = plot(shape_sampled, markerstrokewidth=1, label=nothing, layout=(1, 2))
+plot!(p1, shape, color=nothing, label=nothing, linewidth=2, subplot=1)
+plot!(p1, packed_ic, markerstrokewidth=1, label=nothing, subplot=2)
+plot!(p1, shape, color=nothing, label=nothing, linewidth=2, subplot=2, size=(850, 400))

src/general/semidiscretization.jl

@@ -699,14 +699,14 @@ end
 
 function update_nhs!(neighborhood_search,
                      system::OpenBoundarySPHSystem, neighbor::TotalLagrangianSPHSystem,
-                     u_system, u_neighbor)
+                     u_system, u_neighbor, semi)
     # Don't update. This NHS is never used.
     return neighborhood_search
 end
 
 function update_nhs!(neighborhood_search,
                      system::TotalLagrangianSPHSystem, neighbor::OpenBoundarySPHSystem,
-                     u_system, u_neighbor)
+                     u_system, u_neighbor, semi)
     # Don't update. This NHS is never used.
     return neighborhood_search
 end

src/preprocessing/geometries/polygon.jl

@@ -139,3 +139,21 @@ end
 end
 
 @inline face_normal(edge, geometry::Polygon) = geometry.edge_normals[edge]
+
+# Although "volume" is typically associated with 3D objects, the function name
+# is intentionally kept as `volume` to maintain consistency between 2D and 3D cases.
+# This allows seamless handling of both cases without requiring differentiation
+# between dimensions in the calling code.
+# Note that this function is not part of the public API.
+function volume(polygon::Polygon)
+
+    # Compute the area of the polygon by the shoelace formula.
+    # https://en.wikipedia.org/wiki/Polygon
+    volume = sum(polygon.edge_vertices_ids, init=zero(eltype(polygon))) do edge
+        v1 = polygon.vertices[edge[1]]
+        v2 = polygon.vertices[edge[2]]
+
+        return (v1[1] * v2[2] - v2[1] * v1[2])
+    end
+    return abs(volume) / 2
+end

src/preprocessing/geometries/triangle_mesh.jl

@@ -129,9 +129,19 @@ struct TriangleMesh{NDIMS, ELTYPE}
         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]...)
 
+        for i in eachindex(edge_normals)
+            # Skip zero normals, which would be normalized to `NaN` vectors.
+            # The edge normals are only used for the `SignedDistanceField`, which is
+            # essential for the packing.
+            # Zero normals are caused by exactly or nearly duplicated faces.
+            if !iszero(norm(edge_normals[i]))
+                edge_normals[i] = normalize(edge_normals[i])
+            end
+        end
+
         return new{NDIMS, ELTYPE}(vertices, face_vertices, face_vertices_ids,
                                   face_edges_ids, edge_vertices_ids,
-                                  normalize.(vertex_normals), normalize.(edge_normals),
+                                  normalize.(vertex_normals), edge_normals,
                                   face_normals, min_corner, max_corner)
     end
 end
@@ -232,3 +242,15 @@ function unique_sorted(vertices)
 
     return vertices_sorted[keep]
 end
+
+function volume(mesh::TriangleMesh)
+    volume = sum(mesh.face_vertices) do vertices
+
+        # Formula for the volume of a tetrahedron:
+        # V = (1/6) * |a · (b × c)|, where a, b, and c are vectors defining the tetrahedron.
+        # Reference: https://en.wikipedia.org/wiki/Tetrahedron#Volume
+        return dot(vertices[1], cross(vertices[2], vertices[3])) / 6
+    end
+
+    return volume
+end

src/preprocessing/particle_packing/nhs_faces.jl

@@ -24,6 +24,10 @@ function faces_in_cell(cell, neighborhood_search)
     return PointNeighbors.points_in_cell(cell, neighborhood_search)
 end
 
+@inline function eachneighbor(coords, nhs::TrivialNeighborhoodSearch)
+    return PointNeighbors.eachneighbor(coords, nhs)
+end
+
 @inline function eachneighbor(coords, neighborhood_search::FaceNeighborhoodSearch)
     cell = PointNeighbors.cell_coords(coords, neighborhood_search)
     return neighborhood_search.neighbors[cell]

src/preprocessing/particle_packing/rhs.jl

@@ -1,7 +1,7 @@
 function interact!(dv, v_particle_system, u_particle_system,
                    v_neighbor_system, u_neighbor_system,
-                   system::ParticlePackingSystem, neighbor_system::ParticlePackingSystem,
-                   semi)
+                   system::ParticlePackingSystem{<:Any, false},
+                   neighbor_system::ParticlePackingSystem, semi)
     system_coords = current_coordinates(u_particle_system, system)
     neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)
 
@@ -34,3 +34,10 @@ function interact!(dv, v_particle_system, u_particle_system,
 
     return dv
 end
+
+# Skip for fixed systems
+function interact!(dv, v_particle_system, u_particle_system,
+                   v_neighbor_system, u_neighbor_system,
+                   system::ParticlePackingSystem{<:Any, true}, neighbor_system, semi)
+    return dv
+end

src/preprocessing/particle_packing/signed_distance.jl

@@ -22,65 +22,67 @@ to this surface.
                               a boundary [`ParticlePackingSystem`](@ref).
                               Use the default of `false` when packing without a boundary.
 """
-struct SignedDistanceField{NDIMS, ELTYPE}
-    positions           :: Vector{SVector{NDIMS, ELTYPE}}
-    normals             :: Vector{SVector{NDIMS, ELTYPE}}
-    distances           :: Vector{ELTYPE}
+struct SignedDistanceField{ELTYPE, P, D}
+    positions           :: P
+    normals             :: P
+    distances           :: D
     max_signed_distance :: ELTYPE
     boundary_packing    :: Bool
     particle_spacing    :: ELTYPE
+end
 
-    function SignedDistanceField(geometry, particle_spacing;
-                                 points=nothing,
-                                 max_signed_distance=4 * particle_spacing,
-                                 use_for_boundary_packing=false)
-        NDIMS = ndims(geometry)
-        ELTYPE = eltype(max_signed_distance)
+function SignedDistanceField(geometry, particle_spacing;
+                             points=nothing, neighborhood_search=true,
+                             max_signed_distance=4 * particle_spacing,
+                             use_for_boundary_packing=false)
+    NDIMS = ndims(geometry)
+    ELTYPE = eltype(max_signed_distance)
 
-        sdf_factor = use_for_boundary_packing ? 2 : 1
+    sdf_factor = use_for_boundary_packing ? 2 : 1
 
-        search_radius = sdf_factor * max_signed_distance
+    search_radius = sdf_factor * max_signed_distance
 
+    if neighborhood_search
         nhs = FaceNeighborhoodSearch{NDIMS}(; search_radius)
+    else
+        nhs = TrivialNeighborhoodSearch{NDIMS}(eachpoint=eachface(geometry))
+    end
 
-        initialize!(nhs, geometry)
+    initialize!(nhs, geometry)
 
-        if isnothing(points)
-            min_corner = geometry.min_corner .- search_radius
-            max_corner = geometry.max_corner .+ search_radius
+    if isnothing(points)
+        min_corner = geometry.min_corner .- search_radius
+        max_corner = geometry.max_corner .+ search_radius
 
-            n_particles_per_dimension = Tuple(ceil.(Int,
-                                                    (max_corner .- min_corner) ./
-                                                    particle_spacing))
+        n_particles_per_dimension = Tuple(ceil.(Int,
+                                                (max_corner .- min_corner) ./
+                                                particle_spacing))
 
-            grid = rectangular_shape_coords(particle_spacing, n_particles_per_dimension,
-                                            min_corner; tlsph=true)
+        grid = rectangular_shape_coords(particle_spacing, n_particles_per_dimension,
+                                        min_corner; tlsph=true)
 
-            points = reinterpret(reshape, SVector{NDIMS, ELTYPE}, grid)
-        end
+        points = reinterpret(reshape, SVector{NDIMS, ELTYPE}, grid)
+    end
 
-        positions = copy(points)
+    positions = copy(points)
 
-        # This gives a performance boost for large geometries
-        delete_positions_in_empty_cells!(positions, nhs)
+    # This gives a performance boost for large geometries
+    delete_positions_in_empty_cells!(positions, nhs)
 
-        normals = fill(SVector(ntuple(dim -> Inf, NDIMS)), length(positions))
-        distances = fill(Inf, length(positions))
+    normals = fill(SVector(ntuple(dim -> Inf, NDIMS)), length(positions))
+    distances = fill(Inf, length(positions))
 
-        calculate_signed_distances!(positions, distances, normals,
-                                    geometry, sdf_factor, max_signed_distance, nhs)
+    calculate_signed_distances!(positions, distances, normals,
+                                geometry, sdf_factor, max_signed_distance, nhs)
 
-        return new{NDIMS, ELTYPE}(positions, normals, distances, max_signed_distance,
-                                  use_for_boundary_packing, particle_spacing)
-    end
+    return SignedDistanceField(positions, normals, distances, max_signed_distance,
+                               use_for_boundary_packing, particle_spacing)
 end
 
-@inline Base.ndims(::SignedDistanceField{NDIMS}) where {NDIMS} = NDIMS
-
 function Base.show(io::IO, system::SignedDistanceField)
     @nospecialize system # reduce precompilation time
 
-    print(io, "SignedDistanceField{", ndims(system), "}()")
+    print(io, "SignedDistanceField()")
 end
 
 function Base.show(io::IO, ::MIME"text/plain", system::SignedDistanceField)
@@ -89,7 +91,7 @@ function Base.show(io::IO, ::MIME"text/plain", system::SignedDistanceField)
     if get(io, :compact, false)
         show(io, system)
     else
-        summary_header(io, "SignedDistanceField{$(ndims(system))}")
+        summary_header(io, "SignedDistanceField")
         summary_line(io, "#particles", length(system.distances))
         summary_line(io, "max signed distance", system.max_signed_distance)
         summary_footer(io)
@@ -105,6 +107,8 @@ function trixi2vtk(signed_distance_field::SignedDistanceField;
               filename=filename, output_directory=output_directory)
 end
 
+delete_positions_in_empty_cells!(positions, nhs::TrivialNeighborhoodSearch) = positions
+
 function delete_positions_in_empty_cells!(positions, nhs::FaceNeighborhoodSearch)
     delete_positions = fill(false, length(positions))