From b216b97109f8467151044d8f68c482a0d5dea2fe Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Mon, 28 Dec 2020 11:59:23 +0100
Subject: [PATCH] FETools::get_projection_matrix for simplices

---
 include/deal.II/fe/fe_tools.templates.h       |  19 ++-
 tests/simplex/get_projection_matrix_01.cc     | 139 ++++++++++++++++++
 ...n_matrix_01.with_simplex_support=on.output |   3 +
 3 files changed, 157 insertions(+), 4 deletions(-)
 create mode 100644 tests/simplex/get_projection_matrix_01.cc
 create mode 100644 tests/simplex/get_projection_matrix_01.with_simplex_support=on.output

diff --git a/include/deal.II/fe/fe_tools.templates.h b/include/deal.II/fe/fe_tools.templates.h
index 7833e04cba..0d00c9062a 100644
--- a/include/deal.II/fe/fe_tools.templates.h
+++ b/include/deal.II/fe/fe_tools.templates.h
@@ -1598,23 +1598,34 @@ namespace FETools
     const unsigned int n1 = fe1.n_dofs_per_cell();
     const unsigned int n2 = fe2.n_dofs_per_cell();
 
+    const auto reference_cell_type = fe1.reference_cell_type();
+
+    Assert(fe1.reference_cell_type() == fe2.reference_cell_type(),
+           ExcNotImplemented());
+
     // First, create a local mass matrix for the unit cell
     Triangulation<dim, spacedim> tr;
-    GridGenerator::hyper_cube(tr);
+    ReferenceCell::make_triangulation(reference_cell_type, tr);
+
+    const auto &mapping =
+      ReferenceCell::get_default_linear_mapping<dim, spacedim>(
+        reference_cell_type);
 
     // Choose a Gauss quadrature rule that is exact up to degree 2n-1
     const unsigned int degree =
       std::max(fe1.tensor_degree(), fe2.tensor_degree());
     Assert(degree != numbers::invalid_unsigned_int, ExcNotImplemented());
-    const QGauss<dim> quadrature(degree + 1);
+    const auto quadrature =
+      ReferenceCell::get_gauss_type_quadrature<dim>(reference_cell_type,
+                                                    degree + 1);
 
     // Set up FEValues.
     const UpdateFlags flags =
       update_values | update_quadrature_points | update_JxW_values;
-    FEValues<dim> val1(fe1, quadrature, update_values);
+    FEValues<dim> val1(mapping, fe1, quadrature, update_values);
     val1.reinit(tr.begin_active());
 
-    FEValues<dim> val2(fe2, quadrature, flags);
+    FEValues<dim> val2(mapping, fe2, quadrature, flags);
     val2.reinit(tr.begin_active());
 
     // Integrate and invert mass matrix. This happens in the target space
diff --git a/tests/simplex/get_projection_matrix_01.cc b/tests/simplex/get_projection_matrix_01.cc
new file mode 100644
index 0000000000..095fee6f5e
--- /dev/null
+++ b/tests/simplex/get_projection_matrix_01.cc
@@ -0,0 +1,139 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Test FETools::get_projection_matrix for simplices.
+
+
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/lac/householder.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+class RightHandSideFunction : public Function<dim>
+{
+public:
+  RightHandSideFunction(const unsigned int n_components)
+    : Function<dim>(n_components)
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component = 0) const
+  {
+    if (component == 0)
+      return p[0];
+    else
+      return 0.0;
+  }
+};
+
+template <int dim, int spacedim = dim>
+void
+test()
+{
+  Simplex::FE_P<dim, spacedim> fe_coarse(1);
+  Simplex::FE_P<dim, spacedim> fe_fine(2);
+  MappingFE<dim>               mapping(Simplex::FE_P<dim>(1));
+
+  FullMatrix<double> matrix(fe_fine.n_dofs_per_cell(),
+                            fe_coarse.n_dofs_per_cell());
+  FETools::get_projection_matrix(fe_coarse, fe_fine, matrix);
+
+  const unsigned int n_refinements = 2;
+
+  Triangulation<dim, spacedim> tria;
+  GridGenerator::subdivided_hyper_cube_with_simplices(
+    tria, Utilities::pow(2, n_refinements));
+
+  DoFHandler<dim, spacedim> dof_handler_coarse(tria);
+  dof_handler_coarse.distribute_dofs(fe_coarse);
+
+  DoFHandler<dim, spacedim> dof_handler_fine(tria);
+  dof_handler_fine.distribute_dofs(fe_fine);
+
+  Vector<double> vec_coarse(dof_handler_coarse.n_dofs());
+  Vector<double> vec_fine(dof_handler_fine.n_dofs());
+
+  Vector<double> temp_coarse(fe_coarse.n_dofs_per_cell());
+  Vector<double> temp_fine(fe_fine.n_dofs_per_cell());
+
+  VectorTools::interpolate(mapping,
+                           dof_handler_coarse,
+                           RightHandSideFunction<dim>(1),
+                           vec_coarse);
+
+  for (const auto &cell_coarse : dof_handler_coarse.active_cell_iterators())
+    {
+      cell_coarse->get_dof_values(vec_coarse, temp_coarse);
+
+      matrix.vmult(temp_fine, temp_coarse);
+
+      DoFCellAccessor<dim, spacedim, false>(&tria,
+                                            cell_coarse->level(),
+                                            cell_coarse->index(),
+                                            &dof_handler_fine)
+        .set_dof_values(temp_fine, vec_fine);
+    }
+
+  vec_fine.print(deallog.get_file_stream());
+
+#if false
+  {
+    DataOut<dim> data_out;
+
+    data_out.attach_dof_handler(dof_handler_coarse);
+    data_out.add_data_vector(vec_coarse, "solution");
+
+    data_out.build_patches(mapping, 2);
+
+    std::ofstream output("test_coarse.vtk");
+    data_out.write_vtk(output);
+  }
+
+  {
+    DataOut<dim> data_out;
+
+    data_out.attach_dof_handler(dof_handler_fine);
+    data_out.add_data_vector(vec_fine, "solution");
+
+    data_out.build_patches(mapping, 2);
+
+    std::ofstream output("test_fine.vtk");
+    data_out.write_vtk(output);
+  }
+#endif
+}
+
+int
+main()
+{
+  initlog();
+
+  test<2>();
+  test<3>();
+}
diff --git a/tests/simplex/get_projection_matrix_01.with_simplex_support=on.output b/tests/simplex/get_projection_matrix_01.with_simplex_support=on.output
new file mode 100644
index 0000000000..4ccd48e0a2
--- /dev/null
+++ b/tests/simplex/get_projection_matrix_01.with_simplex_support=on.output
@@ -0,0 +1,3 @@
+
+-2.082e-17 2.500e-01 -2.082e-17 1.250e-01 1.250e-01 0.000e+00 2.500e-01 1.250e-01 2.500e-01 5.000e-01 3.750e-01 3.750e-01 5.000e-01 3.750e-01 5.000e-01 7.500e-01 6.250e-01 6.250e-01 7.500e-01 6.250e-01 7.500e-01 1.000e+00 8.750e-01 8.750e-01 1.000e+00 8.750e-01 1.000e+00 -2.082e-17 1.250e-01 0.000e+00 2.500e-01 1.250e-01 2.500e-01 3.750e-01 5.000e-01 3.750e-01 5.000e-01 6.250e-01 7.500e-01 6.250e-01 7.500e-01 8.750e-01 1.000e+00 8.750e-01 1.000e+00 -2.082e-17 1.250e-01 0.000e+00 2.500e-01 1.250e-01 2.500e-01 3.750e-01 5.000e-01 3.750e-01 5.000e-01 6.250e-01 7.500e-01 6.250e-01 7.500e-01 8.750e-01 1.000e+00 8.750e-01 1.000e+00 -1.110e-16 1.250e-01 0.000e+00 2.500e-01 1.250e-01 2.500e-01 3.750e-01 5.000e-01 3.750e-01 5.000e-01 6.250e-01 7.500e-01 6.250e-01 7.500e-01 8.750e-01 1.000e+00 8.750e-01 1.000e+00 
+-3.331e-16 2.500e-01 -7.772e-16 -7.772e-16 1.250e-01 1.250e-01 -1.890e-16 -5.095e-17 1.250e-01 -3.454e-18 2.500e-01 2.500e-01 2.500e-01 1.250e-01 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 1.250e-01 2.500e-01 -3.331e-16 -1.159e-16 -1.890e-16 1.250e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 3.750e-01 3.750e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 8.750e-01 1.000e+00 1.000e+00 8.750e-01 8.750e-01 8.750e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 -3.331e-16 -7.772e-16 1.250e-01 -1.890e-16 -5.439e-16 1.250e-01 -5.095e-17 -1.890e-16 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 8.750e-01 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 -7.772e-16 1.250e-01 -1.890e-16 -3.454e-18 2.500e-01 2.500e-01 2.500e-01 1.250e-01 1.250e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 -3.331e-16 -1.159e-16 -1.890e-16 1.250e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 8.750e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 -5.274e-16 -2.109e-15 1.250e-01 -1.890e-16 -5.439e-16 1.250e-01 -5.551e-17 -1.890e-16 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 8.750e-01 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 2.500e-01 1.250e-01 2.500e-01 2.500e-01 -7.772e-16 -5.439e-16 1.250e-01 -5.095e-17 -3.331e-16 -1.159e-16 1.250e-01 -1.890e-16 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 3.750e-01 5.000e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 7.500e-01 7.500e-01 6.250e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 8.750e-01 1.000e+00 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 1.000e+00 8.750e-01 -3.454e-18 2.500e-01 1.250e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 -3.331e-16 -1.159e-16 -1.890e-16 1.250e-01 5.000e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 2.500e-01 -7.772e-16 -5.439e-16 1.250e-01 -5.095e-17 -1.890e-16 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 5.000e-01 3.750e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 -3.454e-18 2.500e-01 1.250e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 -5.274e-16 -3.308e-16 -1.890e-16 1.250e-01 5.000e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 -7.772e-16 -5.095e-17 1.250e-01 -3.454e-18 2.500e-01 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 1.250e-01 2.500e-01 -3.331e-16 -1.159e-16 -1.890e-16 1.250e-01 5.000e-01 3.750e-01 5.000e-01 5.000e-01 3.750e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 6.250e-01 7.500e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 1.000e+00 1.000e+00 8.750e-01 8.750e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 2.500e-01 -7.772e-16 -5.439e-16 1.250e-01 -5.095e-17 -1.890e-16 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 5.000e-01 3.750e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 -3.454e-18 2.500e-01 1.250e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 -3.331e-16 -1.159e-16 -1.890e-16 1.250e-01 5.000e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 2.500e-01 -2.109e-15 -5.439e-16 1.250e-01 -5.551e-17 -1.890e-16 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 5.000e-01 3.750e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 2.500e-01 1.250e-01 2.500e-01 2.500e-01 -2.096e-15 -5.439e-16 1.250e-01 -5.095e-17 -8.882e-16 -1.159e-16 1.250e-01 -2.220e-17 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 3.750e-01 5.000e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 7.500e-01 7.500e-01 6.250e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 8.750e-01 1.000e+00 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 1.000e+00 8.750e-01 -3.454e-18 2.500e-01 1.250e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 -8.882e-16 -1.159e-16 -2.220e-17 1.250e-01 5.000e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 2.500e-01 -2.096e-15 -5.439e-16 1.250e-01 -5.095e-17 -2.220e-17 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 5.000e-01 3.750e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 -3.454e-18 2.500e-01 1.250e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 -1.499e-15 -3.308e-16 -2.220e-17 1.250e-01 5.000e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 
-- 
2.39.5