Merge pull request #26073 from recuero/normal_contact_dynamics

Initial work on contact for explicit finite element dynamics

framework/include/base/Assembly.h

@@ -714,6 +714,9 @@ class Assembly
    */
   void reinitNeighborAtPhysical(const Elem * neighbor, const std::vector<Point> & physical_points);
 
+  /**
+   * Reinitializes the neighbor side using reference coordinates.
+   */
   void reinitNeighbor(const Elem * neighbor, const std::vector<Point> & reference_points);
 
   /**

framework/include/constraints/NodeFaceConstraint.h

@@ -115,6 +115,29 @@ class NodeFaceConstraint : public Constraint,
 
   void residualSetup() override;
 
+  /**
+   * @returns material property dependencies
+   */
+  virtual const std::unordered_set<unsigned int> & getMatPropDependencies() const;
+
+  /**
+   * Whether (contact) constraint is of 'explicit dynamics' type.
+   */
+  virtual bool isExplicitConstraint() const { return false; }
+
+  /**
+   * Allows for overwriting boundary variables (explicit dynamics contact).
+   */
+  virtual void overwriteBoundaryVariables(NumericVector<Number> & /*soln*/,
+                                          const Node & /*secondary_node*/) const
+  {
+  }
+
+  /**
+   * @returns IDs of the subdomains that connect to the secondary boundary
+   */
+  std::set<SubdomainID> getSecondaryConnectedBlocks() const;
+
 protected:
   /**
    * Compute the value the secondary node should have at the beginning of a timestep.
@@ -219,20 +242,27 @@ class NodeFaceConstraint : public Constraint,
     return coupledNeighborSecond(var_name, comp);
   }
 
+  /**
+   * Builds the \p _boundary_ids data member and returns it
+   * @returns a set containing the secondary and primary boundary IDs
+   */
+  const std::set<BoundaryID> & buildBoundaryIDs();
+
   /// Boundary ID for the secondary surface
-  unsigned int _secondary;
+  BoundaryID _secondary;
   /// Boundary ID for the primary surface
-  unsigned int _primary;
+  BoundaryID _primary;
 
   MooseVariable & _var;
 
   const MooseArray<Point> & _primary_q_point;
   const QBase * const & _primary_qrule;
 
-public:
+  /// the union of the secondary and primary boundary ids
+  std::set<BoundaryID> _boundary_ids;
+
   PenetrationLocator & _penetration_locator;
 
-protected:
   /// current node being processed
   const Node * const & _current_node;
   const Elem * const & _current_primary;
@@ -287,7 +317,9 @@ class NodeFaceConstraint : public Constraint,
   /// The value of the secondary residual
   Real _secondary_residual;
 
-public:
+  /// An empty material property dependency set for use with \p getMatPropDependencies
+  const std::unordered_set<unsigned int> _empty_mat_prop_deps;
+
   std::vector<dof_id_type> _connected_dof_indices;
 
   /// The Jacobian corresponding to the derivatives of the neighbor/primary residual with respect to
@@ -309,4 +341,12 @@ class NodeFaceConstraint : public Constraint,
   /// overwriting the secondary residual we traditionally want to use a different scaling factor from the
   /// one associated with interior physics
   DenseMatrix<Number> _Ken;
+
+  friend class NonlinearSystemBase;
 };
+
+inline const std::unordered_set<unsigned int> &
+NodeFaceConstraint::getMatPropDependencies() const
+{
+  return _empty_mat_prop_deps;
+}

framework/include/systems/NonlinearSystemBase.h

@@ -44,6 +44,7 @@ class Split;
 class KernelBase;
 class BoundaryCondition;
 class ResidualObject;
+class PenetrationInfo;
 
 // libMesh forward declarations
 namespace libMesh
@@ -361,6 +362,12 @@ class NonlinearSystemBase : public SystemBase, public PerfGraphInterface
 
   virtual void setSolution(const NumericVector<Number> & soln);
 
+  /**
+   * Called from explicit time stepping to overwrite boundary positions (explicit dynamics). This
+   * will close/assemble the passed-in \p soln after overwrite
+   */
+  void overwriteNodeFace(NumericVector<Number> & soln);
+
   /**
    * Update active objects of Warehouses owned by NonlinearSystemBase
    */
@@ -689,6 +696,8 @@ class NonlinearSystemBase : public SystemBase, public PerfGraphInterface
    */
   void setupDM();
 
+  using SystemBase::reinitNodeFace;
+
 protected:
   /**
    * Compute the residual for a given tag
@@ -763,6 +772,15 @@ class NonlinearSystemBase : public SystemBase, public PerfGraphInterface
 
   NumericVector<Number> & solutionInternal() const override { return *_sys.solution; }
 
+  /**
+   * Reinitialize quantities such as variables, residuals, Jacobians, materials for node-face
+   * constraints
+   */
+  void reinitNodeFace(const Node & secondary_node,
+                      const BoundaryID secondary_boundary,
+                      const PenetrationInfo & info,
+                      const bool displaced);
+
   /// solution vector from nonlinear solver
   const NumericVector<Number> * _current_solution;
   /// ghosted form of the residual

framework/src/base/Assembly.C

@@ -1867,8 +1867,8 @@ Assembly::reinitFVFace(const FaceInfo & fi)
   if (_current_qrule_face != qrules(dim).fv_face.get())
   {
     setFaceQRule(qrules(dim).fv_face.get(), dim);
-    // The order of the element that is used for initing here doesn't matter since this will just be
-    // used for constant monomials (which only need a single integration point)
+    // The order of the element that is used for initing here doesn't matter since this will just
+    // be used for constant monomials (which only need a single integration point)
     if (dim == 3)
       _current_qrule_face->init(QUAD4);
     else
@@ -1881,9 +1881,9 @@ Assembly::reinitFVFace(const FaceInfo & fi)
               "Our finite volume quadrature rule should always yield a single point");
 
   // We've initialized the reference points. Now we need to compute the physical location of the
-  // quadrature points. We do not do any FE initialization so we cannot simply copy over FE results
-  // like we do in reinitFEFace. Instead we handle the computation of the physical locations
-  // manually
+  // quadrature points. We do not do any FE initialization so we cannot simply copy over FE
+  // results like we do in reinitFEFace. Instead we handle the computation of the physical
+  // locations manually
   _current_q_points_face.resize(1);
   const auto & ref_points = _current_qrule_face->get_points();
   const auto & ref_point = ref_points[0];
@@ -2344,12 +2344,12 @@ Assembly::reinitLowerDElem(const Elem * elem,
 
   if (pts && !weights)
   {
-    // We only have dummy weights so the JxWs computed during our FE reinits are meaningless and we
-    // cannot use them
+    // We only have dummy weights so the JxWs computed during our FE reinits are meaningless and
+    // we cannot use them
 
     if (_subproblem.getCoordSystem(elem->subdomain_id()) == Moose::CoordinateSystemType::COORD_XYZ)
-      // We are in a Cartesian coordinate system and we can just use the element volume method which
-      // has fast computation for certain element types
+      // We are in a Cartesian coordinate system and we can just use the element volume method
+      // which has fast computation for certain element types
       _current_lower_d_elem_volume = elem->volume();
     else
       // We manually compute the volume taking the curvilinear coordinate transformations into
@@ -2422,8 +2422,10 @@ Assembly::reinitNeighborAtPhysical(const Elem * neighbor,
         "If reinitializing with more than one point, then I am dubious of your use case. Perhaps "
         "you are performing a DG type method and you are reinitializing using points from the "
         "element face. In such a case your neighbor JxW must have its index order 'match' the "
-        "element JxW index order, e.g. imagining a vertical 1D face with two quadrature points, if "
-        "index 0 for elem JxW corresponds to the 'top' quadrature point, then index 0 for neighbor "
+        "element JxW index order, e.g. imagining a vertical 1D face with two quadrature points, "
+        "if "
+        "index 0 for elem JxW corresponds to the 'top' quadrature point, then index 0 for "
+        "neighbor "
         "JxW must also correspond to the 'top' quadrature point. And libMesh/MOOSE has no way to "
         "guarantee that with multiple quadrature points.");
 
@@ -2844,11 +2846,12 @@ Assembly::prepareLowerD()
     for (MooseIndex(_jacobian_block_lower_used) tag = 0; tag < _jacobian_block_lower_used.size();
          tag++)
     {
-      // To cover all possible cases we should have 9 combinations below for every 2-permutation of
-      // Lower,Secondary,Primary. However, 4 cases will in general be covered by calls to prepare()
-      // and prepareNeighbor(). These calls will cover SecondarySecondary (ElementElement),
-      // SecondaryPrimary (ElementNeighbor), PrimarySecondary (NeighborElement), and PrimaryPrimary
-      // (NeighborNeighbor). With these covered we only need to prepare the 5 remaining below
+      // To cover all possible cases we should have 9 combinations below for every 2-permutation
+      // of Lower,Secondary,Primary. However, 4 cases will in general be covered by calls to
+      // prepare() and prepareNeighbor(). These calls will cover SecondarySecondary
+      // (ElementElement), SecondaryPrimary (ElementNeighbor), PrimarySecondary (NeighborElement),
+      // and PrimaryPrimary (NeighborNeighbor). With these covered we only need to prepare the 5
+      // remaining below
 
       // derivatives w.r.t. lower dimensional residuals
       jacobianBlockMortar(Moose::LowerLower, vi, vj, LocalDataKey{}, tag)
@@ -4787,8 +4790,8 @@ Assembly::adCurvatures() const
   _calculate_curvatures = true;
   const Order helper_order = _mesh.hasSecondOrderElements() ? SECOND : FIRST;
   const FEType helper_type(helper_order, LAGRANGE);
-  // Must prerequest the second derivatives. Sadly because there is only one _need_second_derivative
-  // map for both volumetric and face FE objects we must request both here
+  // Must prerequest the second derivatives. Sadly because there is only one
+  // _need_second_derivative map for both volumetric and face FE objects we must request both here
   feSecondPhi<Real>(helper_type);
   feSecondPhiFace<Real>(helper_type);
   return _ad_curvatures;
@@ -4820,8 +4823,8 @@ Assembly::helpersRequestData()
 
   for (unsigned int dim = 0; dim < _mesh_dimension; dim++)
   {
-    // We need these computations in order to compute correct lower-d element volumes in curvilinear
-    // coordinates
+    // We need these computations in order to compute correct lower-d element volumes in
+    // curvilinear coordinates
     _holder_fe_lower_helper[dim]->get_xyz();
     _holder_fe_lower_helper[dim]->get_JxW();
   }
@@ -4868,8 +4871,8 @@ Assembly::havePRefinement(const std::vector<FEFamily> & disable_p_refinement_for
       helper_container[dim] = unique_helper.get();
 
       // If the user did not request the helper type then we should erase it from our FE container
-      // so that they're not penalized (in the "we should be able to do p-refinement sense") for our
-      // perhaps silly helpers
+      // so that they're not penalized (in the "we should be able to do p-refinement sense") for
+      // our perhaps silly helpers
       if (!user_added_helper_type)
       {
         auto & fe_container_dim = libmesh_map_find(fe_container, dim);

framework/src/constraints/NodeFaceConstraint.C

@@ -284,3 +284,17 @@ NodeFaceConstraint::overwriteSecondaryResidual()
 {
   return _overwrite_secondary_residual;
 }
+
+const std::set<BoundaryID> &
+NodeFaceConstraint::buildBoundaryIDs()
+{
+  _boundary_ids.insert(_primary);
+  _boundary_ids.insert(_secondary);
+  return _boundary_ids;
+}
+
+std::set<SubdomainID>
+NodeFaceConstraint::getSecondaryConnectedBlocks() const
+{
+  return _mesh.getBoundaryConnectedBlocks(_secondary);
+}

framework/src/interfaces/Coupleable.C

@@ -21,6 +21,7 @@
 #include "AuxKernel.h"
 #include "ElementUserObject.h"
 #include "NodalUserObject.h"
+#include "NodeFaceConstraint.h"
 
 Coupleable::Coupleable(const MooseObject * moose_object, bool nodal, bool is_fv)
   : _c_parameters(moose_object->parameters()),
@@ -862,12 +863,13 @@ Coupleable::writableVariable(const std::string & var_name, unsigned int comp)
   const auto * aux = dynamic_cast<const AuxKernel *>(this);
   const auto * euo = dynamic_cast<const ElementUserObject *>(this);
   const auto * nuo = dynamic_cast<const NodalUserObject *>(this);
+  const auto * nfc = dynamic_cast<const NodeFaceConstraint *>(this);
 
-  if (!aux && !euo && !nuo)
-    mooseError("writableVariable() can only be called from AuxKernels, ElementUserObjects, or "
-               "NodalUserObjects. '",
+  if (!aux && !euo && !nuo && !nfc)
+    mooseError("writableVariable() can only be called from AuxKernels, ElementUserObjects, "
+               "NodalUserObjects, or NodeFaceConstraints. '",
                _obj->name(),
-               "' is neither of those.");
+               "' is none of those.");
 
   if (aux && !aux->isNodal() && var->isNodal())
     mooseError("The elemental AuxKernel '",
@@ -929,14 +931,23 @@ 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);
-  if (!var->hasBlocks(br->blockIDs()))
+  const auto * nfc = dynamic_cast<const NodeFaceConstraint *>(this);
+
+  if (br && !var->hasBlocks(br->blockIDs()))
     mooseError("The variable '",
                var->name(),
                "' must be defined on all blocks '",
                _obj->name(),
-               "' is defined on");
+               "' is defined on.");
+
+  if (nfc && !var->hasBlocks(nfc->getSecondaryConnectedBlocks()))
+    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>())
@@ -948,6 +959,8 @@ Coupleable::checkWritableVar(MooseWritableVariable * var)
       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(),
@@ -960,6 +973,7 @@ Coupleable::checkWritableVar(MooseWritableVariable * var)
   // 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.
+
   _writable_coupled_variables[_c_tid].insert(var);
 }