fix bugs in cut cell function generate_sampling_points

src/cut_cell_meshes.jl

@@ -151,28 +151,26 @@ function generate_sampling_points(objects, elem, rd, Np_target; N_sampled = 4 *
 
     # map sampled points to the background Cartesian cell
     x_sampled, y_sampled = map_nodes_to_background_cell(elem, r_sampled, s_sampled)
-    is_in_element = fill(true, length(x_sampled))
-    for curve in objects
-        is_in_element .= is_in_element .&& .!(map(x->is_contained(curve, x), zip(x_sampled, y_sampled)))
+    is_in_domain = fill(true, length(x_sampled))
+    for (index, point) in enumerate(zip(x_sampled, y_sampled))
+        is_in_domain[index] = !any(map(obj -> is_contained(obj, point), objects))
     end
-
+    
     # increase number of background points until we are left with `Np_target` sampling points 
-    while sum(is_in_element) < Np_target
+    while sum(is_in_domain) < Np_target
 
-        # check if all the points are in all the objects
-        is_in_element = is_contained.(objects[1], zip(x_sampled, y_sampled)) .== false
-        for _ in 2:length(objects)
-            is_in_element = is_in_element .&& (is_contained.(objects[2], zip(x_sampled, y_sampled)) .== false)
-        end
+        N_sampled *= 2 # double degree of sampling
+        r_sampled, s_sampled = equi_nodes(rd.element_type, N_sampled) # oversampled nodes
+        x_sampled, y_sampled = map_nodes_to_background_cell(elem, r_sampled, s_sampled)
 
-        if sum(is_in_element) < Np_target
-            N_sampled *= 2 # double degree of sampling
-            r_sampled, s_sampled = equi_nodes(rd.element_type, N_sampled) # oversampled nodes
-            x_sampled, y_sampled = map_nodes_to_background_cell(elem, r_sampled, s_sampled)
+        # check if all the points are in all the objects
+        is_in_domain = fill(true, length(x_sampled))
+        for (index, point) in enumerate(zip(x_sampled, y_sampled))
+            is_in_domain[index] = !any(map(obj -> is_contained(obj, point), objects))
         end
     end
 
-    ids_in_element = findall(is_in_element)
+    ids_in_element = findall(is_in_domain)
 
     return x_sampled[ids_in_element], y_sampled[ids_in_element]
 end