Apply constraints

dealii · Oct 31, 2023 · f1701f8 · f1701f8
1 parent 7dc9607
commit f1701f8
Show file tree

Hide file tree

Showing 4 changed files with 70 additions and 3 deletions.
diff --git a/include/deal.II/multigrid/mg_constrained_dofs.h b/include/deal.II/multigrid/mg_constrained_dofs.h
@@ -299,6 +299,7 @@ MGConstrainedDoFs::initialize(
   refinement_edge_indices.resize(nlevels);
   level_constraints.resize(nlevels);
   user_constraints.resize(nlevels);
+
   for (unsigned int l = min_level; l <= max_level; ++l)
     {
       if (use_provided_level_relevant_dofs)

diff --git a/source/lac/affine_constraints.inst.in b/source/lac/affine_constraints.inst.in
@@ -53,6 +53,13 @@ for (S : REAL_AND_COMPLEX_SCALARS; T : DEAL_II_VEC_TEMPLATES)
         LinearAlgebra::distributed::T<S> & vec) const;
   }
 
+for (T : DEAL_II_VEC_TEMPLATES)
+  {
+    template void
+      AffineConstraints<double>::condense<LinearAlgebra::distributed::T<float>>(
+        LinearAlgebra::distributed::T<float> & vec) const;
+  }
+
 
 // ---------------------------------------------------------------------
 //

diff --git a/source/multigrid/mg_transfer_matrix_free.cc b/source/multigrid/mg_transfer_matrix_free.cc
@@ -319,6 +319,23 @@ MGTransferMatrixFree<dim, Number>::prolongate_and_add(
   if (dst_inplace == false)
     dst += dst_vec;
 
+  // apply constraints on the fine vector
+  if (this->mg_constrained_dofs != nullptr &&
+      this->mg_constrained_dofs->get_user_constraint_matrix(to_level)
+          .get_local_lines()
+          .size() > 0)
+    {
+      this->mg_constrained_dofs->get_user_constraint_matrix(to_level).set_zero(
+        dst);
+    }
+
+  if (this->mg_constrained_dofs != nullptr &&
+      this->mg_constrained_dofs->have_boundary_indices())
+    for (const auto i :
+         this->mg_constrained_dofs->get_boundary_indices(to_level))
+      if (dst.locally_owned_elements().is_element(i))
+        dst[i] = 0.0;
+
   if (src_inplace == true)
     src.zero_out_ghost_values();
 }
@@ -562,6 +579,37 @@ MGTransferMatrixFree<dim, Number>::do_restrict_add(
     Utilities::fixed_power<dim>(n_child_dofs_1d);
   constexpr unsigned int three_to_dim = Utilities::pow(3, dim);
 
+  // If we have user defined MG constraints, we must create
+  // a non-const, ghosted version of the source vector to distribute
+  // constraints.
+  const LinearAlgebra::distributed::Vector<Number> *to_use = &src;
+  LinearAlgebra::distributed::Vector<Number>        new_src;
+  if (this->mg_constrained_dofs != nullptr &&
+      this->mg_constrained_dofs->get_user_constraint_matrix(from_level)
+          .get_local_lines()
+          .size() > 0)
+    {
+      LinearAlgebra::distributed::Vector<Number> copy_src(src);
+
+      // zero out dbc
+      if (this->mg_constrained_dofs->have_boundary_indices())
+        for (const auto i :
+             this->mg_constrained_dofs->get_boundary_indices(from_level))
+          if (copy_src.locally_owned_elements().is_element(i))
+            copy_src[i] = 0.0;
+
+      // zero out user constraints
+      this->mg_constrained_dofs->get_user_constraint_matrix(from_level)
+        .set_zero(copy_src);
+
+      // Re-initialize new ghosted vector with correct constraints
+      new_src.reinit(copy_src);
+      new_src = copy_src;
+      new_src.update_ghost_values();
+
+      to_use = &new_src;
+    }
+
   for (unsigned int cell = 0; cell < n_owned_level_cells[from_level - 1];
        cell += vec_size)
     {
@@ -577,7 +625,7 @@ MGTransferMatrixFree<dim, Number>::do_restrict_add(
         for (unsigned int v = 0; v < n_lanes; ++v)
           {
             for (unsigned int i = 0; i < n_child_cell_dofs; ++i)
-              evaluation_data[i][v] = src.local_element(indices[i]);
+              evaluation_data[i][v] = to_use->local_element(indices[i]);
             indices += n_child_cell_dofs;
           }
       }
@@ -672,6 +720,17 @@ MGTransferMatrixFree<dim, Number>::do_restrict_add(
             }
         }
     }
+
+  // condense coarse vector; note dbc and pbc are handled during cell loop
+
+  if (this->mg_constrained_dofs != nullptr &&
+      this->mg_constrained_dofs->get_user_constraint_matrix(from_level - 1)
+          .get_local_lines()
+          .size() > 0)
+    {
+      this->mg_constrained_dofs->get_user_constraint_matrix(from_level - 1)
+        .condense(dst);
+    }
 }
 
 

diff --git a/tests/multigrid-global-coarsening/constraints_cmp_01.output b/tests/multigrid-global-coarsening/constraints_cmp_01.output
@@ -16,12 +16,12 @@ DEAL::
 Process #0
 Local range: [0, 50), global size: 50
 Vector data:
-1.000e+00 -1.000e+00 1.000e+00 0.000e+00 1.000e+00 0.000e+00 1.000e+00 5.000e-01 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 0.000e+00 1.000e+00 5.000e-01 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 0.000e+00 1.000e+00 5.000e-01 1.000e+00 -1.000e+00 1.000e+00 0.000e+00 1.000e+00 0.000e+00 1.000e+00 5.000e-01 
+1.000e+00 0.000e+00 1.000e+00 0.000e+00 1.000e+00 0.000e+00 1.000e+00 5.000e-01 1.000e+00 0.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 0.000e+00 0.000e+00 1.000e+00 5.000e-01 1.000e+00 1.000e+00 1.000e+00 0.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 0.000e+00 0.000e+00 1.000e+00 5.000e-01 1.000e+00 0.000e+00 1.000e+00 0.000e+00 1.000e+00 0.000e+00 1.000e+00 5.000e-01 
 DEAL::
 Process #0
 Local range: [0, 16), global size: 16
 Vector data:
-3.000e+00 0.000e+00 3.000e+00 0.000e+00 3.250e+00 3.250e+00 3.250e+00 3.250e+00 3.000e+00 0.000e+00 3.250e+00 3.250e+00 3.250e+00 3.250e+00 3.000e+00 0.000e+00 
+1.000e+00 0.000e+00 4.000e+00 0.000e+00 3.250e+00 3.250e+00 3.250e+00 3.250e+00 4.000e+00 0.000e+00 3.250e+00 3.250e+00 3.250e+00 3.250e+00 1.000e+00 0.000e+00 
 DEAL::
 Process #0
 Local range: [0, 50), global size: 50