Allow access of exterior dof values for face-centric FEFaceEvaluation

diff --git a/include/deal.II/matrix_free/dof_info.h b/include/deal.II/matrix_free/dof_info.h
index d263fdfeee0932179ffc07e636f43d47a0754ad7..9109769da1c949eac6db5c652352560053b3862c 100644 (file)
--- a/include/deal.II/matrix_free/dof_info.h
+++ b/include/deal.II/matrix_free/dof_info.h
@@ -262,7 +262,8 @@ namespace internal
       template <int length>
       void
       compute_face_index_compression(
-        const std::vector<FaceToCellTopology<length>> &faces);
+        const std::vector<FaceToCellTopology<length>> &faces,
+        bool hold_all_faces_to_owned_cells);
 
       /**
        * This function computes the connectivity of the currently stored

diff --git a/include/deal.II/matrix_free/dof_info.templates.h b/include/deal.II/matrix_free/dof_info.templates.h
index 942aa7df2e6ef7672c593e54ff1cd0e1dfd7dd1b..2c4b2dcc3450ad796841097d3add8dfacb77332e 100644 (file)
--- a/include/deal.II/matrix_free/dof_info.templates.h
+++ b/include/deal.II/matrix_free/dof_info.templates.h
@@ -243,7 +243,8 @@ namespace internal
     template <int length>
     void
     DoFInfo::compute_face_index_compression(
-      const std::vector<FaceToCellTopology<length>> &faces)
+      const std::vector<FaceToCellTopology<length>> &faces,
+      const bool hold_all_faces_to_owned_cells)
     {
       AssertDimension(length, vectorization_length);
 
@@ -256,20 +257,33 @@ namespace internal
       n_vectorization_lanes_filled[dof_access_face_interior].resize(
         faces.size());
 
-      // all interior faces come before the boundary faces
-      unsigned int n_exterior_faces = 0;
-      for (; n_exterior_faces < faces.size(); ++n_exterior_faces)
-        if (faces[n_exterior_faces].cells_exterior[0] ==
+      // all inner faces come before the boundary faces
+      unsigned int n_inner_faces = 0;
+      for (; n_inner_faces < faces.size(); ++n_inner_faces)
+        if (faces[n_inner_faces].cells_exterior[0] ==
+            numbers::invalid_unsigned_int)
+          break;
+
+      // all boundary faces come after the inner faces and before the ghosted
+      // inner faces
+      unsigned int n_boundary_faces = 0;
+      for (; n_inner_faces + n_boundary_faces < faces.size();
+           ++n_boundary_faces)
+        if (faces[n_inner_faces + n_boundary_faces].cells_exterior[0] !=
             numbers::invalid_unsigned_int)
           break;
+
+      const unsigned int size_exterior_faces =
+        hold_all_faces_to_owned_cells ? faces.size() : n_inner_faces;
+
       index_storage_variants[dof_access_face_exterior].resize(
-        n_exterior_faces, IndexStorageVariants::full);
+        size_exterior_faces, IndexStorageVariants::full);
       dof_indices_contiguous[dof_access_face_exterior].resize(
-        n_exterior_faces * length, numbers::invalid_unsigned_int);
+        size_exterior_faces * length, numbers::invalid_unsigned_int);
       dof_indices_interleave_strides[dof_access_face_exterior].resize(
         faces.size() * length, numbers::invalid_unsigned_int);
       n_vectorization_lanes_filled[dof_access_face_exterior].resize(
-        n_exterior_faces);
+        size_exterior_faces);
 
       for (unsigned int face = 0; face < faces.size(); ++face)
         {
@@ -362,7 +376,9 @@ namespace internal
 
           face_computation(dof_access_face_interior,
                            faces[face].cells_interior);
-          if (face < n_exterior_faces)
+          if (face < n_inner_faces ||
+              (hold_all_faces_to_owned_cells &&
+               face >= (n_inner_faces + n_boundary_faces)))
             face_computation(dof_access_face_exterior,
                              faces[face].cells_exterior);
         }

diff --git a/include/deal.II/matrix_free/face_setup_internal.h b/include/deal.II/matrix_free/face_setup_internal.h
index 6a1ef52ffb31a7dea2c92a186acc19e1242d75d8..6fe9850926a620c5abd811c85624e613360e7f6d 100644 (file)
--- a/include/deal.II/matrix_free/face_setup_internal.h
+++ b/include/deal.II/matrix_free/face_setup_internal.h
@@ -852,7 +852,10 @@ namespace internal
                                     else if (face_is_owned[dcell->face(f)
                                                              ->child(c)
                                                              ->index()] ==
-                                             FaceCategory::ghosted)
+                                               FaceCategory::ghosted ||
+                                             face_is_owned[dcell->face(f)
+                                                             ->index()] ==
+                                               FaceCategory::ghosted)
                                       {
                                         inner_ghost_faces.push_back(create_face(
                                           neighbor_face_no,
@@ -863,15 +866,11 @@ namespace internal
                                           is_mixed_mesh));
                                       }
                                     else
-                                      Assert(
-                                        face_is_owned[dcell->face(f)
-                                                        ->index()] ==
-                                            FaceCategory::
-                                              locally_active_done_elsewhere ||
-                                          face_is_owned[dcell->face(f)
-                                                          ->index()] ==
-                                            FaceCategory::ghosted,
-                                        ExcInternalError());
+                                      Assert(face_is_owned[dcell->face(f)
+                                                             ->index()] ==
+                                               FaceCategory::
+                                                 locally_active_done_elsewhere,
+                                             ExcInternalError());
                                   }
                                 else
                                   {

diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
index f7c6d6ff61b3548c95d941111cbac20ea916e911..ee9430614ea6b8b544b83f51fbea3794c0fcb29b 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -2026,7 +2026,8 @@ MatrixFree<dim, Number, VectorizedArrayType>::initialize_indices(
            task_info.refinement_edge_face_partition_data[0]));
 
       for (auto &di : dof_info)
-        di.compute_face_index_compression(face_info.faces);
+        di.compute_face_index_compression(
+          face_info.faces, additional_data.hold_all_faces_to_owned_cells);
 
       // build the inverse map back from the faces array to
       // cell_and_face_to_plain_faces

diff --git a/source/matrix_free/dof_info.cc b/source/matrix_free/dof_info.cc
index 71230646b879693193e5125f18d6ab9c4c86d93b..4804548a36aec3b1ec40fc70f813ddd62cc771d0 100644 (file)
--- a/source/matrix_free/dof_info.cc
+++ b/source/matrix_free/dof_info.cc
@@ -1589,19 +1589,24 @@ namespace internal
 
     template void
     DoFInfo::compute_face_index_compression<1>(
-      const std::vector<FaceToCellTopology<1>> &);
+      const std::vector<FaceToCellTopology<1>> &,
+      bool);
     template void
     DoFInfo::compute_face_index_compression<2>(
-      const std::vector<FaceToCellTopology<2>> &);
+      const std::vector<FaceToCellTopology<2>> &,
+      bool);
     template void
     DoFInfo::compute_face_index_compression<4>(
-      const std::vector<FaceToCellTopology<4>> &);
+      const std::vector<FaceToCellTopology<4>> &,
+      bool);
     template void
     DoFInfo::compute_face_index_compression<8>(
-      const std::vector<FaceToCellTopology<8>> &);
+      const std::vector<FaceToCellTopology<8>> &,
+      bool);
     template void
     DoFInfo::compute_face_index_compression<16>(
-      const std::vector<FaceToCellTopology<16>> &);
+      const std::vector<FaceToCellTopology<16>> &,
+      bool);
 
     template void
     DoFInfo::compute_vector_zero_access_pattern<1>(

diff --git a/tests/matrix_free/matrix_vector_faces_24.cc b/tests/matrix_free/matrix_vector_faces_24.cc
index 20bbe9c8a6ede2fc7a2425dc80a82cd4f82cda07..e96250eaf3ec9a29baf6b71076e12ac59226731f 100644 (file)
--- a/tests/matrix_free/matrix_vector_faces_24.cc
+++ b/tests/matrix_free/matrix_vector_faces_24.cc
@@ -103,4 +103,8 @@ test()
   mf.vmult(out, in);
 
   deallog << "Norm of result:          " << out.l2_norm() << std::endl;
+
+  mf.manual_loop_vmult(out, in);
+
+  deallog << "Norm of result:          " << out.l2_norm() << std::endl;
 }

diff --git a/tests/matrix_free/matrix_vector_faces_24.with_p4est=true.mpirun=2.output b/tests/matrix_free/matrix_vector_faces_24.with_p4est=true.mpirun=2.output
index 71f9902dab3333ddfb1e1da1ddde796426f21b37..ae81bd819f36cd526bda6176ee023b92c69ac6ab 100644 (file)
--- a/tests/matrix_free/matrix_vector_faces_24.with_p4est=true.mpirun=2.output
+++ b/tests/matrix_free/matrix_vector_faces_24.with_p4est=true.mpirun=2.output
@@ -1,5 +1,7 @@
 
 DEAL:2d::Testing FE_DGQ<2>(1)
 DEAL:2d::Norm of result:          4.18
+DEAL:2d::Norm of result:          4.18
 DEAL:3d::Testing FE_DGQ<3>(1)
 DEAL:3d::Norm of result:          3.39
+DEAL:3d::Norm of result:          3.39

diff --git a/tests/matrix_free/matrix_vector_faces_common.h b/tests/matrix_free/matrix_vector_faces_common.h
index 052becdabd753fb5f28c1784a6c1631b4b538903..97fb5ad5de3980037d46c1e08fd14507a5968788 100644 (file)
--- a/tests/matrix_free/matrix_vector_faces_common.h
+++ b/tests/matrix_free/matrix_vector_faces_common.h
@@ -87,6 +87,32 @@ public:
               src);
   }
 
+  void
+  manual_loop_vmult(VectorType &dst, const VectorType &src) const
+  {
+    src.update_ghost_values();
+    dst = 0;
+    local_apply(data, dst, src, std::make_pair(0, data.n_cell_batches()));
+    local_apply_face(data,
+                     dst,
+                     src,
+                     std::make_pair(0, data.n_inner_face_batches()));
+    local_apply_boundary_face(data,
+                              dst,
+                              src,
+                              std::make_pair(data.n_inner_face_batches(),
+                                             data.n_inner_face_batches() +
+                                               data.n_boundary_face_batches()));
+    local_apply_face(data,
+                     dst,
+                     src,
+                     std::make_pair(data.n_inner_face_batches() +
+                                      data.n_boundary_face_batches(),
+                                    data.n_inner_face_batches() +
+                                      data.n_boundary_face_batches() +
+                                      data.n_ghost_inner_face_batches()));
+  }
+
 private:
   void
   local_apply(const MatrixFree<dim, number, VectorizedArrayType> &data,