diff --git a/examples/preprocessing/data/cube.stl b/examples/preprocessing/data/cube.stl
new file mode 100644
index 000000000..84b3e2a48
Binary files /dev/null and b/examples/preprocessing/data/cube.stl differ
diff --git a/examples/preprocessing/data/polyhedron.stl b/examples/preprocessing/data/polyhedron.stl
new file mode 100644
index 000000000..c13bbefce
Binary files /dev/null and b/examples/preprocessing/data/polyhedron.stl differ
diff --git a/examples/preprocessing/data/sphere.stl b/examples/preprocessing/data/sphere.stl
index 02e3b33f4..351ea2404 100644
Binary files a/examples/preprocessing/data/sphere.stl and b/examples/preprocessing/data/sphere.stl differ
diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl
index 817ee4517..df02f1d47 100644
--- a/src/general/initial_condition.jl
+++ b/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`
diff --git a/src/preprocessing/point_in_poly/winding_number_jacobson.jl b/src/preprocessing/point_in_poly/winding_number_jacobson.jl
index 630da58db..63666b4b2 100644
--- a/src/preprocessing/point_in_poly/winding_number_jacobson.jl
+++ b/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
diff --git a/src/setups/complex_shape.jl b/src/setups/complex_shape.jl
index 4b1592ce6..db8eca2c6 100644
--- a/src/setups/complex_shape.jl
+++ b/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)