From: Daniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Date: Tue, 12 Jun 2018 22:48:33 +0000 (+0200)
Subject: Allow using FEFaceEvaluation with degree=-1
X-Git-Tag: v9.1.0-rc1~1046^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a92c2aeff54484d0113083866ae3d856ad7b08b2;p=dealii.git

Allow using FEFaceEvaluation with degree=-1
---

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index a5f081b826..a4762c56ba 100644
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -7124,17 +7124,22 @@ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
                   const bool        evaluate_gradients)
 {
   const unsigned int side = this->face_no % 2;
-  const unsigned int dofs_per_face =
+
+  const unsigned int static_dofs_per_face =
     fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
+                     numbers::invalid_unsigned_int;
+  const unsigned int dofs_per_face =
+    fe_degree > -1 ? static_dofs_per_face :
                      Utilities::pow(this->data->fe_degree + 1, dim - 1);
 
   constexpr unsigned int stack_array_size_threshold = 100;
 
-  VectorizedArray<Number> temp_data[dofs_per_face < stack_array_size_threshold ?
-                                      n_components_ * 2 * dofs_per_face :
-                                      1];
+  VectorizedArray<Number>
+    temp_data[static_dofs_per_face < stack_array_size_threshold ?
+                n_components_ * 2 * dofs_per_face :
+                1];
   VectorizedArray<Number> *__restrict temp1;
-  if (dofs_per_face < stack_array_size_threshold)
+  if (static_dofs_per_face < stack_array_size_threshold)
     temp1 = &temp_data[0];
   else
     temp1 = this->scratch_data;
diff --git a/tests/matrix_free/faces_value_optimization_02.cc b/tests/matrix_free/faces_value_optimization_02.cc
new file mode 100644
index 0000000000..1b82a5c249
--- /dev/null
+++ b/tests/matrix_free/faces_value_optimization_02.cc
@@ -0,0 +1,264 @@
+// ---------------------------------------------------------------------
+//
+// 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// tests matrix-free face evaluation with the option to compress
+// This is the same as face_value_optimization but uses fe_degree=-1
+
+#include <deal.II/base/function.h>
+
+#include <deal.II/fe/fe_dgq.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/constraint_matrix.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include "../tests.h"
+
+std::ofstream logfile("output");
+
+template <int dim, int fe_degree, typename number>
+class MatrixFreeTest
+{
+public:
+  MatrixFreeTest(const MatrixFree<dim, number> &data) : data(data), error(2)
+  {}
+
+  ~MatrixFreeTest()
+  {
+    deallog << "Error between variants: " << error[0] / error[1] << std::endl;
+  }
+
+  void
+  check_error(Vector<number> &src) const
+  {
+    data.loop(&MatrixFreeTest::local_apply,
+              &MatrixFreeTest::local_apply_face,
+              &MatrixFreeTest::local_apply_face,
+              this,
+              src,
+              src);
+  }
+
+private:
+  void
+  local_apply(const MatrixFree<dim, number> &,
+              Vector<number> &,
+              const Vector<number> &,
+              const std::pair<unsigned int, unsigned int> &) const
+  {}
+
+  void
+  local_apply_face(
+    const MatrixFree<dim, number> &data,
+    Vector<number> &,
+    const Vector<number> &                       src,
+    const std::pair<unsigned int, unsigned int> &face_range) const
+  {
+    FEFaceEvaluation<dim, -1, 0, 1, number> ref(data, true);
+    FEFaceEvaluation<dim, -1, 0, 1, number> check(data, true);
+
+    for (unsigned int face = face_range.first; face < face_range.second; face++)
+      {
+        ref.reinit(face);
+        check.reinit(face);
+
+        ref.read_dof_values(src);
+        ref.evaluate(true, false);
+        check.gather_evaluate(src, true, false);
+
+        for (unsigned int q = 0; q < ref.n_q_points; ++q)
+          {
+            VectorizedArray<number> diff =
+              (ref.get_value(q) - check.get_value(q));
+            for (unsigned int v = 0;
+                 v < VectorizedArray<number>::n_array_elements;
+                 ++v)
+              {
+                if (std::abs(diff[v]) > 1e-12)
+                  {
+                    deallog << "Error detected on face" << face << ", v=" << v
+                            << "!" << std::endl;
+                    deallog << "ref: ";
+                    for (unsigned int i = 0; i < ref.dofs_per_cell; ++i)
+                      deallog << ref.get_dof_value(i)[v] << " ";
+                    deallog << std::endl;
+                    deallog << "done: " << check.get_value(q)[v]
+                            << " instead of " << ref.get_value(q)[v]
+                            << std::endl;
+
+                    deallog
+                      << data.get_face_info(face).cells_interior[v] << " "
+                      << (int)data.get_face_info(face).interior_face_no << " "
+                      << (int)data.get_face_info(face).face_orientation << " "
+                      << (int)data.get_face_info(face).subface_index
+                      << std::endl;
+                    deallog << std::endl;
+                  }
+                error[0] += std::abs(diff[v]);
+                error[1] += std::abs(ref.get_value(q)[v]);
+              }
+          }
+      }
+
+    FEFaceEvaluation<dim, fe_degree, fe_degree + 1, 1, number> refr(data,
+                                                                    false);
+    FEFaceEvaluation<dim, fe_degree, fe_degree + 1, 1, number> checkr(data,
+                                                                      false);
+
+    for (unsigned int face = face_range.first;
+         face < std::min(data.n_inner_face_batches(), face_range.second);
+         face++)
+      {
+        refr.reinit(face);
+        checkr.reinit(face);
+
+        refr.read_dof_values(src);
+        refr.evaluate(true, false);
+        checkr.gather_evaluate(src, true, false);
+
+        for (unsigned int q = 0; q < ref.n_q_points; ++q)
+          {
+            VectorizedArray<number> diff =
+              (refr.get_value(q) - checkr.get_value(q));
+            for (unsigned int v = 0;
+                 v < VectorizedArray<number>::n_array_elements;
+                 ++v)
+              {
+                if (std::abs(diff[v]) > 1e-12)
+                  {
+                    deallog << "Error detected on face" << face << ", v=" << v
+                            << "!" << std::endl;
+                    deallog << "ref: ";
+                    for (unsigned int i = 0; i < ref.dofs_per_cell; ++i)
+                      deallog << refr.get_dof_value(i)[v] << " ";
+                    deallog << std::endl;
+                    deallog << "done: " << check.get_value(q)[v]
+                            << " instead of " << ref.get_value(q)[v]
+                            << std::endl;
+
+                    deallog
+                      << data.get_face_info(face).cells_exterior[v] << " "
+                      << (int)data.get_face_info(face).exterior_face_no << " "
+                      << (int)data.get_face_info(face).face_orientation << " "
+                      << (int)data.get_face_info(face).subface_index
+                      << std::endl;
+                    deallog << std::endl;
+                  }
+                error[0] += std::abs(diff[v]);
+                error[1] += std::abs(refr.get_value(q)[v]);
+              }
+          }
+      }
+  }
+
+  const MatrixFree<dim, number> &data;
+  mutable std::vector<double>    error;
+};
+
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_ball(tria);
+  tria.begin_active()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  typename Triangulation<dim>::active_cell_iterator cell, endc;
+  cell = tria.begin_active();
+  endc = tria.end();
+  for (; cell != endc; ++cell)
+    if (cell->center().norm() < 0.5)
+      cell->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  tria.begin(tria.n_levels() - 1)->set_refine_flag();
+  tria.last()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  tria.refine_global(1);
+  cell = tria.begin_active();
+  for (unsigned int i = 0; i < 10 - 3 * dim; ++i)
+    {
+      cell                 = tria.begin_active();
+      endc                 = tria.end();
+      unsigned int counter = 0;
+      for (; cell != endc; ++cell, ++counter)
+        if (counter % (7 - i) == 0)
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+    }
+
+  FE_DGQ<dim>     fe(fe_degree);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+  ConstraintMatrix constraints;
+  constraints.close();
+
+  typedef double number;
+
+  MatrixFree<dim, number> mf_data;
+  {
+    const QGauss<1>                                  quad(fe_degree + 1);
+    typename MatrixFree<dim, number>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+    data.tasks_block_size      = 3;
+    data.mapping_update_flags_inner_faces =
+      (update_gradients | update_JxW_values);
+    data.mapping_update_flags_boundary_faces =
+      (update_gradients | update_JxW_values);
+
+    mf_data.reinit(dof, constraints, quad, data);
+  }
+
+  MatrixFreeTest<dim, fe_degree, number> mf(mf_data);
+
+  Vector<number> in(dof.n_dofs());
+
+  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;
+    }
+
+  mf.check_error(in);
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  deallog << std::setprecision(3);
+
+  {
+    deallog.push("2d");
+    test<2, 1>();
+    test<2, 2>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3, 1>();
+    test<3, 2>();
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/faces_value_optimization_02.output b/tests/matrix_free/faces_value_optimization_02.output
new file mode 100644
index 0000000000..0bb296849a
--- /dev/null
+++ b/tests/matrix_free/faces_value_optimization_02.output
@@ -0,0 +1,5 @@
+
+DEAL:2d::Error between variants: 0
+DEAL:2d::Error between variants: 0
+DEAL:3d::Error between variants: 0
+DEAL:3d::Error between variants: 0