Improve plane interpolation performance

Project.toml

@@ -53,7 +53,7 @@ KernelAbstractions = "0.9"
 MuladdMacro = "0.2"
 OrdinaryDiffEq = "6.91"
 OrdinaryDiffEqCore = "1"
-PointNeighbors = "0.4.8"
+PointNeighbors = "0.4.11"
 Polyester = "0.7.10"
 RecipesBase = "1"
 Reexport = "1"

examples/postprocessing/interpolation_plane.jl

@@ -17,8 +17,8 @@ resolution = 0.005
 # 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)
-original_x = [point[1] for point in original_plane.coord]
-original_y = [point[2] for point in original_plane.coord]
+original_x = original_plane.point_coords[1, :]
+original_y = original_plane.point_coords[2, :]
 original_pressure = original_plane.pressure
 
 # To export the interpolated plane as a VTI file, which can be read by tools like ParaView,
@@ -35,8 +35,8 @@ interpolate_plane_2d_vtk(interpolation_start, interpolation_end, resolution,
 double_smoothing_plane = interpolate_plane_2d(interpolation_start, interpolation_end,
                                               resolution, semi, fluid_system, sol,
                                               smoothing_length=2.0 * smoothing_length)
-double_x = [point[1] for point in double_smoothing_plane.coord]
-double_y = [point[2] for point in double_smoothing_plane.coord]
+double_x = double_smoothing_plane.point_coords[1, :]
+double_y = double_smoothing_plane.point_coords[2, :]
 double_pressure = double_smoothing_plane.pressure
 
 # Plane with half smoothing length.
@@ -47,8 +47,8 @@ double_pressure = double_smoothing_plane.pressure
 half_smoothing_plane = interpolate_plane_2d(interpolation_start, interpolation_end,
                                             resolution, semi, fluid_system, sol,
                                             smoothing_length=0.5 * smoothing_length)
-half_x = [point[1] for point in half_smoothing_plane.coord]
-half_y = [point[2] for point in half_smoothing_plane.coord]
+half_x = half_smoothing_plane.point_coords[1, :]
+half_y = half_smoothing_plane.point_coords[2, :]
 half_pressure = half_smoothing_plane.pressure
 
 scatter1 = Plots.scatter(original_x, original_y, zcolor=original_pressure, marker=:circle,
@@ -82,9 +82,9 @@ resolution = 0.025
 # 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_x = [point[1] for point in original_plane.coord]
-original_y = [point[2] for point in original_plane.coord]
-original_z = [point[3] for point in original_plane.coord]
+original_x = original_plane.point_coords[1, :]
+original_y = original_plane.point_coords[2, :]
+original_z = original_plane.point_coords[3, :]
 original_pressure = original_plane.pressure
 
 scatter_3d = Plots.scatter3d(original_x, original_y, original_z, marker_z=original_pressure,

examples/postprocessing/interpolation_point_line.jl

@@ -9,26 +9,23 @@ trixi_include(@__MODULE__,
 
 position_x = tank_size[1] / 2
 
-# `interpolate_point` can be used to interpolate the properties of the `fluid_system` with the original kernel and `smoothing_length`
-println(interpolate_point([position_x, 0.01], semi, fluid_system, sol))
+# `interpolate_points` can be used to interpolate the properties of the `fluid_system`
+# with the original kernel and `smoothing_length`.
+println(interpolate_points([position_x; 0.01;;], semi, fluid_system, sol))
 # Or with an increased `smoothing_length` smoothing the result
-println(interpolate_point([position_x, 0.01], semi, fluid_system, sol,
-                          smoothing_length=2.0 * smoothing_length))
+println(interpolate_points([position_x; 0.01;;], semi, fluid_system, sol,
+                           smoothing_length=2.0 * smoothing_length))
 
-# A point outside of the domain will result in properties with value 0
-# On the boundary a result can still be obtained
-println(interpolate_point([position_x, 0.0], semi, fluid_system, sol))
+# A point outside of the domain will result in properties with NaN values.
+# On the boundary a result can still be obtained.
+println(interpolate_points([position_x; 0.0;;], semi, fluid_system, sol))
 # Slightly outside of the fluid domain the result is 0
-println(interpolate_point([position_x, -0.01], semi, fluid_system, sol))
+println(interpolate_points([position_x; -0.01;;], semi, fluid_system, sol))
 
 # Multiple points can be interpolated by providing an array
-println(interpolate_point([
-                              [position_x, 0.01],
-                              [position_x, 0.1],
-                              [position_x, 0.0],
-                              [position_x, -0.01],
-                              [position_x, -0.05]
-                          ], semi, fluid_system, sol))
+println(interpolate_points([position_x position_x position_x position_x position_x;
+                            0.01 0.1 0.0 -0.01 -0.05],
+                           semi, fluid_system, sol))
 
 # It is also possible to interpolate along a line
 n_interpolation_points = 10
@@ -40,8 +37,8 @@ result_endpoint = interpolate_line(start_point, end_point, n_interpolation_point
                                    semi, fluid_system, sol, endpoint=false)
 
 # Extracting wall distance for the standard and endpoint cases
-walldistance = [coord[2] for coord in result.coord]
-walldistance_endpoint = [coord[2] for coord in result_endpoint.coord]
+walldistance = result.point_coords[2, :]
+walldistance_endpoint = result_endpoint.point_coords[2, :]
 
 # Instead of using Plots.jl one can also use PythonPlot which uses matplotlib
 # using PythonPlot

src/TrixiParticles.jl

@@ -89,7 +89,7 @@ export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection,
 export nparticles
 export available_data, kinetic_energy, total_mass, max_pressure, min_pressure, avg_pressure,
        max_density, min_density, avg_density
-export interpolate_line, interpolate_point, interpolate_plane_3d, interpolate_plane_2d,
+export interpolate_line, interpolate_points, interpolate_plane_3d, interpolate_plane_2d,
        interpolate_plane_2d_vtk
 export SurfaceTensionAkinci, CohesionForceAkinci, SurfaceTensionMorris,
        SurfaceTensionMomentumMorris