Enable enforcement of non-zero gaps

src/tests/CMakeLists.txt

@@ -33,6 +33,7 @@ set( tribol_tests
      tribol_mortar_wts.cpp
      tribol_nodal_nrmls.cpp
      tribol_quad_integ.cpp
+     tribol_residual_gap.cpp
      tribol_tet_mesh.cpp
      tribol_timestep_vote.cpp
      tribol_twb_integ.cpp

src/tests/tribol_enzyme_element_mortar.cpp

@@ -974,6 +974,43 @@ TEST_F( EnzymeElementMortarTest, NoOverlap )
   FDCheck( x1, x2, n1, p1, x1_stencil, x2_stencil, num_nodes, check_scale, len_collapse_ratio );
 }
 
+TEST_F( EnzymeElementMortarTest, ExactOverlapResidualGap )
+{
+  // clang-format off
+  // Setup two overlapping elements with a gap of 0.1
+  double x1[12] = { 0.0, 1.0, 1.0, 0.0,
+                    0.0, 0.0, 1.0, 1.0,
+                    0.0, 0.0, 0.0, 0.0 };
+  // Define x2 in CW order so ElemReverse makes it CCW
+  double x2[12] = { 0.0, 0.0, 1.0, 1.0,
+                    0.0, 1.0, 1.0, 0.0,
+                    0.1, 0.1, 0.1, 0.1 };
+  double n1[12] = { 0.0, 0.0, 0.0, 0.0,
+                    0.0, 0.0, 0.0, 0.0,
+                    1.0, 1.0, 1.0, 1.0 };
+  double p1[4] = { 1.0, 1.0, 1.0, 1.0 };
+  // clang-format on
+
+  double residual_gap = 0.15;
+
+  double f1[12], f2[12], g1[4];
+  int num_nodes = 4;
+
+  // Call with residual_gap
+  tribol::ComputeMortarForceEnzyme( x1, n1, p1, f1, g1, num_nodes, x2, f2, num_nodes, 1.0e-8, residual_gap );
+
+  // Kinematic gap is 0.1. Effective gap is 0.1 - 0.15 = -0.05.
+  // Integrated gap for unit square (overlap area 1.0) distributed to 4 nodes:
+  // Each node i has integrated gap G_i = \int \phi_i (g_kin - g_r) dA.
+  // For constant g_kin and g_r: G_i = (g_kin - g_r) \int \phi_i dA.
+  // \int \phi_i dA for bilinear quad node is 0.25 (1/4 of element area).
+  // So G_i = -0.05 * 0.25 = -0.0125.
+
+  for ( int i = 0; i < 4; ++i ) {
+    EXPECT_NEAR( g1[i], -0.0125, 1e-12 );
+  }
+}
+
 }  // namespace tribol
 
 //------------------------------------------------------------------------------

src/tests/tribol_finite_diff_energy_mortar.cpp

@@ -400,6 +400,74 @@ TEST( QuadraturePointPenaltyCheck, OpenGapIsInactive )
   EXPECT_EQ( result.energy, 0.0 );
 }
 
+TEST( EnergyMortarResidualGapCheck, AssembledGapIsShiftedByArea )
+{
+  RealT x1[2] = { 0.0, 1.0 };
+  RealT y1[2] = { 0.0, 0.0 };
+  IndexT conn1[2] = { 1, 0 };
+  MeshData mesh1( 0, 1, 2, conn1, LINEAR_EDGE, x1, y1, nullptr, MemorySpace::Host );
+
+  RealT x2[2] = { 0.2, 0.8 };
+  RealT y2[2] = { 0.1, 0.1 };
+  IndexT conn2[2] = { 0, 1 };
+  MeshData mesh2( 1, 1, 2, conn2, LINEAR_EDGE, x2, y2, nullptr, MemorySpace::Host );
+
+  ContactParams params;
+  params.del = 0.1;
+  params.k = 3.0;
+  params.N = 3;
+  params.enzyme_quadrature = true;
+
+  double gap_without_residual[2] = { 0.0, 0.0 };
+  double area_without_residual[2] = { 0.0, 0.0 };
+  EnergyMortarCalculator evaluator_without_residual( params );
+  evaluator_without_residual.compute_gtilde_and_area( InterfacePair( 0, 0 ), mesh1.getView(), mesh2.getView(),
+                                                      gap_without_residual, area_without_residual );
+
+  params.residual_gap = 0.15;
+  double gap_with_residual[2] = { 0.0, 0.0 };
+  double area_with_residual[2] = { 0.0, 0.0 };
+  EnergyMortarCalculator evaluator_with_residual( params );
+  evaluator_with_residual.compute_gtilde_and_area( InterfacePair( 0, 0 ), mesh1.getView(), mesh2.getView(),
+                                                   gap_with_residual, area_with_residual );
+
+  for ( int i = 0; i < 2; ++i ) {
+    EXPECT_NEAR( area_with_residual[i], area_without_residual[i], 1.0e-14 );
+    EXPECT_NEAR( gap_with_residual[i], gap_without_residual[i] - params.residual_gap * area_without_residual[i],
+                 1.0e-14 );
+  }
+}
+
+TEST( EnergyMortarResidualGapCheck, QuadraturePointOpenGapBecomesActive )
+{
+  RealT x1[2] = { 0.0, 1.0 };
+  RealT y1[2] = { 0.0, 0.0 };
+  IndexT conn1[2] = { 1, 0 };
+  MeshData mesh1( 0, 1, 2, conn1, LINEAR_EDGE, x1, y1, nullptr, MemorySpace::Host );
+
+  RealT x2[2] = { 0.2, 0.8 };
+  RealT y2[2] = { 0.1, 0.1 };
+  IndexT conn2[2] = { 0, 1 };
+  MeshData mesh2( 1, 1, 2, conn2, LINEAR_EDGE, x2, y2, nullptr, MemorySpace::Host );
+
+  ContactParams params;
+  params.del = 0.1;
+  params.k = 3.0;
+  params.N = 3;
+  params.enzyme_quadrature = true;
+
+  EnergyMortarCalculator evaluator_without_residual( params );
+  const auto inactive = evaluator_without_residual.compute_quadrature_point_penalty_data(
+      InterfacePair( 0, 0 ), mesh1.getView(), mesh2.getView() );
+  EXPECT_EQ( inactive.energy, 0.0 );
+
+  params.residual_gap = 0.15;
+  EnergyMortarCalculator evaluator_with_residual( params );
+  const auto active = evaluator_with_residual.compute_quadrature_point_penalty_data( InterfacePair( 0, 0 ),
+                                                                                     mesh1.getView(), mesh2.getView() );
+  EXPECT_GT( active.energy, 0.0 );
+}
+
 TEST( QuadraturePointPenaltyCheck, DerivativesMatchFiniteDifference )
 {
   RealT x1[2] = { 0.0, 1.0 };
@@ -417,6 +485,7 @@ TEST( QuadraturePointPenaltyCheck, DerivativesMatchFiniteDifference )
   params.k = 3.0;
   params.N = 3;
   params.enzyme_quadrature = true;
+  params.residual_gap = 0.15;
 
   EnergyMortarCalculator evaluator( params );
   const InterfacePair pair( 0, 0 );
@@ -516,6 +585,7 @@ TEST( GradientCheck, GtildeFDvsAD )
   params_.k = 1.0;                   // penalty stiffness
   params_.N = 3;                     // quadrature points
   params_.enzyme_quadrature = true;  // use the non-Enzyme quadrature path
+  params_.residual_gap = 0.15;
 
   EnergyMortarCalculator evaluator_( params_ );
 
@@ -575,6 +645,7 @@ TEST( HessianCheck, GtildeFDvsAD )
   params_.k = 1.0;
   params_.N = 3;
   params_.enzyme_quadrature = true;
+  params_.residual_gap = 0.15;
 
   EnergyMortarCalculator evaluator_( params_ );
 

src/tests/tribol_mfem_mortar_lm.cpp

@@ -35,7 +35,7 @@
  * @brief This tests the Tribol MFEM interface running a contact patch test.
  *
  */
-class MfemMortarTest : public testing::TestWithParam<std::tuple<int, mfem::Element::Type>> {
+class MfemMortarTest : public testing::TestWithParam<std::tuple<int, mfem::Element::Type, double>> {
  protected:
   tribol::RealT max_disp_;
   void SetUp() override
@@ -46,6 +46,8 @@ class MfemMortarTest : public testing::TestWithParam<std::tuple<int, mfem::Eleme
     // polynomial order of the finite element discretization
     int order = 1;
 
+    double residual_gap = std::get<2>( GetParam() );
+
     // fixed options
     // boundary element attributes of mortar surface
     auto mortar_attrs = std::set<int>( { 4 } );
@@ -144,6 +146,7 @@ class MfemMortarTest : public testing::TestWithParam<std::tuple<int, mfem::Eleme
                                         tribol::SINGLE_MORTAR, tribol::FRICTIONLESS, tribol::LAGRANGE_MULTIPLIER,
                                         tribol::BINNING_GRID );
     tribol::setLagrangeMultiplierOptions( 0, tribol::ImplicitEvalMode::MORTAR_RESIDUAL_JACOBIAN );
+    tribol::setResidualGap( coupling_scheme_id, residual_gap );
 
     // update tribol (compute contact contribution to force and stiffness)
     tribol::updateMfemParallelDecomposition();
@@ -219,14 +222,17 @@ class MfemMortarTest : public testing::TestWithParam<std::tuple<int, mfem::Eleme
 
 TEST_P( MfemMortarTest, check_mortar_displacement )
 {
-  EXPECT_LT( std::abs( max_disp_ - 0.005 ), 1.0e-6 );
+  double residual_gap = std::get<2>( GetParam() );
+  EXPECT_LT( std::abs( max_disp_ - ( 0.005 + residual_gap / 2.0 ) ), 1.0e-6 );
 
   MPI_Barrier( MPI_COMM_WORLD );
 }
 
 INSTANTIATE_TEST_SUITE_P( tribol, MfemMortarTest,
-                          testing::Values( std::make_tuple( 2, mfem::Element::Type::HEXAHEDRON ),
-                                           std::make_tuple( 2, mfem::Element::Type::TETRAHEDRON ) ) );
+                          testing::Values( std::make_tuple( 2, mfem::Element::Type::HEXAHEDRON, 0.0 ),
+                                           std::make_tuple( 2, mfem::Element::Type::HEXAHEDRON, 0.01 ),
+                                           std::make_tuple( 2, mfem::Element::Type::TETRAHEDRON, 0.0 ),
+                                           std::make_tuple( 2, mfem::Element::Type::TETRAHEDRON, 0.01 ) ) );
 
 //------------------------------------------------------------------------------
 #include "axom/slic/core/SimpleLogger.hpp"

src/tests/tribol_residual_gap.cpp

@@ -0,0 +1,179 @@
+// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and
+// other Tribol Project Developers. See the top-level LICENSE file for details.
+//
+// SPDX-License-Identifier: (MIT)
+
+// gtest includes
+#include "gtest/gtest.h"
+
+// Tribol includes
+#include "tribol/interface/tribol.hpp"
+#include "tribol/common/Parameters.hpp"
+#include "tribol/mesh/CouplingScheme.hpp"
+
+using RealT = tribol::RealT;
+
+class ResidualGapTest : public ::testing::Test {
+ protected:
+  void SetUp() override {}
+  void TearDown() override { tribol::finalize(); }
+};
+
+TEST_F( ResidualGapTest, common_plane_residual_gap )
+{
+  // Setup two 3D quads with a 0.1 gap
+  // Quad 1: z = 0, x in [0, 1], y in [0, 1]
+  // Quad 2: z = 0.1, x in [0, 1], y in [0, 1]
+
+  constexpr int numVerts = 4;
+  constexpr int numElems = 1;
+  RealT x1[numVerts] = { 0.0, 1.0, 1.0, 0.0 };
+  RealT y1[numVerts] = { 0.0, 0.0, 1.0, 1.0 };
+  RealT z1[numVerts] = { 0.0, 0.0, 0.0, 0.0 };
+
+  RealT x2[numVerts] = { 0.0, 1.0, 1.0, 0.0 };
+  RealT y2[numVerts] = { 0.0, 0.0, 1.0, 1.0 };
+  RealT z2[numVerts] = { 0.1, 0.1, 0.1, 0.1 };
+
+  // Connectivity for 3D quads
+  tribol::IndexT conn1[numVerts] = { 0, 1, 2, 3 };
+  // Quad 2 normal points -z: 0-3-2-1 is CCW viewed from -z.
+  tribol::IndexT conn2[numVerts] = { 0, 3, 2, 1 };
+
+  tribol::registerMesh( 0, numElems, numVerts, &conn1[0], static_cast<int>( tribol::LINEAR_QUAD ), &x1[0], &y1[0],
+                        &z1[0], tribol::MemorySpace::Host );
+  tribol::registerMesh( 1, numElems, numVerts, &conn2[0], static_cast<int>( tribol::LINEAR_QUAD ), &x2[0], &y2[0],
+                        &z2[0], tribol::MemorySpace::Host );
+
+  RealT fx1[numVerts] = { 0., 0., 0., 0. };
+  RealT fy1[numVerts] = { 0., 0., 0., 0. };
+  RealT fz1[numVerts] = { 0., 0., 0., 0. };
+  RealT fx2[numVerts] = { 0., 0., 0., 0. };
+  RealT fy2[numVerts] = { 0., 0., 0., 0. };
+  RealT fz2[numVerts] = { 0., 0., 0., 0. };
+
+  tribol::registerNodalResponse( 0, &fx1[0], &fy1[0], &fz1[0] );
+  tribol::registerNodalResponse( 1, &fx2[0], &fy2[0], &fz2[0] );
+
+  RealT penalty = 1.0;
+  tribol::setKinematicConstantPenalty( 0, penalty );
+  tribol::setKinematicConstantPenalty( 1, penalty );
+
+  tribol::registerCouplingScheme( 0, 0, 1, tribol::SURFACE_TO_SURFACE, tribol::NO_CASE, tribol::COMMON_PLANE,
+                                  tribol::FRICTIONLESS, tribol::PENALTY, tribol::BINNING_GRID,
+                                  tribol::ExecutionMode::Sequential );
+
+  tribol::setPenaltyOptions( 0, tribol::KINEMATIC, tribol::KINEMATIC_CONSTANT );
+
+  // Default residual_gap is 0.0. Gap is 0.1, so no contact.
+  RealT dt = 1.0;
+  tribol::update( 1, 1.0, dt );
+
+  // there still should be an active pair even with no contact...
+  auto& cs = tribol::CouplingSchemeManager::getInstance().at( 0 );
+  EXPECT_EQ( 1, cs.getNumActivePairs() );
+
+  // ...but the pair should have zero force
+  for ( int i = 0; i < numVerts; ++i ) {
+    EXPECT_NEAR( fz1[i], 0.0, 1.e-10 );
+    EXPECT_NEAR( fz2[i], 0.0, 1.e-10 );
+  }
+
+  // Set residual_gap to 0.15. Now gap (0.1) < residual_gap (0.15), so contact.
+  RealT residual_gap = 0.15;
+  tribol::setResidualGap( 0, residual_gap );
+
+  tribol::update( 2, 2.0, dt );
+
+  EXPECT_EQ( 1, cs.getNumActivePairs() );
+
+  // Check forces.
+  // actual gap g = 0.1.
+  // pressure p = (g - residual_gap) * effective_penalty = (0.1 - 0.15) * 0.5 = -0.025
+  // Force on mesh 1 nodes: repulsion points DOWN (-z).
+  // Total force = p * A = -0.025 * 1.0 = -0.025.
+  // Distributed to 4 nodes: -0.025 / 4 = -0.00625.
+
+  for ( int i = 0; i < numVerts; ++i ) {
+    EXPECT_NEAR( fz1[i], -0.00625, 1.e-10 );
+    EXPECT_NEAR( fz2[i], 0.00625, 1.e-10 );
+  }
+}
+
+TEST_F( ResidualGapTest, mortar_residual_gap )
+{
+  // Setup two 3D quads with a 0.1 gap
+  constexpr int numVerts = 4;
+  constexpr int numElems = 1;
+  RealT x1[numVerts] = { 0.0, 1.0, 1.0, 0.0 };
+  RealT y1[numVerts] = { 0.0, 0.0, 1.0, 1.0 };
+  RealT z1[numVerts] = { 0.0, 0.0, 0.0, 0.0 };
+
+  RealT x2[numVerts] = { 0.0, 1.0, 1.0, 0.0 };
+  RealT y2[numVerts] = { 0.0, 0.0, 1.0, 1.0 };
+  RealT z2[numVerts] = { 0.1, 0.1, 0.1, 0.1 };
+
+  tribol::IndexT conn1[numVerts] = { 0, 1, 2, 3 };
+  tribol::IndexT conn2[numVerts] = { 0, 3, 2, 1 };
+
+  tribol::registerMesh( 0, numElems, numVerts, &conn1[0], (int)( tribol::LINEAR_QUAD ), &x1[0], &y1[0], &z1[0],
+                        tribol::MemorySpace::Host );
+  tribol::registerMesh( 1, numElems, numVerts, &conn2[0], (int)( tribol::LINEAR_QUAD ), &x2[0], &y2[0], &z2[0],
+                        tribol::MemorySpace::Host );
+
+  RealT fx1[numVerts] = { 0., 0., 0., 0. };
+  RealT fy1[numVerts] = { 0., 0., 0., 0. };
+  RealT fz1[numVerts] = { 0., 0., 0., 0. };
+  RealT fx2[numVerts] = { 0., 0., 0., 0. };
+  RealT fy2[numVerts] = { 0., 0., 0., 0. };
+  RealT fz2[numVerts] = { 0., 0., 0., 0. };
+
+  tribol::registerNodalResponse( 0, &fx1[0], &fy1[0], &fz1[0] );
+  tribol::registerNodalResponse( 1, &fx2[0], &fy2[0], &fz2[0] );
+
+  RealT gaps[numVerts] = { 0., 0., 0., 0. };
+  RealT pressures[numVerts] = { 0., 0., 0., 0. };
+  tribol::registerMortarGaps( 1, gaps );
+  tribol::registerMortarPressures( 1, pressures );
+
+  tribol::registerCouplingScheme( 0, 0, 1, tribol::SURFACE_TO_SURFACE, tribol::NO_CASE, tribol::SINGLE_MORTAR,
+                                  tribol::FRICTIONLESS, tribol::LAGRANGE_MULTIPLIER, tribol::BINNING_GRID,
+                                  tribol::ExecutionMode::Sequential );
+
+  tribol::setLagrangeMultiplierOptions( 0, tribol::ImplicitEvalMode::MORTAR_RESIDUAL );
+
+  // Default residual_gap is 0.0. Gap is 0.1.
+  RealT dt = 1.0;
+  tribol::update( 1, 1.0, dt );
+
+  auto& cs = tribol::CouplingSchemeManager::getInstance().at( 0 );
+  EXPECT_EQ( 1, cs.getNumActivePairs() );
+
+  // Nodal gaps in mortar are integrated: 0.1 * 0.25 = 0.025
+  for ( int i = 0; i < numVerts; ++i ) {
+    EXPECT_NEAR( gaps[i], 0.025, 1.e-10 );
+  }
+
+  // Set residual_gap to 0.15.
+  RealT residual_gap = 0.15;
+  tribol::setResidualGap( 0, residual_gap );
+
+  // Clear gaps to avoid accumulation in test
+  for ( int i = 0; i < numVerts; ++i ) gaps[i] = 0.0;
+
+  tribol::update( 2, 2.0, dt );
+
+  // Kinematic gap is (0.1 - 0.15) = -0.05. Integrated: -0.05 * 0.25 = -0.0125
+  for ( int i = 0; i < numVerts; ++i ) {
+    EXPECT_NEAR( gaps[i], -0.0125, 1.e-10 );
+  }
+}
+
+int main( int argc, char* argv[] )
+{
+  int result = 0;
+  ::testing::InitGoogleTest( &argc, argv );
+  axom::slic::SimpleLogger logger;
+  result = RUN_ALL_TESTS();
+  return result;
+}

src/tribol/common/Parameters.hpp

@@ -505,6 +505,8 @@ struct Parameters {
   // Note, auto-contact will require registration of element thicknesses.
   bool auto_contact_check = false;  ///! True if auto-contact checks should be enabled
 
+  RealT residual_gap = 0.0;  ///! User defined residual gap constraint
+
   EnforcementLocation enforcement_location =
       EnforcementLocation::QuadraturePoint;  ///! Defaults to quadrature-point enforcement
 };

src/tribol/geom/CompGeom.hpp

@@ -1034,7 +1034,7 @@ TRIBOL_HOST_DEVICE inline bool CommonPlanePair::exceedsMaxAutoInterpen( const Me
     RealT max_interpen = -1. * params.auto_contact_pen_frac *
                          axom::utilities::min( mesh1.getElementData().m_thickness[faceId1],
                                                mesh2.getElementData().m_thickness[faceId2] );
-    if ( gap < max_interpen ) {
+    if ( ( gap - params.residual_gap ) < max_interpen ) {
       return true;
     }
   }