From 3e82e504cef3a0a13120b726681f14ac24eb9bc9 Mon Sep 17 00:00:00 2001 From: goll Date: Thu, 26 Jan 2012 09:50:36 +0000 Subject: [PATCH] find_active_cell_around_point for the hp case works with a MappingCollection containing only on mapping. git-svn-id: https://svn.dealii.org/trunk@24934 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/source/grid/grid_tools.cc | 88 ++++++++++++++++++------------- 1 file changed, 52 insertions(+), 36 deletions(-) diff --git a/deal.II/source/grid/grid_tools.cc b/deal.II/source/grid/grid_tools.cc index 1cfd13be65..e26e2e7d8d 100644 --- a/deal.II/source/grid/grid_tools.cc +++ b/deal.II/source/grid/grid_tools.cc @@ -840,52 +840,68 @@ namespace GridTools const hp::DoFHandler &container, const Point &p) { + Assert ((mapping.size() == 1) || + (mapping.size() == container.get_fe().size()), + ExcMessage ("Mapping collection needs to have either size 1 " + "or size equal to the number of elements in " + "the FECollection.")); + typedef typename hp::DoFHandler::active_cell_iterator cell_iterator; - // The best distance is set to the - // maximum allowable distance from - // the unit cell; we assume a - // max. deviation of 1e-10 - double best_distance = 1e-10; - int best_level = -1; std::pair > best_cell; - - // Find closest vertex and determine - // all adjacent cells - unsigned int vertex = find_closest_vertex(container, p); - - std::vector adjacent_cells = - find_cells_adjacent_to_vertex(container, vertex); - - typename std::vector::const_iterator - cell = adjacent_cells.begin(), - endc = adjacent_cells.end(); - - for(; cell != endc; ++cell) + //If we have only one element in the MappingCollection, + //we use find_active_cell_around_point using only one + //mapping. + if(mapping.size()==1) + best_cell = find_active_cell_around_point(mapping[0], container, p); + else { - const Point p_cell - = mapping[(*cell)->active_fe_index()].transform_real_to_unit_cell(*cell, p); - // calculate the infinity norm of - // the distance vector to the unit cell. - const double dist = GeometryInfo::distance_to_unit_cell(p_cell); - // We compare if the point is inside the - // unit cell (or at least not too far - // outside). If it is, it is also checked - // that the cell has a more refined state - if (dist < best_distance || - (dist == best_distance && (*cell)->level() > best_level)) - { - best_distance = dist; - best_level = (*cell)->level(); - best_cell = std::make_pair(*cell, p_cell); - } + // The best distance is set to the + // maximum allowable distance from + // the unit cell; we assume a + // max. deviation of 1e-10 + double best_distance = 1e-10; + int best_level = -1; + + + // Find closest vertex and determine + // all adjacent cells + unsigned int vertex = find_closest_vertex(container, p); + + std::vector adjacent_cells = + find_cells_adjacent_to_vertex(container, vertex); + + typename std::vector::const_iterator + cell = adjacent_cells.begin(), + endc = adjacent_cells.end(); + + for(; cell != endc; ++cell) + { + const Point p_cell + = mapping[(*cell)->active_fe_index()].transform_real_to_unit_cell(*cell, p); + + // calculate the infinity norm of + // the distance vector to the unit cell. + const double dist = GeometryInfo::distance_to_unit_cell(p_cell); + + // We compare if the point is inside the + // unit cell (or at least not too far + // outside). If it is, it is also checked + // that the cell has a more refined state + if (dist < best_distance || + (dist == best_distance && (*cell)->level() > best_level)) + { + best_distance = dist; + best_level = (*cell)->level(); + best_cell = std::make_pair(*cell, p_cell); + } } Assert (best_cell.first.state() == IteratorState::valid, ExcPointNotFound(p)); - + } return best_cell; } -- 2.39.5