From c3108d37c619fe52e0113a70cd7a285acaa37d15 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Tue, 22 Mar 2016 15:59:16 +0100 Subject: [PATCH] Implement periodic neighbor functions in TriaAccessor --- include/deal.II/grid/tria_accessor.h | 117 ++++++++++ source/grid/tria_accessor.cc | 321 ++++++++++++++++++++++++++- 2 files changed, 437 insertions(+), 1 deletion(-) diff --git a/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h index da465338c7..96c18fb221 100644 --- a/include/deal.II/grid/tria_accessor.h +++ b/include/deal.II/grid/tria_accessor.h @@ -204,6 +204,13 @@ namespace TriaAccessorExceptions * @ingroup Exceptions */ DeclException0 (ExcFacesHaveNoLevel); + /** + * You are trying to get the periodic neighbor for a face, which does not + * have a periodic neighbor. For more information on this, refer to + * @ref GlossPeriodicConstraints "entry for periodic boundaries". + * @ingroup Exceptions + */ + DeclException0 (ExcNoPeriodicNeighbor); //TODO: Write documentation! /** * @ingroup Exceptions @@ -2486,6 +2493,116 @@ public: */ unsigned int neighbor_face_no (const unsigned int neighbor) const; + /** + * @} + */ + /** + * @name Dealing with periodic neighbor + */ + /** + * @{ + */ + /** + * If the cell has a periodic neighbor at its @c ith face, this function + * returns true, otherwise, the returned value is false. + */ + bool has_periodic_neighbor(const unsigned int i) const; + + /** + * For a cell with its @c ith face at a periodic boundary, + * (see @ref GlossPeriodicConstraints "the entry for periodic baoundaries") + * this function returns an iterator to the cell on the other side + * of the periodic boundary. If there is no periodic boundary at the @c ith + * face, an exception will be thrown. + * In order to check if a cell has periodic neighbor on its @c ith face, + * one should first call the function @c has_periodic_neighbor(). + * The behavior of @c periodic_neighbor() is similar to @c neighbor(), in + * the sense that the returned cell has at most the same level of refinement + * as the current cell. On distributed meshes, by calling + * Triangulation::add_periodicity(), + * we can make sure that the element on the other side of the periodic + * boundary exists in this rank as a ghost cell or a locally owned cell. + */ + TriaIterator > + periodic_neighbor (const unsigned int i) const; + + /** + * Returns an iterator to the periodic neighbor of the cell at a given + * face and subface number. The general guidelines for using this function + * is similar to the function @c neighbor_child_on_subface. The + * implementation of this function is consistent with + * @c periodic_neighbor_of_coarser_periodic_neighbor. For instance, + * assume that we are sitting on a cell named @c cell1, which has a 1 level + * coarser neighbor at its @c ith face. Let us name this coarser neighbor + * @c cell2. Then, by calling + * @c periodic_neighbor_of_coarser_periodic_neighbor, from @c cell1, we get + * a @c face_num and a @c subface_num. Now, if we call + * @c periodic_neighbor_child_on_subface from cell2, with the above face_num + * and subface_num, we get an iterator to @c cell1. + */ + TriaIterator > + periodic_neighbor_child_on_subface (const unsigned int face_no, + const unsigned int subface_no) const; + + /** + * This function is a generalization of + * @c periodic_neighbor_of_periodic_neighbor + * for those cells which have a coarser periodic neighbor. The returned + * pair of numbers can be used in @c periodic_neighbor_child_on_subface + * to get back to the current cell. In other words, the following + * assertion should be true, for a cell with coarser periodic neighbor: + * cell->periodic_neighbor(i)->periodic_neighbor_child_on_subface(face_no, subface_no)==cell + */ + std::pair + periodic_neighbor_of_coarser_periodic_neighbor (const unsigned i) const; + + /** + * This function returns the index of the periodic neighbor. If there is + * no periodic neighbor at the given face, the returned value is -1. + */ + int + periodic_neighbor_index (const unsigned int i) const; + + /** + * This function returns the level of the periodic neighbor. If there is + * no periodic neighbor at the given face, the returned value is -1. + */ + int + periodic_neighbor_level (const unsigned int i) const; + + /** + * For a cell with a periodic neighbor at its @c ith face, this function + * returns the face number of that periodic neighbor such that, the + * current cell is the periodic neighbor of that neighbor. In other words + * the following assertion holds for those cells which have a periodic + * neighbor with the same or a higher level of refinement as the current + * cell: + * @c {cell->periodic_neighbor(i)-> + * periodic_neighbor(cell->periodic_neighbor_of_periodic_neighbor(i))==cell} + * For the cells with a coarser periodic neighbor, one should use + * @c periodic_neighbor_of_coarser_periodic_neighbor and + * @c periodic_neighbor_child_on_subface + * to get back to the current cell. + */ + unsigned int + periodic_neighbor_of_periodic_neighbor (const unsigned int i) const; + + /** + * If a cell has periodic neighbor, this function returns the face number + * of the neighbor, which is connected to this cell. + */ + unsigned int + periodic_neighbor_face_no (const unsigned int i) const; + + /** + * This function returns true if the element on the other side of the + * periodic boundary is coarser and returns false, otherwise. The + * implementation allows this function to work, in the case of + * anisotropic refinement. + */ + bool + periodic_neighbor_is_coarser (const unsigned int i) const; + /** * @} */ diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index c7921ab84d..cb3abc3d0c 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -1,4 +1,4 @@ -// --------------------------------------------------------------------- +// --------------------------------------------------------------------- // // Copyright (C) 1998 - 2015 by the deal.II authors // @@ -1868,6 +1868,325 @@ CellAccessor::neighbor_of_coarser_neighbor (const unsigned int ne +template +bool +CellAccessor::has_periodic_neighbor (const unsigned int i_face) const +{ + /* + * Implementation note: In all of the functions corresponding to periodic faces + * we mainly use the Triangulation::periodic_face_map to find the + * information about periodically connected faces. So, we actually search in + * this std::map and return the cell_face on the other side of periodic boundary. + * For this search process, we have three options: + * + * 1- Using the [] operator of std::map: This option results in a more readalbe + * code, but requires an extra iteration in the map. Becasue when we call [] on + * std::map, with a key which does not exist in the std::map, that key will be + * created and the default value will be returned by []. This is not desirable. + * So, one has to first check if the key exists in the std::map and if it exists, + * then use the [] operator. The existence check is possible using std::map::find() + * or std::map::count(). Using this option will result in two iteration cycles + * through the map. First, existence check, then returning the value. + * + * 2- Using std::map::find(): This option is less readble, but theoretically + * faster. Because, it results in one iteration through std::map object. + * + * We decided to use the 2nd option. + */ + AssertIndexRange (i_face, GeometryInfo::faces_per_cell); + typedef TriaIterator > cell_iterator; + // my_it : is the iterator to the current cell. + cell_iterator my_it(*this); + if (this->tria->periodic_face_map.find(std::pair(my_it, i_face)) + != this->tria->periodic_face_map.end()) + return true; + return false; +} + + + +template +TriaIterator > +CellAccessor:: +periodic_neighbor (const unsigned int i_face) const +{ + /* + * To know, why we are using std::map::find() instead of [] operator, refer + * to the implementation note in has_periodic_neighbor() function. + * + * my_it : the iterator to the current cell. + * my_face_pair : the pair reported by periodic_face_map as its first pair being + * the current cell_face. + */ + AssertIndexRange (i_face, GeometryInfo::faces_per_cell); + typedef TriaIterator > cell_iterator; + typedef std::pair cell_face_pair; + typedef std::pair > oriented_cell_face_pair; + cell_iterator my_it(*this); + typename std::map::const_iterator my_face_pair = + this->tria->periodic_face_map.find(std::pair(my_it, i_face)); + // Assertion is required to check that we are actually on a periodic boundary. + Assert (my_face_pair != this->tria->periodic_face_map.end(), + TriaAccessorExceptions::ExcNoPeriodicNeighbor()); + return my_face_pair->second.first.first; +} + + + +template +TriaIterator > +CellAccessor:: +periodic_neighbor_child_on_subface (const unsigned int i_face, + const unsigned int i_subface) const +{ + /* + * To know, why we are using std::map::find() instead of [] operator, refer + * to the implementation note in has_periodic_neighbor() function. + * + * my_it : the iterator to the current cell. + * my_face_pair : the pair reported by periodic_face_map as its first pair being + * the current cell_face. + * nb_it : the iterator to the neighbor of current cell at i_face. + * face_num_of_nb : the face number of the periodically neighboring face in the + * relevant element. + * nb_parent_face_it: the iterator to the parent face of the periodically + * neighboring face. + */ + AssertIndexRange (i_face, GeometryInfo::faces_per_cell); + typedef TriaIterator > cell_iterator; + typedef std::pair cell_face_pair; + typedef std::pair > oriented_cell_face_pair; + cell_iterator my_it(*this); + typename std::map::const_iterator my_face_pair = + this->tria->periodic_face_map.find(std::pair(my_it, i_face)); + /* + * There should be an assertion, which tells the user that this function should not be + * used for a cell which is not located at a periodic boundary. + */ + Assert (my_face_pair != this->tria->periodic_face_map.end(), + TriaAccessorExceptions::ExcNoPeriodicNeighbor()); + cell_iterator parent_nb_it = my_face_pair->second.first.first; + unsigned int nb_face_num = my_face_pair->second.first.second; + TriaIterator > nb_parent_face_it = parent_nb_it->face(nb_face_num); + /* + * We should check if the parent face of the neighbor has at least the same number of + * children as i_subface. + */ + AssertIndexRange (i_subface, nb_parent_face_it->n_children()); + unsigned int sub_neighbor_num = + GeometryInfo::child_cell_on_face(parent_nb_it->refinement_case(), + nb_face_num, + i_subface, + my_face_pair->second.second[0], + my_face_pair->second.second[1], + my_face_pair->second.second[2], + nb_parent_face_it->refinement_case()); + return parent_nb_it->child(sub_neighbor_num); +} + + + +template +std::pair +CellAccessor:: +periodic_neighbor_of_coarser_periodic_neighbor (const unsigned int i_face) const +{ + /* + * To know, why we are using std::map::find() instead of [] operator, refer + * to the implementation note in has_periodic_neighbor() function. + * + * my_it : the iterator to the current cell. + * my_face_pair : the pair reported by periodic_face_map as its first pair being + * the current cell_face. + * nb_it : the iterator to the periodic neighbor. + * nb_face_pair : the pair reported by periodic_face_map as its first pair being + * the periodic neighbor cell_face. + * p_nb_of_p_nb : the iterator of the periodic neighbor of the periodic neighbor + * of the current cell. + */ + AssertIndexRange (i_face, GeometryInfo::faces_per_cell); + typedef TriaIterator > cell_iterator; + typedef std::pair cell_face_pair; + typedef std::pair > oriented_cell_face_pair; + const int my_face_index = this->face_index(i_face); + cell_iterator my_it(*this); + typename std::map::const_iterator my_face_pair = + this->tria->periodic_face_map.find(std::pair(my_it, i_face)); + /* + * There should be an assertion, which tells the user that this function should not be + * used for a cell which is not located at a periodic boundary. + */ + Assert (my_face_pair != this->tria->periodic_face_map.end(), + TriaAccessorExceptions::ExcNoPeriodicNeighbor()); + cell_iterator nb_it = my_face_pair->second.first.first; + unsigned int face_num_of_nb = my_face_pair->second.first.second; + typename std::map::const_iterator nb_face_pair = + this->tria->periodic_face_map.find(std::pair(nb_it, face_num_of_nb)); + /* + * Since, we store periodic neighbors for every cell (either active or artificial or inactive) + * the nb_face_pair should also be mapped to some cell_face pair. We assert this here. + */ + Assert (nb_face_pair != this->tria->periodic_face_map.end(), + TriaAccessorExceptions::ExcNoPeriodicNeighbor()); + cell_iterator p_nb_of_p_nb = nb_face_pair->second.first.first; + TriaIterator > parent_face_it = p_nb_of_p_nb->face(nb_face_pair->second.first.second); + for (unsigned int i_subface = 0; i_subface < parent_face_it->n_children(); ++i_subface) + if (parent_face_it->child_index(i_subface) == my_face_index) + return (std::pair(face_num_of_nb, i_subface)); + /* + * Obviously, if the execution reaches to this point, some of our assumptions should have + * been false. The most important one is, the user has called this funciton on a face + * which does not have a coarser periodic neighbor. + */ + Assert (false, TriaAccessorExceptions::ExcNeighborIsNotCoarser()); + return std::pair(numbers::invalid_unsigned_int, + numbers::invalid_unsigned_int); +} + + + +template +int CellAccessor:: +periodic_neighbor_index(const unsigned int i_face) const +{ + /* + * To know, why we are using std::map::find() instead of [] operator, refer + * to the implementation note in has_periodic_neighbor() function. + * + * my_it : the iterator to the current cell. + * my_face_pair : the pair reported by periodic_face_map as its first pair being + * the current cell_face. + */ + AssertIndexRange (i_face, GeometryInfo::faces_per_cell); + typedef TriaIterator > cell_iterator; + typedef std::pair cell_face_pair; + typedef std::pair > oriented_cell_face_pair; + cell_iterator my_it(*this); + typename std::map::const_iterator my_face_pair = + this->tria->periodic_face_map.find(std::pair(my_it, i_face)); + Assert (my_face_pair != this->tria->periodic_face_map.end(), + TriaAccessorExceptions::ExcNoPeriodicNeighbor()); + return my_face_pair->second.first.first->index(); +} + + + +template +int CellAccessor:: +periodic_neighbor_level(const unsigned int i_face) const +{ + /* + * To know, why we are using std::map::find() instead of [] operator, refer + * to the implementation note in has_periodic_neighbor() function. + * + * my_it : the iterator to the current cell. + * my_face_pair : the pair reported by periodic_face_map as its first pair being + * the current cell_face. + */ + AssertIndexRange (i_face, GeometryInfo::faces_per_cell); + typedef TriaIterator > cell_iterator; + typedef std::pair cell_face_pair; + typedef std::pair > oriented_cell_face_pair; + cell_iterator my_it(*this); + typename std::map::const_iterator my_face_pair = + this->tria->periodic_face_map.find(std::pair(my_it, i_face)); + Assert (my_face_pair != this->tria->periodic_face_map.end(), + TriaAccessorExceptions::ExcNoPeriodicNeighbor()); + return my_face_pair->second.first.first->level(); +} + + + +template +unsigned int CellAccessor:: +periodic_neighbor_of_periodic_neighbor(const unsigned int i_face) const +{ + return periodic_neighbor_face_no(i_face); +} + + + +template +unsigned int +CellAccessor::periodic_neighbor_face_no (const unsigned int i_face) const +{ + /* + * To know, why we are using std::map::find() instead of [] operator, refer + * to the implementation note in has_periodic_neighbor() function. + * + * my_it : the iterator to the current cell. + * my_face_pair : the pair reported by periodic_face_map as its first pair being + * the current cell_face. + */ + AssertIndexRange (i_face, GeometryInfo::faces_per_cell); + typedef TriaIterator > cell_iterator; + typedef std::pair cell_face_pair; + typedef std::pair > oriented_cell_face_pair; + cell_iterator my_it(*this); + typename std::map::const_iterator my_face_pair = + this->tria->periodic_face_map.find(std::pair(my_it, i_face)); + /* + * There should be an assertion, which tells the user that this function should not be + * called for a cell which is not located at a periodic boundary ! + */ + Assert (my_face_pair != this->tria->periodic_face_map.end(), + TriaAccessorExceptions::ExcNoPeriodicNeighbor()); + return my_face_pair->second.first.second; +} + + + +template +bool +CellAccessor::periodic_neighbor_is_coarser (const unsigned int i_face) const +{ + /* + * To know, why we are using std::map::find() instead of [] operator, refer + * to the implementation note in has_periodic_neighbor() function. + * + * Implementation note: Let p_nb_of_p_nb be the periodic neighbor of the periodic + * neighbor of the current cell. Also, let p_face_of_p_nb_of_p_nb be the periodic + * face of the p_nb_of_p_nb. If p_face_of_p_nb_of_p_nb has children , then the + * periodic neighbor of the current cell is coarser than itself. Although not tested, + * this implementation should work for anisotropic refinement as well. + * + * my_it : the iterator to the current cell. + * my_face_pair : the pair reported by periodic_face_map as its first pair being + * the current cell_face. + * nb_it : the iterator to the periodic neighbor. + * nb_face_pair : the pair reported by periodic_face_map as its first pair being + * the periodic neighbor cell_face. + */ + AssertIndexRange (i_face, GeometryInfo::faces_per_cell); + typedef TriaIterator > cell_iterator; + typedef std::pair cell_face_pair; + typedef std::pair > oriented_cell_face_pair; + cell_iterator my_it(*this); + typename std::map::const_iterator my_face_pair = + this->tria->periodic_face_map.find(std::pair(my_it, i_face)); + /* + * There should be an assertion, which tells the user that this function should not be + * used for a cell which is not located at a periodic boundary. + */ + Assert (my_face_pair != this->tria->periodic_face_map.end(), + TriaAccessorExceptions::ExcNoPeriodicNeighbor()); + cell_iterator nb_it = my_face_pair->second.first.first; + unsigned int face_num_of_nb = my_face_pair->second.first.second; + typename std::map::const_iterator nb_face_pair = + this->tria->periodic_face_map.find(std::pair(nb_it, face_num_of_nb)); + /* + * Since, we store periodic neighbors for every cell (either active or artificial or inactive) + * the nb_face_pair should also be mapped to some cell_face pair. We assert this here. + */ + Assert (nb_face_pair != this->tria->periodic_face_map.end(), + TriaAccessorExceptions::ExcNoPeriodicNeighbor()); + if (nb_face_pair->second.first.first->face(nb_face_pair->second.first.second)->has_children()) + return true; + return false; +} + + + template bool CellAccessor::at_boundary (const unsigned int i) const { -- 2.39.5