Merge pull request #13058 from peterrum/ecl_pbc_orientation

MatrixFree: fix orientation for cell-centric loop and pbc

include/deal.II/matrix_free/evaluation_kernels.h

@@ -3996,6 +3996,8 @@ namespace internal
                               proc.value(
                                 temp1[reorientate(v, i)][v],
                                 vector_ptrs[v][index_array_nodal[v][i]]);
+                            else if (integrate == false)
+                              temp1[i][v] = 0.0;
                         }
                   }
               }

include/deal.II/matrix_free/fe_evaluation.h

@@ -8919,7 +8919,7 @@ FEFaceEvaluation<dim,
                   // invert (see also:
                   // Triangulation::update_periodic_face_map())
                   static const std::array<unsigned int, 8> table{
-                    {0, 1, 0, 3, 6, 5, 4, 7}};
+                    {0, 1, 2, 3, 6, 5, 4, 7}};
 
                   face_orientation = table[face_orientation];
                 }

tests/matrix_free/pbc_orientation_02.cc

@@ -0,0 +1,217 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Like pbc_orientation_02 but with cell-centric loops.
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/operators.h>
+#include <deal.II/matrix_free/tools.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+template <int dim>
+class ExactSolution : public Function<dim>
+{
+public:
+  ExactSolution()
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int /*component*/ = 0) const
+  {
+    // independent of p[2] since we apply PBC in z-direction
+    return p[0] + p[1];
+  }
+};
+
+
+
+void
+generate_grid(Triangulation<3> &triangulation, int orientation)
+{
+  Point<3>              vertices_1[] = {Point<3>(-0., -0., -0.),
+                           Point<3>(+1., -0., -0.),
+                           Point<3>(-0., +1., -0.),
+                           Point<3>(+1., +1., -0.),
+                           Point<3>(-0., -0., +0.5),
+                           Point<3>(+1., -0., +0.5),
+                           Point<3>(-0., +1., +0.5),
+                           Point<3>(+1., +1., +0.5),
+                           Point<3>(-0., -0., +1.),
+                           Point<3>(+1., -0., +1.),
+                           Point<3>(-0., +1., +1.),
+                           Point<3>(+1., +1., +1.)};
+  std::vector<Point<3>> vertices(&vertices_1[0], &vertices_1[12]);
+
+  std::vector<CellData<3>> cells(2, CellData<3>());
+
+  /* cell 0 */
+  int cell_vertices_0[GeometryInfo<3>::vertices_per_cell] = {
+    0, 1, 2, 3, 4, 5, 6, 7};
+
+  /* cell 1 */
+  int cell_vertices_1[8][GeometryInfo<3>::vertices_per_cell] = {
+    {4, 5, 6, 7, 8, 9, 10, 11},
+    {5, 7, 4, 6, 9, 11, 8, 10},
+    {7, 6, 5, 4, 11, 10, 9, 8},
+    {6, 4, 7, 5, 10, 8, 11, 9},
+    {9, 8, 11, 10, 5, 4, 7, 6},
+    {8, 10, 9, 11, 4, 6, 5, 7},
+    {10, 11, 8, 9, 6, 7, 4, 5},
+    {11, 9, 10, 8, 7, 5, 6, 4}};
+
+  for (const unsigned int j : GeometryInfo<3>::vertex_indices())
+    {
+      cells[orientation < 8 ? 0 : 1].vertices[j] = cell_vertices_0[j];
+      cells[orientation < 8 ? 1 : 0].vertices[j] =
+        cell_vertices_1[orientation % 8][j];
+    }
+  cells[0].material_id = 0;
+  cells[1].material_id = 0;
+
+
+  triangulation.create_triangulation(vertices, cells, SubCellData());
+}
+
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+  if (dim == 2)
+    return;
+
+  Triangulation<3> tria;
+  for (unsigned int orientation = 0; orientation < 16; ++orientation)
+    {
+      deallog << "Testing orientation case " << orientation << std::endl;
+      tria.clear();
+      generate_grid(tria, orientation);
+
+      FE_DGQ<3>     fe(fe_degree);
+      DoFHandler<3> dof_handler(tria);
+      dof_handler.distribute_dofs(fe);
+
+      for (auto &cell : tria)
+        for (auto &face : cell.face_iterators())
+          {
+            if (std::abs(face->center()[2] - 0.0) < 10e-6)
+              face->set_boundary_id(4);
+            if (std::abs(face->center()[2] - 1.0) < 10e-6)
+              face->set_boundary_id(5);
+          }
+
+      std::vector<dealii::GridTools::PeriodicFacePair<
+        typename dealii::Triangulation<dim>::cell_iterator>>
+        periodic_faces;
+
+      dealii::GridTools::collect_periodic_faces(tria, 4, 5, 2, periodic_faces);
+
+      tria.add_periodicity(periodic_faces);
+
+      QGauss<1> quadrature(fe_degree + 1);
+      typename MatrixFree<dim, double>::AdditionalData additional_data;
+      additional_data.mapping_update_flags                = update_values;
+      additional_data.mapping_update_flags_inner_faces    = update_values;
+      additional_data.mapping_update_flags_boundary_faces = update_values;
+      additional_data.mapping_update_flags_faces_by_cells = update_values;
+      additional_data.tasks_parallel_scheme =
+        MatrixFree<dim, double>::AdditionalData::none;
+
+      MatrixFreeTools::categorize_by_boundary_ids(
+        dof_handler.get_triangulation(), additional_data);
+
+      MatrixFree<dim, double> data;
+
+      AffineConstraints<double> dummy;
+      dummy.close();
+      data.reinit(dof_handler, dummy, quadrature, additional_data);
+
+      using VectorType = LinearAlgebra::distributed::Vector<double>;
+
+      VectorType vec;
+      data.initialize_dof_vector(vec);
+
+      VectorTools::interpolate(dof_handler, ExactSolution<dim>(), vec);
+
+
+
+      data.template loop_cell_centric<VectorType, VectorType>(
+        [](const auto &data, auto &, const auto &src, const auto cell_range) {
+          FEFaceEvaluation<dim, fe_degree, fe_degree + 1, 1, double> eval_minus(
+            data, true);
+          FEFaceEvaluation<dim, fe_degree, fe_degree + 1, 1, double> eval_plus(
+            data, false);
+
+          for (unsigned int cell = cell_range.first; cell < cell_range.second;
+               ++cell)
+            for (unsigned int face = 0;
+                 face < GeometryInfo<dim>::faces_per_cell;
+                 ++face)
+              {
+                const auto boundary_ids =
+                  data.get_faces_by_cells_boundary_id(cell, face);
+                const auto boundary_id = boundary_ids[0];
+
+                if (boundary_id != numbers::internal_face_boundary_id)
+                  continue;
+
+                eval_minus.reinit(cell, face);
+                eval_minus.gather_evaluate(src, true, false);
+                eval_plus.reinit(cell, face);
+                eval_plus.gather_evaluate(src, true, false);
+
+                for (unsigned int q = 0; q < eval_minus.n_q_points; ++q)
+                  {
+                    const auto u_minus = eval_minus.get_value(q);
+                    const auto u_plus  = eval_plus.get_value(q);
+
+                    for (unsigned int v = 0;
+                         v < VectorizedArray<double>::size();
+                         ++v)
+                      {
+                        Assert(std::abs(u_minus[v] - u_plus[v]) < 1e-10,
+                               ExcMessage("Entries do not match!"));
+                      }
+                  }
+              }
+        },
+        vec,
+        vec);
+    }
+}
+
+int
+main()
+{
+  initlog();
+  test<3, 3>();
+}

tests/matrix_free/pbc_orientation_02.output

@@ -0,0 +1,17 @@
+
+DEAL::Testing orientation case 0
+DEAL::Testing orientation case 1
+DEAL::Testing orientation case 2
+DEAL::Testing orientation case 3
+DEAL::Testing orientation case 4
+DEAL::Testing orientation case 5
+DEAL::Testing orientation case 6
+DEAL::Testing orientation case 7
+DEAL::Testing orientation case 8
+DEAL::Testing orientation case 9
+DEAL::Testing orientation case 10
+DEAL::Testing orientation case 11
+DEAL::Testing orientation case 12
+DEAL::Testing orientation case 13
+DEAL::Testing orientation case 14
+DEAL::Testing orientation case 15