Fix extrude_geometry (#967)

* fix

* fix tests

* fix tests

* fix test

* fix doc tests

---------

Co-authored-by: LasNikas <[email protected]>

Author: LasNikas

src/general/interpolation.jl

@@ -300,12 +300,7 @@ function interpolate_plane_3d(point1, point2, point3, resolution, semi, ref_syst
         throw(ArgumentError("`interpolate_plane_3d` requires a 3D simulation"))
     end
 
-    points_coords, resolution_ = sample_plane((point1, point2, point3), resolution)
-
-    if !isapprox(resolution, resolution_, rtol=5e-2)
-        @info "The desired plane size is not a multiple of the resolution $resolution." *
-              "\nNew resolution is set to $resolution_."
-    end
+    points_coords = sample_face_with_fixed_size((point1, point2, point3), resolution)
 
     # Interpolate using the generated points
     results = interpolate_points(points_coords, semi, ref_system, v_ode, u_ode;
@@ -766,3 +761,49 @@ function divide_by_shepard_coefficient!(field::AbstractArray{T, 3}, shepard_coef
 
     return field
 end
+
+function sample_face_with_fixed_size(face_points::NTuple{3}, resolution)
+    # Verify that points are in 3D space
+    if any(length.(face_points) .!= 3)
+        throw(ArgumentError("all points must be 3D coordinates"))
+    end
+
+    point1_ = SVector{3}(face_points[1])
+    point2_ = SVector{3}(face_points[2])
+    point3_ = SVector{3}(face_points[3])
+
+    # Vectors defining the edges of the parallelogram
+    edge1 = point2_ - point1_
+    edge2 = point3_ - point1_
+
+    # Check if the points are collinear
+    if isapprox(norm(cross(edge1, edge2)), 0; atol=eps())
+        throw(ArgumentError("the vectors `AB` and `AC` must not be collinear"))
+    end
+
+    # Determine the number of points along each edge
+    num_points_edge1 = ceil(Int, norm(edge1) / resolution)
+    num_points_edge2 = ceil(Int, norm(edge2) / resolution)
+
+    coords = zeros(3, (num_points_edge1 + 1) * (num_points_edge2 + 1))
+
+    index = 1
+    for i in 0:num_points_edge1
+        for j in 0:num_points_edge2
+            point_on_plane = point1_ + (i / num_points_edge1) * edge1 +
+                             (j / num_points_edge2) * edge2
+            coords[:, index] = point_on_plane
+            index += 1
+        end
+    end
+
+    resolution_new = min(norm(edge1 / num_points_edge1),
+                         norm(edge2 / num_points_edge2))
+
+    if !isapprox(resolution, resolution_new, rtol=sqrt(eps(resolution)))
+        @info "The desired face size is not a multiple of the resolution $resolution." *
+              "\nNew resolution is set to $resolution_new."
+    end
+
+    return coords
+end

src/setups/extrude_geometry.jl

@@ -64,10 +64,12 @@ direction = [0.0, 0.0, 1.0]
 shape = extrude_geometry(shape; direction, particle_spacing=0.1, n_extrude=4, density=1000.0)
 
 # output
-┌ Info: The desired size is not a multiple of the particle spacing 0.1.
-└ New particle spacing is set to 0.09387239731236392.
-┌ Info: The desired size is not a multiple of the particle spacing 0.1.
-└ New particle spacing is set to 0.09198039027185569.
+┌ Info: The desired line length 1.314213562373095 is not a multiple of the particle spacing 0.1.
+└ New line length is set to 1.3.
+┌ Info: The desired edge 1 length 1.0180339887498948 is not a multiple of the particle spacing 0.1.
+└ New edge 1 length is set to 1.0.
+┌ Info: The desired edge 2 length 0.9198039027185568 is not a multiple of the particle spacing 0.1.
+└ New edge 2 length is set to 0.9.
 ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
 │ InitialCondition{Float64}                                                                        │
 │ ═════════════════════════                                                                        │
@@ -99,16 +101,9 @@ function extrude_geometry(geometry; particle_spacing=-1, direction, n_extrude::I
 
     geometry = shift_plane_corners(geometry, direction_, particle_spacing, place_on_shell)
 
-    face_coords,
-    particle_spacing_ = sample_plane(geometry, particle_spacing;
-                                     place_on_shell=place_on_shell)
+    face_coords = sample_plane(geometry, particle_spacing; place_on_shell=place_on_shell)
 
-    if !isapprox(particle_spacing, particle_spacing_, rtol=5e-2)
-        @info "The desired size is not a multiple of the particle spacing $particle_spacing." *
-              "\nNew particle spacing is set to $particle_spacing_."
-    end
-
-    coords = (face_coords .+ i * particle_spacing_ * direction_ for i in 0:(n_extrude - 1))
+    coords = (face_coords .+ i * particle_spacing * direction_ for i in 0:(n_extrude - 1))
 
     # In this context, `stack` is faster than `hcat(coords...)`
     coordinates = reshape(stack(coords), (NDIMS, size(face_coords, 2) * n_extrude))
@@ -118,7 +113,7 @@ function extrude_geometry(geometry; particle_spacing=-1, direction, n_extrude::I
     end
 
     return InitialCondition(; coordinates, velocity, density, mass, pressure,
-                            particle_spacing=particle_spacing_)
+                            particle_spacing=particle_spacing)
 end
 
 # For corners/endpoints of a plane/line, sample the plane/line with particles.
@@ -132,10 +127,10 @@ function sample_plane(geometry::AbstractMatrix, particle_spacing; place_on_shell
                       fill(place_on_shell ? 0 : particle_spacing / 2, size(geometry, 2))')
 
         # TODO: 2D shapes not only in x-y plane but in any user-defined plane
-        return coords, particle_spacing
+        return coords
     end
 
-    return geometry, particle_spacing
+    return geometry
 end
 
 function sample_plane(shape::InitialCondition, particle_spacing; place_on_shell)
@@ -148,10 +143,10 @@ function sample_plane(shape::InitialCondition, particle_spacing; place_on_shell)
                            size(shape.coordinates, 2))')
 
         # TODO: 2D shapes not only in x-y plane but in any user-defined plane
-        return coords, particle_spacing
+        return coords
     end
 
-    return shape.coordinates, particle_spacing
+    return shape.coordinates
 end
 
 function sample_plane(plane_points, particle_spacing; place_on_shell=nothing)
@@ -166,12 +161,19 @@ function sample_plane(plane_points::NTuple{2}, particle_spacing; place_on_shell=
         throw(ArgumentError("all points must be 2D coordinates"))
     end
 
-    n_points = ceil(Int, norm(plane_points[2] - plane_points[1]) / particle_spacing) + 1
+    p1 = plane_points[1]
+    p2 = plane_points[2]
+
+    line_dir = p2 - p1
+    line_length = norm(line_dir)
 
-    coords = stack(range(plane_points[1], plane_points[2], length=n_points))
-    particle_spacing_new = norm(coords[:, 1] - coords[:, 2])
+    n_particles,
+    new_length = round_n_particles(line_length, particle_spacing,
+                                   "line length")
 
-    return coords, particle_spacing_new
+    coords = stack([p1 + i * particle_spacing * normalize(line_dir) for i in 0:n_particles])
+
+    return coords
 end
 
 function sample_plane(plane_points::NTuple{3}, particle_spacing; place_on_shell=nothing)
@@ -194,25 +196,28 @@ function sample_plane(plane_points::NTuple{3}, particle_spacing; place_on_shell=
     end
 
     # Determine the number of points along each edge
-    num_points_edge1 = ceil(Int, norm(edge1) / particle_spacing)
-    num_points_edge2 = ceil(Int, norm(edge2) / particle_spacing)
-
+    num_points_edge1,
+    new_length = round_n_particles(norm(edge1), particle_spacing,
+                                   "edge 1 length")
+    num_points_edge2,
+    new_length = round_n_particles(norm(edge2), particle_spacing,
+                                   "edge 2 length")
+
+    dir1 = normalize(edge1)
+    dir2 = normalize(edge2)
     coords = zeros(3, (num_points_edge1 + 1) * (num_points_edge2 + 1))
 
     index = 1
     for i in 0:num_points_edge1
         for j in 0:num_points_edge2
-            point_on_plane = point1_ + (i / num_points_edge1) * edge1 +
-                             (j / num_points_edge2) * edge2
+            point_on_plane = point1_ + i * particle_spacing * dir1 +
+                             j * particle_spacing * dir2
             coords[:, index] = point_on_plane
             index += 1
         end
     end
 
-    particle_spacing_new = min(norm(edge1 / num_points_edge1),
-                               norm(edge2 / num_points_edge2))
-
-    return coords, particle_spacing_new
+    return coords
 end
 
 # Shift corners of the plane/line inwards by half a particle spacing with `place_on_shell=false`

test/examples/examples.jl

@@ -184,7 +184,8 @@
                 r"WARNING: using deprecated binding PlotUtils.*\n",
                 r"WARNING: Makie.* is deprecated.*\n",
                 r"  likely near none:1\n",
-                r", use .* instead.\n"
+                r", use .* instead.\n",
+                r"┌ Info: The desired face size is not a multiple of the resolution [0-9.]+\.\n└ New resolution is set to [0-9.]+\.\n"
             ]
             @test sol.retcode == ReturnCode.Success
         end

test/schemes/boundary/open_boundary/characteristic_variables.jl

@@ -2,7 +2,7 @@
     particle_spacing = 0.1
 
     # Number of boundary particles in the influence of fluid particles
-    influenced_particles = [20, 52, 26]
+    influenced_particles = [20, 50, 26]
 
     open_boundary_layers = 8
     sound_speed = 20.0