Add set operations for geometries

src/general/initial_condition.jl

@@ -259,6 +259,30 @@ function Base.setdiff(initial_condition::InitialCondition, initial_conditions...
     return setdiff(result, Base.tail(initial_conditions)...)
 end
 
+function Base.setdiff(initial_condition::InitialCondition,
+                      geometries::Union{Polygon, TriangleMesh}...)
+    geometry = first(geometries)
+
+    if ndims(geometry) != ndims(initial_condition)
+        throw(ArgumentError("all passed geometries must have the same dimensionality as the initial condition"))
+    end
+
+    delete_indices, _ = WindingNumberJacobson(; geometry)(geometry,
+                                                          initial_condition.coordinates)
+
+    coordinates = initial_condition.coordinates[:, .!delete_indices]
+    velocity = initial_condition.velocity[:, .!delete_indices]
+    mass = initial_condition.mass[.!delete_indices]
+    density = initial_condition.density[.!delete_indices]
+    pressure = initial_condition.pressure[.!delete_indices]
+
+    result = InitialCondition{ndims(initial_condition)}(coordinates, velocity, mass,
+                                                        density, pressure,
+                                                        initial_condition.particle_spacing)
+
+    return setdiff(result, Base.tail(geometries)...)
+end
+
 Base.setdiff(initial_condition::InitialCondition) = initial_condition
 
 function Base.intersect(initial_condition::InitialCondition, initial_conditions...)
@@ -288,6 +312,30 @@ function Base.intersect(initial_condition::InitialCondition, initial_conditions.
     return intersect(result, Base.tail(initial_conditions)...)
 end
 
+function Base.intersect(initial_condition::InitialCondition,
+                        geometries::Union{Polygon, TriangleMesh}...)
+    geometry = first(geometries)
+
+    if ndims(geometry) != ndims(initial_condition)
+        throw(ArgumentError("all passed geometries must have the same dimensionality as the initial condition"))
+    end
+
+    keep_indices, _ = WindingNumberJacobson(; geometry)(geometry,
+                                                        initial_condition.coordinates)
+
+    coordinates = initial_condition.coordinates[:, keep_indices]
+    velocity = initial_condition.velocity[:, keep_indices]
+    mass = initial_condition.mass[keep_indices]
+    density = initial_condition.density[keep_indices]
+    pressure = initial_condition.pressure[keep_indices]
+
+    result = InitialCondition{ndims(initial_condition)}(coordinates, velocity, mass,
+                                                        density, pressure,
+                                                        initial_condition.particle_spacing)
+
+    return intersect(result, Base.tail(geometries)...)
+end
+
 Base.intersect(initial_condition::InitialCondition) = initial_condition
 
 # Find particles in `coords1` that are closer than `max_distance` to any particle in `coords2`

src/preprocessing/point_in_poly/winding_number_jacobson.jl

@@ -69,7 +69,7 @@ struct WindingNumberJacobson{ELTYPE, W}
     winding               :: W
 
     function WindingNumberJacobson(; geometry=nothing, winding_number_factor=sqrt(eps()),
-                                   hierarchical_winding=false)
+                                   hierarchical_winding=true)
         if hierarchical_winding && geometry isa Nothing
             throw(ArgumentError("`geometry` must be of type `Polygon` (2D) or `TriangleMesh` (3D) when using hierarchical winding"))
         end

src/setups/complex_shape.jl

@@ -44,7 +44,7 @@ For more information about the method see [`WindingNumberJacobson`](@ref) or [`W
 function ComplexShape(geometry::Union{TriangleMesh, Polygon}; particle_spacing, density,
                       pressure=0.0, mass=nothing, velocity=zeros(ndims(geometry)),
                       point_in_geometry_algorithm=WindingNumberJacobson(; geometry,
-                                                                        hierarchical_winding=false,
+                                                                        hierarchical_winding=true,
                                                                         winding_number_factor=sqrt(eps())),
                       store_winding_number=false, grid_offset::Real=0.0,
                       max_nparticles=10^7, pad_initial_particle_grid=2particle_spacing)