MGTransferGlobalCoarsening: test operator with ones on operator

diff --git a/tests/multigrid-global-coarsening/multigrid_a_01.cc b/tests/multigrid-global-coarsening/multigrid_a_01.cc
index f1df401b39c8f98f2854f1bde2e200c04a813863..d7bf93e9c43bfcbcb9df2fd9a63b58e28e3fcf39 100644 (file)
--- a/tests/multigrid-global-coarsening/multigrid_a_01.cc
+++ b/tests/multigrid-global-coarsening/multigrid_a_01.cc
@@ -25,7 +25,8 @@ template <int dim, typename Number = double>
 void
 test(const unsigned int n_refinements,
      const unsigned int fe_degree_fine,
-     const bool         do_simplex_mesh)
+     const bool         do_simplex_mesh,
+     const bool         ones_on_diagonal)
 {
   using VectorType = LinearAlgebra::distributed::Vector<Number>;
 
@@ -89,7 +90,12 @@ test(const unsigned int n_refinements,
       constraint.close();
 
       // set up operator
-      op.reinit(*mapping, dof_handler, *quad, constraint);
+      op.reinit(*mapping,
+                dof_handler,
+                *quad,
+                constraint,
+                numbers::invalid_unsigned_int,
+                ones_on_diagonal);
     }
 
   // set up transfer operator
@@ -161,11 +167,17 @@ main(int argc, char **argv)
 
   deallog.precision(8);
 
-  for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements)
-    for (unsigned int degree = 1; degree <= 4; ++degree)
-      test<2>(n_refinements, degree, false /*quadrilateral*/);
-
-  for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements)
-    for (unsigned int degree = 1; degree <= 2; ++degree)
-      test<2>(n_refinements, degree, true /*triangle*/);
+  for (const auto ones_on_diagonal : {false, true})
+    {
+      for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements)
+        for (unsigned int degree = 1; degree <= 4; ++degree)
+          test<2>(n_refinements,
+                  degree,
+                  false /*quadrilateral*/,
+                  ones_on_diagonal);
+
+      for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements)
+        for (unsigned int degree = 1; degree <= 2; ++degree)
+          test<2>(n_refinements, degree, true /*triangle*/, ones_on_diagonal);
+    }
 }

diff --git a/tests/multigrid-global-coarsening/multigrid_a_01.mpirun=1.with_trilinos=true.output b/tests/multigrid-global-coarsening/multigrid_a_01.mpirun=1.with_trilinos=true.output
index 07875b2dbeac1d95689033c025389a3d72d6b584..5a1eb699ca74c12d9c27b5ad9ca46df19064d6ea 100644 (file)
--- a/tests/multigrid-global-coarsening/multigrid_a_01.mpirun=1.with_trilinos=true.output
+++ b/tests/multigrid-global-coarsening/multigrid_a_01.mpirun=1.with_trilinos=true.output
@@ -17,3 +17,21 @@ DEAL:0::2 1 3 tri  3
 DEAL:0::2 2 3 tri  3
 DEAL:0::2 1 4 tri  3
 DEAL:0::2 2 4 tri  3
+DEAL:0::2 1 2 quad 2
+DEAL:0::2 2 2 quad 3
+DEAL:0::2 3 2 quad 3
+DEAL:0::2 4 2 quad 3
+DEAL:0::2 1 3 quad 3
+DEAL:0::2 2 3 quad 3
+DEAL:0::2 3 3 quad 3
+DEAL:0::2 4 3 quad 3
+DEAL:0::2 1 4 quad 3
+DEAL:0::2 2 4 quad 3
+DEAL:0::2 3 4 quad 3
+DEAL:0::2 4 4 quad 3
+DEAL:0::2 1 2 tri  3
+DEAL:0::2 2 2 tri  3
+DEAL:0::2 1 3 tri  3
+DEAL:0::2 2 3 tri  3
+DEAL:0::2 1 4 tri  3
+DEAL:0::2 2 4 tri  3

diff --git a/tests/multigrid-global-coarsening/multigrid_p_01.cc b/tests/multigrid-global-coarsening/multigrid_p_01.cc
index bf55c9ba8e760022efa377ebaf9849a3b97e2bb2..5ca1e98eebdca9dd40e4426c442ee836cb454dd6 100644 (file)
--- a/tests/multigrid-global-coarsening/multigrid_p_01.cc
+++ b/tests/multigrid-global-coarsening/multigrid_p_01.cc
@@ -26,7 +26,8 @@ void
 test(const unsigned int n_refinements,
      const unsigned int fe_degree_fine,
      const bool         do_simplex_mesh,
-     const unsigned int mesh_type)
+     const unsigned int mesh_type,
+     const bool         ones_on_diagonal)
 {
   using VectorType = LinearAlgebra::distributed::Vector<Number>;
 
@@ -117,7 +118,12 @@ test(const unsigned int n_refinements,
       constraint.close();
 
       // set up operator
-      op.reinit(*mapping, dof_handler, *quad, constraint);
+      op.reinit(*mapping,
+                dof_handler,
+                *quad,
+                constraint,
+                numbers::invalid_unsigned_int,
+                ones_on_diagonal);
     }
 
   // set up transfer operator
@@ -193,14 +199,26 @@ main(int argc, char **argv)
 
   deallog.precision(8);
 
-  for (unsigned int mesh_type = 0; mesh_type < 2; ++mesh_type)
-    {
-      for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements)
-        for (unsigned int degree = 2; degree <= 4; ++degree)
-          test<2>(n_refinements, degree, false /*quadrilateral*/, mesh_type);
-
-      for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements)
-        for (unsigned int degree = 2; degree <= 2; ++degree)
-          test<2>(n_refinements, degree, true /*triangle*/, mesh_type);
-    }
+  for (const auto ones_on_diagonal : {false, true})
+    for (unsigned int mesh_type = 0; mesh_type < 2; ++mesh_type)
+      {
+        for (unsigned int n_refinements = 2; n_refinements <= 4;
+             ++n_refinements)
+          for (unsigned int degree = 2; degree <= 4; ++degree)
+            test<2>(n_refinements,
+                    degree,
+                    false /*quadrilateral*/,
+                    mesh_type,
+                    ones_on_diagonal);
+
+
+        for (unsigned int n_refinements = 2; n_refinements <= 4;
+             ++n_refinements)
+          for (unsigned int degree = 2; degree <= 2; ++degree)
+            test<2>(n_refinements,
+                    degree,
+                    true /*triangle*/,
+                    mesh_type,
+                    ones_on_diagonal);
+      }
 }

diff --git a/tests/multigrid-global-coarsening/multigrid_p_01.mpirun=1.with_trilinos=true.output b/tests/multigrid-global-coarsening/multigrid_p_01.mpirun=1.with_trilinos=true.output
index f2980eab6f2767892d9e186e695fc5b61a2e630d..4637dd6777a702f548bd882af2b1c2c56d42d5f4 100644 (file)
--- a/tests/multigrid-global-coarsening/multigrid_p_01.mpirun=1.with_trilinos=true.output
+++ b/tests/multigrid-global-coarsening/multigrid_p_01.mpirun=1.with_trilinos=true.output
@@ -77,3 +77,81 @@ DEAL:0::Norm interpolated solution on level 1: 0.41924192
 DEAL:0::2 2 4 tri  4
 DEAL:0::Norm interpolated solution on level 0: 0.40718520
 DEAL:0::Norm interpolated solution on level 1: 0.80533381
+DEAL:0::2 2 2 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.33877875
+DEAL:0::Norm interpolated solution on level 1: 0.66014529
+DEAL:0::2 3 2 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.33872299
+DEAL:0::Norm interpolated solution on level 1: 0.99016130
+DEAL:0::2 4 2 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.33877890
+DEAL:0::Norm interpolated solution on level 1: 0.66010393
+DEAL:0::Norm interpolated solution on level 2: 1.3201957
+DEAL:0::2 2 3 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.66459543
+DEAL:0::Norm interpolated solution on level 1: 1.3203627
+DEAL:0::2 3 3 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.66458798
+DEAL:0::Norm interpolated solution on level 1: 1.9805368
+DEAL:0::2 4 3 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.66459725
+DEAL:0::Norm interpolated solution on level 1: 1.3203574
+DEAL:0::Norm interpolated solution on level 2: 2.6407133
+DEAL:0::2 2 4 quad 3
+DEAL:0::Norm interpolated solution on level 0: 1.3225826
+DEAL:0::Norm interpolated solution on level 1: 2.6407347
+DEAL:0::2 3 4 quad 3
+DEAL:0::Norm interpolated solution on level 0: 1.3225816
+DEAL:0::Norm interpolated solution on level 1: 3.9611011
+DEAL:0::2 4 4 quad 3
+DEAL:0::Norm interpolated solution on level 0: 1.3225830
+DEAL:0::Norm interpolated solution on level 1: 2.6407340
+DEAL:0::Norm interpolated solution on level 2: 5.2814679
+DEAL:0::2 2 2 tri  3
+DEAL:0::Norm interpolated solution on level 0: 0.33806168
+DEAL:0::Norm interpolated solution on level 1: 0.66011744
+DEAL:0::2 2 3 tri  3
+DEAL:0::Norm interpolated solution on level 0: 0.66445481
+DEAL:0::Norm interpolated solution on level 1: 1.3203595
+DEAL:0::2 2 4 tri  3
+DEAL:0::Norm interpolated solution on level 0: 1.3225585
+DEAL:0::Norm interpolated solution on level 1: 2.6407344
+DEAL:0::2 2 2 quad 2
+DEAL:0::Norm interpolated solution on level 0: 0.13157424
+DEAL:0::Norm interpolated solution on level 1: 0.24449632
+DEAL:0::2 3 2 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.13144458
+DEAL:0::Norm interpolated solution on level 1: 0.35587970
+DEAL:0::2 4 2 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.13139421
+DEAL:0::Norm interpolated solution on level 1: 0.24404614
+DEAL:0::Norm interpolated solution on level 2: 0.46558561
+DEAL:0::2 2 3 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.23622379
+DEAL:0::Norm interpolated solution on level 1: 0.45479680
+DEAL:0::2 3 3 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.23657515
+DEAL:0::Norm interpolated solution on level 1: 0.66773399
+DEAL:0::2 4 3 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.23652211
+DEAL:0::Norm interpolated solution on level 1: 0.45456995
+DEAL:0::Norm interpolated solution on level 2: 0.87733284
+DEAL:0::2 2 4 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.43119420
+DEAL:0::Norm interpolated solution on level 1: 0.85199621
+DEAL:0::2 3 4 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.43119824
+DEAL:0::Norm interpolated solution on level 1: 1.2615628
+DEAL:0::2 4 4 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.43118187
+DEAL:0::Norm interpolated solution on level 1: 0.85199072
+DEAL:0::Norm interpolated solution on level 2: 1.6669458
+DEAL:0::2 2 2 tri  3
+DEAL:0::Norm interpolated solution on level 0: 0.13575160
+DEAL:0::Norm interpolated solution on level 1: 0.24280231
+DEAL:0::2 2 3 tri  3
+DEAL:0::Norm interpolated solution on level 0: 0.21915908
+DEAL:0::Norm interpolated solution on level 1: 0.41924192
+DEAL:0::2 2 4 tri  4
+DEAL:0::Norm interpolated solution on level 0: 0.40718520
+DEAL:0::Norm interpolated solution on level 1: 0.80533381

diff --git a/tests/multigrid-global-coarsening/multigrid_util.h b/tests/multigrid-global-coarsening/multigrid_util.h
index af2b57b6d8a42e223a174416bec92899fa2f62ad..a0a65e1db2ef91f511725c8ae1c82bc792f4367f 100644 (file)
--- a/tests/multigrid-global-coarsening/multigrid_util.h
+++ b/tests/multigrid-global-coarsening/multigrid_util.h
@@ -75,7 +75,8 @@ public:
          const DoFHandler<dim>           &dof_handler,
          const Quadrature<dim>           &quad,
          const AffineConstraints<number> &constraints,
-         const unsigned int mg_level = numbers::invalid_unsigned_int)
+         const unsigned int mg_level         = numbers::invalid_unsigned_int,
+         const bool         ones_on_diagonal = false)
   {
     // Clear internal data structures (if operator is reused).
     this->system_matrix.clear();
@@ -92,6 +93,12 @@ public:
     data.mg_level             = mg_level;
 
     matrix_free.reinit(mapping, dof_handler, constraints, quad, data);
+
+    constrained_indices.clear();
+
+    if (ones_on_diagonal)
+      for (auto i : this->matrix_free.get_constrained_dofs())
+        constrained_indices.push_back(i);
   }
 
   virtual types::global_dof_index
@@ -122,6 +129,10 @@ public:
   {
     this->matrix_free.cell_loop(
       &Operator::do_cell_integral_range, this, dst, src, true);
+
+    for (unsigned int i = 0; i < constrained_indices.size(); ++i)
+      dst.local_element(constrained_indices[i]) =
+        src.local_element(constrained_indices[i]);
   }
 
   void
@@ -269,6 +280,8 @@ private:
   AffineConstraints<number> constraints;
 
   mutable TrilinosWrappers::SparseMatrix system_matrix;
+
+  std::vector<unsigned int> constrained_indices;
 };
 
 struct GMGParameters