From 58b69618d3a6270385161aa3942822290d036232 Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 11 Jan 2007 03:39:42 +0000 Subject: [PATCH] Implement find_active_cell_around_point also for hp cells and mapping collections. git-svn-id: https://svn.dealii.org/trunk@14294 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/grid/grid_tools.h | 28 +++++++- deal.II/deal.II/source/grid/grid_tools.cc | 81 +++++++++++++++++++++-- 2 files changed, 101 insertions(+), 8 deletions(-) diff --git a/deal.II/deal.II/include/grid/grid_tools.h b/deal.II/deal.II/include/grid/grid_tools.h index 72de802a0b..99936e5a2d 100644 --- a/deal.II/deal.II/include/grid/grid_tools.h +++ b/deal.II/deal.II/include/grid/grid_tools.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors +// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -25,6 +25,15 @@ DEAL_II_NAMESPACE_OPEN +template class DoFHandler; +template class Mapping; +namespace hp +{ + template class DoFHandler; + template class MappingCollection; +} + + /** * This class is a collection of algorithms working on triangulations, * such as shifting or rotating triangulations, but also finding a @@ -373,6 +382,23 @@ class GridTools const Container &container, const Point &p); + /** + * A version of the previous function + * where we use that mapping on a given + * cell that corresponds to the active + * finite element index of that + * cell. This is obviously only useful + * for hp problems, since the active + * finite element index for all other DoF + * handlers is always zero. + */ + template + static + std::pair::active_cell_iterator, Point > + find_active_cell_around_point (const hp::MappingCollection &mapping, + const hp::DoFHandler &container, + const Point &p); + /** * Use the METIS partitioner to generate * a partitioning of the active cells diff --git a/deal.II/deal.II/source/grid/grid_tools.cc b/deal.II/deal.II/source/grid/grid_tools.cc index fb0171b79e..b076738097 100644 --- a/deal.II/deal.II/source/grid/grid_tools.cc +++ b/deal.II/deal.II/source/grid/grid_tools.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors +// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -25,6 +25,7 @@ #include #include #include +#include #include #include @@ -685,7 +686,8 @@ typename Container::active_cell_iterator GridTools::find_active_cell_around_point (const Container &container, const Point &p) { - return find_active_cell_around_point(StaticMappingQ1::mapping, container, p).first; + return find_active_cell_around_point(StaticMappingQ1::mapping, + container, p).first; } @@ -718,18 +720,18 @@ GridTools::find_active_cell_around_point (const Mapping &mapping, for(; cell != endc; ++cell) { - Point p_cell = mapping.transform_real_to_unit_cell(*cell, p); + const Point p_cell = mapping.transform_real_to_unit_cell(*cell, p); // calculate the infinity norm of // the distance vector to the unit cell. - double dist = GeometryInfo::distance_to_unit_cell(p_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)) + if (dist < best_distance || + (dist == best_distance && (*cell)->level() > best_level)) { best_distance = dist; best_level = (*cell)->level(); @@ -737,12 +739,71 @@ GridTools::find_active_cell_around_point (const Mapping &mapping, } } - Assert(best_cell.first.state() == IteratorState::valid, ExcPointNotFound(p)); + Assert (best_cell.first.state() == IteratorState::valid, + ExcPointNotFound(p)); return best_cell; } + +template +std::pair::active_cell_iterator, Point > +GridTools::find_active_cell_around_point (const hp::MappingCollection &mapping, + const hp::DoFHandler &container, + const Point &p) +{ + 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) + { + 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; +} + + + template void GridTools:: @@ -1188,6 +1249,12 @@ GridTools::find_active_cell_around_point (const Mapping &, const hp::DoFHandler &, const Point &); +template +std::pair::active_cell_iterator, Point > +GridTools::find_active_cell_around_point (const hp::MappingCollection &, + const hp::DoFHandler &, + const Point &); + template void GridTools::partition_triangulation (const unsigned int, -- 2.39.5