From 5581e81db6242033725a5fbb18d57efd746e4b18 Mon Sep 17 00:00:00 2001 From: Richard Schussnig Date: Mon, 22 May 2023 13:40:48 +0200 Subject: [PATCH] fix bug in find_active_cell_around_point and add test --- source/grid/grid_tools_dof_handlers.cc | 22 +- .../search_adjacent_cells.cc | 282 ++++++++++++++++++ .../search_adjacent_cells.mpirun=1.output | 6 + 3 files changed, 304 insertions(+), 6 deletions(-) create mode 100644 tests/remote_point_evaluation/search_adjacent_cells.cc create mode 100644 tests/remote_point_evaluation/search_adjacent_cells.mpirun=1.output diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc index 8ada7e4b1b..08453c597e 100644 --- a/source/grid/grid_tools_dof_handlers.cc +++ b/source/grid/grid_tools_dof_handlers.cc @@ -495,6 +495,22 @@ namespace GridTools { if ((*cell)->is_artificial() == false) { + // marked_vertices are used to filter cell candidates + if (marked_vertices.size() > 0) + { + bool any_vertex_marked = false; + for (const auto &v : (*cell)->vertex_indices()) + { + if (marked_vertices[(*cell)->vertex_index(v)]) + { + any_vertex_marked = true; + break; + } + } + if (!any_vertex_marked) + continue; + } + try { const Point p_cell = @@ -540,12 +556,6 @@ namespace GridTools // update the number of cells searched cells_searched += adjacent_cells.size(); - // if the user provided a custom mask for vertices, - // terminate the search without trying to expand the search - // to all cells of the triangulation, as done below. - if (marked_vertices.size() > 0) - cells_searched = n_active_cells; - // if we have not found the cell in // question and have not yet searched every // cell, we expand our search to diff --git a/tests/remote_point_evaluation/search_adjacent_cells.cc b/tests/remote_point_evaluation/search_adjacent_cells.cc new file mode 100644 index 0000000000..9ea900755d --- /dev/null +++ b/tests/remote_point_evaluation/search_adjacent_cells.cc @@ -0,0 +1,282 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include + +#include + +#include + +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include + +#include +#include + +#include "../tests.h" + +template +double +displacement(const Point & point, + const unsigned int component, + const double factor) +{ + if (component == 0) + return (factor * std::pow(point[1], 2) * std::pow(point[2], 2)); + else + return 0.0; +} + +template +void +construct_triangulation(Triangulation &tria) +{ + double const L_F = 0.7; + double const B_F = 1.0; + double const H_F = 0.5; + + double const T_S = 0.05; + double const B_S = 0.6; + double const H_S = 0.4; + + double const L_IN = 0.6; + + unsigned int const N_CELLS_X_OUTFLOW = 1; + unsigned int const N_CELLS_Y_LOWER = 2; + unsigned int const N_CELLS_Z_MIDDLE = 2; + + std::vector> tria_vec; + tria_vec.resize(4); + + dealii::GridGenerator::subdivided_hyper_rectangle( + tria_vec[0], + std::vector({N_CELLS_X_OUTFLOW, 1, N_CELLS_Z_MIDDLE}), + dealii::Point<3>(L_IN + T_S, H_S - H_F / 2.0, -B_S / 2.0), + dealii::Point<3>(L_F, H_F / 2.0, B_S / 2.0)); + + dealii::GridGenerator::subdivided_hyper_rectangle( + tria_vec[1], + std::vector( + {N_CELLS_X_OUTFLOW, N_CELLS_Y_LOWER, N_CELLS_Z_MIDDLE}), + dealii::Point<3>(L_IN + T_S, -H_F / 2.0, -B_S / 2.0), + dealii::Point<3>(L_F, H_S - H_F / 2.0, B_S / 2.0)); + + dealii::GridGenerator::subdivided_hyper_rectangle( + tria_vec[2], + std::vector({N_CELLS_X_OUTFLOW, N_CELLS_Y_LOWER, 1}), + dealii::Point<3>(L_IN + T_S, -H_F / 2.0, -B_F / 2.0), + dealii::Point<3>(L_F, H_S - H_F / 2.0, -B_S / 2.0)); + + dealii::GridGenerator::subdivided_hyper_rectangle( + tria_vec[3], + std::vector({1, N_CELLS_Y_LOWER, 1}), + dealii::Point<3>(L_IN, -H_F / 2.0, -B_F / 2.0), + dealii::Point<3>(L_IN + T_S, H_S - H_F / 2.0, -B_S / 2.0)); + + std::vector const *> tria_vec_ptr(tria_vec.size()); + for (unsigned int i = 0; i < tria_vec.size(); ++i) + tria_vec_ptr[i] = &tria_vec[i]; + + dealii::GridGenerator::merge_triangulations(tria_vec_ptr, tria, 1.e-10); +} + +template +void +test(const unsigned int mapping_degree, + const double mapping_factor, + const double tolerance, + const unsigned int n_points) +{ + Triangulation tria; + construct_triangulation(tria); + + FESystem fe(FE_Q(mapping_degree), dim); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + Quadrature quadrature(fe.base_element(0).get_unit_support_points()); + FEValues fe_values(fe, + quadrature, + update_quadrature_points | update_values); + MappingQCache mapping(mapping_degree); + mapping.initialize( + tria, + [&](const typename dealii::Triangulation::cell_iterator &cell_tria) + -> std::vector> { + std::vector> grid_coordinates(quadrature.size()); + + fe_values.reinit(cell_tria); + // extract displacement and add to original position + for (unsigned int i = 0; i < grid_coordinates.size(); ++i) + { + grid_coordinates[i] = fe_values.quadrature_point(i); + grid_coordinates[i][0] += + displacement(grid_coordinates[i], 0, mapping_factor); + } + return grid_coordinates; + }); + + // initialize points with noise + std::vector> points(n_points * n_points); + unsigned int point_idx = 0; + for (unsigned int i = 0; i < n_points; ++i) + { + for (unsigned int j = 0; j < n_points; ++j) + { + points[point_idx][0] = 0.65; + points[point_idx][1] = -0.25 + 0.4 * (1.0 / (n_points - 1) * j); + points[point_idx][2] = -0.3 + 0.6 * (1.0 / (n_points - 1) * i); + + points[point_idx][0] += + displacement(points[point_idx], 0, mapping_factor) + + (tolerance * 0.1) * (-0.5 + (rand() / (double)RAND_MAX)); + + point_idx += 1; + } + } + + // initialize RPE without any marked points + dealii::Utilities::MPI::RemotePointEvaluation rpe(tolerance); + rpe.reinit(points, tria, mapping); + + unsigned int n_points_not_found_rpe = 0; + std::vector> points_not_found; + if (!rpe.all_points_found()) + { + // get vector of points not found + points_not_found.reserve(points.size()); + + for (unsigned int i = 0; i < points.size(); ++i) + { + if (!rpe.point_found(i)) + { + n_points_not_found_rpe += 1; + points_not_found.push_back(points[i]); + } + } + } + std::cout << "points not found (no marked points) : " + << n_points_not_found_rpe << "\n"; + + // initialize RPE with all points marked + std::vector marked_vertices(tria.n_vertices(), true); + dealii::Utilities::MPI::RemotePointEvaluation rpe2( + tolerance, false, 0, [marked_vertices]() { return marked_vertices; }); + + rpe2.reinit(points, tria, mapping); + + unsigned int n_points_not_found_rpe2 = 0; + std::vector> points_not_found2; + if (!rpe2.all_points_found()) + { + // get vector of points not found + std::vector> points_not_found2; + points_not_found2.reserve(points.size()); + + for (unsigned int i = 0; i < points.size(); ++i) + { + if (!rpe2.point_found(i)) + { + n_points_not_found_rpe2 += 1; + points_not_found2.push_back(points[i]); + } + } + } + std::cout << "points not found (all points marked) : " + << n_points_not_found_rpe2 << "\n"; + + // output in case of failure + if (n_points_not_found_rpe > 0 || n_points_not_found_rpe2 > 0) + { + // output triangulation + DataOut data_out; + DataOutBase::VtkFlags flags; + flags.write_higher_order_cells = true; + data_out.set_flags(flags); + data_out.attach_triangulation(tria); + data_out.build_patches(mapping, + mapping_degree, + DataOut::curved_inner_cells); + std::ofstream stream("grid.vtu"); + data_out.write_vtu(stream); + + // enclosing dummy triangulation for point plots + dealii::BoundingBox bounding_box(points); + auto const boundary_points = + bounding_box.create_extended_relative(1e-3).get_boundary_points(); + + dealii::Triangulation particle_dummy_tria; + dealii::GridGenerator::hyper_rectangle(particle_dummy_tria, + boundary_points.first, + boundary_points.second); + + dealii::MappingQGeneric particle_dummy_mapping( + 1 /* mapping_degree */); + + // output all points + { + dealii::Particles::ParticleHandler particle_handler( + particle_dummy_tria, particle_dummy_mapping); + particle_handler.insert_particles(points); + dealii::Particles::DataOut particle_output; + particle_output.build_patches(particle_handler); + std::ofstream filestream("all_points.vtu"); + particle_output.write_vtu(filestream); + } + + // output points not found (no vertices marked) + { + dealii::Particles::ParticleHandler particle_handler( + particle_dummy_tria, particle_dummy_mapping); + particle_handler.insert_particles(points_not_found); + dealii::Particles::DataOut particle_output; + particle_output.build_patches(particle_handler); + std::ofstream filestream("points_not_found.vtu"); + particle_output.write_vtu(filestream); + } + + // output points not found + { + dealii::Particles::ParticleHandler particle_handler( + particle_dummy_tria, particle_dummy_mapping); + particle_handler.insert_particles(points_not_found2); + dealii::Particles::DataOut particle_output; + particle_output.build_patches(particle_handler); + std::ofstream filestream("points_not_found2.vtu"); + particle_output.write_vtu(filestream); + } + } +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + deallog << std::setprecision(8) << std::fixed; + + test<3>(2, 100.0, 1e-5, 100); + test<3>(3, 100.0, 1e-5, 100); + test<3>(4, 100.0, 1e-5, 100); +} diff --git a/tests/remote_point_evaluation/search_adjacent_cells.mpirun=1.output b/tests/remote_point_evaluation/search_adjacent_cells.mpirun=1.output new file mode 100644 index 0000000000..7a35e449d0 --- /dev/null +++ b/tests/remote_point_evaluation/search_adjacent_cells.mpirun=1.output @@ -0,0 +1,6 @@ +points not found (no marked points) : 0 +points not found (all points marked) : 0 +points not found (no marked points) : 0 +points not found (all points marked) : 0 +points not found (no marked points) : 0 +points not found (all points marked) : 0 \ No newline at end of file -- 2.39.5