From 9766b0e990cbd27dd4df0fa7f1b8b5f2fa2d63a4 Mon Sep 17 00:00:00 2001 From: wolf Date: Wed, 2 Feb 2000 12:02:30 +0000 Subject: [PATCH] Implement and use CellAccessor::neighbor_of_neighbor git-svn-id: https://svn.dealii.org/trunk@2325 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/grid/geometry_info.h | 9 +++ deal.II/deal.II/include/grid/tria_accessor.h | 36 ++++++++- deal.II/deal.II/source/grid/geometry_info.cc | 13 ++++ deal.II/deal.II/source/grid/tria.cc | 76 +++++++++---------- deal.II/deal.II/source/grid/tria_accessor.cc | 52 +++++++++++++ .../source/numerics/error_estimator.cc | 19 ++--- 6 files changed, 149 insertions(+), 56 deletions(-) diff --git a/deal.II/deal.II/include/grid/geometry_info.h b/deal.II/deal.II/include/grid/geometry_info.h index eec8024104..f79fc7993e 100644 --- a/deal.II/deal.II/include/grid/geometry_info.h +++ b/deal.II/deal.II/include/grid/geometry_info.h @@ -127,6 +127,15 @@ struct GeometryInfo static const unsigned int hexes_per_cell = (2*GeometryInfo::hexes_per_cell + GeometryInfo::quads_per_cell); + /** + * List of numbers which is + * denote which face is opposite + * to a given face. In 1d, this + * list is #{1,0}#, in 2d #{2, 3, 0, 1}#, + * in 3d #{1, 0, 4, 5, 2, 3}#. + */ + static const unsigned int opposite_face[faces_per_cell]; + /** * This field store which child cells * are adjacent to a certain face of diff --git a/deal.II/deal.II/include/grid/tria_accessor.h b/deal.II/deal.II/include/grid/tria_accessor.h index 6d87bb6143..29a92ffbeb 100644 --- a/deal.II/deal.II/include/grid/tria_accessor.h +++ b/deal.II/deal.II/include/grid/tria_accessor.h @@ -207,6 +207,10 @@ class TriaAccessor * Exception */ DeclException0 (ExcCantCompareIterators); + /** + * Exception + */ + DeclException0 (ExcNeighborIsCoarser); /*@}*/ protected: @@ -1711,10 +1715,11 @@ class TriaObjectAccessor<3, dim> : public TriaAccessor * the possibility to check whether they are at the boundary etc. This class * offers access to all this data. * - * @author Wolfgang Bangerth, 1998 + * @author Wolfgang Bangerth, 1998, 1999, 2000 */ template -class CellAccessor : public TriaObjectAccessor { +class CellAccessor : public TriaObjectAccessor +{ public: /** * Propagate the AccessorData type @@ -1761,6 +1766,33 @@ class CellAccessor : public TriaObjectAccessor { void set_neighbor (const unsigned int i, const TriaIterator > &pointer) const; + /** + * Return the how-many'th + * neighbor this cell is of + * #cell->neighbor(neighbor)#, + * i.e. return the number #n# + * such that + * #cell->neighbor(neighbor)->neighbor(n)==cell#. This + * function is the right one if + * you want to know how to get + * back from a neighbor to the + * present cell. + * + * Note that this operation is + * only useful if the neighbor is + * not on a coarser level than + * the present cell + * (i.e. #cell->neighbor(neighbor)->level()# + * needs to be equal to + * #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. + */ + unsigned int neighbor_of_neighbor (const unsigned int neighbor) const; + /** * Return whether the #i#th vertex or * face (depending on the dimension) is diff --git a/deal.II/deal.II/source/grid/geometry_info.cc b/deal.II/deal.II/source/grid/geometry_info.cc index b68cc09f4f..a5cfea3080 100644 --- a/deal.II/deal.II/source/grid/geometry_info.cc +++ b/deal.II/deal.II/source/grid/geometry_info.cc @@ -12,4 +12,17 @@ //const unsigned int GeometryInfo::children_per_cell; +template <> +const unsigned int GeometryInfo<1>::opposite_face[GeometryInfo<1>::faces_per_cell] += { 0, 1 }; + + +template <> +const unsigned int GeometryInfo<2>::opposite_face[GeometryInfo<2>::faces_per_cell] += { 2, 3, 0, 1 }; + + +template <> +const unsigned int GeometryInfo<3>::opposite_face[GeometryInfo<3>::faces_per_cell] += { 1, 0, 4, 5, 2, 3 }; diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index 4834e6bbcd..dae91a404f 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -4242,43 +4242,41 @@ void Triangulation<2>::execute_refinement () { // two children which // we can use. cell_iterator neighbor = cell->neighbor(nb); - for (unsigned int nb_nb=0; nb_nb<4; ++nb_nb) - if (neighbor->neighbor(nb_nb)==cell) // this cell is the nb_nb-th // neighbor or neighbor(nb) - { - neighbors_neighbor[2*nb] = neighbors_neighbor[2*nb+1] = nb_nb; - // vertex 1 of child 0 - // is always the interior - // one - new_vertices[2*nb+1] = neighbor->line(nb_nb) - ->child(0)->vertex_index(1); - - if (nb < 2) - { - new_lines[2*nb] = neighbor->line(nb_nb)->child(0); - new_lines[2*nb+1]= neighbor->line(nb_nb)->child(1); - } else { - // lines 2 and 3 have - // opposite sense - new_lines[2*nb] = neighbor->line(nb_nb)->child(1); - new_lines[2*nb+1]= neighbor->line(nb_nb)->child(0); - }; - - // finally find out which - // are the two neighbor - // subcells, adjacent to - // the two sublines - const int child_mapping[4][2] = {{0,1},{1,2},{3,2},{0,3}}; - if (nb < 2) - { - neighbors[2*nb] = neighbor->child(child_mapping[nb_nb][0]); - neighbors[2*nb+1]= neighbor->child(child_mapping[nb_nb][1]); - } else { - neighbors[2*nb] = neighbor->child(child_mapping[nb_nb][1]); - neighbors[2*nb+1]= neighbor->child(child_mapping[nb_nb][0]); - }; - }; + const unsigned int nb_nb = cell->neighbor_of_neighbor (nb); + + neighbors_neighbor[2*nb] = neighbors_neighbor[2*nb+1] = nb_nb; + // vertex 1 of child 0 + // is always the interior + // one + new_vertices[2*nb+1] = neighbor->line(nb_nb) + ->child(0)->vertex_index(1); + + if (nb < 2) + { + new_lines[2*nb] = neighbor->line(nb_nb)->child(0); + new_lines[2*nb+1]= neighbor->line(nb_nb)->child(1); + } else { + // lines 2 and 3 have + // opposite sense + new_lines[2*nb] = neighbor->line(nb_nb)->child(1); + new_lines[2*nb+1]= neighbor->line(nb_nb)->child(0); + }; + + // finally find out which + // are the two neighbor + // subcells, adjacent to + // the two sublines + static const unsigned int child_mapping[4][2] = {{0,1},{1,2},{3,2},{0,3}}; + if (nb < 2) + { + neighbors[2*nb] = neighbor->child(child_mapping[nb_nb][0]); + neighbors[2*nb+1]= neighbor->child(child_mapping[nb_nb][1]); + } else { + neighbors[2*nb] = neighbor->child(child_mapping[nb_nb][1]); + neighbors[2*nb+1]= neighbor->child(child_mapping[nb_nb][0]); + }; } else @@ -5560,7 +5558,7 @@ void Triangulation<3>::execute_refinement () { for (unsigned int face=0; face::faces_per_cell; ++face) { - cell_iterator neighbor = hex->neighbor(face); + const cell_iterator neighbor = hex->neighbor(face); // if no neighbor if (neighbor.state() != valid) @@ -5604,11 +5602,7 @@ void Triangulation<3>::execute_refinement () { // of the neighbor // adjacent to which // the present cell is - unsigned int nb_nb; - for (nb_nb=0; nb_nb::faces_per_cell; - ++nb_nb) - if (neighbor->neighbor(nb_nb) == hex) - break; + const unsigned int nb_nb = hex->neighbor_of_neighbor(face); Assert (nb_nb::faces_per_cell, ExcInternalError()); diff --git a/deal.II/deal.II/source/grid/tria_accessor.cc b/deal.II/deal.II/source/grid/tria_accessor.cc index 301e7093b4..ae058d1562 100644 --- a/deal.II/deal.II/source/grid/tria_accessor.cc +++ b/deal.II/deal.II/source/grid/tria_accessor.cc @@ -1599,6 +1599,58 @@ void CellAccessor::set_neighbor (const unsigned int i, +template +unsigned int CellAccessor::neighbor_of_neighbor (const unsigned int neighbor) const +{ + // make sure that the neighbor is + // not on a coarser level + Assert (neighbor_level(neighbor) == present_level, + ExcNeighborIsCoarser()); + Assert (neighbor < GeometryInfo::faces_per_cell, + ExcInvalidNeighbor(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 neighbor_guess + = GeometryInfo::opposite_face[neighbor]; + + if ((neighbor_cell->neighbor_index (neighbor_guess) == present_index) && + (neighbor_cell->neighbor_level (neighbor_guess) == present_level)) + return neighbor_guess; + else + // if the guess was false, then + // we need to loop over all + // neighbors and find the number + // the hard way + { + for (unsigned int face=0; face::faces_per_cell; ++face) + if ((neighbor_cell->neighbor_index (face) == present_index) && + (neighbor_cell->neighbor_level (face) == present_level)) + return face; + + // we should never get here, + // since then we did not find + // our way back... + Assert (false, ExcInternalError()); + return static_cast(-1); + }; +}; + + + template bool CellAccessor::at_boundary (const unsigned int i) const { Assert (used(), ExcCellNotUsed()); diff --git a/deal.II/deal.II/source/numerics/error_estimator.cc b/deal.II/deal.II/source/numerics/error_estimator.cc index f4ca76d363..7f08d1532f 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -578,18 +578,14 @@ integrate_over_regular_face (Data &data, // of gradient across this face { Assert (cell->neighbor(face_no).state() == valid, - ExcInternalError()); - unsigned int neighbor_neighbor; - DoFHandler::active_cell_iterator neighbor = cell->neighbor(face_no); + ExcInternalError()); + + const DoFHandler::active_cell_iterator neighbor = cell->neighbor(face_no); // find which number the current // face has relative to the neighboring // cell - for (neighbor_neighbor=0; neighbor_neighbor::faces_per_cell; - ++neighbor_neighbor) - if (neighbor->neighbor(neighbor_neighbor) == cell) - break; - + const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor (face_no); Assert (neighbor_neighbor::faces_per_cell, ExcInternalError()); // get restriction of finite element @@ -720,6 +716,7 @@ integrate_over_irregular_face (Data &data, FESubfaceValues &fe_subface_values) { const DoFHandler::cell_iterator neighbor = cell->neighbor(face_no); + Assert (neighbor.state() == valid, ExcInternalError()); Assert (neighbor->has_children(), ExcInternalError()); // set up a vector of the gradients @@ -736,11 +733,7 @@ integrate_over_irregular_face (Data &data, // store which number #cell# has in the // list of neighbors of #neighbor# - unsigned int neighbor_neighbor; - for (neighbor_neighbor=0; neighbor_neighbor::faces_per_cell; - ++neighbor_neighbor) - if (neighbor->neighbor(neighbor_neighbor) == cell) - break; + const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor (face_no); Assert (neighbor_neighbor::faces_per_cell, ExcInternalError()); // loop over all subfaces -- 2.39.5