]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add internal function for identification of faces. 6210/head
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 16 Apr 2018 10:00:27 +0000 (12:00 +0200)
committerKatharina Kormann <katharina.kormann@tum.de>
Fri, 20 Apr 2018 20:09:38 +0000 (22:09 +0200)
include/deal.II/matrix_free/face_setup_internal.h [new file with mode: 0644]

diff --git a/include/deal.II/matrix_free/face_setup_internal.h b/include/deal.II/matrix_free/face_setup_internal.h
new file mode 100644 (file)
index 0000000..368364c
--- /dev/null
@@ -0,0 +1,1068 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii_face_setup_internal_h
+#define dealii_face_setup_internal_h
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/memory_consumption.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/matrix_free/task_info.h>
+#include <deal.II/matrix_free/face_info.h>
+
+#include <fstream>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace internal
+{
+  namespace MatrixFreeFunctions
+  {
+    /**
+     * A struct that is used to represent a collection of faces of a process
+     * with one of its neighbor within the setup done in struct FaceInfo.
+     *
+     * @author Katharina Kormann, Martin Kronbichler, 2018
+     */
+    struct FaceIdentifier
+    {
+      FaceIdentifier()
+        :
+        n_hanging_faces_smaller_subdomain (0),
+        n_hanging_faces_larger_subdomain  (0)
+      {}
+
+      std::vector<std::pair<CellId,CellId> > shared_faces;
+      unsigned int n_hanging_faces_smaller_subdomain;
+      unsigned int n_hanging_faces_larger_subdomain;
+    };
+
+
+
+    /**
+     * A struct that extracts the faces relevant to a given set of cells,
+     * including the assignment of which of the two neighboring processors at
+     * a subdomain boundary with MPI should do the integration (from both
+     * sides). This data structure is used for the setup of the connectivity
+     * between faces and cells and for identification of the dof indices to be
+     * used for face integrals.
+     *
+     * @author Katharina Kormann, Martin Kronbichler, 2018
+     */
+    template <int dim>
+    struct FaceSetup
+    {
+      FaceSetup();
+
+      /**
+       * Perform the initial detection of faces before reading the indices on
+       * the cells. This does not add the faces yet but only decides on
+       * whether some of the faces should be considered for processing
+       * locally.
+       */
+      template <typename MFAddData>
+      void initialize(const dealii::Triangulation<dim> &triangulation,
+                      const MFAddData                  &additional_data,
+                      std::vector<std::pair<unsigned int,unsigned int> > &cell_levels);
+
+      /**
+       * Upon completion of the dof indices, this function extracts the
+       * information relevant for FaceToCellTopology and categorizes the faces
+       * into interior faces, boundary faces, and ghost faces (not processed
+       * locally but adjacent to some of the cells present locally).
+       */
+      void
+      generate_faces(const dealii::Triangulation<dim> &triangulation,
+                     const std::vector<std::pair<unsigned int,unsigned int> > &cell_levels,
+                     TaskInfo                         &task_info);
+
+      /**
+       * Fills the information about the cell, the face number, and numbers
+       * within the plain array representation in MatrixFree into
+       * FaceToCellTopology (without vectorization, which is something applied
+       * later).
+       */
+      FaceToCellTopology<1>
+      create_face(const unsigned int                                        face_no,
+                  const typename dealii::Triangulation<dim>::cell_iterator &cell,
+                  const unsigned int                                        number_cell_interior,
+                  const typename dealii::Triangulation<dim>::cell_iterator &neighbor,
+                  const unsigned int                                        number_cell_exterior);
+
+      bool use_active_cells;
+
+      /**
+       * A type that categorizes faces in the first initialize() function such
+       * that we can later get their correct value in generate_faces().
+       */
+      enum class FaceCategory : char
+      {
+        locally_active_at_boundary,
+        locally_active_done_here,
+        locally_active_done_elsewhere,
+        ghosted,
+        multigrid_refinement_edge
+      };
+
+      std::vector<FaceCategory> face_is_owned;
+      std::vector<bool> at_processor_boundary;
+      std::vector<unsigned int> cells_close_to_boundary;
+      std::vector<FaceToCellTopology<1> > inner_faces;
+      std::vector<FaceToCellTopology<1> > boundary_faces;
+      std::vector<FaceToCellTopology<1> > inner_ghost_faces;
+      std::vector<FaceToCellTopology<1> > refinement_edge_faces;
+    };
+
+
+
+    /**
+     * Actually form the batches for vectorized execution of face integrals.
+     */
+    template <int vectorization_width>
+    void
+    collect_faces_vectorization
+    (const std::vector<FaceToCellTopology<1> > &faces_in,
+     const std::vector<bool>                   &hard_vectorization_boundary,
+     std::vector<unsigned int>                 &face_partition_data,
+     std::vector<FaceToCellTopology<vectorization_width> > &faces_out);
+
+
+
+    /* -------------------------------------------------------------------- */
+
+#ifndef DOXYGEN
+
+    template <int dim>
+    FaceSetup<dim>::FaceSetup()
+      :
+      use_active_cells (true)
+    {}
+
+
+
+    template <int dim>
+    template <typename MFAddData>
+    void
+    FaceSetup<dim>::initialize(const dealii::Triangulation<dim> &triangulation,
+                               const MFAddData                  &additional_data,
+                               std::vector<std::pair<unsigned int,unsigned int> > &cell_levels)
+    {
+      use_active_cells = additional_data.level_mg_handler == numbers::invalid_unsigned_int;
+
+#ifdef DEBUG
+      // safety check
+      if (use_active_cells)
+        for (unsigned int i=0; i<cell_levels.size(); ++i)
+          {
+            typename dealii::Triangulation<dim>::cell_iterator
+            dcell(&triangulation, cell_levels[i].first,
+                  cell_levels[i].second);
+            Assert(dcell->active(), ExcInternalError());
+          }
+#endif
+
+      // step 1: add ghost cells for those cells that we identify as
+      // interesting
+
+      at_processor_boundary.resize(cell_levels.size(), false);
+      cells_close_to_boundary.clear();
+      face_is_owned.resize(dim > 1 ? triangulation.n_raw_faces() :
+                           triangulation.n_vertices(), FaceCategory::locally_active_done_elsewhere);
+
+      // go through the mesh and divide the faces on the processor
+      // boundaries as evenly as possible between the processors
+      std::map<types::subdomain_id, FaceIdentifier> inner_faces_at_proc_boundary;
+      if (dynamic_cast<const parallel::Triangulation<dim>*>(&triangulation))
+        {
+          const types::subdomain_id my_domain = triangulation.locally_owned_subdomain();
+          for (unsigned int i=0; i<cell_levels.size(); ++i)
+            {
+              if (i>0 && cell_levels[i] == cell_levels[i-1])
+                continue;
+              typename dealii::Triangulation<dim>::cell_iterator
+              dcell(&triangulation, cell_levels[i].first,
+                    cell_levels[i].second);
+              for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+                {
+                  if (dcell->at_boundary(f) && !dcell->has_periodic_neighbor(f))
+                    continue;
+                  typename dealii::Triangulation<dim>::cell_iterator
+                  neighbor = dcell->neighbor_or_periodic_neighbor(f);
+
+                  // faces at hanging nodes are always treated by the processor
+                  // who owns the element on the fine side. but we need to count
+                  // the number of inner faces in order to balance the remaining
+                  // faces properly
+                  const CellId id_mine = dcell->id();
+                  if (use_active_cells && neighbor->has_children())
+                    for (unsigned int c=0; c<dcell->face(f)->n_children(); ++c)
+                      {
+                        typename dealii::Triangulation<dim>::cell_iterator
+                        neighbor_c = dcell->at_boundary(f) ?
+                                     dcell->periodic_neighbor_child_on_subface(f, c) :
+                                     dcell->neighbor_child_on_subface(f, c);
+                        const types::subdomain_id neigh_domain = neighbor_c->subdomain_id();
+                        if (my_domain < neigh_domain)
+                          inner_faces_at_proc_boundary[neigh_domain].n_hanging_faces_larger_subdomain++;
+                        else if (my_domain > neigh_domain)
+                          inner_faces_at_proc_boundary[neigh_domain].n_hanging_faces_smaller_subdomain++;
+                      }
+                  else
+                    {
+                      const types::subdomain_id neigh_domain = use_active_cells ?
+                                                               neighbor->subdomain_id() :
+                                                               neighbor->level_subdomain_id();
+                      if (neighbor->level() < dcell->level() &&
+                          use_active_cells)
+                        {
+                          if (my_domain < neigh_domain)
+                            inner_faces_at_proc_boundary[neigh_domain].n_hanging_faces_smaller_subdomain++;
+                          else if (my_domain > neigh_domain)
+                            inner_faces_at_proc_boundary[neigh_domain].n_hanging_faces_larger_subdomain++;
+                        }
+                      else if (neighbor->level() == dcell->level() &&
+                               my_domain != neigh_domain)
+                        {
+                          // always list the cell whose owner has the lower
+                          // subdomain id first. this applies to both processors
+                          // involved, so both processors will generate the same
+                          // list that we will later order
+                          const CellId id_neigh = neighbor->id();
+                          inner_faces_at_proc_boundary[neigh_domain].shared_faces.
+                          push_back(my_domain < neigh_domain
+                                    ?
+                                    std::make_pair(id_mine, id_neigh)
+                                    :
+                                    std::make_pair(id_neigh, id_mine));
+                        }
+                    }
+                }
+            }
+
+          // sort the cell ids related to each neighboring processor. This
+          // algorithm is symmetric so every processor combination should
+          // arrive here and no deadlock should be possible
+          for (std::map<types::subdomain_id, FaceIdentifier>::iterator it
+               = inner_faces_at_proc_boundary.begin();
+               it != inner_faces_at_proc_boundary.end(); ++it)
+            {
+              Assert(it->first != my_domain,
+                     ExcInternalError("Should not send info to myself"));
+              std::sort(it->second.shared_faces.begin(), it->second.shared_faces.end());
+              it->second.shared_faces.erase(std::unique(it->second.shared_faces.begin(),
+                                                        it->second.shared_faces.end()),
+                                            it->second.shared_faces.end());
+
+              // safety check: both involved processors should see the same list
+              // because the pattern of ghosting is symmetric. We test this by
+              // looking at the length of the lists of faces
+#if defined(DEAL_II_WITH_MPI) && defined(DEBUG)
+              MPI_Comm comm = MPI_COMM_SELF;
+              if (const parallel::Triangulation<dim> *ptria =
+                    dynamic_cast<const parallel::Triangulation<dim> *>(&triangulation))
+                comm = ptria->get_communicator();
+
+              MPI_Status status;
+              unsigned int mysize = it->second.shared_faces.size();
+              unsigned int othersize = numbers::invalid_unsigned_int;
+              MPI_Sendrecv(&mysize, 1, MPI_UNSIGNED, it->first, 600+my_domain,
+                           &othersize, 1, MPI_UNSIGNED, it->first, 600+it->first,
+                           comm, &status);
+              AssertDimension(mysize, othersize);
+              mysize = it->second.n_hanging_faces_smaller_subdomain;
+              MPI_Sendrecv(&mysize, 1, MPI_UNSIGNED, it->first, 700+my_domain,
+                           &othersize, 1, MPI_UNSIGNED, it->first, 700+it->first,
+                           comm, &status);
+              AssertDimension(mysize, othersize);
+              mysize = it->second.n_hanging_faces_larger_subdomain;
+              MPI_Sendrecv(&mysize, 1, MPI_UNSIGNED, it->first, 800+my_domain,
+                           &othersize, 1, MPI_UNSIGNED, it->first, 800+it->first,
+                           comm, &status);
+              AssertDimension(mysize, othersize);
+#endif
+
+              // Arrange the face "ownership" such that cells that are access
+              // by more than one face (think of a cell in a corner) get
+              // ghosted. This arrangement has the advantage that we need to
+              // send less data because the same data is used twice. The
+              // strategy applied here is to ensure the same order of face
+              // pairs on both processors that share some faces, and make the
+              // same decision on both sides.
+
+              // Create a vector with cell ids sorted over the processor with
+              // the larger rank. In the code below we need to be able to
+              // identify the same cell once for the processor with higher
+              // rank and once for the processor with the lower rank. The
+              // format for the processor with the higher rank is already
+              // contained in `shared_faces`, whereas we need a copy that we
+              // sort differently for the other way around.
+              std::vector<std::tuple<CellId,CellId,unsigned int> > other_range(it->second.shared_faces.size());
+              for (unsigned int i=0; i<other_range.size(); ++i)
+                other_range[i] = std::make_tuple(it->second.shared_faces[i].second,
+                                                 it->second.shared_faces[i].first,
+                                                 i);
+              std::sort(other_range.begin(), other_range.end());
+
+              // the vector 'assignment' sets whether a particular cell
+              // appears more often and acts as a pre-selection of the rank. A
+              // value of 1 means that the process with the higher rank gets
+              // those faces, a value -1 means that the process with the lower
+              // rank gets it, whereas a value 0 means that the decision can
+              // be made in an arbitrary way.
+              unsigned int n_faces_lower_proc = 0, n_faces_higher_proc = 0;
+              std::vector<char> assignment(other_range.size(), 0);
+              if (it->second.shared_faces.size() > 0)
+                {
+                  // identify faces that go to the processor with the higher
+                  // rank
+                  unsigned int count = 0;
+                  for (unsigned int i=1; i<it->second.shared_faces.size(); ++i)
+                    if (it->second.shared_faces[i].first ==
+                        it->second.shared_faces[i-1-count].first)
+                      ++count;
+                    else
+                      {
+                        AssertThrow(count < 2*dim, ExcInternalError());
+                        if (count > 0)
+                          {
+                            for (unsigned int k=0; k<=count; ++k)
+                              assignment[i-1-k] = 1;
+                            n_faces_higher_proc += count+1;
+                          }
+                        count = 0;
+                      }
+
+                  // identify faces that definitely go to the processor with
+                  // the lower rank - this must use the sorting of CellId
+                  // variables from the processor with the higher rank, i.e.,
+                  // other_range rather than `shared_faces`.
+                  count = 0;
+                  for (unsigned int i=1; i<other_range.size(); ++i)
+                    if (std::get<0>(other_range[i]) ==
+                        std::get<0>(other_range[i-1-count]))
+                      ++count;
+                    else
+                      {
+                        AssertThrow(count < 2*dim, ExcInternalError());
+                        if (count > 0)
+                          {
+                            for (unsigned int k=0; k<=count; ++k)
+                              {
+                                Assert(it->second.shared_faces[std::get<2>(other_range[i-1])].second ==
+                                       it->second.shared_faces[std::get<2>(other_range[i-1-k])].second,
+                                       ExcInternalError());
+                                // only assign to -1 if higher rank was not
+                                // yet set
+                                if (assignment[std::get<2>(other_range[i-1-k])] == 0)
+                                  {
+                                    assignment[std::get<2>(other_range[i-1-k])] = -1;
+                                    ++n_faces_lower_proc;
+                                  }
+                              }
+                          }
+                        count = 0;
+                      }
+                }
+
+
+              // divide the faces evenly between the two processors. the
+              // processor with small rank takes the first half, the processor
+              // with larger rank the second half. Adjust for the hanging
+              // faces that always get assigned to one side, and the faces we
+              // have already assigned due to the criterion above
+              n_faces_lower_proc += it->second.n_hanging_faces_smaller_subdomain;
+              n_faces_higher_proc += it->second.n_hanging_faces_larger_subdomain;
+              const unsigned int n_total_faces_at_proc_boundary =
+                (it->second.shared_faces.size() +
+                 it->second.n_hanging_faces_smaller_subdomain +
+                 it->second.n_hanging_faces_larger_subdomain);
+              unsigned int split_index = n_total_faces_at_proc_boundary/2;
+              if (split_index < n_faces_lower_proc)
+                split_index = 0;
+              else if (split_index < n_total_faces_at_proc_boundary -
+                       n_faces_higher_proc)
+                split_index -= n_faces_lower_proc;
+              else
+                split_index = n_total_faces_at_proc_boundary - n_faces_higher_proc - n_faces_lower_proc;
+
+              // make sure the splitting is consistent between both sides
+#if defined(DEAL_II_WITH_MPI) && defined(DEBUG)
+              MPI_Sendrecv(&split_index, 1, MPI_UNSIGNED, it->first, 900+my_domain,
+                           &othersize, 1, MPI_UNSIGNED, it->first, 900+it->first,
+                           comm, &status);
+              AssertDimension(split_index, othersize);
+              MPI_Sendrecv(&n_faces_lower_proc, 1, MPI_UNSIGNED, it->first, 1000+my_domain,
+                           &othersize, 1, MPI_UNSIGNED, it->first, 1000+it->first,
+                           comm, &status);
+              AssertDimension(n_faces_lower_proc, othersize);
+              MPI_Sendrecv(&n_faces_higher_proc, 1, MPI_UNSIGNED, it->first, 1100+my_domain,
+                           &othersize, 1, MPI_UNSIGNED, it->first, 1100+it->first,
+                           comm, &status);
+              AssertDimension(n_faces_higher_proc, othersize);
+#endif
+
+              // collect the faces on both sides
+              std::vector<std::pair<CellId,CellId> > owned_faces_lower, owned_faces_higher;
+              for (unsigned int i=0; i<assignment.size(); ++i)
+                if (assignment[i] < 0)
+                  owned_faces_lower.push_back(it->second.shared_faces[i]);
+                else if (assignment[i] > 0)
+                  owned_faces_higher.push_back(it->second.shared_faces[i]);
+              AssertIndexRange(split_index, it->second.shared_faces.size() + 1 -
+                               owned_faces_lower.size() - owned_faces_higher.size());
+
+              unsigned int i=0, c=0;
+              for (; i<assignment.size() && c < split_index; ++i)
+                if (assignment[i] == 0)
+                  {
+                    owned_faces_lower.push_back(it->second.shared_faces[i]);
+                    ++c;
+                  }
+              for ( ; i<assignment.size(); ++i)
+                if (assignment[i] == 0)
+                  {
+                    owned_faces_higher.push_back(it->second.shared_faces[i]);
+                  }
+
+#ifdef DEBUG
+              // check consistency of faces on both sides
+              std::vector<std::pair<CellId,CellId> > check_faces;
+              check_faces.insert(check_faces.end(), owned_faces_lower.begin(),
+                                 owned_faces_lower.end());
+              check_faces.insert(check_faces.end(), owned_faces_higher.begin(),
+                                 owned_faces_higher.end());
+              std::sort(check_faces.begin(), check_faces.end());
+              AssertDimension(check_faces.size(), it->second.shared_faces.size());
+              for (unsigned int i=0; i<check_faces.size(); ++i)
+                Assert(check_faces[i] == it->second.shared_faces[i], ExcInternalError());
+#endif
+
+              // now only set half of the faces as the ones to keep
+              if (my_domain < it->first)
+                it->second.shared_faces.swap(owned_faces_lower);
+              else
+                it->second.shared_faces.swap(owned_faces_higher);
+
+              std::sort(it->second.shared_faces.begin(), it->second.shared_faces.end());
+            }
+        }
+
+      // fill in the additional cells that we need access to via ghosting to
+      // cell_levels
+      std::set<std::pair<unsigned int, unsigned int> > ghost_cells;
+      for (unsigned int i=0; i<cell_levels.size(); ++i)
+        {
+          typename dealii::Triangulation<dim>::cell_iterator
+          dcell(&triangulation, cell_levels[i].first,
+                cell_levels[i].second);
+          if (use_active_cells)
+            Assert(dcell->active(), ExcNotImplemented());
+          for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+            {
+              if (dcell->at_boundary(f) && !dcell->has_periodic_neighbor(f))
+                face_is_owned[dcell->face(f)->index()] = FaceCategory::locally_active_at_boundary;
+
+              // treat boundaries of cells of different refinement level
+              // inside the domain in case of multigrid separately
+              else if ((dcell->at_boundary(f) == false ||
+                        dcell->has_periodic_neighbor(f)) &&
+                       additional_data.level_mg_handler != numbers::invalid_unsigned_int
+                       &&
+                       dcell->neighbor_or_periodic_neighbor(f)->level() < dcell->level())
+                {
+                  face_is_owned[dcell->face(f)->index()] = FaceCategory::multigrid_refinement_edge;
+                }
+              else
+                {
+                  typename dealii::Triangulation<dim>::cell_iterator neighbor =
+                    dcell->neighbor_or_periodic_neighbor(f);
+
+                  // neighbor is refined -> face will be treated by neighbor
+                  if (use_active_cells && neighbor->has_children() &&
+                      additional_data.hold_all_faces_to_owned_cells == false)
+                    continue;
+
+                  bool add_to_ghost = false;
+                  const types::subdomain_id
+                  id1 = use_active_cells ? dcell->subdomain_id() : dcell->level_subdomain_id(),
+                  id2 = use_active_cells ? (neighbor->has_children() ?
+                                            dcell->neighbor_child_on_subface(f,0)->subdomain_id() :
+                                            neighbor->subdomain_id()) :
+                        neighbor->level_subdomain_id();
+
+                  // Check whether the current face should be processed
+                  // locally (instead of being processed from the other
+                  // side). We process a face locally when we are more refined
+                  // (in the active cell case) or when the face is listed in
+                  // the `shared_faces` data structure that we built above.
+                  if ((id1 == id2 && (use_active_cells == false ||
+                                      neighbor->has_children() == false)) ||
+                      dcell->level() > neighbor->level() ||
+                      std::binary_search(inner_faces_at_proc_boundary[id2].shared_faces.begin(),
+                                         inner_faces_at_proc_boundary[id2].shared_faces.end(),
+                                         std::make_pair(id1<id2 ? dcell->id() : neighbor->id(),
+                                                        id1<id2 ? neighbor->id() : dcell->id())))
+                    {
+                      face_is_owned[dcell->face(f)->index()] = FaceCategory::locally_active_done_here;
+                      if (dcell->level() == neighbor->level())
+                        face_is_owned[neighbor->face(dcell->has_periodic_neighbor(f) ?
+                                                     dcell->periodic_neighbor_face_no(f) :
+                                                     dcell->neighbor_face_no(f))->index()]
+                          = FaceCategory::locally_active_done_here;
+
+                      // If neighbor is a ghost element (i.e. dcell->subdomain_id
+                      // ! dcell->neighbor(f)->subdomain_id()), we need to add its
+                      // index into cell level list.
+                      if (use_active_cells)
+                        add_to_ghost = (dcell->subdomain_id() !=
+                                        neighbor->subdomain_id());
+                      else
+                        add_to_ghost = (dcell->level_subdomain_id() !=
+                                        neighbor->level_subdomain_id());
+                    }
+                  else if (additional_data.hold_all_faces_to_owned_cells == false)
+                    {
+                      // mark the cell to be close to the boundary
+                      cells_close_to_boundary.emplace_back(i);
+                    }
+                  else
+                    {
+                      // add all cells to ghost layer...
+                      face_is_owned[dcell->face(f)->index()] = FaceCategory::ghosted;
+                      if (use_active_cells)
+                        {
+                          if (neighbor->has_children())
+                            for (unsigned int s=0; s<dcell->face(f)->n_children(); ++s)
+                              if (dcell->at_boundary(f))
+                                {
+                                  if (dcell->periodic_neighbor_child_on_subface(f,s)->subdomain_id() !=
+                                      dcell->subdomain_id())
+                                    add_to_ghost = true;
+                                }
+                              else
+                                {
+                                  if (dcell->neighbor_child_on_subface(f,s)->subdomain_id() !=
+                                      dcell->subdomain_id())
+                                    add_to_ghost = true;
+                                }
+                          else
+                            add_to_ghost = (dcell->subdomain_id() !=
+                                            neighbor->subdomain_id());
+                        }
+                      else
+                        add_to_ghost = (dcell->level_subdomain_id() !=
+                                        neighbor->level_subdomain_id());
+                    }
+
+                  if (add_to_ghost)
+                    {
+                      ghost_cells.insert
+                      (std::pair<unsigned int,unsigned int>
+                       (neighbor->level(), neighbor->index()));
+                      at_processor_boundary[i] = true;
+                    }
+                }
+            }
+        }
+
+      // step 2: append the ghost cells at the end of the locally owned
+      // cells
+      for (std::set<std::pair<unsigned int, unsigned int> >::iterator
+           it=ghost_cells.begin(); it!=ghost_cells.end(); ++it)
+        cell_levels.push_back(*it);
+
+      // step 3: clean up the cells close to the boundary
+      std::sort(cells_close_to_boundary.begin(), cells_close_to_boundary.end());
+      cells_close_to_boundary.erase(std::unique(cells_close_to_boundary.begin(),
+                                                cells_close_to_boundary.end()),
+                                    cells_close_to_boundary.end());
+      std::vector<unsigned int> final_cells;
+      final_cells.reserve(cells_close_to_boundary.size());
+      for (unsigned int i=0; i<cells_close_to_boundary.size(); ++i)
+        if (at_processor_boundary[cells_close_to_boundary[i]] == false)
+          final_cells.push_back(cells_close_to_boundary[i]);
+      cells_close_to_boundary = std::move(final_cells);
+    }
+
+
+
+    template <int dim>
+    void
+    FaceSetup<dim>
+    ::generate_faces(const dealii::Triangulation<dim> &triangulation,
+                     const std::vector<std::pair<unsigned int,unsigned int> > &cell_levels,
+                     TaskInfo &task_info)
+    {
+      // step 1: create the inverse map between cell iterators and the
+      // cell_level_index field
+      std::map<std::pair<unsigned int,unsigned int>, unsigned int>
+      map_to_vectorized;
+      for (unsigned int cell=0; cell<cell_levels.size(); ++cell)
+        if (cell == 0 || cell_levels[cell] != cell_levels[cell-1])
+          {
+            typename dealii::Triangulation<dim>::cell_iterator dcell
+            (&triangulation, cell_levels[cell].first, cell_levels[cell].second);
+            std::pair<unsigned int,unsigned int> level_index(dcell->level(),
+                                                             dcell->index());
+            map_to_vectorized[level_index] = cell;
+          }
+
+      // step 2: fill the information about inner faces and boundary faces
+      const unsigned int vectorization_length = task_info.vectorization_length;
+      task_info.face_partition_data.resize(task_info.cell_partition_data.size()-1, 0);
+      task_info.boundary_partition_data.resize(task_info.cell_partition_data.size()-1, 0);
+      std::vector<unsigned char> face_visited(face_is_owned.size(), 0);
+      for (unsigned int partition=0; partition<task_info.cell_partition_data.size()-2; ++partition)
+        {
+          unsigned int boundary_counter = 0;
+          unsigned int inner_counter = 0;
+          for (unsigned int cell=
+                 task_info.cell_partition_data[partition]*vectorization_length;
+               cell<task_info.cell_partition_data[partition+1]*vectorization_length;
+               ++cell)
+            if (cell == 0 || cell_levels[cell] != cell_levels[cell-1])
+              {
+                typename dealii::Triangulation<dim>::cell_iterator dcell
+                (&triangulation, cell_levels[cell].first, cell_levels[cell].second);
+                for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+                  {
+                    // boundary face
+                    if (face_is_owned[dcell->face(f)->index()] == FaceCategory::locally_active_at_boundary)
+                      {
+                        Assert(dcell->at_boundary(f), ExcInternalError());
+                        ++boundary_counter;
+                        FaceToCellTopology<1> info;
+                        info.cells_interior[0] = cell;
+                        info.cells_exterior[0] = numbers::invalid_unsigned_int;
+                        info.interior_face_no = f;
+                        info.exterior_face_no = dcell->face(f)->boundary_id();
+                        info.subface_index = GeometryInfo<dim>::max_children_per_cell;
+                        info.face_orientation = 0;
+                        boundary_faces.push_back(info);
+
+                        face_visited[dcell->face(f)->index()]++;
+                      }
+                    // interior face, including faces over periodic boundaries
+                    else
+                      {
+                        typename dealii::Triangulation<dim>::cell_iterator neighbor
+                          = dcell->neighbor_or_periodic_neighbor(f);
+                        if (use_active_cells && neighbor->has_children())
+                          {
+                            for (unsigned int c=0; c<dcell->face(f)->n_children(); ++c)
+                              {
+                                typename dealii::Triangulation<dim>::cell_iterator
+                                neighbor_c = dcell->at_boundary(f) ?
+                                             dcell->periodic_neighbor_child_on_subface(f, c) :
+                                             dcell->neighbor_child_on_subface(f, c);
+                                const types::subdomain_id neigh_domain = neighbor_c->subdomain_id();
+                                const unsigned int neighbor_face_no =
+                                  dcell->has_periodic_neighbor(f) ?
+                                  dcell->periodic_neighbor_face_no(f) :
+                                  dcell->neighbor_face_no(f);
+                                if (neigh_domain != dcell->subdomain_id() ||
+                                    face_visited[dcell->face(f)->child(c)->index()] == 1)
+                                  {
+                                    std::pair<unsigned int,unsigned int>
+                                    level_index(neighbor_c->level(), neighbor_c->index());
+                                    if (face_is_owned[dcell->face(f)->child(c)->index()]==
+                                        FaceCategory::locally_active_done_here)
+                                      {
+                                        ++inner_counter;
+                                        inner_faces.
+                                        push_back(create_face(neighbor_face_no, neighbor_c,
+                                                              map_to_vectorized[level_index],
+                                                              dcell, cell));
+
+                                      }
+                                    else if (face_is_owned[dcell->face(f)->child(c)->index()] ==
+                                             FaceCategory::ghosted)
+                                      {
+                                        inner_ghost_faces.
+                                        push_back(create_face(neighbor_face_no, neighbor_c,
+                                                              map_to_vectorized[level_index],
+                                                              dcell, cell));
+                                      }
+                                    else
+                                      Assert(face_is_owned[dcell->face(f)->index()] ==
+                                             FaceCategory::locally_active_done_elsewhere,
+                                             ExcInternalError());
+                                  }
+                                else
+                                  {
+                                    face_visited[dcell->face(f)->child(c)->index()] = 1;
+                                  }
+                              }
+                          }
+                        else
+                          {
+                            const types::subdomain_id my_domain = use_active_cells ?
+                                                                  dcell->subdomain_id() : dcell->level_subdomain_id();
+                            const types::subdomain_id neigh_domain = use_active_cells ?
+                                                                     neighbor->subdomain_id() : neighbor->level_subdomain_id();
+                            if (neigh_domain != my_domain ||
+                                face_visited[dcell->face(f)->index()] == 1)
+                              {
+                                std::pair<unsigned int,unsigned int>
+                                level_index(neighbor->level(), neighbor->index());
+                                if (face_is_owned[dcell->face(f)->index()]==
+                                    FaceCategory::locally_active_done_here)
+                                  {
+                                    ++inner_counter;
+                                    inner_faces.
+                                    push_back(create_face(f, dcell, cell, neighbor,
+                                                          map_to_vectorized[level_index]));
+                                  }
+                                else if (face_is_owned[dcell->face(f)->index()] ==
+                                         FaceCategory::ghosted)
+                                  {
+                                    inner_ghost_faces.
+                                    push_back(create_face(f, dcell, cell, neighbor,
+                                                          map_to_vectorized[level_index]));
+                                  }
+                              }
+                            else
+                              {
+                                face_visited[dcell->face(f)->index()] = 1;
+                                if (dcell->has_periodic_neighbor(f))
+                                  face_visited[neighbor->face(dcell->periodic_neighbor_face_no(f))->index()] = 1;
+                              }
+                            if (face_is_owned[dcell->face(f)->index()] ==
+                                FaceCategory::multigrid_refinement_edge)
+                              {
+                                refinement_edge_faces.
+                                push_back(create_face(f, dcell, cell, neighbor,
+                                                      refinement_edge_faces.size()));
+                              }
+                          }
+                      }
+                  }
+              }
+          task_info.face_partition_data[partition+1]=
+            task_info.face_partition_data[partition] + inner_counter;
+          task_info.boundary_partition_data[partition+1] =
+            task_info.boundary_partition_data[partition] + boundary_counter;
+        }
+      task_info.ghost_face_partition_data.resize(2);
+      task_info.ghost_face_partition_data[0] = 0;
+      task_info.ghost_face_partition_data[1] = inner_ghost_faces.size();
+      task_info.refinement_edge_face_partition_data.resize(2);
+      task_info.refinement_edge_face_partition_data[0] = 0;
+      task_info.refinement_edge_face_partition_data[1] = refinement_edge_faces.size();
+    }
+
+
+
+    template <int dim>
+    FaceToCellTopology<1>
+    FaceSetup<dim>
+    ::create_face(const unsigned int                                        face_no,
+                  const typename dealii::Triangulation<dim>::cell_iterator &cell,
+                  const unsigned int                                        number_cell_interior,
+                  const typename dealii::Triangulation<dim>::cell_iterator &neighbor,
+                  const unsigned int                                        number_cell_exterior)
+    {
+      FaceToCellTopology<1> info;
+      info.cells_interior[0] = number_cell_interior;
+      info.cells_exterior[0] = number_cell_exterior;
+      info.interior_face_no = face_no;
+      if (cell->has_periodic_neighbor(face_no))
+        info.exterior_face_no = cell->periodic_neighbor_face_no(face_no);
+      else
+        info.exterior_face_no = cell->neighbor_face_no(face_no);
+
+      info.subface_index = GeometryInfo<dim>::max_children_per_cell;
+      Assert(neighbor->level() <= cell->level(),
+             ExcInternalError());
+      if (cell->level() > neighbor->level())
+        {
+          if (cell->has_periodic_neighbor(face_no))
+            info.subface_index =
+              cell->periodic_neighbor_of_coarser_periodic_neighbor(face_no).second;
+          else
+            info.subface_index =
+              cell->neighbor_of_coarser_neighbor(face_no).second;
+        }
+
+      info.face_orientation = 0;
+      unsigned int  left_face_orientation =
+        !cell->face_orientation(face_no) + 2 * cell->face_flip(face_no) +
+        4 * cell->face_rotation(face_no);
+      unsigned int right_face_orientation =
+        !neighbor->face_orientation(info.exterior_face_no) +
+        2 * neighbor->face_flip(info.exterior_face_no) +
+        4 * neighbor->face_rotation(info.exterior_face_no);
+      if (left_face_orientation != 0)
+        {
+          info.face_orientation = 8 + left_face_orientation;
+          Assert(right_face_orientation == 0,
+                 ExcMessage("Face seems to be wrongly oriented from both sides"));
+        }
+      else
+        info.face_orientation = right_face_orientation;
+      return info;
+    }
+
+
+
+    /**
+     * This simple comparison for collect_faces_vectorization() identifies
+     * faces of the same type, i.e., where all of the interior and exterior
+     * face number, subface index and orientation are the same. This is used
+     * to batch similar faces together for vectorization.
+     */
+    bool compare_faces_for_vectorization
+    (const FaceToCellTopology<1> &face1,
+     const FaceToCellTopology<1> &face2)
+    {
+      if (face1.interior_face_no != face2.interior_face_no)
+        return false;
+      if (face1.exterior_face_no != face2.exterior_face_no)
+        return false;
+      if (face1.subface_index != face2.subface_index)
+        return false;
+      if (face1.face_orientation != face2.face_orientation)
+        return false;
+      return true;
+    }
+
+
+
+    /**
+     * This comparator is used within collect_faces_vectorization() to create
+     * a sorting of FaceToCellTopology objects based on their
+     * identifiers. This is used to obtain a good data locality when
+     * processing the face integrals.
+     */
+    template <int length>
+    struct FaceComparator
+    {
+      bool operator() (const FaceToCellTopology<length> &face1,
+                       const FaceToCellTopology<length> &face2)
+      {
+        for (unsigned int i=0; i<length; ++i)
+          if (face1.cells_interior[i] < face2.cells_interior[i])
+            return true;
+          else if (face1.cells_interior[i] > face2.cells_interior[i])
+            return false;
+        for (unsigned int i=0; i<length; ++i)
+          if (face1.cells_exterior[i] < face2.cells_exterior[i])
+            return true;
+          else if (face1.cells_exterior[i] > face2.cells_exterior[i])
+            return false;
+        if (face1.interior_face_no < face2.interior_face_no)
+          return true;
+        else if (face1.interior_face_no > face2.interior_face_no)
+          return false;
+        if (face1.exterior_face_no < face2.exterior_face_no)
+          return true;
+        else if (face1.exterior_face_no > face2.exterior_face_no)
+          return false;
+
+        // we do not need to check for subface_index and orientation because
+        // those cannot be different if when all the other values are the
+        // same.
+        AssertDimension(face1.subface_index, face2.subface_index);
+        AssertDimension(face1.face_orientation, face2.face_orientation);
+
+        return false;
+      }
+    };
+
+
+
+    template <int vectorization_width>
+    void
+    collect_faces_vectorization
+    (const std::vector<FaceToCellTopology<1> > &faces_in,
+     const std::vector<bool>                   &hard_vectorization_boundary,
+     std::vector<unsigned int>                 &face_partition_data,
+     std::vector<FaceToCellTopology<vectorization_width> > &faces_out)
+    {
+      FaceToCellTopology<vectorization_width> macro_face;
+      std::vector<std::vector<unsigned int> > faces_type;
+
+      unsigned int face_start = face_partition_data[0],
+                   face_end = face_partition_data[0];
+
+      face_partition_data[0] = faces_out.size();
+      for ( unsigned int partition=0; partition<face_partition_data.size()-1; ++partition)
+        {
+          std::vector<std::vector<unsigned int> > new_faces_type;
+
+          // start with the end point for the last partition
+          face_start = face_end;
+          face_end = face_partition_data[partition+1];
+
+          // set the partitioner to the new vectorized lengths
+          face_partition_data[partition+1] = face_partition_data[partition];
+
+          // loop over the faces in the current partition and reorder according to
+          // the face type
+          for (unsigned int face=face_start; face<face_end; ++face)
+            {
+              for (unsigned int type=0; type<faces_type.size(); ++type)
+                {
+                  // Compare current face with first face of type type
+                  if ( compare_faces_for_vectorization(faces_in[face],
+                                                       faces_in[faces_type[type][0]]) )
+                    {
+                      faces_type[type].push_back(face);
+                      goto face_found;
+                    }
+                }
+              faces_type.push_back(std::vector<unsigned int>(1, face));
+face_found :
+              {}
+            }
+
+          // insert new faces in sorted list to get good data locality
+          std::set<FaceToCellTopology<vectorization_width>,
+              FaceComparator<vectorization_width> > new_faces;
+          for (unsigned int type=0; type<faces_type.size(); ++type)
+            {
+              macro_face.interior_face_no = faces_in[faces_type[type][0]].interior_face_no;
+              macro_face.exterior_face_no = faces_in[faces_type[type][0]].exterior_face_no;
+              macro_face.subface_index = faces_in[faces_type[type][0]].subface_index;
+              macro_face.face_orientation = faces_in[faces_type[type][0]].face_orientation;
+              unsigned int no_faces = faces_type[type].size();
+              std::vector<unsigned char> touched(no_faces, 0);
+
+              // do two passes through the data. The first is to identify similar
+              // faces within the same index range as the cells which will allow
+              // for vectorized read operations, the second picks up all the rest
+              unsigned int n_vectorized = 0;
+              for (unsigned int f=0; f<no_faces; ++f)
+                if (faces_in[faces_type[type][f]].cells_interior[0] % vectorization_width == 0)
+                  {
+                    bool is_contiguous = true;
+                    if (f + vectorization_width > no_faces)
+                      is_contiguous = false;
+                    else
+                      for (unsigned int v=1; v<vectorization_width; ++v)
+                        if (faces_in[faces_type[type][f+v]].cells_interior[0] !=
+                            faces_in[faces_type[type][f]].cells_interior[0]+v)
+                          is_contiguous = false;
+                    if (is_contiguous)
+                      {
+                        AssertIndexRange(f, faces_type[type].size()-vectorization_width+1);
+                        for (unsigned int v=0; v<vectorization_width; ++v)
+                          {
+                            macro_face.cells_interior[v] = faces_in[faces_type[type][f+v]].cells_interior[0];
+                            macro_face.cells_exterior[v] = faces_in[faces_type[type][f+v]].cells_exterior[0];
+                            touched[f+v] = 1;
+                          }
+                        new_faces.insert(macro_face);
+                        f += vectorization_width-1;
+                        n_vectorized += vectorization_width;
+                      }
+                  }
+
+              std::vector<unsigned int> untouched;
+              untouched.reserve(no_faces - n_vectorized);
+              for (unsigned int f=0; f<no_faces; ++f)
+                if (touched[f] == 0)
+                  untouched.push_back(f);
+              unsigned int v=0;
+              for (auto f : untouched)
+                {
+                  macro_face.cells_interior[v] = faces_in[faces_type[type][f]].cells_interior[0];
+                  macro_face.cells_exterior[v] = faces_in[faces_type[type][f]].cells_exterior[0];
+                  ++v;
+                  if (v == vectorization_width)
+                    {
+                      new_faces.insert(macro_face);
+                      v = 0;
+                    }
+                }
+              if (v > 0 && v < vectorization_width)
+                {
+                  // must add non-filled face
+                  if (hard_vectorization_boundary[partition+1] ||
+                      partition == face_partition_data.size()-2)
+                    {
+                      for ( ; v < vectorization_width; ++v)
+                        {
+                          // Dummy cell, not used
+                          macro_face.cells_interior[v] = numbers::invalid_unsigned_int;
+                          macro_face.cells_exterior[v] = numbers::invalid_unsigned_int;
+                        }
+                      new_faces.insert(macro_face);
+                    }
+                  else
+                    {
+                      // postpone to the next partition
+                      std::vector<unsigned int> untreated(v);
+                      for (unsigned int f=0; f<v; ++f)
+                        untreated[f] = faces_type[type][*(untouched.end()-1-f)];
+                      new_faces_type.push_back(untreated);
+                    }
+                }
+            }
+
+          // insert sorted list to vector of faces
+          for (auto it = new_faces.begin(); it != new_faces.end(); ++it)
+            faces_out.push_back(*it);
+          face_partition_data[partition+1] += new_faces.size();
+
+          // set the faces that were left over to faces_type for the next round
+          faces_type = std::move(new_faces_type);
+        }
+
+#ifdef DEBUG
+      // final safety checks
+      for (unsigned int i=0; i<faces_type.size(); ++i)
+        AssertDimension(faces_type[i].size(), 0U);
+
+      AssertDimension(faces_out.size(), face_partition_data.back());
+      unsigned int nfaces = 0;
+      for (unsigned int i=face_partition_data[0]; i<face_partition_data.back(); ++i)
+        for (unsigned int v=0; v<vectorization_width; ++v)
+          nfaces += (faces_out[i].cells_interior[v] != numbers::invalid_unsigned_int);
+      AssertDimension(nfaces, faces_in.size());
+
+      std::vector<std::pair<unsigned int,unsigned int>> in_faces, out_faces;
+      for (unsigned int i=0; i<faces_in.size(); ++i)
+        in_faces.emplace_back(faces_in[i].cells_interior[0],
+                              faces_in[i].cells_exterior[0]);
+      for (unsigned int i=face_partition_data[0]; i<face_partition_data.back(); ++i)
+        for (unsigned int v=0; v<vectorization_width &&
+             faces_out[i].cells_interior[v] != numbers::invalid_unsigned_int; ++v)
+          out_faces.emplace_back(faces_out[i].cells_interior[v],
+                                 faces_out[i].cells_exterior[v]);
+      std::sort(in_faces.begin(), in_faces.end());
+      std::sort(out_faces.begin(), out_faces.end());
+      AssertDimension(in_faces.size(), out_faces.size());
+      for (unsigned int i=0; i<in_faces.size(); ++i)
+        {
+          AssertDimension(in_faces[i].first, out_faces[i].first);
+          AssertDimension(in_faces[i].second, out_faces[i].second);
+        }
+#endif
+    }
+
+#endif // ifndef DOXYGEN
+
+  }
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#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.