Implement balance-momentum equations for explicit contact (idaholab#2…

…5666)

Kinematic constraint works for primary surface being rigid and works
also in parallel.

Requires verification and adding documentation.

framework/include/constraints/NodeFaceConstraint.h

@@ -117,6 +117,11 @@ class NodeFaceConstraint : public Constraint,
 
   virtual const std::set<unsigned int> & getMatPropDependencies() const;
 
+  /**
+   * Virtual method to avoid initial setups.
+   */
+  virtual bool isExplicitConstraint() const { return false; }
+
 protected:
   /**
    * Compute the value the secondary node should have at the beginning of a timestep.
@@ -221,7 +226,7 @@ class NodeFaceConstraint : public Constraint,
     return coupledNeighborSecond(var_name, comp);
   }
 
-  const std::set<BoundaryID> & getBoundaryIDs();
+  const std::set<BoundaryID> & buildBoundaryIDs();
 
   /// Boundary ID for the secondary surface
   unsigned int _secondary;
@@ -237,6 +242,8 @@ class NodeFaceConstraint : public Constraint,
   std::set<BoundaryID> _boundary_ids;
 
 public:
+  const std::set<SubdomainID> getSecondaryConnectecBlocks() const;
+
   PenetrationLocator & _penetration_locator;
 
 protected:

framework/include/systems/NonlinearSystemBase.h

@@ -978,4 +978,7 @@ class NonlinearSystemBase : public SystemBase, public PerfGraphInterface
 
   /// The number of scaling groups
   std::size_t _num_scaling_groups;
+
+  /// Mutex for auxiliary variables
+  static Threads::spin_mutex writable_variable_mutex;
 };

framework/src/constraints/NodeFaceConstraint.C

@@ -286,9 +286,15 @@ NodeFaceConstraint::overwriteSecondaryResidual()
 }
 
 const std::set<BoundaryID> &
-NodeFaceConstraint::getBoundaryIDs()
+NodeFaceConstraint::buildBoundaryIDs()
 {
   _boundary_ids.insert(_primary);
   _boundary_ids.insert(_secondary);
   return _boundary_ids;
 }
+
+const std::set<SubdomainID>
+NodeFaceConstraint::getSecondaryConnectecBlocks() const
+{
+  return _mesh.getBoundaryConnectedBlocks(_secondary);
+}

framework/src/interfaces/Coupleable.C

@@ -931,44 +931,51 @@ Coupleable::writableCoupledValue(const std::string & var_name, unsigned int comp
 void
 Coupleable::checkWritableVar(MooseWritableVariable * var)
 {
-  // check block restrictions for compatibility
+  // check domain restrictions for compatibility
   const auto * br = dynamic_cast<const BlockRestrictable *>(this);
+  const auto * nfc = dynamic_cast<const NodeFaceConstraint *>(this);
 
-  if (!br)
-    mooseWarning("This object '", _obj->name(), "'is not of the block restrictable type.");
-  else
-  {
-    if (!var->hasBlocks(br->blockIDs()))
-      mooseError("The variable '",
+  //   mooseWarning("This object '", _obj->name(), "'is not of the block restrictable type.");
+
+  if (br && !var->hasBlocks(br->blockIDs()))
+    mooseError("The variable '",
+               var->name(),
+               "' must be defined on all blocks '",
+               _obj->name(),
+               "' is defined on.");
+
+  if (nfc && !var->hasBlocks(nfc->getSecondaryConnectecBlocks()))
+    mooseError("The variable '",
+               var->name(),
+               " must be defined on all blocks '",
+               _obj->name(),
+               "'s secondary surface is defined on.");
+
+  // make sure only one object can access a variable
+  for (const auto & ci : _obj->getMooseApp().getInterfaceObjects<Coupleable>())
+    if (ci != this && ci->_writable_coupled_variables[_c_tid].count(var))
+    {
+      // if both this and ci are block restrictable then we check if the block restrictions
+      // are not overlapping. If they don't we permit the call.
+      const auto * br_other = dynamic_cast<const BlockRestrictable *>(ci);
+      if (br && br_other && br->blockRestricted() && br_other->blockRestricted() &&
+          !MooseUtils::setsIntersect(br->blockIDs(), br_other->blockIDs()))
+        continue;
+      else if (nfc)
+        continue;
+
+      mooseError("'",
+                 ci->_obj->name(),
+                 "' already obtained a writable reference to '",
                  var->name(),
-                 "' must be defined on all blocks '",
-                 _obj->name(),
-                 "' is defined on");
+                 "'. Only one object can obtain such a reference per variable and subdomain in a "
+                 "simulation.");
+    }
 
-    // make sure only one object can access a variable
-    for (const auto & ci : _obj->getMooseApp().getInterfaceObjects<Coupleable>())
-      if (ci != this && ci->_writable_coupled_variables[_c_tid].count(var))
-      {
-        // if both this and ci are block restrictable then we check if the block restrictions
-        // are not overlapping. If they don't we permit the call.
-        const auto * br_other = dynamic_cast<const BlockRestrictable *>(ci);
-        if (br && br_other && br->blockRestricted() && br_other->blockRestricted() &&
-            !MooseUtils::setsIntersect(br->blockIDs(), br_other->blockIDs()))
-          continue;
-
-        mooseError("'",
-                   ci->_obj->name(),
-                   "' already obtained a writable reference to '",
-                   var->name(),
-                   "'. Only one object can obtain such a reference per variable and subdomain in a "
-                   "simulation.");
-      }
-  }
   // var is unique across threads, so we could forego having a separate set per thread, but we
   // need quick access to the list of all variables that need to be inserted into the solution
   // vector by a given thread.
 
-  Moose::out << "Inserting variable... " << var->name() << " in checkWritableVar. \n";
   _writable_coupled_variables[_c_tid].insert(var);
 }
 

framework/src/systems/NonlinearSystemBase.C

@@ -1143,6 +1143,8 @@ NonlinearSystemBase::setConstraintSecondaryValues(NumericVector<Number> & soluti
 
             for (const auto & nfc : constraints)
             {
+              if (nfc->isExplicitConstraint())
+                continue;
               // Return if this constraint does not correspond to the primary-secondary pair
               // prepared by the outer loops.
               // This continue statement is required when, e.g. one secondary surface constrains
@@ -1162,10 +1164,9 @@ NonlinearSystemBase::setConstraintSecondaryValues(NumericVector<Number> & soluti
                 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
                 for (auto * var : nfc->getWritableCoupledVariables())
                 {
-                  var->insert(_fe_problem.getAuxiliarySystem().solution());
+                  if (var->isNodalDefined())
+                    var->insert(_fe_problem.getAuxiliarySystem().solution());
                 }
-                _fe_problem.getAuxiliarySystem().solution().close();
-                _fe_problem.getAuxiliarySystem().system().update();
               }
             }
           }
@@ -1350,16 +1351,23 @@ NonlinearSystemBase::constraintResiduals(NumericVector<Number> & residual, bool
                 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
                 for (auto * var : nfc->getWritableCoupledVariables())
                 {
-                  var->insert(_fe_problem.getAuxiliarySystem().solution());
+                  if (var->isNodalDefined())
+                    var->insert(_fe_problem.getAuxiliarySystem().solution());
                 }
-                _fe_problem.getAuxiliarySystem().solution().close();
-                _fe_problem.getAuxiliarySystem().system().update();
               }
             }
           }
         }
       }
     }
+    const bool has_explicit_contact(true);
+    if (has_explicit_contact)
+    {
+      _fe_problem.getAuxiliarySystem().solution().close();
+      _fe_problem.getAuxiliarySystem().system().update();
+      solutionOld().close();
+    }
+
     if (_assemble_constraints_separately)
     {
       // Make sure that secondary contribution to primary are assembled, and ghosts have been
@@ -2214,6 +2222,8 @@ NonlinearSystemBase::constraintJacobians(bool displaced)
 
             for (const auto & nfc : constraints)
             {
+              if (nfc->isExplicitConstraint())
+                continue;
               // Return if this constraint does not correspond to the primary-secondary pair
               // prepared by the outer loops.
               // This continue statement is required when, e.g. one secondary surface constrains

framework/src/variables/MooseVariableData.C

@@ -1305,25 +1305,14 @@ MooseVariableData<OutputType>::setDofValue(const OutputData & value, unsigned in
   dof_values[index] = value;
   _has_dof_values = true;
 
-  Moose::out << "MooseVariableData<OutputType>::setDofValue: " << value << " and index: " << index
-             << "\n";
   auto & u = _vector_tag_u[_solution_tag];
   for (unsigned int qp = 0; qp < u.size(); qp++)
   {
-    Moose::out << "Making that modification\n";
     u[qp] = (*_phi)[0][qp] * dof_values[0];
-    Moose::out << " dof_values[0] is: " << dof_values[0] << "\n";
-    Moose::out << "First u of " << qp << " is: " << u[qp] << "\n";
 
     for (unsigned int i = 1; i < dof_values.size(); i++)
-    {
       u[qp] += (*_phi)[i][qp] * dof_values[i];
-      Moose::out << " dof_values[i] is: " << dof_values[i] << "\n";
-      Moose::out << "Second u of " << qp << " is: " << u[qp] << "\n";
-    }
   }
-
-  Moose::out << "Final is (of 0): " << u[0] << "\n";
 }
 
 template <typename OutputType>

modules/contact/doc/content/source/actions/ExplicitDynamicsContactAction.md

@@ -0,0 +1,7 @@
+# Contact Action
+
+!syntax description /Contact/ExplicitDynamicsContactAction
+
+## Description
+
+ExplicitDynamicsContactAction is a MOOSE action that constructs objects needed for mechanical contact enforcement for explicit dynamics finite element simulations.

modules/contact/doc/content/source/userobjects/NodalDensity.md

@@ -0,0 +1,13 @@
+# NodalDensity
+
+## Description
+
+The `NodalDensity` object computes a nodally assigned value of the
+density material property. It is used by explicit dynamics contact.
+
+
+!syntax parameters /UserObjects/NodalDensity
+
+!syntax inputs /UserObjects/NodalDensity
+
+!syntax children /UserObjects/NodalDensity

modules/contact/doc/content/source/userobjects/NodalWaveSpeed.md

@@ -0,0 +1,13 @@
+# NodalWaveSpeed
+
+## Description
+
+The `NodalWaveSpeed` object computes a nodally assigned value of the
+wave speed material property.
+
+
+!syntax parameters /UserObjects/NodalWaveSpeed
+
+!syntax inputs /UserObjects/NodalWaveSpeed
+
+!syntax children /UserObjects/NodalWaveSpeed

modules/contact/explicit_dynamics/first_test.i

@@ -47,6 +47,7 @@
     type = MeshCollectionGenerator
     inputs = ' block_one_id block_two_id'
   []
+  allow_renumbering = false
 []
 
 [Variables]
@@ -71,27 +72,9 @@
   []
   [accel_z]
   []
-  [stress_zz]
-  []
-  [strain_zz]
-  []
 []
 
 [AuxKernels]
-  # [stress_zz]
-  #   type = RankTwoAux
-  #   rank_two_tensor = stress
-  #   index_i = 2
-  #   index_j = 2
-  #   variable = stress_zz
-  # []
-  # [strain_zz]
-  #   type = RankTwoAux
-  #   rank_two_tensor = elastic_strain
-  #   index_i = 2
-  #   index_j = 2
-  #   variable = strain_zz
-  # []
   [accel_x]
     type = TestNewmarkTI
     variable = accel_x
@@ -207,24 +190,11 @@
 [ExplicitDynamicsContact]
   [my_contact]
     model = frictionless
-    # formulation = penalty
     primary = base_front
     secondary = ball_back
-    # penalty = 1e+05
-    # normalize_penalty = true
   []
 []
 
-# [Controls]
-#   [mycontrol]
-#     type = TimePeriod
-#     disable_objects = 'BCs/z_bot'
-#     start_time = 1.0e-3
-#     end_time = 1.0e9
-#     execute_on = 'INITIAL TIMESTEP_END'
-#   []
-# []
-
 [Materials]
   [elasticity_tensor_block_one]
     type = ComputeIsotropicElasticityTensor
@@ -269,16 +239,6 @@
 []
 
 [Postprocessors]
-  [accel_58z]
-    type = NodalVariableValue
-    nodeid = 1
-    variable = accel_z
-  []
-  [vel_58z]
-    type = NodalVariableValue
-    nodeid = 1
-    variable = vel_z
-  []
   [disp_58z]
     type = NodalVariableValue
     nodeid = 1
@@ -296,7 +256,6 @@
 []
 
 [Outputs]
-  # interval = 50
   exodus = true
   csv = true
 []