From 2e04675459a4f0915cfb235e2861cbd6080fbc2f Mon Sep 17 00:00:00 2001 From: Niklas Fehn Date: Wed, 17 Jun 2020 09:56:51 +0200 Subject: [PATCH] fix bug in fast version of GridTools::find_active_cell_around_point() --- source/grid/grid_tools.cc | 45 ++++++++++++++++++++++----------------- 1 file changed, 25 insertions(+), 20 deletions(-) diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 535e1d4392..502111f4d0 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -1682,27 +1682,32 @@ namespace GridTools { auto cell = vertex_to_cells[closest_vertex_index].begin(); std::advance(cell, neighbor_permutation[i]); - const Point p_unit = - mapping.transform_real_to_unit_cell(*cell, p); - if (GeometryInfo::is_inside_unit_cell(p_unit, tolerance)) + if (!(*cell)->is_artificial()) { - cell_and_position.first = *cell; - cell_and_position.second = p_unit; - found_cell = true; - approx_cell = false; - break; - } - // The point is not inside this cell: checking how far outside - // it is and whether we want to use this cell as a backup if we - // can't find a cell within which the point lies. - const double dist = - GeometryInfo::distance_to_unit_cell(p_unit); - if (dist < best_distance) - { - best_distance = dist; - cell_and_position_approx.first = *cell; - cell_and_position_approx.second = p_unit; - approx_cell = true; + const Point p_unit = + mapping.transform_real_to_unit_cell(*cell, p); + if (GeometryInfo::is_inside_unit_cell(p_unit, + tolerance)) + { + cell_and_position.first = *cell; + cell_and_position.second = p_unit; + found_cell = true; + approx_cell = false; + break; + } + // The point is not inside this cell: checking how far + // outside it is and whether we want to use this cell as a + // backup if we can't find a cell within which the point + // lies. + const double dist = + GeometryInfo::distance_to_unit_cell(p_unit); + if (dist < best_distance) + { + best_distance = dist; + cell_and_position_approx.first = *cell; + cell_and_position_approx.second = p_unit; + approx_cell = true; + } } } catch (typename Mapping::ExcTransformationFailed &) -- 2.39.5