DiffSolver::set_exact_constraint_enforcement

Author: roystgnr

include/solvers/diff_solver.h

@@ -285,8 +285,36 @@ class DiffSolver : public ReferenceCountedObject<DiffSolver>,
    */
   std::unique_ptr<LinearSolutionMonitor> linear_solution_monitor;
 
+  /**
+   * Enable (or disable; it is \p true by default) exact enforcement
+   * of constraints at the solver level, correcting any constrained
+   * DoF coefficients in \p current_local_solution as well as applying
+   * nonlinear residual and Jacobian terms based on constraint
+   * equations.
+   *
+   * This is probably only safe to disable if user code is setting
+   * nonlinear residual and Jacobian terms based on constraint
+   * equations at an element-by-element level, by combining the
+   * \p asymmetric_constraint_rows option with the
+   * \p residual_constrain_element_vector processing option in
+   * \p DofMap.
+   */
+  virtual void set_exact_constraint_enforcement(bool enable) {
+    _exact_constraint_enforcement = enable;
+  }
+
+  bool exact_constraint_enforcement() {
+    return _exact_constraint_enforcement;
+  }
+
 protected:
 
+  /**
+   * Whether we should enforce exact constraints globally during a
+   * solve.
+   */
+  bool _exact_constraint_enforcement;
+
   /**
    * The largest solution norm which the DiffSolver has yet seen will be stored
    * here, to be used for stopping criteria based on relative_step_tolerance

src/solvers/diff_solver.C

@@ -41,6 +41,7 @@ DiffSolver::DiffSolver (sys_type & s) :
   relative_step_tolerance(0.),
   initial_linear_tolerance(1e-12),
   minimum_linear_tolerance(TOLERANCE*TOLERANCE),
+  _exact_constraint_enforcement(true),
   max_solution_norm(0.),
   max_residual_norm(0.),
   _outer_iterations(0),

src/solvers/newton_solver.C

@@ -295,7 +295,8 @@ unsigned int NewtonSolver::solve()
   rhs.close();
 
 #ifdef LIBMESH_ENABLE_CONSTRAINTS
-  _system.get_dof_map().enforce_constraints_exactly(_system);
+  if (this->_exact_constraint_enforcement)
+    _system.get_dof_map().enforce_constraints_exactly(_system);
 #endif
 
   SparseMatrix<Number> & matrix = *(_system.matrix);
@@ -426,8 +427,9 @@ unsigned int NewtonSolver::solve()
       _system.update ();
       // The linear solver may not have fit our constraints exactly
 #ifdef LIBMESH_ENABLE_CONSTRAINTS
-      _system.get_dof_map().enforce_constraints_exactly
-        (_system, &linear_solution, /* homogeneous = */ true);
+      if (this->_exact_constraint_enforcement)
+        _system.get_dof_map().enforce_constraints_exactly
+          (_system, &linear_solution, /* homogeneous = */ true);
 #endif
 
       const unsigned int linear_steps = rval.first;
@@ -553,7 +555,8 @@ unsigned int NewtonSolver::solve()
 
   // The linear solver may not have fit our constraints exactly
 #ifdef LIBMESH_ENABLE_CONSTRAINTS
-  _system.get_dof_map().enforce_constraints_exactly(_system);
+  if (this->_exact_constraint_enforcement)
+    _system.get_dof_map().enforce_constraints_exactly(_system);
 #endif
 
   // We may need to localize a parallel solution

src/solvers/petsc_diff_solver.C

@@ -115,7 +115,8 @@ extern "C"
     sys.update();
 
     // We may need to correct a non-conforming solution
-    sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
+    if (solver.exact_constraint_enforcement())
+      sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
 
     // Do DiffSystem assembly
     sys.assembly(true, false);
@@ -168,7 +169,8 @@ extern "C"
     sys.update();
 
     // We may need to correct a non-conforming solution
-    sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
+    if (solver.exact_constraint_enforcement())
+      sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
 
     // Do DiffSystem assembly
     sys.assembly(false, true);
@@ -318,7 +320,8 @@ unsigned int PetscDiffSolver::solve()
   LIBMESH_CHKERR(ierr);
 
 #ifdef LIBMESH_ENABLE_CONSTRAINTS
-  _system.get_dof_map().enforce_constraints_exactly(_system);
+  if (this->_exact_constraint_enforcement)
+    _system.get_dof_map().enforce_constraints_exactly(_system);
 #endif
 
   SNESConvergedReason reason;