Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Particle packing #529

Open
wants to merge 322 commits into
base: main
Choose a base branch
from
Open

Conversation

LasNikas
Copy link
Collaborator

@LasNikas LasNikas commented May 15, 2024

based on #439 and #436

@svchb
Copy link
Collaborator

svchb commented Nov 13, 2024

@LasNikas Please fix the tests first than I will review this PR.

Copy link
Member

@efaulhaber efaulhaber left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Review of the first part. Remaining: Examples, NHS, signed distance and following.

Comment on lines +265 to +270
Pages = [joinpath("preprocessing", "particle_packing", "system.jl")]
```

```@autodocs
Modules = [TrixiParticles]
Pages = [joinpath("preprocessing", "particle_packing", "signed_distance.jl")]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not add the whole particle packing directory?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See e.g. the setups:

```@autodocs
Modules = [TrixiParticles]
Pages = map(file -> joinpath("setups", file), readdir(joinpath("..", "src", "setups")))
```

examples/preprocessing/complex_shape_3d.jl Outdated Show resolved Hide resolved
examples/preprocessing/data/sharp_triangle.asc Outdated Show resolved Hide resolved
src/TrixiParticles.jl Outdated Show resolved Hide resolved
src/general/initial_condition.jl Outdated Show resolved Hide resolved
src/general/initial_condition.jl Outdated Show resolved Hide resolved
Comment on lines +85 to +93
vertex_normals = Vector{SVector{2, eltype(geometry)}}()
vertices = Vector{SVector{2, eltype(geometry)}}()

for edge in eachface(geometry)
push!(vertex_normals, geometry.vertex_normals[edge][1])
push!(vertex_normals, geometry.vertex_normals[edge][2])
push!(vertices, geometry.vertices[geometry.edge_vertices_ids[edge][1]])
push!(vertices, geometry.vertices[geometry.edge_vertices_ids[edge][2]])
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why did this change?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This version is more robust. consider a deleteat! on geometry.

src/preprocessing/particle_packing/rhs.jl Outdated Show resolved Hide resolved
src/preprocessing/particle_packing/signed_distance.jl Outdated Show resolved Hide resolved
src/preprocessing/particle_packing/nhs_faces.jl Outdated Show resolved Hide resolved
src/preprocessing/particle_packing/nhs_faces.jl Outdated Show resolved Hide resolved
Comment on lines 166 to 226
function triangle_plane_intersection(point, plane_normal, min_corner, max_corner, cell_size)
dirx = SVector(cell_size[1], zero(eltype(point)), zero(eltype(point)))
diry = SVector(zero(eltype(point)), cell_size[2], zero(eltype(point)))
dirz = SVector(zero(eltype(point)), zero(eltype(point)), cell_size[3])

cell_vertex_1 = min_corner
cell_vertex_2 = min_corner + dirx

pos_diff_1 = cell_vertex_1 - point
pos_diff_2 = cell_vertex_2 - point

# Corners: bottom north-west and bottom north-east
dot1 = sign(dot(pos_diff_1, plane_normal))
dot2 = sign(dot(pos_diff_2, plane_normal))
!(dot1 == dot2) && return true

cell_vertex_3 = min_corner + diry
pos_diff_3 = cell_vertex_3 - point

# Corners: bottom north-east and top north-west
dot3 = sign(dot(pos_diff_3, plane_normal))
!(dot2 == dot3) && return true

cell_vertex_4 = min_corner + dirz
pos_diff_4 = cell_vertex_4 - point

# Corners: top norht-west and bottom south-west
dot4 = sign(dot(pos_diff_4, plane_normal))
!(dot3 == dot4) && return true

cell_vertex_5 = max_corner
pos_diff_5 = cell_vertex_5 - point

# Corners: bottom south-west and top south-east
dot5 = sign(dot(pos_diff_5, plane_normal))
!(dot4 == dot5) && return true

cell_vertex_6 = max_corner - dirx
pos_diff_6 = cell_vertex_6 - point

# Corners: top south-east and top south-west
dot6 = sign(dot(pos_diff_6, plane_normal))
!(dot5 == dot6) && return true

cell_vertex_7 = max_corner - diry
pos_diff_7 = cell_vertex_7 - point

# Corners: top south-west and bottom south-east
dot7 = sign(dot(pos_diff_7, plane_normal))
!(dot6 == dot7) && return true

cell_vertex_8 = max_corner - dirz
pos_diff_8 = cell_vertex_8 - point

# Corners: bottom south-east and top north-east
dot8 = sign(dot(pos_diff_8, plane_normal))
!(dot7 == dot8) && return true

# All edge vertices are on one side of the plane
return false
end
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
function triangle_plane_intersection(point, plane_normal, min_corner, max_corner, cell_size)
dirx = SVector(cell_size[1], zero(eltype(point)), zero(eltype(point)))
diry = SVector(zero(eltype(point)), cell_size[2], zero(eltype(point)))
dirz = SVector(zero(eltype(point)), zero(eltype(point)), cell_size[3])
cell_vertex_1 = min_corner
cell_vertex_2 = min_corner + dirx
pos_diff_1 = cell_vertex_1 - point
pos_diff_2 = cell_vertex_2 - point
# Corners: bottom north-west and bottom north-east
dot1 = sign(dot(pos_diff_1, plane_normal))
dot2 = sign(dot(pos_diff_2, plane_normal))
!(dot1 == dot2) && return true
cell_vertex_3 = min_corner + diry
pos_diff_3 = cell_vertex_3 - point
# Corners: bottom north-east and top north-west
dot3 = sign(dot(pos_diff_3, plane_normal))
!(dot2 == dot3) && return true
cell_vertex_4 = min_corner + dirz
pos_diff_4 = cell_vertex_4 - point
# Corners: top norht-west and bottom south-west
dot4 = sign(dot(pos_diff_4, plane_normal))
!(dot3 == dot4) && return true
cell_vertex_5 = max_corner
pos_diff_5 = cell_vertex_5 - point
# Corners: bottom south-west and top south-east
dot5 = sign(dot(pos_diff_5, plane_normal))
!(dot4 == dot5) && return true
cell_vertex_6 = max_corner - dirx
pos_diff_6 = cell_vertex_6 - point
# Corners: top south-east and top south-west
dot6 = sign(dot(pos_diff_6, plane_normal))
!(dot5 == dot6) && return true
cell_vertex_7 = max_corner - diry
pos_diff_7 = cell_vertex_7 - point
# Corners: top south-west and bottom south-east
dot7 = sign(dot(pos_diff_7, plane_normal))
!(dot6 == dot7) && return true
cell_vertex_8 = max_corner - dirz
pos_diff_8 = cell_vertex_8 - point
# Corners: bottom south-east and top north-east
dot8 = sign(dot(pos_diff_8, plane_normal))
!(dot7 == dot8) && return true
# All edge vertices are on one side of the plane
return false
end
function triangle_plane_intersection(point, plane_normal, min_corner, max_corner, cell_size)
# Generate all eight corners of the cell (axis-aligned bounding box).
# Each corner is calculated by adding a combination of cell_size components
# to the min_corner, based on the indices i, j, k, which take values 0 or 1.
# This creates all possible combinations, resulting in the eight corners of the cell.
corners = [
min_corner + SVector(i * cell_size[1], j * cell_size[2], k * cell_size[3])
for i in (0, 1), j in (0, 1), k in (0, 1)
]
# For each corner, compute the vector from a known point on the plane to the corner.
# Then, calculate the dot product of this vector with the plane's normal vector.
# The dot product indicates the relative position of the corner to the plane:
# - Positive value: corner is on the same side as the normal vector points.
# - Zero: corner lies exactly on the plane.
# - Negative value: corner is on the opposite side from where the normal vector points.
# The sign() function converts the dot product into -1, 0, or 1 accordingly.
signs = [sign(dot(corner - point, plane_normal)) for corner in corners]
# Determine if the plane intersects the cell by checking the signs.
# - If all signs are the same (all corners on one side of the plane), there's no intersection.
# - If signs vary (some corners on different sides), the plane intersects the cell.
# If maximum(signs) > minimum(signs), it means there are both positive and negative signs.
return maximum(signs) > minimum(signs)
end

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for this. Your version is definitely more readable. I'm wondering if my version is slightly faster as we don't have to check each combination in best case.

However, I prefer your readable version. What's @efaulhaber opinion on this?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please just test what the performance is like compared to the current version.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's not that easy. I mean it's case dependent. In best case (first condition is true), it should be faster. But in worst case it should be same performance, no?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes but I would just compare it for a large example and check if the impact is actually measurable.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the performance difference is negligible, that's why I also prefer your version.

src/preprocessing/particle_packing/nhs_faces.jl Outdated Show resolved Hide resolved
Comment on lines +17 to +33
- `boundary`: Geometry returned by [`load_geometry`](@ref).
- `background_pressure`: Constant background pressure to physically pack the particles.
A large `background_pressure` can cause high accelerations which requires a properly adjusted time-step criterion.
- `tlsph`: With the [`TotalLagrangianSPHSystem`](@ref), particles need to be placed
on the boundary of the shape and not one particle radius away, as for fluids.
When `tlsph=true`, particles will be placed on the boundary of the shape.
- `is_boundary`: When `is_boundary=true`, boundary particles will be sampled and packed in an offset surface of the `boundary`.
The thickness of the boundary is specified by passing the [`SignedDistanceField`](@ref) of `boundary` with:
- `use_for_boundary_packing=true`
- `max_signed_distance=boundary_thickness`
- `signed_distance_field`: To constrain particles onto the surface, the information about
the signed distance from a particle to a face is required.
The precalculated signed distances will be interpolated to each particle during the packing procedure.
- `smoothing_kernel`: Smoothing kernel to be used for this system.
See [Smoothing Kernels](@ref smoothing_kernel).
- `smoothing_length`: Smoothing length to be used for this system.
See [Smoothing Kernels](@ref smoothing_kernel).
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please align

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it is aligned?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

we have usually also aligned the ":"

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ah didn't know that. But it seems that we haven't done that consistently.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

like so?
image

src/preprocessing/particle_packing/system.jl Outdated Show resolved Hide resolved
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request high priority
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants