]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add a helper function GridTools::collect_periodic_cell_pairs for periodic
authormaier <maier@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 21 May 2012 22:21:25 +0000 (22:21 +0000)
committermaier <maier@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 21 May 2012 22:21:25 +0000 (22:21 +0000)
boundary conditions

git-svn-id: https://svn.dealii.org/trunk@25529 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/grid/grid_tools.h
deal.II/source/grid/grid_tools.cc
deal.II/source/grid/grid_tools.inst.in

index 9358faf344b44b996216fce8b264e3c825fa5572..2bd4b431efea7855be842c84a70deb0638174a13 100644 (file)
@@ -901,6 +901,74 @@ namespace GridTools
                          const std::set<types::boundary_id_t> &boundary_ids
                          = std::set<types::boundary_id_t>());
 
+
+                                   /**
+                                    * This function loops over all cells
+                                    * from @p begin to @p end and collects a
+                                    * set of periodic cell pairs for
+                                    * @p direction:
+                                    *
+                                    * Given a @p direction,
+                                    * define a 'left' boundary as all
+                                    * boundary faces belonging to
+                                    * @p boundary_component with corresponding
+                                    * unit normal (of the @ref
+                                    * GlossReferenceCell "reference cell") in
+                                    * negative @p direction, i.e. all boundary
+                                    * faces with
+                                    * <tt>face(2*direction)->at_boundary()==true</tt>.
+                                    * Analogously, a 'right' boundary
+                                    * consisting of all faces of @p
+                                    * boundary_component with unit normal
+                                    * in positive @p direction,
+                                    * <tt>face(2*direction+1)->at_boundary()==true</tt>.
+                                    *
+                                    * This function tries to match
+                                    * all cells with faces belonging to the
+                                    * left' boundary with the cells with faces
+                                    * belonging to the 'right' boundary
+                                    * with the help of an
+                                    * orthogonal equality relation:
+                                    * Two cells do match if the vertices of
+                                    * the respective boundary faces
+                                    * can be transformed into each other by
+                                    * parallel translation in @p direction.
+                                    *
+                                    * The @p offset is a vector tangential to
+                                    * the faces that is added to the location
+                                    * of vertices of the 'left' boundary when
+                                    * attempting to match them to the
+                                    * corresponding vertices of the 'right'
+                                    * boundary. This can be used to implement
+                                    * conditions such as $u(0,y)=u(1,y+1)$.
+                                    */
+  template<typename CellIterator>
+  std::map<CellIterator, CellIterator>
+  collect_periodic_cell_pairs (const CellIterator                          &begin,
+                               const typename identity<CellIterator>::type &end,
+                               const types::boundary_id_t                  boundary_component,
+                               int                                         direction,
+                               const dealii::Tensor<1,CellIterator::AccessorType::space_dimension>
+                                 &offset = dealii::Tensor<1,CellIterator::AccessorType::space_dimension>());
+
+
+                                   /**
+                                    * Same as above but matches all active
+                                    * cells, i.e. this function calls the
+                                    * above function with
+                                    * @p dof_handler.begin_active() and
+                                    * @p dof_handler.end() as the first two
+                                    * arguments.
+                                    */
+  template<typename DH>
+  std::map<typename DH::cell_iterator, typename DH::cell_iterator>
+  collect_periodic_cell_pairs (const DH                   &dof_handler,
+                               const types::boundary_id_t boundary_component,
+                               int                        direction,
+                               const dealii::Tensor<1,DH::space_dimension>
+                                 &offset = Tensor<1,DH::space_dimension>());
+
+
                                  /**
                                   * Exception
                                   */
index 8c6cea3de26ade3ae70fb1622945ea6643ba7267..25ffb101e44f001c378dc3eec656a280ec678c78 100644 (file)
@@ -2009,7 +2009,139 @@ namespace GridTools
     return surface_to_volume_mapping;
   }
 
-}
+                                   // Internally used in
+                                   // collect_periodic_cell_pairs.
+  template<typename CellIterator>
+  inline bool orthogonal_equality (const CellIterator &cell1,
+                                   const CellIterator &cell2,
+                                   const int          direction,
+                                   const dealii::Tensor<1,CellIterator::AccessorType::space_dimension>
+                                                      &offset)
+  {
+                                   // An orthogonal equality test for
+                                   // cells:
+                                   // Two cells are equal if the
+                                   // corresponding vertices of the
+                                   // boundary faces can be transformed
+                                   // into each other parallel to the
+                                   // given direction.
+                                   // It is assumed that
+                                   // cell1->face(2*direction) and
+                                   // cell2->face(2*direction+1) are the
+                                   // corresponding boundary faces.
+                                   // The optional argument offset will
+                                   // be added to each vertex of cell1
+
+    static const int dim = CellIterator::AccessorType::dimension;
+    static const int space_dim = CellIterator::AccessorType::space_dimension;
+    Assert(dim == space_dim,
+           ExcNotImplemented());
+                                   // ... otherwise we would need two
+                                   // directions: One for selecting the
+                                   // faces, thus living in dim,
+                                   // the other direction for selecting the
+                                   // direction for the orthogonal
+                                   // projection, thus living in
+                                   // space_dim.
+
+    for (int i = 0; i < space_dim; ++i) {
+                                   // Only compare
+                                   // coordinate-components != direction:
+      if (i == direction)
+        continue;
+
+      for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
+        if (fabs (  cell1->face(2*direction)->vertex(v)(i)
+                  + offset[i]
+                  - cell2->face(2*direction+1)->vertex(v)(i))
+                    > 1.e-10)
+          return false;
+    }
+    return true;
+  }
+
+
+  template<typename CellIterator>
+  std::map<CellIterator, CellIterator>
+  collect_periodic_cell_pairs (const CellIterator                          &begin,
+                               const typename identity<CellIterator>::type &end,
+                               const types::boundary_id_t                  boundary_component,
+                               int                                         direction,
+                               const dealii::Tensor<1,CellIterator::AccessorType::space_dimension>
+                                 &offset = dealii::Tensor<1,CellIterator::AccessorType::space_dimension>())
+  {
+    static const int space_dim = CellIterator::AccessorType::space_dimension;
+    Assert (0<=direction && direction<space_dim,
+            ExcIndexRange (direction, 0, space_dim));
+
+    std::set<CellIterator> cells1;
+    std::set<CellIterator> cells2;
+
+                                   // collect boundary cells, i.e. cells
+                                   // with face(2*direction), resp.
+                                   // face(2*direction+1) being on the
+                                   // boundary and having the correct
+                                   // color:
+    CellIterator cell = begin;
+    for (; cell!= end; ++cell) {
+      if (cell->face(2*direction)->at_boundary() &&
+          cell->face(2*direction)->boundary_indicator() == boundary_component) {
+        cells1.insert(cell);
+      }
+      if (cell->face(2*direction+1)->at_boundary() &&
+          cell->face(2*direction+1)->boundary_indicator() == boundary_component) {
+        cells2.insert(cell);
+      }
+    }
+
+    Assert (cells1.size() == cells2.size(),
+            ExcMessage ("Unmatched cells on periodic boundaries"));
+
+    std::map<CellIterator, CellIterator> matched_cells;
+
+                                   // Match with a complexity of O(n^2).
+                                   // This could be improved...
+    typedef typename std::set<CellIterator>::const_iterator SetIterator;
+    for (SetIterator it1 = cells1.begin(); it1 != cells1.end(); ++it1) {
+      for (SetIterator it2 = cells2.begin(); it2 != cells2.end(); ++it2) {
+        if (orthogonal_equality(*it1, *it2, direction, offset)) {
+                                   // We have a match, so insert the
+                                   // matching pairs and remove the matched
+                                   // cell in cells2 to speed up the
+                                   // matching:
+          matched_cells[*it1] = *it2;
+          cells2.erase(it2);
+          break;
+        }
+      }
+    }
+
+    Assert (matched_cells.size() == cells1.size() &&
+            cells2.size() == 0,
+            ExcMessage ("Unmatched cells on periodic boundaries"));
+
+    return matched_cells;
+  }
+
+
+  template<typename DH>
+  std::map<typename DH::cell_iterator, typename DH::cell_iterator>
+  collect_periodic_cell_pairs (const DH                   &dof_handler,
+                               const types::boundary_id_t boundary_component,
+                               int                        direction,
+                               const dealii::Tensor<1,DH::space_dimension>
+                                 &offset = Tensor<1,DH::space_dimension>())
+  {
+    return collect_periodic_cell_pairs<typename DH::cell_iterator>
+                                      (dof_handler.begin_active(),
+                                       dof_handler.end(),
+                                       boundary_component,
+                                       direction,
+                                       offset);
+  }
+
+
+} /* namespace GridTools */
 
 
 // explicit instantiations
index e9da14bd40f5440bafdabd77c29d71a82122af3d..63a370af004be7f17fb81d639c72106f52e2bd5c 100644 (file)
@@ -176,3 +176,30 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
 #endif
 
   }
+
+
+for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS)
+{
+#if deal_II_dimension <= deal_II_space_dimension
+  namespace GridTools \{
+
+    template
+    std::map<X::cell_iterator, X::cell_iterator>
+    collect_periodic_cell_pairs(const X::cell_iterator &,
+                                const X::cell_iterator &,
+                                const types::boundary_id_t,
+                                int,
+                                const Tensor<1,deal_II_space_dimension> &);
+
+    template
+    std::map<X::cell_iterator, X::cell_iterator>
+    collect_periodic_cell_pairs (const X&,
+                                 const types::boundary_id_t,
+                                 int,
+                                 const Tensor<1,deal_II_space_dimension> &);
+
+
+  \}
+  #endif
+}
+

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.