From: Peter Munch <peterrmuench@gmail.com>
Date: Tue, 1 Sep 2020 10:48:10 +0000 (+0200)
Subject: Work on interpolation from cell to face quadrature points
X-Git-Tag: v9.3.0-rc1~1137^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F10870%2Fhead;p=dealii.git

Work on interpolation from cell to face quadrature points
---

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index 4451027de4..aaad644911 100644
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -1686,24 +1686,77 @@ namespace internal
                 Number *                                      output,
                 const bool                                    do_gradients,
                 const unsigned int                            face_no)
+    {
+      Assert(static_cast<unsigned int>(fe_degree) ==
+                 data.data.front().fe_degree ||
+               fe_degree == -1,
+             ExcInternalError());
+
+      interpolate_generic<do_evaluate, add_into_output>(
+        input,
+        output,
+        do_gradients,
+        face_no,
+        data.data.front().fe_degree + 1,
+        data.data.front().shape_data_on_face,
+        data.dofs_per_component_on_cell,
+        2 * data.dofs_per_component_on_face);
+    }
+
+    /**
+     * Interpolate the values on the cell quadrature points onto a face.
+     */
+    template <bool do_evaluate, bool add_into_output>
+    static void
+    interpolate_quadrature(const MatrixFreeFunctions::ShapeInfo<Number> &data,
+                           const Number *                                input,
+                           Number *                                      output,
+                           const bool         do_gradients,
+                           const unsigned int face_no)
+    {
+      Assert(static_cast<unsigned int>(fe_degree + 1) ==
+                 data.data.front().quadrature.size() ||
+               fe_degree == -1,
+             ExcInternalError());
+
+      interpolate_generic<do_evaluate, add_into_output>(
+        input,
+        output,
+        do_gradients,
+        face_no,
+        data.data.front().quadrature.size(),
+        data.data.front().quadrature_data_on_face,
+        data.n_q_points,
+        data.n_q_points_face);
+    }
+
+  private:
+    template <bool do_evaluate, bool add_into_output>
+    static void
+    interpolate_generic(const Number *               input,
+                        Number *                     output,
+                        const bool                   do_gradients,
+                        const unsigned int           face_no,
+                        const unsigned int           n_points_1d,
+                        const AlignedVector<Number> *shape_data,
+                        const unsigned int           dofs_per_component_on_cell,
+                        const unsigned int           dofs_per_component_on_face)
     {
       internal::EvaluatorTensorProduct<internal::evaluate_general,
                                        dim,
                                        fe_degree + 1,
                                        0,
                                        Number>
-        evalf(data.data.front().shape_data_on_face[face_no % 2],
+        evalf(shape_data[face_no % 2],
               AlignedVector<Number>(),
               AlignedVector<Number>(),
-              data.data.front().fe_degree + 1,
+              n_points_1d,
               0);
 
-      const unsigned int in_stride = do_evaluate ?
-                                       data.dofs_per_component_on_cell :
-                                       2 * data.dofs_per_component_on_face;
-      const unsigned int out_stride = do_evaluate ?
-                                        2 * data.dofs_per_component_on_face :
-                                        data.dofs_per_component_on_cell;
+      const unsigned int in_stride =
+        do_evaluate ? dofs_per_component_on_cell : dofs_per_component_on_face;
+      const unsigned int out_stride =
+        do_evaluate ? dofs_per_component_on_face : dofs_per_component_on_cell;
       const unsigned int face_direction = face_no / 2;
       for (unsigned int c = 0; c < n_components; c++)
         {
diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h
index f20abbc98d..b30412de08 100644
--- a/include/deal.II/matrix_free/shape_info.h
+++ b/include/deal.II/matrix_free/shape_info.h
@@ -229,6 +229,18 @@ namespace internal
        */
       AlignedVector<Number> shape_data_on_face[2];
 
+      /**
+       * Collects all data of 1D nodal shape values (defined by the Lagrange
+       * polynomials in the points of the quadrature rule) evaluated at the
+       * point 0 and 1 (the vertices) in one data structure.
+       *
+       * This data structure can be used to interpolate from the cell to the
+       * face quadrature points.
+       *
+       * @note In contrast to shape_data_on_face, only the vales are evaluated.
+       */
+      AlignedVector<Number> quadrature_data_on_face[2];
+
       /**
        * Stores one-dimensional values of shape functions on subface. Since
        * there are two subfaces, store two variants.
diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index 8b12453bab..41a90beb54 100644
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -112,8 +112,10 @@ namespace internal
         univariate_shape_data.shape_gradients_collocation;
       auto &shape_hessians_collocation =
         univariate_shape_data.shape_hessians_collocation;
-      auto &inverse_shape_values  = univariate_shape_data.inverse_shape_values;
-      auto &shape_data_on_face    = univariate_shape_data.shape_data_on_face;
+      auto &inverse_shape_values = univariate_shape_data.inverse_shape_values;
+      auto &shape_data_on_face   = univariate_shape_data.shape_data_on_face;
+      auto &quadrature_data_on_face =
+        univariate_shape_data.quadrature_data_on_face;
       auto &values_within_subface = univariate_shape_data.values_within_subface;
       auto &gradients_within_subface =
         univariate_shape_data.gradients_within_subface;
@@ -298,6 +300,23 @@ namespace internal
             fe->shape_grad_grad(my_i, q_point)[0][0];
         }
 
+      if (n_q_points_1d < 200)
+        {
+          quadrature_data_on_face[0].resize(quad.size() * 3);
+          quadrature_data_on_face[1].resize(quad.size() * 3);
+
+          dealii::FE_DGQArbitraryNodes<1> fe_quad(quad);
+
+          for (unsigned int i = 0; i < quad.size(); ++i)
+            {
+              Point<1> q_point;
+              q_point[0]                    = 0;
+              quadrature_data_on_face[0][i] = fe_quad.shape_value(i, q_point);
+              q_point[0]                    = 1;
+              quadrature_data_on_face[1][i] = fe_quad.shape_value(i, q_point);
+            }
+        }
+
       // get gradient and Hessian transformation matrix for the polynomial
       // space associated with the quadrature rule (collocation space). We
       // need to avoid the case with more than a few hundreds of quadrature
@@ -772,6 +791,8 @@ namespace internal
         {
           memory +=
             MemoryConsumption::memory_consumption(shape_data_on_face[i]);
+          memory +=
+            MemoryConsumption::memory_consumption(quadrature_data_on_face[i]);
           memory +=
             MemoryConsumption::memory_consumption(values_within_subface[i]);
           memory +=
diff --git a/tests/matrix_free/interpolate_quadrature_01.cc b/tests/matrix_free/interpolate_quadrature_01.cc
new file mode 100644
index 0000000000..8110ae3fa4
--- /dev/null
+++ b/tests/matrix_free/interpolate_quadrature_01.cc
@@ -0,0 +1,157 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+
+// Tests internal::FEFaceNormalEvaluationImpl::interpolate_quadrature()
+// by comparing the results of FEFaceEvaluation::gather_evaluate().
+
+template <int dim>
+class ExactSolution : public Function<dim>
+{
+public:
+  ExactSolution()
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int /*component*/ = 0) const
+  {
+    return p[0] + p[1];
+  }
+};
+
+template <int dim,
+          int fe_degree,
+          int n_points                 = fe_degree + 1,
+          typename Number              = double,
+          typename VectorizedArrayType = VectorizedArray<Number>>
+void
+test(const unsigned int n_refinements = 1)
+{
+  using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria, 0, 1, true);
+
+  tria.refine_global(n_refinements);
+
+  FE_DGQ<dim>     fe(fe_degree);
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  MappingQ<dim> mapping(1);
+  QGauss<1>     quad(n_points);
+
+  AffineConstraints<Number> constraint;
+
+  using MF = MatrixFree<dim, Number, VectorizedArrayType>;
+
+  typename MF::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.hold_all_faces_to_owned_cells       = true;
+
+  MF matrix_free;
+  matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data);
+
+  VectorType src, dst;
+
+  matrix_free.initialize_dof_vector(src);
+  matrix_free.initialize_dof_vector(dst);
+
+  VectorTools::interpolate(dof_handler, ExactSolution<dim>(), src);
+
+  FEEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType> phi(
+    matrix_free);
+  FEFaceEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType>
+    phi_m(matrix_free, true);
+
+  matrix_free.template loop_cell_centric<VectorType, VectorType>(
+    [&](const auto &, auto &, const auto &src, const auto range) {
+      for (unsigned int cell = range.first; cell < range.second; ++cell)
+        {
+          phi.reinit(cell);
+
+          phi.gather_evaluate(src, true, false);
+
+          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
+               face++)
+            {
+              phi_m.reinit(cell, face);
+              phi_m.gather_evaluate(src, true, false);
+
+              AlignedVector<VectorizedArrayType> temp(
+                phi_m.static_dofs_per_cell);
+
+              internal::FEFaceNormalEvaluationImpl<dim,
+                                                   n_points - 1,
+                                                   1,
+                                                   VectorizedArrayType>::
+                template interpolate_quadrature<true, false>(
+                  matrix_free.get_shape_info(),
+                  phi.begin_values(),
+                  temp.data(),
+                  false,
+                  face);
+
+
+              for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+                {
+                  const auto u_cell = temp[q];
+                  const auto u_face = phi_m.get_value(q);
+
+                  for (unsigned int v = 0; v < VectorizedArray<double>::size();
+                       ++v)
+                    {
+                      Assert(std::abs(u_cell[v] - u_face[v]) < 1e-10,
+                             ExcMessage("Entries do not match!"));
+                    }
+                }
+            }
+        }
+    },
+    dst,
+    src);
+}
+
+int
+main()
+{
+  initlog();
+  test<2, 1, 2, double, VectorizedArray<double>>();
+  test<3, 1, 2, double, VectorizedArray<double>>();
+
+  deallog << "OK!" << std::endl;
+}
diff --git a/tests/matrix_free/interpolate_quadrature_01.output b/tests/matrix_free/interpolate_quadrature_01.output
new file mode 100644
index 0000000000..5cfb783b8f
--- /dev/null
+++ b/tests/matrix_free/interpolate_quadrature_01.output
@@ -0,0 +1,2 @@
+
+DEAL::OK!