DiffSystem::set_constrain_in_solver()

This should let us play with the option easily in FEMSystem codes, which
should be good for verification purposes.

Author: roystgnr

include/systems/diff_system.h

@@ -378,6 +378,16 @@ class DifferentiableSystem : public ImplicitSystem,
    */
   bool print_element_jacobians;
 
+  /**
+   * set_constrain_in_solver to false to apply constraints only via
+   * residual terms in the systems to be solved.
+   */
+  virtual void set_constrain_in_solver(bool enable);
+
+  virtual bool get_constrain_in_solver() {
+    return _constrain_in_solver;
+  }
+
 protected:
 
   /**
@@ -412,6 +422,12 @@ class DifferentiableSystem : public ImplicitSystem,
   void add_dot_var_dirichlet_bcs( unsigned int var_idx, unsigned int dot_var_idx);
 #endif
 
+  /**
+   * _constrain_in_solver defaults to true; if false then we apply
+   * constraints only via residual terms in the systems to be solved.
+   */
+  bool _constrain_in_solver;
+
 private:
   /**
    * Stack of pointers to objects to use for physics assembly

src/systems/diff_system.C

@@ -48,6 +48,7 @@ DifferentiableSystem::DifferentiableSystem(EquationSystems & es,
   print_element_solutions(false),
   print_element_residuals(false),
   print_element_jacobians(false),
+  _constrain_in_solver(true),
   _diff_physics(),
   _diff_qoi()
 {
@@ -417,6 +418,11 @@ void DifferentiableSystem::pop_physics ()
 }
 
 
+void DifferentiableSystem::set_constrain_in_solver(bool enable)
+{
+  _constrain_in_solver = enable;
+  this->time_solver->diff_solver()->set_exact_constraint_enforcement(enable);
+}
 
 
 } // namespace libMesh

src/systems/fem_system.C

@@ -266,20 +266,22 @@ void add_element_system(FEMSystem & _sys,
       libMesh::out.precision(old_precision);
     }
 
-  // We turn off the asymmetric constraint application;
-  // enforce_constraints_exactly() should be called in the solver
+  // We turn off the asymmetric constraint application iff we expect
+  // enforce_constraints_exactly() to be called in the solver
+  const bool constrain_in_solver = _sys.get_constrain_in_solver();
+
   if (_get_residual && _get_jacobian)
     {
       if (_constrain_heterogeneously)
         _sys.get_dof_map().heterogenously_constrain_element_matrix_and_vector
           (_femcontext.get_elem_jacobian(),
            _femcontext.get_elem_residual(),
-           _femcontext.get_dof_indices(), false);
+           _femcontext.get_dof_indices(), !constrain_in_solver);
       else if (!_no_constraints)
         _sys.get_dof_map().constrain_element_matrix_and_vector
           (_femcontext.get_elem_jacobian(),
            _femcontext.get_elem_residual(),
-           _femcontext.get_dof_indices(), false);
+           _femcontext.get_dof_indices(), !constrain_in_solver);
       // Do nothing if (_no_constraints)
     }
   else if (_get_residual)
@@ -288,10 +290,11 @@ void add_element_system(FEMSystem & _sys,
         _sys.get_dof_map().heterogenously_constrain_element_vector
           (_femcontext.get_elem_jacobian(),
            _femcontext.get_elem_residual(),
-           _femcontext.get_dof_indices(), false);
+           _femcontext.get_dof_indices(), !constrain_in_solver);
       else if (!_no_constraints)
         _sys.get_dof_map().constrain_element_vector
-          (_femcontext.get_elem_residual(), _femcontext.get_dof_indices(), false);
+          (_femcontext.get_elem_residual(),
+           _femcontext.get_dof_indices(), !constrain_in_solver);
       // Do nothing if (_no_constraints)
     }
   else if (_get_jacobian)
@@ -302,8 +305,14 @@ void add_element_system(FEMSystem & _sys,
       if (!_no_constraints)
         _sys.get_dof_map().constrain_element_matrix (_femcontext.get_elem_jacobian(),
                                                      _femcontext.get_dof_indices(),
-                                                     false);
+                                                     !constrain_in_solver);
     }
+
+  if (_get_residual && !_no_constraints && !constrain_in_solver)
+    _sys.get_dof_map().residual_constrain_element_vector
+      (_femcontext.get_elem_residual(),
+       _femcontext.get_dof_indices(),
+       *_sys.current_local_solution);
 #else
   libmesh_ignore(_constrain_heterogeneously, _no_constraints);
 #endif // #ifdef LIBMESH_ENABLE_CONSTRAINTS