Merge pull request #12768 from kronbichler/mf_no_inner_faces

MatrixFree: Do not build inner face data if not requested

include/deal.II/matrix_free/face_setup_internal.h

@@ -82,6 +82,7 @@ namespace internal
         const dealii::Triangulation<dim> &triangulation,
         const unsigned int                mg_level,
         const bool                        hold_all_faces_to_owned_cells,
+        const bool                        build_inner_faces,
         std::vector<std::pair<unsigned int, unsigned int>> &cell_levels);
 
       /**
@@ -165,6 +166,7 @@ namespace internal
       const dealii::Triangulation<dim> &triangulation,
       const unsigned int                mg_level,
       const bool                        hold_all_faces_to_owned_cells,
+      const bool                        build_inner_faces,
       std::vector<std::pair<unsigned int, unsigned int>> &cell_levels)
     {
       use_active_cells = mg_level == numbers::invalid_unsigned_int;
@@ -583,6 +585,8 @@ namespace internal
               if (dcell->at_boundary(f) && !dcell->has_periodic_neighbor(f))
                 face_is_owned[dcell->face(f)->index()] =
                   FaceCategory::locally_active_at_boundary;
+              else if (!build_inner_faces)
+                continue;
 
               // treat boundaries of cells of different refinement level
               // inside the domain in case of multigrid separately

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

@@ -1029,6 +1029,7 @@ namespace internal
     const std::vector<unsigned int> &cell_vectorization_category,
     const bool                       cell_vectorization_categories_strict,
     const bool                       do_face_integrals,
+    const bool                       build_inner_faces,
     const bool                       overlap_communication_computation,
     MatrixFreeFunctions::TaskInfo &  task_info,
     std::vector<std::pair<unsigned int, unsigned int>> &cell_level_index,
@@ -1041,6 +1042,7 @@ namespace internal
       face_setup.initialize(dof_handler[0]->get_triangulation(),
                             mg_level,
                             hold_all_faces_to_owned_cells,
+                            build_inner_faces,
                             cell_level_index);
 
     const unsigned int n_dof_handlers = dof_handler.size();
@@ -1740,6 +1742,7 @@ MatrixFree<dim, Number, VectorizedArrayType>::initialize_indices(
     additional_data.cell_vectorization_category,
     additional_data.cell_vectorization_categories_strict,
     do_face_integrals,
+    additional_data.mapping_update_flags_inner_faces != update_default,
     additional_data.overlap_communication_computation,
     task_info,
     cell_level_index,

tests/matrix_free/matrix_vector_faces_33.cc

@@ -0,0 +1,126 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 - 2020 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// tests matrix-free face evaluation when only boundary faces but not inner
+// faces are enabled. Otherwise the same test as matrix_vector_faces_05 (FE_Q,
+// Laplacian, weak imposition of Dirichlet boundary condition)
+
+#include <deal.II/base/function.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include "../tests.h"
+
+#include "matrix_vector_faces_common.h"
+
+template <int dim, int fe_degree>
+void
+test()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(5 - dim);
+
+  FE_Q<dim>       fe(fe_degree);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+  AffineConstraints<double> constraints;
+  constraints.close();
+
+  deallog << "Testing " << dof.get_fe().get_name();
+  deallog << std::endl;
+  // std::cout << "Number of cells: " <<
+  // dof.get_triangulation().n_active_cells() << std::endl; std::cout << "Number
+  // of degrees of freedom: " << dof.n_dofs() << std::endl; std::cout << "Number
+  // of constraints: " << constraints.n_constraints() << std::endl;
+
+  MappingQ<dim> mapping(dof.get_fe().degree + 1);
+
+  Vector<double> in(dof.n_dofs()), out(dof.n_dofs());
+  Vector<double> out_dist(out);
+
+  // Set random seed for reproducibility
+  Testing::srand(42);
+  for (unsigned int i = 0; i < dof.n_dofs(); ++i)
+    {
+      if (constraints.is_constrained(i))
+        continue;
+      const double entry = Testing::rand() / (double)RAND_MAX;
+      in(i)              = entry;
+    }
+
+  constexpr unsigned int n_q_points_1d = fe_degree + 1;
+
+  // assemble sparse matrix with MeshWorker
+  SparsityPattern      sparsity;
+  SparseMatrix<double> matrix;
+  {
+    DynamicSparsityPattern d_sparsity(dof.n_dofs());
+    DoFTools::make_flux_sparsity_pattern(dof, d_sparsity);
+    sparsity.copy_from(d_sparsity);
+  }
+  matrix.reinit(sparsity);
+  MeshWorker::IntegrationInfoBox<dim> info_box;
+  UpdateFlags                         update_flags =
+    update_values | update_gradients | update_jacobians;
+  info_box.add_update_flags_all(update_flags);
+  info_box.initialize_gauss_quadrature(n_q_points_1d,
+                                       n_q_points_1d,
+                                       n_q_points_1d);
+  info_box.initialize(dof.get_fe(), mapping);
+
+  MeshWorker::DoFInfo<dim> dof_info(dof);
+
+  MeshWorker::Assembler::MatrixSimple<SparseMatrix<double>> assembler;
+  assembler.initialize(matrix);
+
+  MatrixIntegrator<dim> integrator;
+  MeshWorker::integration_loop<dim, dim>(
+    dof.begin_active(), dof.end(), dof_info, info_box, integrator, assembler);
+
+  matrix.vmult(out, in);
+
+  // zero constrained dofs
+  for (unsigned int i = 0; i < dof.n_dofs(); ++i)
+    if (constraints.is_constrained(i))
+      out(i) = 0;
+
+  MatrixFree<dim, double> mf_data;
+  const QGauss<1>         quad(n_q_points_1d > 0 ? n_q_points_1d :
+                                                   dof.get_fe().degree + 1);
+  typename MatrixFree<dim, double>::AdditionalData data;
+  data.tasks_parallel_scheme = MatrixFree<dim, double>::AdditionalData::none;
+  data.tasks_block_size      = 3;
+  data.mapping_update_flags_boundary_faces =
+    (update_gradients | update_JxW_values);
+
+  mf_data.reinit(mapping, dof, constraints, quad, data);
+
+  // Check that only boundary faces are set up as requested
+  Assert(mf_data.n_inner_face_batches() == 0, ExcInternalError());
+  Assert(mf_data.n_boundary_face_batches() > 0, ExcInternalError());
+
+  MatrixFreeTest<dim, fe_degree, n_q_points_1d, double, Vector<double>, 1> mf(
+    mf_data);
+  mf.vmult(out_dist, in);
+
+  out_dist -= out;
+  const double diff_norm = out_dist.linfty_norm() / out.linfty_norm();
+  deallog << "Norm of difference:          " << diff_norm << std::endl;
+
+  deallog << std::endl;
+}

tests/matrix_free/matrix_vector_faces_33.output

@@ -0,0 +1,13 @@
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Norm of difference:          0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(2)
+DEAL:2d::Norm of difference:          0
+DEAL:2d::
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Norm of difference:          0
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Norm of difference:          0
+DEAL:3d::