Fix convert_particle! (#926)

* fix

* fix face normal

* fix again

* avoid signbit

* fix unit test

* Update src/preprocessing/geometries/geometries.jl

Co-authored-by: Erik Faulhaber <[email protected]>

---------

Co-authored-by: LasNikas <[email protected]>
Co-authored-by: Erik Faulhaber <[email protected]>

Author: 3 people

src/preprocessing/geometries/geometries.jl

@@ -102,8 +102,13 @@ function planar_geometry_to_face(planar_geometry::TriangleMesh)
         throw(ArgumentError("geometry is not planar"))
     end
 
+    # The `face_normal` computed above might not be exactly orthogonal to the plane
+    # spanned by `edge1` and `edge2`, but this is important for some computations later.
+    computed_face_normal = SVector(Tuple(normalize(cross(edge1, edge2))))
+    computed_face_normal *= sign(dot(computed_face_normal, face_normal))
+
     return (; face=(face_vertices[:, 1], face_vertices[:, 2], face_vertices[:, 3]),
-            face_normal=face_normal)
+            face_normal=computed_face_normal)
 end
 
 # According to:

src/schemes/boundary/open_boundary/boundary_zones.jl

@@ -327,7 +327,7 @@ function set_up_boundary_zone(boundary_face, face_normal, density, particle_spac
     # First vector of `spanning_vectors` is normal to the boundary face.
     dot_face_normal = dot(normalize(spanning_set[:, 1]), face_normal)
 
-    if !isapprox(abs(dot_face_normal), 1, rtol=1e-5)
+    if !isapprox(abs(dot_face_normal), 1)
         throw(ArgumentError("`face_normal` is not normal to the boundary face"))
     end
 
@@ -426,6 +426,12 @@ function update_boundary_zone_indices!(system, u, boundary_zones, semi)
                 system.boundary_zone_indices[particle] = zone_id
             end
         end
+
+        # This only occurs if `face_normal` is not exactly normal to the `boundary_face`.
+        # In such cases, particles that are actually outside the simulation domain (outflow particles)
+        # may be incorrectly kept active as inflow particles and therefore cannot be assigned to any boundary zone.
+        # This issue should not occur and has been fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/926 .
+        @assert system.boundary_zone_indices[particle] != 0 "No boundary zone found for active buffer particle"
     end
 
     return system

src/schemes/boundary/open_boundary/system.jl

@@ -301,11 +301,16 @@ end
 @inline function convert_particle!(system::OpenBoundarySystem, fluid_system,
                                    boundary_zone, particle, particle_new,
                                    v, u, v_fluid, u_fluid)
+    # Position relative to the origin of the transition face
     relative_position = current_coords(u, system, particle) - boundary_zone.zone_origin
 
     # Check if particle is in- or outside the fluid domain.
     # `face_normal` is always pointing into the fluid domain.
-    if signbit(dot(relative_position, boundary_zone.face_normal))
+    # Since this function is called for a particle that left the boundary zone,
+    # it is sufficient to check if the dot product between the relative position and the face normal is negative
+    # to determine if it exited the boundary zone through the free surface (outflow).
+    if dot(relative_position, boundary_zone.face_normal) < 0
+        # Particle is outside the fluid domain
         deactivate_particle!(system, particle, u)
 
         return system

test/preprocessing/geometries/geometries.jl

@@ -138,7 +138,7 @@
         expected_face = ([-0.10239515072676975, 0.2644994251485518, -0.36036119092034713],
                          [0.3064669575380171, 0.2392044626289733, -0.10866880239395837],
                          [-0.02275190052262935, 0.299506937268509, -0.034649329562556])
-        ecpected_normal = [0.14372397390844055, 0.979596249614303, -0.14047991694743392]
+        ecpected_normal = [0.1475390528366601, 0.9789123605574772, -0.1412898376948919]
 
         @test any(isapprox.(face, expected_face))
         @test isapprox(face_normal, ecpected_normal)