From 188b386e0c2e3fbe2e2648cab86d8b65550b6157 Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Mon, 21 Dec 2020 13:30:43 +0100
Subject: [PATCH] MatrixFreeTools::compute_diagonal() for simplex grid

---
 tests/matrix_free/compute_diagonal_05.cc      | 93 +++++++++++++++++++
 ...ue.with_simplex_support=on.mpirun=1.output | 31 +++++++
 2 files changed, 124 insertions(+)
 create mode 100644 tests/matrix_free/compute_diagonal_05.cc
 create mode 100644 tests/matrix_free/compute_diagonal_05.with_mpi=true.with_simplex_support=on.mpirun=1.output

diff --git a/tests/matrix_free/compute_diagonal_05.cc b/tests/matrix_free/compute_diagonal_05.cc
new file mode 100644
index 0000000000..79a5ec2fbe
--- /dev/null
+++ b/tests/matrix_free/compute_diagonal_05.cc
@@ -0,0 +1,93 @@
+// ---------------------------------------------------------------------
+//
+// 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 MatrixFreeTools::compute_diagonal() for a Laplace operator on a simplex
+// mesh.
+
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include "compute_diagonal_util.h"
+
+template <int dim,
+          int n_components             = 1,
+          typename Number              = double,
+          typename VectorizedArrayType = VectorizedArray<Number>>
+void
+test()
+{
+  Triangulation<dim> tria;
+
+  GridGenerator::subdivided_hyper_cube_with_simplices(tria, 2);
+
+  const unsigned int fe_degree = 2;
+  const unsigned int n_points  = 3;
+
+  const Simplex::FE_P<dim> fe_q(fe_degree);
+  const FESystem<dim>      fe(fe_q, n_components);
+
+  // setup dof-handlers
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  AffineConstraints<Number> constraint;
+
+  MappingFE<dim> mapping(Simplex::FE_P<dim>{1});
+
+  VectorTools::interpolate_boundary_values(mapping,
+                                           dof_handler,
+                                           0,
+                                           Functions::ZeroFunction<dim>(
+                                             n_components),
+                                           constraint);
+
+  constraint.close();
+
+  typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
+    additional_data;
+  additional_data.mapping_update_flags = update_values | update_gradients;
+  Simplex::QGauss<dim> quad(fe_degree + 1);
+
+  MatrixFree<dim, Number, VectorizedArrayType> matrix_free;
+  matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data);
+
+
+  Test<dim, -1, 0, n_components, Number, VectorizedArrayType> test(
+    matrix_free,
+    constraint,
+    [](FEEvaluation<dim, -1, 0, n_components, Number, VectorizedArrayType>
+         &phi) {
+      phi.evaluate(false, true);
+      for (unsigned int q = 0; q < phi.n_q_points; ++q)
+        phi.submit_gradient(phi.get_gradient(q), q);
+      phi.integrate(false, true);
+    });
+
+  test.do_test();
+}
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    all;
+
+  test<2, 1>(); // scalar
+  test<2, 2>(); // vector
+}
diff --git a/tests/matrix_free/compute_diagonal_05.with_mpi=true.with_simplex_support=on.mpirun=1.output b/tests/matrix_free/compute_diagonal_05.with_mpi=true.with_simplex_support=on.mpirun=1.output
new file mode 100644
index 0000000000..c034e69e6a
--- /dev/null
+++ b/tests/matrix_free/compute_diagonal_05.with_mpi=true.with_simplex_support=on.mpirun=1.output
@@ -0,0 +1,31 @@
+
+Process #0
+Local range: [0, 25), global size: 25
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 5.333e+00 0.000e+00 4.000e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 5.333e+00 0.000e+00 5.333e+00 0.000e+00 0.000e+00 5.333e+00 0.000e+00 0.000e+00 0.000e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:0::15.6063
+Process #0
+Local range: [0, 25), global size: 25
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 5.333e+00 0.000e+00 4.000e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 5.333e+00 0.000e+00 5.333e+00 0.000e+00 0.000e+00 5.333e+00 0.000e+00 0.000e+00 0.000e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:0::15.6063
+Process #0
+Local range: [0, 25), global size: 25
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 5.333e+00 0.000e+00 4.000e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 5.333e+00 0.000e+00 5.333e+00 0.000e+00 0.000e+00 5.333e+00 0.000e+00 0.000e+00 0.000e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:0::15.6063
+Process #0
+Local range: [0, 50), global size: 50
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 4.000e+00 4.000e+00 5.333e+00 5.333e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 5.333e+00 5.333e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:0::22.0706
+Process #0
+Local range: [0, 50), global size: 50
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 4.000e+00 4.000e+00 5.333e+00 5.333e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 5.333e+00 5.333e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:0::22.0706
+Process #0
+Local range: [0, 50), global size: 50
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 4.000e+00 4.000e+00 5.333e+00 5.333e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 5.333e+00 5.333e+00 5.333e+00 5.333e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:0::22.0706
-- 
2.39.5