API for consistent constraints via residuals

This is a relatively clean way to handle drift away from constraints at
an elemental rather than a global level.

The distinction between how we want to do heterogeneous constraints in
linear vs in nonlinear systems gets a little tricky.

Author: roystgnr

include/base/dof_map.h

@@ -1249,6 +1249,10 @@ class DofMap : public ReferenceCountedObject<DofMap>,
    * appropriate for finding linearized updates to a solution in which
    * heterogeneous constraints are already satisfied.
    *
+   * Note the sign difference from the nonlinear heterogeneous constraint
+   * method: Solving u:=K\f has the opposite sign convention from
+   * u:=u_in-J\r, and we apply heterogeneous constraints accordingly.
+   *
    * By default, the constraints for the primal solution of this
    * system are used.  If a non-negative \p qoi_index is passed in,
    * then the constraints for the corresponding adjoint solution are
@@ -1272,7 +1276,67 @@ class DofMap : public ReferenceCountedObject<DofMap>,
       (matrix, rhs, elem_dofs, asymmetric_constraint_rows, qoi_index);
   }
 
+  /**
+   * Constrains the element Jacobian and residual.  The element
+   * Jacobian is square, and the elem_dofs should correspond to the
+   * global DOF indices of both the rows and columns of the element
+   * matrix.
+   *
+   * The residual-constraining version of this method creates linear
+   * systems in which heterogeneously constrained degrees of freedom
+   * create non-zero residual terms when not at their correct offset
+   * values, as would be appropriate for finding a solution to a
+   * nonlinear problem in a quasi-Newton solve.
+   *
+   * Note the sign difference from the linear heterogeneous constraint
+   * method: Solving u:=u_in-J\r has the opposite sign convention from
+   * u:=K\f, and we apply heterogeneous constraints accordingly.
+   *
+   * The \p solution vector passed in should be a serialized or
+   * ghosted primal solution
+   */
+  void heterogeneously_constrain_element_jacobian_and_residual (DenseMatrix<Number> & matrix,
+                                                                DenseVector<Number> & rhs,
+                                                                std::vector<dof_id_type> & elem_dofs,
+                                                                NumericVector<Number> & solution_local) const;
+
+  /**
+   * Constrains the element residual.  The element Jacobian is square,
+   * and the elem_dofs should correspond to the global DOF indices of
+   * both the rows and columns of the element matrix.
+   *
+   * The residual-constraining version of this method creates linear
+   * systems in which heterogeneously constrained degrees of freedom
+   * create non-zero residual terms when not at their correct offset
+   * values, as would be appropriate for finding a solution to a
+   * nonlinear problem in a quasi-Newton solve.
+   *
+   * The \p solution vector passed in should be a serialized or
+   * ghosted primal solution
+   */
+  void heterogeneously_constrain_element_residual (DenseVector<Number> & rhs,
+                                                   std::vector<dof_id_type> & elem_dofs,
+                                                   NumericVector<Number> & solution_local) const;
+
 
+  /**
+   * Constrains the element residual.  The element Jacobian is square,
+   * and the elem_dofs should correspond to the global DOF indices of
+   * both the rows and columns of the element matrix, and the dof
+   * constraint should not include any heterogeneous terms.
+   *
+   * The residual-constraining version of this method creates linear
+   * systems in which heterogeneously constrained degrees of freedom
+   * create non-zero residual terms when not at their correct offset
+   * values, as would be appropriate for finding a solution to a
+   * nonlinear problem in a quasi-Newton solve.
+   *
+   * The \p solution vector passed in should be a serialized or
+   * ghosted primal solution
+   */
+  void constrain_element_residual (DenseVector<Number> & rhs,
+                                   std::vector<dof_id_type> & elem_dofs,
+                                   NumericVector<Number> & solution_local) const;
 
   /**
    * Constrains a dyadic element matrix B = v w'.  This method

src/base/dof_map_constraints.C

@@ -2622,6 +2622,207 @@ void DofMap::heterogeneously_constrain_element_matrix_and_vector (DenseMatrix<Nu
 }
 
 
+void DofMap::heterogeneously_constrain_element_jacobian_and_residual
+  (DenseMatrix<Number> & matrix,
+   DenseVector<Number> & rhs,
+   std::vector<dof_id_type> & elem_dofs,
+   NumericVector<Number> & solution_local) const
+{
+  libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
+  libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
+  libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
+
+  libmesh_assert (solution_local.type() == SERIAL ||
+                  solution_local.type() == GHOSTED);
+
+  // check for easy return
+  if (this->_dof_constraints.empty())
+    return;
+
+  // The constrained matrix is built up as C^T K C.
+  // The constrained RHS is built up as C^T F
+  // Asymmetric residual terms are added if we do not have x = Cx+h
+  DenseMatrix<Number> C;
+  DenseVector<Number> H;
+
+  this->build_constraint_matrix_and_vector (C, H, elem_dofs);
+
+  LOG_SCOPE("hetero_cnstrn_elem_jac_res()", "DofMap");
+
+  // It is possible that the matrix is not constrained at all.
+  if ((C.m() != matrix.m()) ||
+      (C.n() != elem_dofs.size()))
+    return;
+
+  // Compute the matrix-vector product C^T F
+  DenseVector<Number> old_rhs(rhs);
+  C.vector_mult_transpose(rhs, old_rhs);
+
+  // Compute the matrix-matrix-matrix product C^T K C
+  matrix.left_multiply_transpose  (C);
+  matrix.right_multiply (C);
+
+  libmesh_assert_equal_to (matrix.m(), matrix.n());
+  libmesh_assert_equal_to (matrix.m(), elem_dofs.size());
+  libmesh_assert_equal_to (matrix.n(), elem_dofs.size());
+
+  for (unsigned int i=0,
+       n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
+       i != n_elem_dofs; i++)
+    {
+      const dof_id_type dof_id = elem_dofs[i];
+
+      const DofConstraints::const_iterator
+        pos = _dof_constraints.find(dof_id);
+
+      if (pos != _dof_constraints.end())
+        {
+          for (auto j : make_range(matrix.n()))
+            matrix(i,j) = 0.;
+
+          // If the DOF is constrained
+          matrix(i,i) = 1.;
+
+          // This will put a nonsymmetric entry in the constraint
+          // row to ensure that the linear system produces the
+          // correct value for the constrained DOF.
+          const DofConstraintRow & constraint_row = pos->second;
+
+          for (const auto & item : constraint_row)
+            for (unsigned int j=0; j != n_elem_dofs; j++)
+              if (elem_dofs[j] == item.first)
+                matrix(i,j) = -item.second;
+
+          const DofConstraintValueMap::const_iterator valpos =
+            _primal_constraint_values.find(dof_id);
+
+          Number & rhs_val = rhs(i);
+          rhs_val = (valpos == _primal_constraint_values.end()) ?
+            0 : -valpos->second;
+          for (const auto & [constraining_dof, coef] : constraint_row)
+            rhs_val -= coef * solution_local(constraining_dof);
+          rhs_val += solution_local(dof_id);
+        }
+    }
+}
+
+
+void DofMap::heterogeneously_constrain_element_residual
+  (DenseVector<Number> & rhs,
+   std::vector<dof_id_type> & elem_dofs,
+   NumericVector<Number> & solution_local) const
+{
+  libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
+
+  libmesh_assert (solution_local.type() == SERIAL ||
+                  solution_local.type() == GHOSTED);
+
+  // check for easy return
+  if (this->_dof_constraints.empty())
+    return;
+
+  // The constrained RHS is built up as C^T F
+  // Asymmetric residual terms are added if we do not have x = Cx+h
+  DenseMatrix<Number> C;
+  DenseVector<Number> H;
+
+  this->build_constraint_matrix_and_vector (C, H, elem_dofs);
+
+  LOG_SCOPE("hetero_cnstrn_elem_res()", "DofMap");
+
+  // It is possible that the element is not constrained at all.
+  if ((C.m() != rhs.size()) ||
+      (C.n() != elem_dofs.size()))
+    return;
+
+  // Compute the matrix-vector product C^T F
+  DenseVector<Number> old_rhs(rhs);
+  C.vector_mult_transpose(rhs, old_rhs);
+
+  for (unsigned int i=0,
+       n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
+       i != n_elem_dofs; i++)
+    {
+      const dof_id_type dof_id = elem_dofs[i];
+
+      const DofConstraints::const_iterator
+        pos = _dof_constraints.find(dof_id);
+
+      if (pos != _dof_constraints.end())
+        {
+          // This will put a nonsymmetric entry in the constraint
+          // row to ensure that the linear system produces the
+          // correct value for the constrained DOF.
+          const DofConstraintRow & constraint_row = pos->second;
+
+          const DofConstraintValueMap::const_iterator valpos =
+            _primal_constraint_values.find(dof_id);
+
+          Number & rhs_val = rhs(i);
+          rhs_val = (valpos == _primal_constraint_values.end()) ?
+            0 : -valpos->second;
+          for (const auto & [constraining_dof, coef] : constraint_row)
+            rhs_val -= coef * solution_local(constraining_dof);
+          rhs_val += solution_local(dof_id);
+        }
+    }
+}
+
+
+void DofMap::constrain_element_residual
+  (DenseVector<Number> & rhs,
+   std::vector<dof_id_type> & elem_dofs,
+   NumericVector<Number> & solution_local) const
+{
+  libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
+
+  libmesh_assert (solution_local.type() == SERIAL ||
+                  solution_local.type() == GHOSTED);
+
+  // check for easy return
+  if (this->_dof_constraints.empty())
+    return;
+
+  // The constrained RHS is built up as C^T F
+  DenseMatrix<Number> C;
+
+  this->build_constraint_matrix (C, elem_dofs);
+
+  LOG_SCOPE("cnstrn_elem_residual()", "DofMap");
+
+  // It is possible that the matrix is not constrained at all.
+  if (C.n() != elem_dofs.size())
+    return;
+
+  // Compute the matrix-vector product C^T F
+  DenseVector<Number> old_rhs(rhs);
+  C.vector_mult_transpose(rhs, old_rhs);
+
+  for (unsigned int i=0,
+       n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
+       i != n_elem_dofs; i++)
+    {
+      const dof_id_type dof_id = elem_dofs[i];
+
+      const DofConstraints::const_iterator
+        pos = _dof_constraints.find(dof_id);
+
+      if (pos != _dof_constraints.end())
+        {
+          // This will put a nonsymmetric entry in the constraint
+          // row to ensure that the linear system produces the
+          // correct value for the constrained DOF.
+          const DofConstraintRow & constraint_row = pos->second;
+
+          Number & rhs_val = rhs(i);
+          rhs_val = 0;
+          for (const auto & [constraining_dof, coef] : constraint_row)
+            rhs_val -= coef * solution_local(constraining_dof);
+          rhs_val += solution_local(dof_id);
+        }
+    }
+}
+
 
 void DofMap::heterogeneously_constrain_element_vector (const DenseMatrix<Number> & matrix,
                                                        DenseVector<Number> & rhs,

src/solvers/newton_solver.C

@@ -83,7 +83,7 @@ Real NewtonSolver::line_search(Real tol,
       _system.update();
 
       // Check residual with fractional Newton step
-      _system.assembly (true, false);
+      _system.assembly(true, false, !this->_exact_constraint_enforcement);
 
       rhs.close();
       current_residual = rhs.l2_norm();
@@ -191,7 +191,7 @@ Real NewtonSolver::line_search(Real tol,
 
       // We may need to localize a parallel solution
       _system.update();
-      _system.assembly (true, false);
+      _system.assembly(true, false, !this->_exact_constraint_enforcement);
 
       rhs.close();
       Real fu = current_residual = rhs.l2_norm();
@@ -317,7 +317,7 @@ unsigned int NewtonSolver::solve()
       if (verbose)
         libMesh::out << "Assembling the System" << std::endl;
 
-      _system.assembly(true, true);
+      _system.assembly(true, true, !this->_exact_constraint_enforcement);
       rhs.close();
       Real current_residual = rhs.l2_norm();
 
@@ -472,7 +472,7 @@ unsigned int NewtonSolver::solve()
           !continue_after_max_iterations)
         {
           _system.update ();
-          _system.assembly(true, false);
+          _system.assembly(true, false, !this->_exact_constraint_enforcement);
 
           rhs.close();
           current_residual = rhs.l2_norm();

src/solvers/petsc_diff_solver.C

@@ -119,7 +119,7 @@ extern "C"
       sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
 
     // Do DiffSystem assembly
-    sys.assembly(true, false);
+    sys.assembly(true, false, !solver.exact_constraint_enforcement());
     R_system.close();
 
     // Swap back
@@ -173,7 +173,7 @@ extern "C"
       sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
 
     // Do DiffSystem assembly
-    sys.assembly(false, true);
+    sys.assembly(false, true, !solver.exact_constraint_enforcement());
     J_system.close();
 
     // Swap back