From: Ralf Hartmann Date: Wed, 21 Nov 2001 10:55:21 +0000 (+0000) Subject: New CellAccessor::neighbor_of_coarser_neighbor function. X-Git-Tag: v8.0.0~18601 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c8b2c3c245fba107087e587a01c6bd49cf5b8f4a;p=dealii.git New CellAccessor::neighbor_of_coarser_neighbor function. git-svn-id: https://svn.dealii.org/trunk@5230 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/grid/tria_accessor.h b/deal.II/deal.II/include/grid/tria_accessor.h index 2027172307..e040cd209e 100644 --- a/deal.II/deal.II/include/grid/tria_accessor.h +++ b/deal.II/deal.II/include/grid/tria_accessor.h @@ -18,6 +18,7 @@ #include #include + template class Point; template class Triangulation; @@ -27,6 +28,12 @@ template class TriaActiveIterator; class Line; class Quad; class Hexahedron; +namespace std +{ + template + struct pair; +} + // note: the file tria_accessor.templates.h is included at the end of @@ -221,6 +228,10 @@ class TriaAccessor * Exception */ DeclException0 (ExcNeighborIsCoarser); + /** + * Exception + */ + DeclException0 (ExcNeighborIsNotCoarser); /*@}*/ protected: @@ -1857,15 +1868,25 @@ class CellAccessor : public TriaObjectAccessor * the present cell * (i.e. @p{cell->neighbor(neighbor)->level()} * needs to be equal to - * @p{cell->level()}, since - * otherwise the neighbors of the - * neighbor cell are on a coarser - * level than the present one and - * you can't get back from there - * to this cell. + * @p{cell->level()}. Use the + * @p{neighbor_of_coarser_neighbor} + * function in that case. */ unsigned int neighbor_of_neighbor (const unsigned int neighbor) const; + /** + * This function is a + * generalization of the + * @p{neighbor_of_neighbor} + * function for the case of a + * coarser neighbor. It returns a + * pair of numbers, face_no and + * subface_no, with the following + * property: + * @p{cell->neighbor(neighbor)->face(face_no)->child(subface_no)==cell}. + */ + std::pair neighbor_of_coarser_neighbor (const unsigned int neighbor) const; + /** * Return whether the @p{i}th * vertex or face (depending on diff --git a/deal.II/deal.II/source/grid/tria_accessor.cc b/deal.II/deal.II/source/grid/tria_accessor.cc index 9701c1dd08..11d45037db 100644 --- a/deal.II/deal.II/source/grid/tria_accessor.cc +++ b/deal.II/deal.II/source/grid/tria_accessor.cc @@ -1806,6 +1806,71 @@ unsigned int CellAccessor::neighbor_of_neighbor (const unsigned int neighbo +template +pair +CellAccessor::neighbor_of_coarser_neighbor (const unsigned int neighbor) const +{ + // make sure that the neighbor is + // on a coarser level + Assert (neighbor_level(neighbor) < present_level, + typename TriaAccessor::ExcNeighborIsNotCoarser()); + Assert (neighbor < GeometryInfo::faces_per_cell, + typename TriaAccessor::ExcInvalidNeighbor(neighbor)); + + const TriaIterator > this_face=face(neighbor); + const TriaIterator > neighbor_cell = this->neighbor(neighbor); + + // usually, on regular patches of + // the grid, this cell is just on + // the opposite side of the + // neighbor that the neighbor is of + // this cell. for example in 2d, if + // we want to know the + // neighbor_of_neighbor if + // neighbor==1 (the right + // neighbor), then we will get 3 + // (the left neighbor) in most + // cases. look up this relationship + // in the table provided by + // GeometryInfo and try it + const unsigned int face_no_guess + = GeometryInfo::opposite_face[neighbor]; + + const TriaIterator > face_guess + =neighbor_cell->face(face_no_guess); + + if (face_guess->has_children()) + for (unsigned int subface_no=0; subface_no::subfaces_per_face; ++subface_no) + if (face_guess->child(subface_no)==this_face) + return pair (face_no_guess, subface_no); + + // if the guess was false, then + // we need to loop over all faces + // and subfaces and find the + // number the hard way + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + if (face_no!=face_no_guess) + { + const TriaIterator > face + =neighbor_cell->face(face_no); + if (face->has_children()) + for (unsigned int subface_no=0; subface_no::subfaces_per_face; ++subface_no) + if (face->child(subface_no)==this_face) + return pair (face_no, subface_no); + } + } + + // we should never get here, + // since then we did not find + // our way back... + Assert (false, ExcInternalError()); + return pair (static_cast(-1), + static_cast(-1)); +}; + + + template bool CellAccessor::at_boundary (const unsigned int i) const {