From e355f32f875ea9c12b2228dc40e7b5c56f61748f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 22 Feb 2023 00:22:58 -0600 Subject: [PATCH] fix bugs in cut cell function `generate_sampling_points` --- src/cut_cell_meshes.jl | 28 +++++++++++++--------------- 1 file changed, 13 insertions(+), 15 deletions(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index e326cb2a..8002dfd9 100644 --- a/src/cut_cell_meshes.jl +++ b/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