add test

diff --git a/tests/matrix_free/point_evaluation_29.cc b/tests/matrix_free/point_evaluation_29.cc
new file mode 100644 (file)
index 0000000..a0bc335
--- /dev/null
+++ b/tests/matrix_free/point_evaluation_29.cc
@@ -0,0 +1,221 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2023 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.
+ *
+ * ---------------------------------------------------------------------
+ */
+
+// Test FCL integration of FEFacePointEvaluation against old version.
+
+#include <deal.II/base/mpi.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/fe_point_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/non_matching/mapping_info.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <../tests.h>
+
+using namespace dealii;
+
+template <int dim, typename Number>
+void
+do_test(unsigned int degree)
+{
+  using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(2);
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(FESystem<dim>(FE_DGQ<dim>(degree), dim + 1));
+
+  AffineConstraints<Number> constraints;
+  constraints.close();
+
+  typename MatrixFree<dim, Number>::AdditionalData data;
+  data.mapping_update_flags                = update_values;
+  data.mapping_update_flags_inner_faces    = update_values;
+  data.mapping_update_flags_boundary_faces = update_values;
+
+  MappingQ1<dim> mapping;
+
+  MatrixFree<dim, Number> matrix_free;
+  matrix_free.reinit(
+    mapping, dof_handler, constraints, QGauss<dim>(degree + 1), data);
+
+
+  // setup NM::MappingInfo
+  constexpr unsigned int n_lanes = VectorizedArray<Number>::size();
+
+  std::vector<Quadrature<dim - 1>> quadrature_vector;
+  quadrature_vector.reserve((matrix_free.n_inner_face_batches() +
+                             matrix_free.n_boundary_face_batches()) *
+                            n_lanes);
+
+  std::vector<std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>>
+    vector_face_accessors;
+  vector_face_accessors.reserve((matrix_free.n_inner_face_batches() +
+                                 matrix_free.n_boundary_face_batches()) *
+                                n_lanes);
+
+  // fill container for inner face batches
+  unsigned int face_batch = 0;
+  for (unsigned int face_batch = 0;
+       face_batch < matrix_free.n_inner_face_batches() +
+                      matrix_free.n_boundary_face_batches();
+       ++face_batch)
+    {
+      for (unsigned int v = 0; v < n_lanes; ++v)
+        {
+          if (v < matrix_free.n_active_entries_per_face_batch(face_batch))
+            vector_face_accessors.push_back(
+              matrix_free.get_face_iterator(face_batch, v));
+          else
+            vector_face_accessors.push_back(
+              matrix_free.get_face_iterator(face_batch, 0));
+
+          quadrature_vector.emplace_back(Quadrature<dim - 1>(degree + 1));
+        }
+    }
+
+  NonMatching::MappingInfo<dim, dim, Number> nm_mapping_info(
+    mapping, update_values | update_JxW_values);
+  nm_mapping_info.reinit_faces(vector_face_accessors, quadrature_vector);
+
+  VectorType dst_1;
+  matrix_free.initialize_dof_vector(dst_1);
+  VectorType dst_2;
+  matrix_free.initialize_dof_vector(dst_2);
+  VectorType src;
+  matrix_free.initialize_dof_vector(src);
+
+
+  class InitialSolution : public Function<dim>
+  {
+  public:
+    InitialSolution()
+      : Function<dim>(dim + 1)
+    {}
+
+    double
+    value(const Point<dim> &p, const unsigned int comp) const final
+    {
+      if (comp < dim)
+        return p[comp] * p[comp] * (comp + 1);
+      else
+        return p[0] * p[0] * (comp + 1);
+    }
+  };
+
+  VectorTools::interpolate(*matrix_free.get_mapping_info().mapping,
+                           matrix_free.get_dof_handler(),
+                           InitialSolution(),
+                           src);
+
+  const auto boundary_function_1 = [&](const auto &data,
+                                       auto &      dst,
+                                       const auto &src,
+                                       const auto  face_range) {
+    FEFacePointEvaluation<dim, dim, dim, Number> phi_pnt(
+      nm_mapping_info, data.get_dof_handler().get_fe(), true, 1);
+    AlignedVector<Number> buffer(data.get_dof_handler().get_fe().dofs_per_cell);
+
+    for (unsigned int face = face_range.first; face < face_range.second; ++face)
+      {
+        for (unsigned int v = 0;
+             v < matrix_free.n_active_entries_per_face_batch(face);
+             ++v)
+          {
+            const auto &[cell, f] =
+              matrix_free.get_face_iterator(face, v, true);
+
+            phi_pnt.reinit(face * n_lanes + v);
+
+            cell->get_dof_values(src, buffer.begin(), buffer.end());
+            phi_pnt.evaluate(buffer, EvaluationFlags::values);
+
+            for (unsigned int q : phi_pnt.quadrature_point_indices())
+              phi_pnt.submit_value(phi_pnt.get_value(q), q);
+
+            phi_pnt.integrate(buffer, EvaluationFlags::values);
+
+            cell->distribute_local_to_global(buffer.begin(), buffer.end(), dst);
+          }
+      }
+  };
+
+  const auto boundary_function_2 =
+    [&](const auto &data, auto &dst, const auto &src, const auto face_range) {
+      FEFaceEvaluation<dim, -1, 0, dim, Number> phi(matrix_free, true, 0, 0, 1);
+      FEFacePointEvaluation<dim, dim, dim, Number> phi_pnt(
+        nm_mapping_info, data.get_dof_handler().get_fe(), true, 1);
+
+      for (unsigned int face = face_range.first; face < face_range.second;
+           ++face)
+        {
+          phi.reinit(face);
+          phi.read_dof_values(src);
+          phi.project_to_face(EvaluationFlags::values);
+          for (unsigned int v = 0;
+               v < matrix_free.n_active_entries_per_face_batch(face);
+               ++v)
+            {
+              phi_pnt.reinit(face * n_lanes + v);
+              phi_pnt.evaluate_in_face(&phi.get_scratch_data().begin()[0][v],
+                                       EvaluationFlags::values);
+
+              for (unsigned int q : phi_pnt.quadrature_point_indices())
+                phi_pnt.submit_value(phi_pnt.get_value(q), q);
+
+              phi_pnt.integrate_in_face(&phi.get_scratch_data().begin()[0][v],
+                                        EvaluationFlags::values);
+            }
+          phi.collect_from_face(EvaluationFlags::values);
+          phi.distribute_local_to_global(dst);
+        }
+    };
+
+
+  matrix_free.template loop<VectorType, VectorType>(
+    {}, {}, boundary_function_1, dst_1, src, true);
+
+  matrix_free.template loop<VectorType, VectorType>(
+    {}, {}, boundary_function_2, dst_2, src, true);
+
+  Assert(std::abs(dst_1.l2_norm() - dst_2.l2_norm()) < 1.0e-12);
+}
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+  MPILogInitAll                    all;
+
+  do_test<2, double>(2);
+
+  deallog << "OK" << std::endl;
+
+  return 0;
+}

diff --git a/tests/matrix_free/point_evaluation_29.with_p4est=true.mpirun=1.output b/tests/matrix_free/point_evaluation_29.with_p4est=true.mpirun=1.output
new file mode 100644 (file)
index 0000000..be8d055
--- /dev/null
+++ b/tests/matrix_free/point_evaluation_29.with_p4est=true.mpirun=1.output
@@ -0,0 +1,2 @@
+
+DEAL:0::OK