Fix hanging node constraints

diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h
index 21531c3785bf71ef8d6a24999ffdf4b6f591fa9a..a6fd34cde20cbc17b0073a97f1d15294e7ffab5f 100644 (file)
--- a/include/deal.II/matrix_free/tools.h
+++ b/include/deal.II/matrix_free/tools.h
@@ -25,6 +25,7 @@
 #include <deal.II/matrix_free/portable_matrix_free.h>
 #include <deal.II/matrix_free/vector_access_internal.h>
 
+
 DEAL_II_NAMESPACE_OPEN
 
 /**
@@ -1375,7 +1376,7 @@ namespace MatrixFreeTools
     {}
 
     KOKKOS_FUNCTION void
-    operator()(const unsigned int /*cell*/,
+    operator()(const unsigned int                                      cell,
                const typename Portable::MatrixFree<dim, Number>::Data *gpu_data,
                Portable::SharedData<dim, Number> *shared_data,
                const Number *,
@@ -1384,15 +1385,31 @@ namespace MatrixFreeTools
       Portable::FEEvaluation<dim, fe_degree, fe_degree + 1, 1, Number> fe_eval(
         gpu_data, shared_data);
       constexpr int dofs_per_cell = decltype(fe_eval)::tensor_dofs_per_cell;
-      double        diagonal[dofs_per_cell] = {};
+      Number        diagonal[dofs_per_cell] = {};
       for (unsigned int i = 0; i < dofs_per_cell; ++i)
         {
           Kokkos::parallel_for(
             Kokkos::TeamThreadRange(shared_data->team_member, dofs_per_cell),
             [&](int j) { fe_eval.submit_dof_value(i == j ? 1 : 0, j); });
+
+          Portable::internal::
+            resolve_hanging_nodes<dim, fe_degree, false, Number>(
+              shared_data->team_member,
+              gpu_data->constraint_weights,
+              gpu_data->constraint_mask(cell),
+              Kokkos::subview(shared_data->values, Kokkos::ALL, 0));
+
           fe_eval.evaluate(m_evaluation_flags);
           fe_eval.apply_for_each_quad_point(m_quad_operation);
           fe_eval.integrate(m_integration_flags);
+
+          Portable::internal::
+            resolve_hanging_nodes<dim, fe_degree, true, Number>(
+              shared_data->team_member,
+              gpu_data->constraint_weights,
+              gpu_data->constraint_mask(cell),
+              Kokkos::subview(shared_data->values, Kokkos::ALL, 0));
+
           Kokkos::single(Kokkos::PerTeam(shared_data->team_member),
                          [&] { diagonal[i] = fe_eval.get_dof_value(i); });
         }
@@ -1401,7 +1418,27 @@ namespace MatrixFreeTools
         for (unsigned int i = 0; i < dofs_per_cell; ++i)
           fe_eval.submit_dof_value(diagonal[i], i);
       });
-      fe_eval.distribute_local_to_global(dst);
+
+      // We need to do the same as distribute_local_to_global but without
+      // constraints since we have already taken care of them earlier
+      if (gpu_data->use_coloring)
+        {
+          Kokkos::parallel_for(Kokkos::TeamThreadRange(shared_data->team_member,
+                                                       dofs_per_cell),
+                               [&](const int &i) {
+                                 dst[gpu_data->local_to_global(cell, i)] +=
+                                   shared_data->values(i, 0);
+                               });
+        }
+      else
+        {
+          Kokkos::parallel_for(
+            Kokkos::TeamThreadRange(shared_data->team_member, dofs_per_cell),
+            [&](const int &i) {
+              Kokkos::atomic_add(&dst[gpu_data->local_to_global(cell, i)],
+                                 shared_data->values(i, 0));
+            });
+        }
     };
 
     static constexpr unsigned int n_local_dofs = QuadOperation::n_local_dofs;

diff --git a/tests/matrix_free_kokkos/compute_diagonal_01.cc b/tests/matrix_free_kokkos/compute_diagonal_01.cc
index 2865a5069acd231e90a101469fb337eb10fd2eeb..6df68126c443532a86451ff4dba077ae26948aa5 100644 (file)
--- a/tests/matrix_free_kokkos/compute_diagonal_01.cc
+++ b/tests/matrix_free_kokkos/compute_diagonal_01.cc
@@ -30,12 +30,10 @@ test()
   GridGenerator::hyper_cube(tria, -numbers::PI / 2, numbers::PI / 2);
 
   tria.refine_global(2);
-  // FIXME
-  // for (auto &cell : tria.active_cell_iterators())
-  //  if (cell->is_active() && cell->is_locally_owned() &&
-  //      cell->center()[0] < 0.0)
-  // tria.begin_active()->set_refine_flag();
-  // tria.execute_coarsening_and_refinement();
+  for (auto &cell : tria.active_cell_iterators())
+    if (cell->is_active() && cell->is_locally_owned() &&
+        cell->center()[0] < 0.0)
+      tria.execute_coarsening_and_refinement();
 
   const FE_Q<dim>     fe_q(fe_degree);
   const FESystem<dim> fe(fe_q, n_components);
@@ -54,8 +52,6 @@ test()
                                            constraint);
   constraint.close();
 
-  // constraint.print(std::cout);
-
   typename Portable::MatrixFree<dim, Number>::AdditionalData additional_data;
   additional_data.mapping_update_flags = update_values | update_gradients;
 

diff --git a/tests/matrix_free_kokkos/compute_diagonal_util.h b/tests/matrix_free_kokkos/compute_diagonal_util.h
index ae71b3f1c48d42e1a804b972dc81566594974e06..12017ebdc6e066e71fd37e0a98c3a968fe9aea98 100644 (file)
--- a/tests/matrix_free_kokkos/compute_diagonal_util.h
+++ b/tests/matrix_free_kokkos/compute_diagonal_util.h
@@ -120,7 +120,6 @@ public:
               Assert(diagonal_global_host[i] > 0,
                      ExcMessage("Diagonal non-positive at position " +
                                 std::to_string(i)));
-              // std::cout << i << ": " << diagonal_global_host[i] << '\n';
             }
         }
 
@@ -151,12 +150,9 @@ public:
           {
             if (!constraints.is_constrained(i))
               {
-                if (std::abs(A_ref(i, i) - diagonal_global_host(i)) > 1e-6)
-                  // std::cout << i << ": " << A_ref(i,i) << " should be " <<
-                  // diagonal_global_host(i) << '\n';
-                  Assert(std::abs(A_ref(i, i) - diagonal_global_host(i)) < 1e-6,
-                         ExcMessage("Wrong diagonal entry at position " +
-                                    std::to_string(i)));
+                Assert(std::abs(A_ref(i, i) - diagonal_global_host(i)) < 1e-6,
+                       ExcMessage("Wrong diagonal entry at position " +
+                                  std::to_string(i)));
               }
           }
       }