Plane interpolation

examples/fluid/rectangular_tank_3d.jl

@@ -1,9 +1,8 @@
 using TrixiParticles
-using Plots
 
 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "rectangular_tank_2d.jl"),
-              fluid_particle_spacing=0.05,
-              initial_fluid_size=(2.0, 1.0, 0.9),
-              tank_size=(2.0, 1.0, 1.2),
-              acceleration=(0.0, 0.0, -9.81),
-              smoothing_kernel=SchoenbergCubicSplineKernel{3}(), tspan=(0.0, 1.0))
+              fluid_particle_spacing=0.05, initial_fluid_size=(1.0, 1.0, 0.9),
+              tank_size=(1.0, 1.0, 1.2), acceleration=(0.0, 0.0, -9.81),
+              smoothing_kernel=SchoenbergCubicSplineKernel{3}(), tspan=(0.0, 1.0),
+              maxiters=10^5, fluid_density_calculator=ContinuityDensity(),
+              clip_negative_pressure=false)

examples/postprocessing/interpolation_plane.jl

@@ -15,16 +15,15 @@ resolution = 0.005
 # We can interpolate a plane by providing the lower left and top right coordinates of a plane.
 # Per default the same `smoothing_length` will be used as during the simulation.
 original_plane = interpolate_plane_2d(interpolation_start, interpolation_end, resolution,
-                                      semi,
-                                      fluid_system, sol)
+                                      semi, fluid_system, sol)
 original_x = [point[1] for point in original_plane.coord]
 original_y = [point[2] for point in original_plane.coord]
 original_density = original_plane.density
 
 # Plane with double smoothing length
-# Utilizing a higher `smoothing_length` in SPH interpolation enhances the level of smoothing,
-# thereby reducing the visibility of disturbances. Concurrently, it also increases the distance
-# from free surfaces where the fluid undergoes cutoff. This adjustment in `smoothing_length`
+# Utilizing a higher `smoothing_length` in SPH interpolation increases the amount of smoothing,
+# thereby reducing the visibility of disturbances. It also increases the distance
+# from free surfaces where the fluid is cut off. This adjustment in `smoothing_length`
 # can affect both the accuracy and smoothness of the interpolated results.
 double_smoothing_plane = interpolate_plane_2d(interpolation_start, interpolation_end,
                                               resolution, semi, fluid_system, sol,
@@ -34,7 +33,7 @@ double_y = [point[2] for point in double_smoothing_plane.coord]
 double_density = double_smoothing_plane.density
 
 # Plane with half smoothing length
-# Employing a lower `smoothing_length` in SPH interpolation reduces the level of smoothing,
+# Employing a lower `smoothing_length` in SPH interpolation reduces the amount of smoothing,
 # consequently increasing the visibility of disturbances. Simultaneously, it allows for a more
 # precise cutoff of the fluid at free surfaces. This change in `smoothing_length` can impact the
 # balance between the detail of disturbances captured and the precision of fluid representation near surfaces.
@@ -47,27 +46,25 @@ half_y = [point[2] for point in half_smoothing_plane.coord]
 half_density = half_smoothing_plane.density
 
 scatter1 = scatter(original_x, original_y, zcolor=original_density, marker=:circle,
-                   markersize=2, markercolor=:viridis, markerstrokewidth=0,
-                   clim=(1000, 1010), colorbar=true)
+                   markersize=2, markercolor=:viridis, markerstrokewidth=0)
 scatter2 = scatter(double_x, double_y, zcolor=double_density, marker=:circle, markersize=2,
-                   markercolor=:viridis, markerstrokewidth=0, clim=(1000, 1010),
-                   colorbar=true)
+                   markercolor=:viridis, markerstrokewidth=0)
 scatter3 = scatter(half_x, half_y, zcolor=half_density, marker=:circle, markersize=2,
-                   markercolor=:viridis, markerstrokewidth=0, clim=(1000, 1010),
-                   colorbar=true)
+                   markercolor=:viridis, markerstrokewidth=0)
 
 plot1 = plot(scatter1, xlabel="X Coordinate", ylabel="Y Coordinate",
              title="Density Distribution", colorbar_title="Density", ylim=(0.0, 1.0),
-             legend=false)
+             legend=false, clim=(1000, 1010), colorbar=true)
 plot2 = plot(scatter2, xlabel="X Coordinate", ylabel="Y Coordinate",
              title="Density with 2x Smoothing Length", colorbar_title="Density",
-             ylim=(0.0, 1.0), legend=false)
+             ylim=(0.0, 1.0), legend=false, clim=(1000, 1010), colorbar=true)
 plot3 = plot(scatter3, xlabel="X Coordinate", ylabel="Y Coordinate",
              title="Density with 0.5x Smoothing Length", colorbar_title="Density",
-             ylim=(0.0, 1.0), legend=false)
+             ylim=(0.0, 1.0), legend=false, clim=(1000, 1010), colorbar=true)
 
 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "rectangular_tank_3d.jl"),
-              tspan=(0.0, 0.1))
+              tspan=(0.0, 0.01), initial_fluid_size=(2.0, 1.0, 0.9),
+              tank_size=(2.0, 1.0, 1.2))
 
 # Interpolation parameters
 p1 = [0.0, 0.0, 0.0]
@@ -84,19 +81,33 @@ original_z = [point[3] for point in original_plane.coord]
 original_density = original_plane.density
 
 scatter_3d = scatter3d(original_x, original_y, original_z, marker_z=original_density,
-                       color=:viridis, legend=false)
+                       color=:viridis, markerstrokewidth=0)
 
 plot_3d = plot(scatter_3d, xlabel="X", ylabel="Y", zlabel="Z",
-               title="3D Scatter Plot with Density Coloring", legend=false)
+               title="3D Scatter Plot with Density Coloring", legend=false,
+               clim=(1000, 1010), colorbar=false)
+
+# Interpolation parameters
+p1 = [0.5, 0.0, 0.0]
+p2 = [0.5, 1.0, 0.0]
+p3 = [0.5, 0.0, 1.0]
+resolution = 0.01
+
+# We can also interpolate a 3D plane but in this case we have to provide 3 points instead!
+original_plane = interpolate_plane_3d(p1, p2, p3, resolution, semi,
+                                      fluid_system, sol)
+original_y = [point[2] for point in original_plane.coord]
+original_z = [point[3] for point in original_plane.coord]
+original_density = original_plane.density
 
 # by ignoring the z coordinate we can also plot this into a 2D plane
-scatter_3d_in_2d = scatter(original_x, original_y, zcolor=original_density,
+scatter_3d_in_2d = scatter(original_y, original_z, zcolor=original_density,
                            marker=:circle, markersize=4,
-                           markercolor=:viridis, markerstrokewidth=0, clim=(1000, 1010),
-                           colorbar=true)
+                           markercolor=:viridis, markerstrokewidth=0)
 
-plot_3d_in_2d = plot(scatter_3d_in_2d, xlabel="X", ylabel="Y", zlabel="Z",
-                     title="3D in 2D Scatter Plot", legend=false)
+plot_3d_in_2d = plot(scatter_3d_in_2d, xlabel="Y", ylabel="Z",
+                     title="3D in 2D Scatter Plot", legend=false, clim=(1000, 1010),
+                     colorbar=true)
 
 combined_plot = plot(plot1, plot2, plot3, plot_3d, plot_3d_in_2d, layout=(3, 2),
-                     size=(1000, 1500), margin=5mm)
+                     size=(1000, 1500), margin=3mm)

src/general/interpolation.jl

@@ -1,19 +1,19 @@
 using LinearAlgebra
 
 @doc raw"""
-    interpolate_plane_2d(lower_left, top_right, resolution, semi, ref_system, sol;
-                      smoothing_length=ref_system.smoothing_length, cut_off_bnd=true)
+    interpolate_plane_2d(min_corner, max_corner, resolution, semi, ref_system, sol;
+                         smoothing_length=ref_system.smoothing_length, cut_off_bnd=true)
 
-Interpolates properties across a plane or a volume in a TrixiParticles simulation.
+Interpolates properties along a plane in a TrixiParticles simulation.
 The region for interpolation is defined by its lower left and top right corners,
 with a specified resolution determining the density of the interpolation points.
 
 The function generates a grid of points within the defined region,
-    spaced uniformly according to the given resolution.
+spaced uniformly according to the given resolution.
 
 # Arguments
-- `lower_left`: The lower left corner of the interpolation region.
-- `top_right`: The top right corner of the interpolation region.
+- `min_corner`: The lower left corner of the interpolation region.
+- `max_corner`: The top right corner of the interpolation region.
 - `resolution`: The distance between adjacent interpolation points in the grid.
 - `semi`: The semidiscretization used for the simulation.
 - `ref_system`: The reference system for the interpolation.
@@ -44,12 +44,11 @@ function interpolate_plane_2d(min_corner, max_corner, resolution, semi, ref_syst
                               cut_off_bnd=true)
     dims = length(min_corner)
     if dims != 2 || length(max_corner) != 2
-        error("Function is intended for 2D coordinates only")
+        throw(ArgumentError("Function is intended for 2D coordinates only"))
     end
 
-    # Check that lower_left is indeed lower and to the left of top_right
     if any(min_corner .> max_corner)
-        error("lower_left should be lower and to the left of top_right in all dimensions")
+        throw(ArgumentError("`min_corner` should be smaller than `max_corner` in every dimension"))
     end
 
     # Calculate the number of points in each dimension based on the resolution
@@ -79,7 +78,7 @@ end
     interpolate_plane_3d(point1, point2, point3, resolution, semi, ref_system, sol;
                          smoothing_length=ref_system.smoothing_length, cut_off_bnd=true)
 
-Interpolates properties across a plane in a 3D space in a TrixiParticles simulation.
+Interpolates properties along a plane in a 3D space in a TrixiParticles simulation.
 The plane for interpolation is defined by three points in 3D space,
 with a specified resolution determining the density of the interpolation points.
 
@@ -113,43 +112,49 @@ The function generates a grid of points on a parallelogram within the plane defi
 # Interpolating across a plane defined by points [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], and [0.0, 1.0, 0.0]
 # with a resolution of 0.1
 results = interpolate_plane_3d([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], 0.1, semi, ref_system, sol)
+```
 """
 function interpolate_plane_3d(point1, point2, point3, resolution, semi, ref_system, sol;
                               smoothing_length=ref_system.smoothing_length,
                               cut_off_bnd=true)
     # Verify that points are in 3D space
     if length(point1) != 3 || length(point2) != 3 || length(point3) != 3
-        error("All points must be 3D coordinates")
+        throw(ArgumentError("all points must be 3D coordinates"))
     end
 
+    point1_ = SVector{3}(point1)
+    point2_ = SVector{3}(point2)
+    point3_ = SVector{3}(point3)
+
     # Vectors defining the edges of the parallelogram
-    edge1 = point2 - point1
-    edge2 = point3 - point1
+    edge1 = point2_ - point1_
+    edge2 = point3_ - point1_
 
     # Check if the points are collinear
     if norm(cross(edge1, edge2)) == 0
-        error("The points must not be collinear")
+        throw(ArgumentError("the points 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)
 
     # Create a set of points on the plane
-    points_coords = []
+    points_coords = Vector{SVector{3, Float64}}(undef,
+                                                (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
-            push!(points_coords, point_on_plane)
+            points_coords[index] = point_on_plane
+            index += 1
         end
     end
 
-    # Convert points_coords to an array of SVectors
-    svector_points = [SVector{3}(pt...) for pt in points_coords]
-
     # Interpolate using the generated points
-    results = interpolate_point(svector_points, semi, ref_system, sol,
+    results = interpolate_point(points_coords, semi, ref_system, sol,
                                 smoothing_length=smoothing_length,
                                 cut_off_bnd=cut_off_bnd)
 

test/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"

test/examples/examples.jl

@@ -18,6 +18,22 @@
             @test sol.retcode == ReturnCode.Success
         end
 
+        @trixi_testset "fluid/rectangular_tank_3d.jl" begin
+            @test_nowarn trixi_include(@__MODULE__,
+                                       joinpath(examples_dir(), "fluid",
+                                                "rectangular_tank_3d.jl"), tspan=(0.0, 0.1))
+            @test sol.retcode == ReturnCode.Success
+        end
+
+        @trixi_testset "fluid/rectangular_tank_3d.jl with SummationDensity" begin
+            @test_nowarn trixi_include(@__MODULE__,
+                                       joinpath(examples_dir(), "fluid",
+                                                "rectangular_tank_3d.jl"), tspan=(0.0, 0.1),
+                                       fluid_density_calculator=SummationDensity(),
+                                       clip_negative_pressure=true)
+            @test sol.retcode == ReturnCode.Success
+        end
+
         @trixi_testset "fluid/rectangular_tank_edac_2d.jl" begin
             @test_nowarn trixi_include(@__MODULE__,
                                        joinpath(examples_dir(), "fluid",
@@ -136,4 +152,20 @@
                                                 "n_body_benchmark_reference_faster.jl"))
         end
     end
+
+    @testset verbose=true "Postprocessing" begin
+        @trixi_testset "postprocessing/interpolation_plane.jl" begin
+            @test_nowarn trixi_include(@__MODULE__,
+                                       joinpath(examples_dir(), "postprocessing",
+                                                "interpolation_plane.jl"),
+                                       tspan=(0.0, 0.01))
+            @test sol.retcode == ReturnCode.Success
+        end
+        @trixi_testset "postprocessing/interpolation_point_line.jl" begin
+            @test_nowarn trixi_include(@__MODULE__,
+                                       joinpath(examples_dir(), "postprocessing",
+                                                "interpolation_point_line.jl"))
+            @test sol.retcode == ReturnCode.Success
+        end
+    end
 end