]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add test: matrix_free/compute_diagonal_13
authorNils Much <nils.much@tum.de>
Tue, 8 Apr 2025 16:32:50 +0000 (18:32 +0200)
committerNils Much <nils.much@tum.de>
Thu, 10 Apr 2025 14:26:13 +0000 (16:26 +0200)
tests/matrix_free/compute_diagonal_13.cc [new file with mode: 0644]
tests/matrix_free/compute_diagonal_13.with_p4est=true.mpirun=1.output [new file with mode: 0644]
tests/matrix_free/compute_diagonal_util.h

diff --git a/tests/matrix_free/compute_diagonal_13.cc b/tests/matrix_free/compute_diagonal_13.cc
new file mode 100644 (file)
index 0000000..80d271f
--- /dev/null
@@ -0,0 +1,112 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2020 - 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+// Same as compute_diagonal_01 but with two different DoFHandlers attached to
+// MatrixFree, one of which is in hp-mode
+
+#include "compute_diagonal_util.h"
+
+template <int dim,
+          int fe_degree                = 1,
+          typename Number              = double,
+          typename VectorizedArrayType = VectorizedArray<Number>>
+void
+test()
+{
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+
+  GridGenerator::hyper_cube(tria, -numbers::PI / 2, numbers::PI / 2);
+
+  tria.refine_global();
+  for (auto &cell : tria.active_cell_iterators())
+    if (cell->is_active() && cell->is_locally_owned() &&
+        cell->center()[0] < 0.0)
+      cell->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+
+  FE_Q<dim> fe0(fe_degree);
+  FE_Q<dim> fe1(fe_degree + 1);
+
+  hp::FECollection<dim> fe_collection;
+  fe_collection.push_back(fe0);
+  fe_collection.push_back(fe1);
+
+  DoFHandler<dim> dof0(tria);
+  dof0.distribute_dofs(fe0);
+
+  DoFHandler<dim> dof1(tria);
+  unsigned int    counter = 0;
+  for (const auto &cell : dof1.active_cell_iterators())
+    cell->set_active_fe_index(++counter % 2);
+  dof1.distribute_dofs(fe_collection);
+
+  std::vector<const DoFHandler<dim> *> dof(2);
+  dof[0] = &dof0;
+  dof[1] = &dof1;
+
+  AffineConstraints<Number> constraint0;
+  DoFTools::make_hanging_node_constraints(*dof[0], constraint0);
+  VectorTools::interpolate_boundary_values(
+    *dof[0], 0, Functions::ZeroFunction<dim, Number>(), constraint0);
+  constraint0.close();
+
+  AffineConstraints<Number> constraint1;
+  DoFTools::make_hanging_node_constraints(*dof[1], constraint1);
+  VectorTools::interpolate_boundary_values(
+    *dof[1], 0, Functions::ZeroFunction<dim, Number>(), constraint1);
+  constraint1.close();
+
+  std::vector<const AffineConstraints<Number> *> constraints(2);
+  constraints[0] = &constraint0;
+  constraints[1] = &constraint1;
+
+
+  typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
+    additional_data;
+  additional_data.mapping_update_flags = update_values | update_gradients;
+
+  MappingQ<dim> mapping(1);
+
+  std::vector<Quadrature<dim>> quad;
+  for (unsigned int no = 0; no < 2; ++no)
+    quad.push_back(QGauss<dim>(fe_degree + 1 + no));
+
+  MatrixFree<dim, Number, VectorizedArrayType> matrix_free;
+  matrix_free.reinit(mapping, dof, constraints, quad, additional_data);
+
+
+  Test<dim, -1, 0, 1, Number, VectorizedArrayType> test(
+    matrix_free,
+    *constraints[0],
+    [](FEEvaluation<dim, -1, 0, 1, Number, VectorizedArrayType> &phi) {
+      phi.evaluate(EvaluationFlags::gradients);
+      for (unsigned int q = 0; q < phi.n_q_points; ++q)
+        phi.submit_gradient(phi.get_gradient(q), q);
+      phi.integrate(EvaluationFlags::gradients);
+    },
+    0,
+    0);
+
+  test.do_test();
+}
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    all;
+
+  test<2>(); // scalar
+}
diff --git a/tests/matrix_free/compute_diagonal_13.with_p4est=true.mpirun=1.output b/tests/matrix_free/compute_diagonal_13.with_p4est=true.mpirun=1.output
new file mode 100644 (file)
index 0000000..3dece03
--- /dev/null
@@ -0,0 +1,16 @@
+
+Process #0
+Local range: [0, 18), global size: 18
+Vector data:
+0.000e+00 0.000e+00 3.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:0::5.50757
+Process #0
+Local range: [0, 18), global size: 18
+Vector data:
+0.000e+00 0.000e+00 3.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:0::5.50757
+Process #0
+Local range: [0, 18), global size: 18
+Vector data:
+0.000e+00 0.000e+00 3.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:0::5.50757
index 615d87016d995c30a67738842f986fa1dddefbc5..563e1743680ddc6b01c9888a0352ff4a3d4ddb0d 100644 (file)
@@ -65,10 +65,14 @@ public:
                                              n_components,
                                              Number,
                                              VectorizedArrayType> &)>
-         &cell_operation)
+                         &cell_operation,
+       const unsigned int dof_no  = 0,
+       const unsigned int quad_no = 0)
     : matrix_free(matrix_free)
     , constraints(constraints)
     , cell_operation(cell_operation)
+    , dof_no(dof_no)
+    , quad_no(quad_no)
   {}
 
   void
@@ -84,12 +88,12 @@ public:
                              (Utilities::MPI::n_mpi_processes(
                                 matrix_free.get_task_info().communicator) == 1);
 
+    const auto &dof_handler = matrix_free.get_dof_handler(dof_no);
+
     if (test_matrix)
       {
-        DynamicSparsityPattern dsp(matrix_free.get_dof_handler().n_dofs());
-        DoFTools::make_sparsity_pattern(matrix_free.get_dof_handler(),
-                                        dsp,
-                                        constraints);
+        DynamicSparsityPattern dsp(dof_handler.n_dofs());
+        DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints);
         sparsity_pattern.copy_from(dsp);
         A1.reinit(sparsity_pattern);
         A2.reinit(sparsity_pattern);
@@ -99,16 +103,18 @@ public:
     double error_local_1, error_local_2, error_global;
 
     {
-      matrix_free.initialize_dof_vector(diagonal_global);
+      matrix_free.initialize_dof_vector(diagonal_global, dof_no);
       MatrixFreeTools::compute_diagonal<dim,
                                         fe_degree,
                                         n_points,
                                         n_components,
                                         Number,
                                         VectorizedArrayType>(
-        matrix_free, diagonal_global, [&](auto &phi) {
-          this->cell_operation(phi);
-        });
+        matrix_free,
+        diagonal_global,
+        [&](auto &phi) { this->cell_operation(phi); },
+        dof_no,
+        quad_no);
 
       diagonal_global.print(deallog.get_file_stream());
       error_local_1 = diagonal_global.l2_norm();
@@ -117,11 +123,13 @@ public:
 
     {
       VectorType diagonal_global;
-      matrix_free.initialize_dof_vector(diagonal_global);
+      matrix_free.initialize_dof_vector(diagonal_global, dof_no);
       MatrixFreeTools::compute_diagonal(matrix_free,
                                         diagonal_global,
                                         &Test::cell_function,
-                                        this);
+                                        this,
+                                        dof_no,
+                                        quad_no);
 
       diagonal_global.print(deallog.get_file_stream());
       error_local_2 = diagonal_global.l2_norm();
@@ -139,24 +147,32 @@ public:
                                         Number,
                                         VectorizedArrayType,
                                         SparseMatrix<Number>>(
-          matrix_free, constraints, A1, [&](auto &phi) {
-            this->cell_operation(phi);
-          });
+          matrix_free,
+          constraints,
+          A1,
+          [&](auto &phi) { this->cell_operation(phi); },
+          dof_no,
+          quad_no);
       }
 
     if (test_matrix)
       {
-        MatrixFreeTools::compute_matrix(
-          matrix_free, constraints, A2, &Test::cell_function, this);
+        MatrixFreeTools::compute_matrix(matrix_free,
+                                        constraints,
+                                        A2,
+                                        &Test::cell_function,
+                                        this,
+                                        dof_no,
+                                        quad_no);
       }
 
     // compute diagonal globally
     {
       VectorType src, temp;
 
-      matrix_free.initialize_dof_vector(src);
-      matrix_free.initialize_dof_vector(diagonal_global_reference);
-      matrix_free.initialize_dof_vector(temp);
+      matrix_free.initialize_dof_vector(src, dof_no);
+      matrix_free.initialize_dof_vector(diagonal_global_reference, dof_no);
+      matrix_free.initialize_dof_vector(temp, dof_no);
 
       for (unsigned int i = 0; i < src.size(); ++i)
         {
@@ -223,7 +239,7 @@ public:
                  n_components,
                  Number,
                  VectorizedArrayType>
-      phi(data, pair);
+      phi(data, pair, dof_no, quad_no);
     for (auto cell = pair.first; cell < pair.second; ++cell)
       {
         phi.reinit(cell);
@@ -252,5 +268,7 @@ public:
                                         n_components,
                                         Number,
                                         VectorizedArrayType> &)>
-    cell_operation;
+                     cell_operation;
+  const unsigned int dof_no;
+  const unsigned int quad_no;
 };

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.