* @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
*/
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<CellAccessor<dim, spacedim> >
+ 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<CellAccessor<dim, spacedim> >
+ 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<unsigned int, unsigned int>
+ 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;
+
/**
* @}
*/
-// ---------------------------------------------------------------------
+// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2015 by the deal.II authors
//
+template <int dim, int spacedim>
+bool
+CellAccessor<dim, spacedim>::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<dim>::faces_per_cell);
+ typedef TriaIterator<CellAccessor<dim, spacedim> > 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<cell_iterator, unsigned int>(my_it, i_face))
+ != this->tria->periodic_face_map.end())
+ return true;
+ return false;
+}
+
+
+
+template <int dim, int spacedim>
+TriaIterator<CellAccessor<dim,spacedim> >
+CellAccessor<dim, spacedim>::
+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<dim>::faces_per_cell);
+ typedef TriaIterator<CellAccessor<dim, spacedim> > cell_iterator;
+ typedef std::pair<cell_iterator, unsigned int> cell_face_pair;
+ typedef std::pair<cell_face_pair, std::bitset<3> > oriented_cell_face_pair;
+ cell_iterator my_it(*this);
+ typename std::map<const cell_face_pair, oriented_cell_face_pair>::const_iterator my_face_pair =
+ this->tria->periodic_face_map.find(std::pair<cell_iterator, unsigned int>(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 <int dim, int spacedim>
+TriaIterator<CellAccessor<dim, spacedim> >
+CellAccessor<dim, spacedim>::
+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<dim>::faces_per_cell);
+ typedef TriaIterator<CellAccessor<dim, spacedim> > cell_iterator;
+ typedef std::pair<cell_iterator, unsigned int> cell_face_pair;
+ typedef std::pair<cell_face_pair, std::bitset<3> > oriented_cell_face_pair;
+ cell_iterator my_it(*this);
+ typename std::map<const cell_face_pair, oriented_cell_face_pair>::const_iterator my_face_pair =
+ this->tria->periodic_face_map.find(std::pair<cell_iterator, unsigned int>(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<TriaAccessor<dim - 1, dim, spacedim> > 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<dim>::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 <int dim, int spacedim>
+std::pair<unsigned int, unsigned int>
+CellAccessor<dim, spacedim>::
+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<dim>::faces_per_cell);
+ typedef TriaIterator<CellAccessor<dim, spacedim> > cell_iterator;
+ typedef std::pair<cell_iterator, unsigned int> cell_face_pair;
+ typedef std::pair<cell_face_pair, std::bitset<3> > oriented_cell_face_pair;
+ const int my_face_index = this->face_index(i_face);
+ cell_iterator my_it(*this);
+ typename std::map<const cell_face_pair, oriented_cell_face_pair>::const_iterator my_face_pair =
+ this->tria->periodic_face_map.find(std::pair<cell_iterator, unsigned int>(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 cell_face_pair, oriented_cell_face_pair>::const_iterator nb_face_pair =
+ this->tria->periodic_face_map.find(std::pair<cell_iterator, unsigned int>(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<TriaAccessor<dim - 1, dim, spacedim> > 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<unsigned int, unsigned int>(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<unsigned int, unsigned int>(numbers::invalid_unsigned_int,
+ numbers::invalid_unsigned_int);
+}
+
+
+
+template <int dim, int spacedim>
+int CellAccessor<dim, spacedim>::
+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<dim>::faces_per_cell);
+ typedef TriaIterator<CellAccessor<dim, spacedim> > cell_iterator;
+ typedef std::pair<cell_iterator, unsigned int> cell_face_pair;
+ typedef std::pair<cell_face_pair, std::bitset<3> > oriented_cell_face_pair;
+ cell_iterator my_it(*this);
+ typename std::map<const cell_face_pair, oriented_cell_face_pair>::const_iterator my_face_pair =
+ this->tria->periodic_face_map.find(std::pair<cell_iterator, unsigned int>(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 dim, int spacedim>
+int CellAccessor<dim, spacedim>::
+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<dim>::faces_per_cell);
+ typedef TriaIterator<CellAccessor<dim, spacedim> > cell_iterator;
+ typedef std::pair<cell_iterator, unsigned int> cell_face_pair;
+ typedef std::pair<cell_face_pair, std::bitset<3> > oriented_cell_face_pair;
+ cell_iterator my_it(*this);
+ typename std::map<const cell_face_pair, oriented_cell_face_pair>::const_iterator my_face_pair =
+ this->tria->periodic_face_map.find(std::pair<cell_iterator, unsigned int>(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 <int dim, int spacedim>
+unsigned int CellAccessor<dim, spacedim>::
+periodic_neighbor_of_periodic_neighbor(const unsigned int i_face) const
+{
+ return periodic_neighbor_face_no(i_face);
+}
+
+
+
+template <int dim, int spacedim>
+unsigned int
+CellAccessor<dim, spacedim>::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<dim>::faces_per_cell);
+ typedef TriaIterator<CellAccessor<dim, spacedim> > cell_iterator;
+ typedef std::pair<cell_iterator, unsigned int> cell_face_pair;
+ typedef std::pair<cell_face_pair, std::bitset<3> > oriented_cell_face_pair;
+ cell_iterator my_it(*this);
+ typename std::map<const cell_face_pair, oriented_cell_face_pair>::const_iterator my_face_pair =
+ this->tria->periodic_face_map.find(std::pair<cell_iterator, unsigned int>(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 <int dim, int spacedim>
+bool
+CellAccessor<dim, spacedim>::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<dim>::faces_per_cell);
+ typedef TriaIterator<CellAccessor<dim, spacedim> > cell_iterator;
+ typedef std::pair<cell_iterator, unsigned int> cell_face_pair;
+ typedef std::pair<cell_face_pair, std::bitset<3> > oriented_cell_face_pair;
+ cell_iterator my_it(*this);
+ typename std::map<const cell_face_pair, oriented_cell_face_pair>::const_iterator my_face_pair =
+ this->tria->periodic_face_map.find(std::pair<cell_iterator, unsigned int>(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 cell_face_pair, oriented_cell_face_pair>::const_iterator nb_face_pair =
+ this->tria->periodic_face_map.find(std::pair<cell_iterator, unsigned int>(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 <int dim, int spacedim>
bool CellAccessor<dim, spacedim>::at_boundary (const unsigned int i) const
{