Add get_face_iterator

Conflicts:
	include/deal.II/matrix_free/matrix_free.templates.h

diff --git a/doc/news/changes/minor/20191016MartinKronbichlerPeterMunchMichalWichrowski b/doc/news/changes/minor/20191016MartinKronbichlerPeterMunchMichalWichrowski
new file mode 100644 (file)
index 0000000..8c3fc31
--- /dev/null
+++ b/doc/news/changes/minor/20191016MartinKronbichlerPeterMunchMichalWichrowski
@@ -0,0 +1,4 @@
+New: Provide cell iterators to interior and exterior cells of
+a macro face through the function MatrixFree::get_face_iterator().
+<br>
+(Martin Kronbichler, Peter Munch, Michal Wichrowski, 2019/10/16)

diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h
index e70bdfbaaad3ecb140954e7fd37d37629587afa7..6377fa2ee9779774ad1a549e2965a691ffabe6a8 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.h
+++ b/include/deal.II/matrix_free/matrix_free.h
@@ -1464,6 +1464,24 @@ public:
   get_cell_level_and_index(const unsigned int macro_cell_number,
                            const unsigned int vector_number) const;
 
+  /**
+   * Return the cell iterator in deal.II speak to a interior/exterior cell of
+   * given face  in the renumbering of this structure. The second element
+   * of the pair is the face number so that the face iterator can be accessed:
+   * pair.first()->face(pair.second() );
+   *
+   * Note that the face iterators in deal.II go through cells differently to
+   * what the face/boundary loop of this class does. This is because several
+   * faces are worked on together (vectorization), and since faces with neighbor
+   * cells on different MPI processors need to be accessed at a certain time
+   * when accessing remote data and overlapping communication with computation.
+   */
+  std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
+  get_face_iterator(const unsigned int face_batch_number,
+                    const unsigned int vector_number,
+                    const bool         interior     = true,
+                    const unsigned int fe_component = 0) const;
+
   /**
    * This returns the cell iterator in deal.II speak to a given cell in the
    * renumbering of this structure. This function returns an exception in case

diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
index 75ab6f44ab085d6b3b2d63f486e288f999e676b5..387594387b72174488f7d9b64e51ae8f0584abb8 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -206,6 +206,59 @@ MatrixFree<dim, Number, VectorizedArrayType>::get_cell_level_and_index(
 
 
 
+template <int dim, typename Number, typename VectorizedArrayType>
+std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
+MatrixFree<dim, Number, VectorizedArrayType>::get_face_iterator(
+  const unsigned int face_batch_number,
+  const unsigned int vector_number,
+  const bool         interior,
+  const unsigned int fe_component) const
+{
+  AssertIndexRange(fe_component, dof_handlers.n_dof_handlers);
+  if (interior)
+    {
+      AssertIndexRange(face_batch_number, n_ghost_inner_face_batches());
+    }
+  else
+    {
+      AssertIndexRange(face_batch_number, n_inner_face_batches());
+    }
+
+  AssertIndexRange(vector_number,
+                   n_active_entries_per_face_batch(face_batch_number));
+
+  const DoFHandler<dim> *dofh = nullptr;
+  if (dof_handlers.active_dof_handler == DoFHandlers::usual)
+    {
+      AssertDimension(dof_handlers.dof_handler.size(),
+                      dof_handlers.n_dof_handlers);
+      dofh = dof_handlers.dof_handler[fe_component];
+    }
+  else
+    {
+      Assert(false,
+             ExcMessage("Cannot return DoFHandler<dim>::cell_iterator "
+                        "for underlying DoFHandler!"));
+    }
+
+  const internal::MatrixFreeFunctions::FaceToCellTopology<
+    VectorizedArrayType::n_array_elements>
+    face2cell_info = get_face_info(face_batch_number);
+
+  const unsigned int cell_index =
+    interior ? face2cell_info.cells_interior[vector_number] :
+               face2cell_info.cells_exterior[vector_number];
+
+  std::pair<unsigned int, unsigned int> index = cell_level_index[cell_index];
+
+  return {typename DoFHandler<dim>::cell_iterator(
+            &dofh->get_triangulation(), index.first, index.second, dofh),
+          interior ? face2cell_info.interior_face_no :
+                     face2cell_info.exterior_face_no};
+}
+
+
+
 template <int dim, typename Number, typename VectorizedArrayType>
 typename hp::DoFHandler<dim>::active_cell_iterator
 MatrixFree<dim, Number, VectorizedArrayType>::get_hp_cell_iterator(