]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement periodic neighbor functions in TriaAccessor
authorDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Tue, 22 Mar 2016 14:59:16 +0000 (15:59 +0100)
committerDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Tue, 29 Mar 2016 11:24:09 +0000 (13:24 +0200)
include/deal.II/grid/tria_accessor.h
source/grid/tria_accessor.cc

index da465338c75c9f04675ff48e41b6d29316187188..96c18fb2212a73abfe4a80702a8e91708262360a 100644 (file)
@@ -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<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;
+
   /**
    * @}
    */
index c7921ab84d90454d6a8c7fb593cbc900816ea47c..cb3abc3d0c274890d54e45ff8390a9eac04dd119 100644 (file)
@@ -1,4 +1,4 @@
-// ---------------------------------------------------------------------
+// ---------------------------------------------------------------------
 //
 // Copyright (C) 1998 - 2015 by the deal.II authors
 //
@@ -1868,6 +1868,325 @@ CellAccessor<dim, spacedim>::neighbor_of_coarser_neighbor (const unsigned int ne
 
 
 
+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
 {

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.