Only update RedecompMesh in MFEM interface after sufficient displacement

src/tests/CMakeLists.txt

@@ -122,6 +122,7 @@ if ( BUILD_REDECOMP AND TRIBOL_USE_MPI )
       tribol_mfem_common_plane.cpp
       tribol_mfem_mortar_lm.cpp
       tribol_proximity_check.cpp
+      tribol_redecomp_tol.cpp
       )
 
   set(combined_test_depends tribol gtest)

src/tests/tribol_redecomp_tol.cpp

@@ -0,0 +1,139 @@
+// 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)
+
+#include <set>
+
+#include <gtest/gtest.h>
+
+// Tribol includes
+#include "tribol/config.hpp"
+#include "tribol/common/Parameters.hpp"
+#include "tribol/interface/tribol.hpp"
+#include "tribol/interface/mfem_tribol.hpp"
+#include "tribol/mesh/CouplingScheme.hpp"
+#include "tribol/mesh/MfemData.hpp"
+
+// Shared includes
+#include "shared/mesh/MeshBuilder.hpp"
+
+#ifdef TRIBOL_USE_UMPIRE
+// Umpire includes
+#include "umpire/ResourceManager.hpp"
+#endif
+
+// MFEM includes
+#include "mfem.hpp"
+
+/**
+ * @brief This tests the Tribol MFEM interface's ability to skip expensive RedecompMesh
+ * recreation when nodal displacement change is below a certain threshold.
+ */
+class MfemRedecompSkipTest : public testing::Test {
+ protected:
+  void SetUp() override {}
+};
+
+TEST_F( MfemRedecompSkipTest, test_skip_redecomp )
+{
+  tribol::ExecutionMode exec_mode = tribol::ExecutionMode::Sequential;
+
+  // 1. Initialize simplified problem: Two unit cubes separated by a small gap
+  // Mesh is refined once, so element size is 0.5.
+  // clang-format off
+  double initial_sep = 0.1;
+  mfem::ParMesh mesh = shared::ParMeshBuilder( MPI_COMM_WORLD, shared::MeshBuilder::Unify( {
+    shared::MeshBuilder::CubeMesh( 1, 1, 1 ),
+    shared::MeshBuilder::CubeMesh( 1, 1, 1 )
+      .translate( { 0.0, 0.0, 1.0 + initial_sep } )
+      .updateAttrib( 1, 2 )
+      .updateBdrAttrib( 1, 7 ) // Bottom of top cube -> 7
+      .updateBdrAttrib( 6, 8 ) // set the top surface to boundary attribute 8 (prevent multiple 6)
+  } ).refine( 1 ) );
+  // clang-format on
+
+  auto fe_coll = mfem::H1_FECollection( 1, mesh.SpaceDimension() );
+  auto par_fe_space = mfem::ParFiniteElementSpace( &mesh, &fe_coll, mesh.SpaceDimension() );
+  auto coords = mfem::ParGridFunction( &par_fe_space );
+  mesh.GetNodes( coords );
+
+  int coupling_scheme_id = 0;
+  // Top of bottom cube: 6, Bottom of top cube: 7
+  tribol::registerMfemCouplingScheme( coupling_scheme_id, 0, 1, mesh, coords, { 6 }, { 7 }, tribol::SURFACE_TO_SURFACE,
+                                      tribol::NO_CASE, tribol::COMMON_PLANE, tribol::FRICTIONLESS, tribol::PENALTY,
+                                      tribol::BINNING_BVH, exec_mode );
+
+  // Set threshold explicitly to 0.2
+  tribol::setMfemRedecompTriggerDisplacement( coupling_scheme_id, 0.2 );
+  tribol::setMfemKinematicConstantPenalty( coupling_scheme_id, 1.0e5, 1.0e5 );
+
+  // 2. Initial redecomp
+  tribol::updateMfemParallelDecomposition();
+  auto* mfem_data = tribol::CouplingSchemeManager::getInstance().at( coupling_scheme_id ).getMfemMeshData();
+  const mfem::Mesh* redecomp_mesh_ptr_1 = &mfem_data->GetRedecompMesh();
+
+  // Verify we can compute forces (should be zero as gap > 0)
+  double dt = 0.1;
+  tribol::update( 0, 0.0, dt );
+  mfem::Vector force_1( coords.Size() );
+  force_1 = 0.0;
+  tribol::getMfemResponse( coupling_scheme_id, force_1 );
+
+  // 3. Small shift (0.1 < 0.2)
+  // Shift ALL nodes by 0.1. This preserves relative gap, so forces should remain same (zero).
+  coords += 0.1;
+
+  tribol::updateMfemParallelDecomposition();
+  const mfem::Mesh* redecomp_mesh_ptr_2 = &mfem_data->GetRedecompMesh();
+
+  // Expect no rebuild
+  EXPECT_EQ( redecomp_mesh_ptr_1, redecomp_mesh_ptr_2 );
+
+  tribol::update( 1, 0.1, dt );
+  mfem::Vector force_2( coords.Size() );
+  force_2 = 0.0;
+  tribol::getMfemResponse( coupling_scheme_id, force_2 );
+  // Forces should still be zero (invariant)
+  EXPECT_LT( force_2.Norml2(), 1.0e-10 );
+
+  // 4. Large shift (accumulated 0.3 > 0.2)
+  // Shift ALL nodes by another 0.2 (total 0.3)
+  coords += 0.2;
+
+  tribol::updateMfemParallelDecomposition();
+  const mfem::Mesh* redecomp_mesh_ptr_3 = &mfem_data->GetRedecompMesh();
+
+  // Expect rebuild
+  EXPECT_NE( redecomp_mesh_ptr_1, redecomp_mesh_ptr_3 );
+
+  tribol::update( 2, 0.2, dt );
+  mfem::Vector force_3( coords.Size() );
+  force_3 = 0.0;
+  tribol::getMfemResponse( coupling_scheme_id, force_3 );
+  // Forces should still be zero (invariant)
+  EXPECT_LT( force_3.Norml2(), 1.0e-10 );
+
+  tribol::finalize();
+}
+
+int main( int argc, char* argv[] )
+{
+  int result = 0;
+
+  MPI_Init( &argc, &argv );
+
+  ::testing::InitGoogleTest( &argc, argv );
+
+#ifdef TRIBOL_USE_UMPIRE
+  umpire::ResourceManager::getInstance();  // initialize umpire's ResouceManager
+#endif
+
+  axom::slic::SimpleLogger logger;
+
+  result = RUN_ALL_TESTS();
+
+  MPI_Finalize();
+
+  return result;
+}

src/tribol/common/BasicTypes.hpp

@@ -12,6 +12,9 @@
 // C includes
 #include <cstddef>
 
+// C++ includes
+#include <type_traits>
+
 // MPI includes
 #ifdef TRIBOL_USE_MPI
 #include <mpi.h>
@@ -20,6 +23,9 @@
 // Axom includes
 #include "axom/core/Types.hpp"
 
+// MFEM includes
+#include "mfem.hpp"
+
 namespace tribol {
 
 #ifdef TRIBOL_USE_MPI
@@ -53,6 +59,9 @@ using RealT = double;
 
 #endif
 
+// mfem's real_t should match ours
+static_assert( std::is_same_v<RealT, mfem::real_t>, "tribol::RealT and mfem::real_t are required to match" );
+
 #define TRIBOL_UNUSED_VAR AXOM_UNUSED_VAR
 #define TRIBOL_UNUSED_PARAM AXOM_UNUSED_PARAM
 

src/tribol/interface/mfem_tribol.cpp

@@ -119,6 +119,20 @@ void setMfemLORFactor( IndexT cs_id, int lor_factor )
   cs->getMfemMeshData()->SetLORFactor( lor_factor );
 }
 
+void setMfemRedecompTriggerDisplacement( IndexT cs_id, RealT val )
+{
+  auto cs = CouplingSchemeManager::getInstance().findData( cs_id );
+  SLIC_ERROR_ROOT_IF(
+      !cs, axom::fmt::format( "Coupling scheme cs_id={0} does not exist. Call tribol::registerMfemCouplingScheme() "
+                              "to create a coupling scheme with this cs_id.",
+                              cs_id ) );
+  SLIC_ERROR_ROOT_IF(
+      !cs->hasMfemData(),
+      "Coupling scheme does not contain MFEM data. "
+      "Create the coupling scheme using registerMfemCouplingScheme() to set the trigger displacement." );
+  cs->getMfemMeshData()->SetRedecompTriggerDisplacement( val );
+}
+
 void setMfemKinematicConstantPenalty( IndexT cs_id, RealT mesh1_penalty, RealT mesh2_penalty )
 {
   auto cs = CouplingSchemeManager::getInstance().findData( cs_id );
@@ -369,7 +383,7 @@ mfem::ParGridFunction& getMfemPressure( IndexT cs_id )
   return cs->getMfemSubmeshData()->GetSubmeshPressure();
 }
 
-void updateMfemParallelDecomposition( int n_ranks )
+void updateMfemParallelDecomposition( int n_ranks, bool force_new_redecomp )
 {
   for ( auto& cs_pair : CouplingSchemeManager::getInstance() ) {
     auto& cs = cs_pair.second;
@@ -386,9 +400,9 @@ void updateMfemParallelDecomposition( int n_ranks )
       if ( mfem_data->GetLORFactor() > 1 ) {
         effective_binning_proximity *= static_cast<RealT>( mfem_data->GetLORFactor() );
       }
-      // creates a new redecomp mesh based on updated coordinates and updates transfer operators and displacement,
-      // velocity, and response grid functions based on new redecomp mesh
-      mfem_data->UpdateMfemMeshData( effective_binning_proximity, n_ranks );
+      // creates a new redecomp mesh based on updated coordinates (if criteria is met) and updates transfer operators
+      // and displacement, velocity, and response grid functions based on new redecomp mesh
+      auto new_redecomp = mfem_data->UpdateMfemMeshData( effective_binning_proximity, n_ranks, force_new_redecomp );
       auto coord_ptrs = mfem_data->GetRedecompCoordsPtrs();
 
       registerMesh( mesh_ids[0], mfem_data->GetMesh1NE(), mfem_data->GetNV(), mfem_data->GetMesh1Conn(),
@@ -417,12 +431,12 @@ void updateMfemParallelDecomposition( int n_ranks )
         auto submesh_data = cs.getMfemSubmeshData();
         // updates submesh-native grid functions and transfer operators on
         // the new redecomp mesh
-        submesh_data->UpdateMfemSubmeshData( mfem_data->GetRedecompMesh() );
+        submesh_data->UpdateMfemSubmeshData( mfem_data->GetRedecompMesh(), new_redecomp );
         auto g_ptrs = submesh_data->GetRedecompGapPtrs();
         registerMortarGaps( mesh_ids[1], g_ptrs[0] );
         auto p_ptrs = submesh_data->GetRedecompPressurePtrs();
         registerMortarPressures( mesh_ids[1], p_ptrs[0] );
-        if ( cs.hasMfemJacobianData() ) {
+        if ( cs.hasMfemJacobianData() && new_redecomp ) {
           // updates Jacobian transfer operator for new redecomp mesh
           cs.getMfemJacobianData()->UpdateJacobianXfer();
         }

src/tribol/interface/mfem_tribol.hpp

@@ -85,6 +85,16 @@ void registerMfemCouplingScheme( IndexT cs_id, int mesh_id_1, int mesh_id_2, con
  */
 void setMfemLORFactor( IndexT cs_id, int lor_factor );
 
+/**
+ * @brief Sets the displacement threshold for triggering a new parallel decomposition
+ *
+ * @pre Coupling scheme cs_id must be registered using registerMfemCouplingScheme()
+ *
+ * @param [in] cs_id The ID of the coupling scheme
+ * @param [in] val Threshold value
+ */
+void setMfemRedecompTriggerDisplacement( IndexT cs_id, RealT val );
+
 /**
  * @brief Clears existing penalty data and sets kinematic constant penalty
  *
@@ -289,10 +299,12 @@ mfem::ParGridFunction& getMfemPressure( IndexT cs_id );
  * @brief Updates mesh parallel decomposition and related grid functions/Jacobian when coordinates are updated
  *
  * @param n_ranks Number of ranks in the parallel decomposition; automatically determine when set to 0 (default)
+ * @param force_new_redecomp If true, construct a new RedecompMesh even if displacement threshold is not met (default =
+ * false)
  *
  * @pre Coupling schemes must be registered using registerMfemCouplingScheme()
  */
-void updateMfemParallelDecomposition( int n_ranks = 0 );
+void updateMfemParallelDecomposition( int n_ranks = 0, bool force_new_redecomp = false );
 
 /**
  * @brief Create VisIt output of the parallel repartitioned RedecompMesh