Re-enable a mix of nodal and node-face constraints

framework/include/constraints/NodalConstraint.h

@@ -46,7 +46,7 @@ class NodalConstraint : public Constraint,
   {
     mooseError("NodalConstraint do not need computeResidual()");
   }
-  virtual void computeResidual(NumericVector<Number> & residual);
+  virtual void computeResidual(const NumericVector<Number> & residual);
 
   /**
    * Computes the jacobian for the current element.
@@ -55,7 +55,7 @@ class NodalConstraint : public Constraint,
   {
     mooseError("NodalConstraint do not need computeJacobian()");
   }
-  virtual void computeJacobian(SparseMatrix<Number> & jacobian);
+  virtual void computeJacobian(const SparseMatrix<Number> & jacobian);
 
   /**
    * The variable number that this object operates on.

framework/include/systems/NonlinearSystemBase.h

@@ -823,7 +823,14 @@ class NonlinearSystemBase : public SolverSystem, public PerfGraphInterface
    * Enforce nodal constraints
    */
   void enforceNodalConstraintsResidual(NumericVector<Number> & residual);
-  void enforceNodalConstraintsJacobian();
+
+  /**
+   * Enforce nodal constraints in the Jacobian
+   * @param jacobian The Jacobian to read from while constructing the Jacobians corresponding to the
+   * nodal constraints
+   * @returns Whether there were active nodal constraints
+   */
+  bool enforceNodalConstraintsJacobian(const SparseMatrix<Number> & jacobian);
 
   /**
    * Do mortar constraint residual/jacobian computations

framework/src/constraints/NodalConstraint.C

@@ -57,7 +57,7 @@ NodalConstraint::NodalConstraint(const InputParameters & parameters)
 }
 
 void
-NodalConstraint::computeResidual(NumericVector<Number> & residual)
+NodalConstraint::computeResidual(const NumericVector<Number> & residual)
 {
   if ((_weights.size() == 0) && (_primary_node_vector.size() == 1))
     _weights.push_back(1.0);
@@ -99,7 +99,7 @@ NodalConstraint::computeResidual(NumericVector<Number> & residual)
 }
 
 void
-NodalConstraint::computeJacobian(SparseMatrix<Number> & jacobian)
+NodalConstraint::computeJacobian(const SparseMatrix<Number> & jacobian)
 {
   if ((_weights.size() == 0) && (_primary_node_vector.size() == 1))
     _weights.push_back(1.0);

framework/src/systems/NonlinearSystemBase.C

@@ -1126,13 +1126,12 @@ NonlinearSystemBase::enforceNodalConstraintsResidual(NumericVector<Number> & res
   }
 }
 
-void
-NonlinearSystemBase::enforceNodalConstraintsJacobian()
+bool
+NonlinearSystemBase::enforceNodalConstraintsJacobian(const SparseMatrix<Number> & jacobian_to_view)
 {
   if (!hasMatrix(systemMatrixTag()))
     mooseError(" A system matrix is required");
 
-  auto & jacobian = getMatrix(systemMatrixTag());
   THREAD_ID tid = 0; // constraints are going to be done single-threaded
 
   if (_constraints.hasActiveNodalConstraints())
@@ -1147,11 +1146,15 @@ NonlinearSystemBase::enforceNodalConstraintsJacobian()
       {
         _fe_problem.reinitNodes(primary_node_ids, tid);
         _fe_problem.reinitNodesNeighbor(secondary_node_ids, tid);
-        nc->computeJacobian(jacobian);
+        nc->computeJacobian(jacobian_to_view);
       }
     }
     _fe_problem.addCachedJacobian(tid);
+
+    return true;
   }
+  else
+    return false;
 }
 
 void
@@ -3052,32 +3055,39 @@ NonlinearSystemBase::computeJacobianInternal(const std::set<TagID> & tags)
       // Some constraints need to be able to read values from the Jacobian, which requires that it
       // be closed/assembled
       auto & system_matrix = getMatrix(systemMatrixTag());
-#if PETSC_RELEASE_GREATER_EQUALS(3, 23, 0)
-      SparseMatrix<Number> * view_jac_ptr;
       std::unique_ptr<SparseMatrix<Number>> hash_copy;
-      if (system_matrix.use_hash_table())
+      const SparseMatrix<Number> * view_jac_ptr;
+      auto make_readable_jacobian = [&]()
       {
-        hash_copy = libMesh::cast_ref<PetscMatrix<Number> &>(system_matrix).copy_from_hash();
-        view_jac_ptr = hash_copy.get();
-      }
-      else
-        view_jac_ptr = &system_matrix;
-      auto & jacobian_to_view = *view_jac_ptr;
+#if PETSC_RELEASE_GREATER_EQUALS(3, 23, 0)
+        if (system_matrix.use_hash_table())
+        {
+          hash_copy = libMesh::cast_ref<PetscMatrix<Number> &>(system_matrix).copy_from_hash();
+          view_jac_ptr = hash_copy.get();
+        }
+        else
+          view_jac_ptr = &system_matrix;
 #else
-      auto & jacobian_to_view = system_matrix;
+        view_jac_ptr = &system_matrix;
 #endif
-      if (&jacobian_to_view == &system_matrix)
-        system_matrix.close();
+        if (view_jac_ptr == &system_matrix)
+          system_matrix.close();
+      };
+
+      make_readable_jacobian();
 
       // Nodal Constraints
-      enforceNodalConstraintsJacobian();
+      const bool had_nodal_constraints = enforceNodalConstraintsJacobian(*view_jac_ptr);
+      if (had_nodal_constraints)
+        // We have to make a new readable Jacobian
+        make_readable_jacobian();
 
       // Undisplaced Constraints
-      constraintJacobians(jacobian_to_view, false);
+      constraintJacobians(*view_jac_ptr, false);
 
       // Displaced Constraints
       if (_fe_problem.getDisplacedProblem())
-        constraintJacobians(jacobian_to_view, true);
+        constraintJacobians(*view_jac_ptr, true);
     }
   }
   PARALLEL_CATCH;

modules/solid_mechanics/include/constraints/NodalFrictionalConstraint.h

@@ -20,8 +20,8 @@ class NodalFrictionalConstraint : public NodalConstraint
 
   virtual void meshChanged() override;
 
-  virtual void computeResidual(NumericVector<Number> & residual) override;
-  virtual void computeJacobian(SparseMatrix<Number> & jacobian) override;
+  virtual void computeResidual(const NumericVector<Number> & residual) override;
+  virtual void computeJacobian(const SparseMatrix<Number> & jacobian) override;
   using NodalConstraint::computeJacobian;
   using NodalConstraint::computeResidual;
 

modules/solid_mechanics/include/constraints/NodalStickConstraint.h

@@ -20,8 +20,8 @@ class NodalStickConstraint : public NodalConstraint
 
   virtual void meshChanged() override;
 
-  virtual void computeResidual(NumericVector<Number> & residual) override;
-  virtual void computeJacobian(SparseMatrix<Number> & jacobian) override;
+  virtual void computeResidual(const NumericVector<Number> & residual) override;
+  virtual void computeJacobian(const SparseMatrix<Number> & jacobian) override;
   using NodalConstraint::computeJacobian;
   using NodalConstraint::computeResidual;
 

modules/solid_mechanics/src/constraints/NodalFrictionalConstraint.C

@@ -214,7 +214,7 @@ NodalFrictionalConstraint::updateConstrainedNodes()
 }
 
 void
-NodalFrictionalConstraint::computeResidual(NumericVector<Number> &
+NodalFrictionalConstraint::computeResidual(const NumericVector<Number> &
                                            /*residual*/)
 {
   const auto & primarydof = _var.dofIndices();
@@ -261,7 +261,7 @@ NodalFrictionalConstraint::computeQpResidual(Moose::ConstraintType type)
 }
 
 void
-NodalFrictionalConstraint::computeJacobian(SparseMatrix<Number> & /*jacobian*/)
+NodalFrictionalConstraint::computeJacobian(const SparseMatrix<Number> & /*jacobian*/)
 {
   // Calculate Jacobian entries and cache those entries along with the row and column indices
   std::vector<dof_id_type> secondarydof = _var.dofIndicesNeighbor();

modules/solid_mechanics/src/constraints/NodalStickConstraint.C

@@ -192,7 +192,7 @@ NodalStickConstraint::updateConstrainedNodes()
 }
 
 void
-NodalStickConstraint::computeJacobian(SparseMatrix<Number> & jacobian)
+NodalStickConstraint::computeJacobian(const SparseMatrix<Number> & jacobian)
 {
   // Calculate Jacobian enteries and cache those entries along with the row and column indices
   std::vector<dof_id_type> secondarydof = _var.dofIndicesNeighbor();
@@ -236,7 +236,7 @@ NodalStickConstraint::computeJacobian(SparseMatrix<Number> & jacobian)
 }
 
 void
-NodalStickConstraint::computeResidual(NumericVector<Number> & residual)
+NodalStickConstraint::computeResidual(const NumericVector<Number> & residual)
 {
   std::vector<dof_id_type> primarydof = _var.dofIndices();
   std::vector<dof_id_type> secondarydof = _var.dofIndicesNeighbor();