Skip to content

Commit 17fb2fe

Browse files
svchbsvchb
committed
more files
1 parent bf7c911 commit 17fb2fe

File tree

4 files changed

+460
-146
lines changed

4 files changed

+460
-146
lines changed

‎examples/preprocessing/complex_shape_2d.jl‎

Lines changed: 99 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -1,30 +1,109 @@
1+
# ==========================================================================================
2+
# 2D Complex Shape Sampling and Winding Number Visualization
3+
#
4+
# This example demonstrates how to:
5+
# 1. Load a 2D geometry from an ASCII file (e.g., a curve).
6+
# 2. Sample particles within this complex geometry using the `ComplexShape` functionality.
7+
# 3. Utilize the Winding Number algorithm to determine if points are inside or outside.
8+
# 4. Visualize the sampled particles and the winding number field.
9+
#
10+
# The example uses an "inverted_open_curve" geometry, where standard inside/outside
11+
# definitions might be ambiguous without a robust point-in-polygon test like winding numbers.
12+
# ==========================================================================================
13+
114
using TrixiParticles
2-
using Plots
15+
using Plots # For visualizing the winding number field
16+
using Plots.PlotMeasures # For plot margins, if needed
17+
18+
# ------------------------------------------------------------------------------
19+
# Parameters
20+
# ------------------------------------------------------------------------------
21+
# Particle spacing for sampling the complex shape
22+
particle_spacing = 0.05 # meters
23+
24+
# Geometry file details
25+
geometry_filename_stem = "inverted_open_curve"
26+
# Assuming the data file is in "examples/preprocessing/data/" relative to TrixiParticles.jl root
27+
geometry_file_path = joinpath(pkgdir(TrixiParticles), "examples", "preprocessing", "data",
28+
geometry_filename_stem * ".asc")
29+
30+
# Parameters for the Winding Number algorithm
31+
# `winding_number_factor` helps classify points near the boundary.
32+
# `hierarchical_winding` uses a multi-level approach for robustness and efficiency.
33+
winding_number_threshold_factor = 0.4 # Points with winding number > factor are considered inside
34+
use_hierarchical_winding = true
35+
36+
# ------------------------------------------------------------------------------
37+
# Load Geometry and Sample Particles
38+
# ------------------------------------------------------------------------------
39+
# Load the 2D geometry from the specified file.
40+
# `geometry` will be a `BoundaryPath` object or similar.
41+
println("Loading 2D geometry from: $geometry_file_path")
42+
geometry = load_geometry(geometry_file_path)
43+
44+
# Optional: Export the loaded geometry to a VTK file for inspection.
45+
# This helps verify that the geometry was loaded correctly.
46+
output_vtk_geometry_filename = "out/$(geometry_filename_stem)_boundary.vtk"
47+
trixi2vtk(geometry, filename=output_vtk_geometry_filename)
48+
println("Exported loaded geometry to $output_vtk_geometry_filename")
49+
50+
# Define the point-in-geometry algorithm using Winding Numbers.
51+
# This algorithm determines which grid points (for sampling) fall inside the geometry.
52+
point_in_geometry_test = WindingNumberJacobson(geometry=geometry,
53+
winding_number_factor=winding_number_threshold_factor,
54+
hierarchical_winding=use_hierarchical_winding)
355

4-
particle_spacing = 0.05
56+
# Sample particles within the complex shape.
57+
# `ComplexShape` generates an `InitialCondition` object containing particle positions and properties.
58+
# `store_winding_number=true` stores the calculated winding number for each sampled point.
59+
println("Sampling particles within the complex shape...")
60+
sampled_shape_data = ComplexShape(geometry;
61+
particle_spacing=particle_spacing,
62+
density=1.0, # Dummy density for visualization
63+
store_winding_number=true,
64+
point_in_geometry_algorithm=point_in_geometry_test)
565

6-
filename = "inverted_open_curve"
7-
file = joinpath("examples", "preprocessing", "data", filename * ".asc")
66+
# The actual particle data is in `sampled_shape_data.initial_condition`.
67+
initial_condition_particles = sampled_shape_data.initial_condition
868

9-
geometry = load_geometry(file)
69+
# Export the sampled particles (initial condition) to a VTK file.
70+
output_vtk_particles_filename = "out/$(geometry_filename_stem)_sampled_particles.vtp"
71+
trixi2vtk(initial_condition_particles, filename=output_vtk_particles_filename)
72+
println("Exported sampled particles to $output_vtk_particles_filename")
1073

11-
trixi2vtk(geometry)
74+
# ------------------------------------------------------------------------------
75+
# Visualize Winding Numbers (Optional)
76+
# ------------------------------------------------------------------------------
77+
# `sampled_shape_data.grid` contains all grid points considered during sampling.
78+
# `sampled_shape_data.winding_numbers` contains the winding number for each of these grid points.
1279

13-
point_in_geometry_algorithm = WindingNumberJacobson(; geometry,
14-
winding_number_factor=0.4,
15-
hierarchical_winding=true)
80+
# Extract coordinates of all grid points and their winding numbers.
81+
grid_point_coordinates = stack(sampled_shape_data.grid) # Converts vector of SVectors to matrix
82+
grid_winding_numbers = sampled_shape_data.winding_numbers
1683

17-
# Returns `InitialCondition`
18-
shape_sampled = ComplexShape(geometry; particle_spacing, density=1.0,
19-
store_winding_number=true,
20-
point_in_geometry_algorithm)
84+
# Create an `InitialCondition` object for plotting all grid points colored by winding number.
85+
# This helps visualize the winding number field.
86+
grid_visualization_ic = InitialCondition(coordinates=grid_point_coordinates,
87+
density=1.0, # Dummy density
88+
particle_spacing=particle_spacing) # For marker size in plot
2189

22-
trixi2vtk(shape_sampled.initial_condition)
90+
println("Plotting the winding number field for all considered grid points...")
91+
winding_number_plot = plot(grid_visualization_ic,
92+
zcolor=grid_winding_numbers, # Color by winding number
93+
markersize=2, markerstrokewidth=0,
94+
colorbar_title="Winding Number",
95+
title="Winding Number Field of '$geometry_filename_stem'",
96+
aspect_ratio=:equal)
97+
# Overlay the original geometry boundary for context
98+
plot!(winding_number_plot, geometry, linecolor=:black, linewidth=1.5, label="Geometry Boundary")
2399

24-
coordinates = stack(shape_sampled.grid)
25-
# trixi2vtk(shape_sampled.signed_distance_field)
26-
# trixi2vtk(coordinates, w=shape_sampled.winding_numbers)
100+
display(winding_number_plot)
101+
println("Complex shape 2D example finished. Check 'out/' directory for VTK files and plot window.")
27102

28-
# Plot the winding number field
29-
plot(InitialCondition(; coordinates, density=1.0, particle_spacing),
30-
zcolor=shape_sampled.winding_numbers)
103+
# Additional advanced visualization (commented out, requires specific understanding):
104+
# - `sampled_shape_data.signed_distance_field`: If `create_signed_distance_field=true` was used.
105+
# trixi2vtk(sampled_shape_data.signed_distance_field)
106+
# - Winding numbers for only the *sampled* (inside) particles:
107+
# To get this, one would filter `sampled_shape_data.grid` based on where particles were placed
108+
# in `initial_condition_particles` and then map the winding numbers.
109+
# The current `sampled_shape_data.winding_numbers` is for all grid points.

‎examples/preprocessing/complex_shape_3d.jl‎

Lines changed: 91 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -1,21 +1,98 @@
1+
# ==========================================================================================
2+
# 3D Complex Shape Sampling (e.g., from STL)
3+
#
4+
# This example demonstrates how to:
5+
# 1. Load a 3D geometry from an STL file (e.g., a sphere).
6+
# 2. Sample particles either as a fluid volume within this geometry or as a boundary layer.
7+
# 3. Optionally create a Signed Distance Field (SDF) from the geometry.
8+
# 4. Export the results to VTK files for visualization.
9+
#
10+
# The Winding Number algorithm is typically used for robust point-in-volume tests for 3D.
11+
# ==========================================================================================
12+
113
using TrixiParticles
214

3-
particle_spacing = 0.05
15+
# ------------------------------------------------------------------------------
16+
# Parameters
17+
# ------------------------------------------------------------------------------
18+
# Particle spacing for sampling
19+
particle_spacing = 0.05 # meters
20+
21+
# Geometry file details
22+
geometry_filename_stem = "sphere" # Example: a sphere STL
23+
geometry_file_path = joinpath(pkgdir(TrixiParticles), "examples", "preprocessing", "data",
24+
geometry_filename_stem * ".stl")
25+
26+
# Parameters for the Winding Number algorithm (if used for fluid sampling)
27+
# For 3D, the winding number concept is more complex. A common threshold is 0.5.
28+
# `hierarchical_winding` is generally recommended for 3D.
29+
winding_number_threshold_factor_3d = 0.5 # Typical threshold for 3D closed surfaces
30+
use_hierarchical_winding_3d = true
31+
32+
# Sampling options for ComplexShape
33+
sample_as_fluid_volume = true # If true, sample particles filling the geometry interior.
34+
sample_as_boundary_layer = false # If true, sample particles on/near the geometry surface.
35+
# Note: `sample_boundary` in `ComplexShape` refers to sampling boundary *particles*,
36+
# not just evaluating the geometry boundary.
37+
# If `sample_as_fluid_volume` is true, `sample_boundary=false` is common.
38+
# If only a boundary layer is needed, `sample_as_fluid_volume=false` and `sample_boundary=true` might be used.
39+
40+
# Thickness of the boundary layer if `sample_as_boundary_layer` or `sample_boundary=true` is active.
41+
boundary_sampling_thickness = 5 * particle_spacing
42+
43+
# Option to create a Signed Distance Field (SDF) from the geometry.
44+
generate_signed_distance_field = true
45+
46+
# ------------------------------------------------------------------------------
47+
# Load Geometry and Sample Particles
48+
# ------------------------------------------------------------------------------
49+
println("Loading 3D geometry from: $geometry_file_path")
50+
geometry = load_geometry(geometry_file_path) # `geometry` will be a `TriangulatedSurface` or similar.
51+
52+
# Optional: Export the loaded geometry for verification.
53+
output_vtk_geometry_filename_3d = "out/$(geometry_filename_stem)_3d_boundary.vtp"
54+
trixi2vtk(geometry, filename=output_vtk_geometry_filename_3d)
55+
println("Exported loaded 3D geometry to $output_vtk_geometry_filename_3d")
56+
57+
# Define the point-in-geometry algorithm (used if `sample_as_fluid_volume=true`).
58+
point_in_geometry_test_3d = WindingNumberJacobson(geometry=geometry,
59+
winding_number_factor=winding_number_threshold_factor_3d,
60+
hierarchical_winding=use_hierarchical_winding_3d)
461

5-
filename = "sphere"
6-
file = joinpath("examples", "preprocessing", "data", filename * ".stl")
62+
# Sample particles using ComplexShape.
63+
# The behavior depends on `sample_boundary` and how the algorithm is used.
64+
# Original example implies fluid sampling (`sample_boundary=false`).
65+
println("Sampling particles for the 3D complex shape...")
66+
sampled_shape_data_3d = ComplexShape(geometry;
67+
particle_spacing=particle_spacing,
68+
density=1.0, # Dummy density
69+
# `boundary_thickness` is used if `sample_boundary=true`
70+
boundary_thickness=boundary_sampling_thickness,
71+
# `sample_boundary=false` means we primarily want interior/fluid particles
72+
sample_boundary=sample_as_boundary_layer,
73+
# Algorithm to determine if grid points are inside (for fluid)
74+
point_in_geometry_algorithm=point_in_geometry_test_3d,
75+
create_signed_distance_field=generate_signed_distance_field)
776

8-
geometry = load_geometry(file)
77+
# Export the sampled fluid particles (initial condition).
78+
# These are the particles inside the geometry if `sample_boundary=false` was effectively used
79+
# for fluid volume sampling.
80+
output_vtk_fluid_particles_3d = "out/$(geometry_filename_stem)_3d_fluid_particles.vtp"
81+
trixi2vtk(sampled_shape_data_3d.initial_condition, filename=output_vtk_fluid_particles_3d)
82+
println("Exported sampled 3D fluid particles to $output_vtk_fluid_particles_3d")
983

10-
point_in_geometry_algorithm = WindingNumberJacobson(; geometry,
11-
# winding_number_factor=0.4,
12-
hierarchical_winding=true)
84+
# If a boundary layer was sampled (`sample_boundary=true`), it might be in a different part
85+
# of `sampled_shape_data_3d` or combined, depending on `ComplexShape` internal logic.
86+
# The original example implies the main output `initial_condition` is the fluid.
1387

14-
# Returns `InitialCondition`
15-
shape_sampled = ComplexShape(geometry; particle_spacing, density=1.0,
16-
boundary_thickness=5 * particle_spacing,
17-
create_signed_distance_field=true,
18-
sample_boundary=false, point_in_geometry_algorithm)
88+
# Export the Signed Distance Field (SDF) if generated.
89+
# The SDF is a grid where each point stores its shortest distance to the geometry surface.
90+
if generate_signed_distance_field && !isnothing(sampled_shape_data_3d.signed_distance_field)
91+
output_vtk_sdf_3d = "out/$(geometry_filename_stem)_3d_sdf.vti" # SDF is usually a grid (VTI)
92+
trixi2vtk(sampled_shape_data_3d.signed_distance_field, filename=output_vtk_sdf_3d)
93+
println("Exported 3D Signed Distance Field to $output_vtk_sdf_3d")
94+
else
95+
println("Signed Distance Field was not generated or is not available.")
96+
end
1997

20-
trixi2vtk(shape_sampled.initial_condition, filename="initial_condition_fluid")
21-
trixi2vtk(shape_sampled.signed_distance_field)
98+
println("Complex shape 3D example finished. Check 'out/' directory for VTK files.")

0 commit comments

Comments
 (0)