Energy mortar gap enforcement at quadrature points

src/redecomp/transfer/MatrixTransfer.cpp

@@ -13,6 +13,8 @@
 #include "shared/math/ParSparseMat.hpp"
 #include "shared/infrastructure/Profiling.hpp"
 
+#include <vector>
+
 namespace redecomp {
 
 void MatrixTransfer::validateTransferInputs( const axom::Array<int>& test_elem_idx,
@@ -302,9 +304,9 @@ shared::ParSparseMat MatrixTransfer::TransferToParallel( const axom::Array<int>&
   }
 
   auto num_rows = parent_test_fes_.GetVSize();
-  int* I_ptr = new int[num_rows + 1];
-  HYPRE_BigInt* J_ptr = new HYPRE_BigInt[num_unique_nonzeros];
-  double* data_ptr = new double[num_unique_nonzeros];
+  std::vector<int> I_ptr( num_rows + 1, 0 );
+  std::vector<HYPRE_BigInt> J_ptr( num_unique_nonzeros );
+  std::vector<double> data_ptr( num_unique_nonzeros );
 
   // Initialize I_ptr with zeros
   for ( int i = 0; i <= num_rows; ++i ) {
@@ -337,7 +339,7 @@ shared::ParSparseMat MatrixTransfer::TransferToParallel( const axom::Array<int>&
 
   // Construct rectangular HypreParMatrix
   shared::ParSparseMat J_full( getMPIUtility().MPIComm(), num_rows, parent_test_fes_.GlobalVSize(),
-                               parent_trial_fes_.GlobalVSize(), I_ptr, J_ptr, data_ptr,
+                               parent_trial_fes_.GlobalVSize(), I_ptr.data(), J_ptr.data(), data_ptr.data(),
                                parent_test_fes_.GetDofOffsets(), parent_trial_fes_.GetDofOffsets() );
 
   if ( !parallel_assemble ) {

src/shared/math/ParSparseMat.cpp

@@ -97,10 +97,11 @@ void ParSparseMatView::eliminateRows( const mfem::Array<int>& rows )
   invokeHypreMethod<MemorySpace::Host>( [&]() { mat_->EliminateRows( rows ); } );
 }
 
-void ParSparseMatView::eliminateRowsCols( const mfem::Array<int>& rows_cols )
+ParSparseMat ParSparseMatView::eliminateRowsCols( const mfem::Array<int>& rows_cols )
 {
   ensureHostMemory( *this );
-  invokeHypreMethod<MemorySpace::Host>( [&]() { mat_->EliminateRowsCols( rows_cols ); } );
+  return ParSparseMat(
+      createHypreParMatrix<MemorySpace::Host>( [&]() { return mat_->EliminateRowsCols( rows_cols ); } ) );
 }
 
 ParSparseMat ParSparseMatView::eliminateCols( const mfem::Array<int>& cols )

src/shared/math/ParSparseMat.hpp

@@ -137,8 +137,9 @@ class ParSparseMatView {
    * @brief Eliminates chosen rows and columns from the matrix
    *
    * @param rows_cols Array of rows/columns to eliminate
+   * @return Matrix containing the eliminated entries
    */
-  void eliminateRowsCols( const mfem::Array<int>& rows_cols );
+  ParSparseMat eliminateRowsCols( const mfem::Array<int>& rows_cols );
 
   /**
    * @brief Eliminates chosen columns from the matrix

src/tests/shared_par_sparse_mat.cpp

@@ -349,12 +349,12 @@ TEST_F( ParSparseMatTest, Elimination )
 
   // Eliminate row 0 (globally)
   // Determine if I own row 0
-  mfem::Array<int> rows_to_elim;
+  mfem::Array<int> indices_to_elim;
   int row_starts_idx = HYPRE_AssumedPartitionCheck() ? 0 : rank;
   if ( row_starts[row_starts_idx] == 0 ) {
-    rows_to_elim.Append( 0 );
+    indices_to_elim.Append( 0 );
   }
-  A.eliminateRows( rows_to_elim );
+  A.eliminateRows( indices_to_elim );
 
   // Check if row 0 is identity (or zero with diagonal 1)
   // Diagonal matrix means we can just check multiplication
@@ -377,6 +377,15 @@ TEST_F( ParSparseMatTest, Elimination )
     }
   }
 
+  A = shared::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 3.0 );
+  shared::ParSparseMat Ae = A.eliminateRowsCols( indices_to_elim );
+  x.fill( 1.0 );
+  y = A * x;
+  shared::ParVector ye = Ae * x;
+  for ( int i = 0; i < y.size(); ++i ) {
+    EXPECT_NEAR( y[i] + ye[i], 3.0, 1e-12 );
+  }
+
   A = shared::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 3.0 );
   int num_procs;
   MPI_Comm_size( MPI_COMM_WORLD, &num_procs );
@@ -385,7 +394,7 @@ TEST_F( ParSparseMatTest, Elimination )
   auto last_local_col = A.width() - 1;
   mfem::Array<int> cols_to_elim( { last_local_col } );
 
-  shared::ParSparseMat Ae = A.eliminateCols( cols_to_elim );
+  Ae = A.eliminateCols( cols_to_elim );
 
   // Now check A * e_last = 0
   // Create vector with 1 at last_local_col, 0 elsewhere
@@ -396,7 +405,7 @@ TEST_F( ParSparseMatTest, Elimination )
   EXPECT_NEAR( y_last[last_local_col], 0.0, 1e-12 );
 
   // Check Ae * e_last = original value
-  shared::ParVector ye = Ae * x_last;
+  ye = Ae * x_last;
   double expected_val = 3.0;
 
   EXPECT_NEAR( ye[last_local_col], expected_val, 1e-12 );

src/tests/tribol_energy_mortar_patch.cpp

@@ -39,7 +39,7 @@
  *   u_y(x,y) = eps_yy * y
  *   u_x(x,y) = eps_xx * x
  */
-class MfemMortarEnergyPatchTest : public testing::TestWithParam<std::tuple<int>> {
+class MfemMortarEnergyPatchTest : public testing::TestWithParam<std::tuple<int, tribol::EnforcementLocation>> {
  protected:
   tribol::RealT max_disp_;
   double l2_err_vec_;
@@ -56,6 +56,7 @@ class MfemMortarEnergyPatchTest : public testing::TestWithParam<std::tuple<int>>
   void SetUp() override
   {
     int ref_levels = std::get<0>( GetParam() );
+    tribol::EnforcementLocation enforcement_location = std::get<1>( GetParam() );
     int order = 1;
 
     auto mortar_attrs = std::set<int>( { 5 } );
@@ -175,7 +176,8 @@ class MfemMortarEnergyPatchTest : public testing::TestWithParam<std::tuple<int>>
     tribol::registerMfemCouplingScheme( cs_id, mesh1_id, mesh2_id, mesh, coords, mortar_attrs, nonmortar_attrs,
                                         tribol::SURFACE_TO_SURFACE, tribol::NO_SLIDING, tribol::ENERGY_MORTAR,
                                         tribol::FRICTIONLESS, tribol::PENALTY, tribol::BINNING_GRID );
-    tribol::setMfemKinematicConstantPenalty( cs_id, 10000.0, 10000.0 );
+    tribol::setEnforcementLocation( cs_id, enforcement_location );
+    tribol::setMfemKinematicConstantPenalty( cs_id, 100.0, 100.0 );
 
     mfem::Vector X( par_fe_space.GetTrueVSize() );
     X = 0.0;
@@ -356,7 +358,10 @@ TEST_P( MfemMortarEnergyPatchTest, check_patch_test )
   MPI_Barrier( MPI_COMM_WORLD );
 }
 
-INSTANTIATE_TEST_SUITE_P( tribol, MfemMortarEnergyPatchTest, testing::Values( std::make_tuple( 2 ) ) );
+INSTANTIATE_TEST_SUITE_P( tribol, MfemMortarEnergyPatchTest,
+                          testing::Combine( testing::Values( 2 ),
+                                            testing::Values( tribol::EnforcementLocation::Nodal,
+                                                             tribol::EnforcementLocation::QuadraturePoint ) ) );
 
 //------------------------------------------------------------------------------
 #include "axom/slic/core/SimpleLogger.hpp"

src/tests/tribol_finite_diff_energy_mortar.cpp

@@ -372,6 +372,122 @@ FiniteDiffResult EnergyMortarCalculator::validate_hessian( const InterfacePair&
   return result;
 }
 
+TEST( QuadraturePointPenaltyCheck, OpenGapIsInactive )
+{
+  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( params );
+  const auto result =
+      evaluator.compute_quadrature_point_penalty_data( InterfacePair( 0, 0 ), mesh1.getView(), mesh2.getView() );
+  EXPECT_EQ( result.energy, 0.0 );
+}
+
+TEST( QuadraturePointPenaltyCheck, DerivativesMatchFiniteDifference )
+{
+  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( params );
+  const InterfacePair pair( 0, 0 );
+  const auto analytical = evaluator.compute_quadrature_point_penalty_data( pair, mesh1.getView(), mesh2.getView() );
+  ASSERT_GT( analytical.energy, 0.0 );
+
+  const std::array<RealT, 2> x1_orig{ x1[0], x1[1] };
+  const std::array<RealT, 2> y1_orig{ y1[0], y1[1] };
+  const std::array<RealT, 2> x2_orig{ x2[0], x2[1] };
+  const std::array<RealT, 2> y2_orig{ y2[0], y2[1] };
+
+  auto restore = [&]() {
+    x1[0] = x1_orig[0];
+    x1[1] = x1_orig[1];
+    y1[0] = y1_orig[0];
+    y1[1] = y1_orig[1];
+    x2[0] = x2_orig[0];
+    x2[1] = x2_orig[1];
+    y2[0] = y2_orig[0];
+    y2[1] = y2_orig[1];
+    mesh1.setPosition( x1, y1, nullptr );
+    mesh2.setPosition( x2, y2, nullptr );
+  };
+
+  auto perturb = [&]( int dof, double delta ) {
+    if ( dof < 4 ) {
+      const int endpoint = dof / 2;
+      const int component = dof % 2;
+      const int node = conn1[endpoint];
+      ( component == 0 ? x1[node] : y1[node] ) += delta;
+      mesh1.setPosition( x1, y1, nullptr );
+    } else {
+      const int endpoint = ( dof - 4 ) / 2;
+      const int component = ( dof - 4 ) % 2;
+      const int node = conn2[endpoint];
+      ( component == 0 ? x2[node] : y2[node] ) += delta;
+      mesh2.setPosition( x2, y2, nullptr );
+    }
+  };
+
+  const double gradient_eps = 1.0e-7;
+  const double gradient_tol = 1.0e-6;
+  for ( int dof = 0; dof < 8; ++dof ) {
+    restore();
+    perturb( dof, gradient_eps );
+    const double energy_plus =
+        evaluator.compute_quadrature_point_penalty_energy( pair, mesh1.getView(), mesh2.getView() );
+    restore();
+    perturb( dof, -gradient_eps );
+    const double energy_minus =
+        evaluator.compute_quadrature_point_penalty_energy( pair, mesh1.getView(), mesh2.getView() );
+    const double fd_force = ( energy_plus - energy_minus ) / ( 2.0 * gradient_eps );
+    EXPECT_NEAR( fd_force, analytical.force[dof], gradient_tol ) << "force mismatch at dof " << dof;
+  }
+
+  const double hessian_eps = 1.0e-6;
+  const double hessian_tol = 1.0e-4;
+  for ( int col = 0; col < 8; ++col ) {
+    restore();
+    perturb( col, hessian_eps );
+    const auto force_plus =
+        evaluator.compute_quadrature_point_penalty_data( pair, mesh1.getView(), mesh2.getView() ).force;
+    restore();
+    perturb( col, -hessian_eps );
+    const auto force_minus =
+        evaluator.compute_quadrature_point_penalty_data( pair, mesh1.getView(), mesh2.getView() ).force;
+    for ( int row = 0; row < 8; ++row ) {
+      const double fd_stiffness = ( force_plus[row] - force_minus[row] ) / ( 2.0 * hessian_eps );
+      EXPECT_NEAR( fd_stiffness, analytical.stiffness[row * 8 + col], hessian_tol )
+          << "stiffness mismatch at row " << row << ", col " << col;
+    }
+  }
+  restore();
+}
+
 TEST( GradientCheck, GtildeFDvsAD )
 {
   // ── Geometry: two facing LINEAR_EDGE segments ────────────────────────────
@@ -486,4 +602,4 @@ TEST( HessianCheck, GtildeFDvsAD )
   }
 }
 
-}  // namespace tribol
+}  // namespace tribol

src/tests/tribol_mfem_jacobian.cpp

@@ -1053,6 +1053,7 @@ TEST_F( MfemJacobianTest, mfem_penalty_jacobian_retrieval )
   tribol::registerMfemCouplingScheme( cs_id, mesh1_id, mesh2_id, mesh, coords, mortar_attrs, nonmortar_attrs,
                                       tribol::SURFACE_TO_SURFACE, tribol::NO_CASE, tribol::ENERGY_MORTAR,
                                       tribol::FRICTIONLESS, tribol::PENALTY, tribol::BINNING_GRID );
+  tribol::setEnforcementLocation( cs_id, tribol::EnforcementLocation::QuadraturePoint );
 
   tribol::setPenaltyOptions( cs_id, tribol::KINEMATIC, tribol::KINEMATIC_CONSTANT );
   tribol::setMfemKinematicConstantPenalty( cs_id, 1.0, 1.0 );

src/tribol/common/Parameters.hpp

@@ -155,6 +155,15 @@ enum EnforcementMethod
   NUM_ENFORCEMENT_METHODS
 };
 
+/*!
+ * \brief Enumerates where contact constraints are enforced
+ */
+enum class EnforcementLocation
+{
+  Nodal,           ///! Enforce contact using assembled nodal gaps
+  QuadraturePoint  ///! Enforce contact independently at quadrature points
+};
+
 /*!
  * \brief Enumerates the available spatial binning methods
  */
@@ -495,6 +504,9 @@ struct Parameters {
   // constituent face elements, then we don't consider the face-pair a contact candidate.
   // Note, auto-contact will require registration of element thicknesses.
   bool auto_contact_check = false;  ///! True if auto-contact checks should be enabled
+
+  EnforcementLocation enforcement_location =
+      EnforcementLocation::QuadraturePoint;  ///! Defaults to quadrature-point enforcement
 };
 
 }  // namespace tribol

src/tribol/interface/mfem_tribol.cpp

@@ -14,7 +14,6 @@
 
 #include "tribol.hpp"
 #include "tribol/mesh/CouplingScheme.hpp"
-#include "tribol/physics/ContactFormulationFactory.hpp"
 
 namespace tribol {
 
@@ -293,6 +292,7 @@ void registerMfemCouplingScheme( IndexT cs_id, int mesh_id_1, int mesh_id_2, con
   auto& cs = CouplingSchemeManager::getInstance().at( cs_id );
   cs.setMPIComm( mesh.GetComm() );
   if ( contact_method == ENERGY_MORTAR && enforcement_method == LAGRANGE_MULTIPLIER ) {
+    cs.getParameters().enforcement_location = EnforcementLocation::Nodal;
     SLIC_WARNING_ROOT(
         "ENERGY_MORTAR with Lagrange multiplier enforcement is experimental, has no testing, and has "
         "no support from Tribol developers." );
@@ -338,7 +338,7 @@ void registerMfemCouplingScheme( IndexT cs_id, int mesh_id_1, int mesh_id_2, con
   }
   cs.setMfemMeshData( std::move( mfem_data ) );
   if ( contact_method == ENERGY_MORTAR ) {
-    cs.setContactFormulation( createContactFormulation( &cs ) );
+    cs.updateContactFormulation();
   }
 }
 

src/tribol/interface/tribol.cpp

@@ -26,6 +26,7 @@
 #include "axom/slic.hpp"
 
 // C/C++ includes
+
 #include <string>
 #include <unordered_map>
 #include <fstream>
@@ -170,6 +171,20 @@ void setTimestepPenFrac( IndexT cs_id, RealT frac )
 
 }  // end setTimestepPenFrac()
 
+//------------------------------------------------------------------------------
+void setEnforcementLocation( IndexT cs_id, EnforcementLocation location )
+{
+  auto cs = CouplingSchemeManager::getInstance().findData( cs_id );
+
+  SLIC_ERROR_ROOT_IF( !cs, "tribol::setEnforcementLocation(): call tribol::registerCouplingScheme() "
+                               << "prior to calling this routine." );
+
+  cs->getParameters().enforcement_location = location;
+
+  // Automatically rebuild the formulation to reflect the new setting
+  cs->updateContactFormulation();
+}
+
 //------------------------------------------------------------------------------
 void setTimestepScale( IndexT cs_id, RealT scale )
 {

src/tribol/interface/tribol.hpp

@@ -131,6 +131,16 @@ void setAutoContactPenScale( IndexT cs_id, RealT scale );
  */
 void setTimestepPenFrac( IndexT cs_id, RealT frac );
 
+/*!
+ * \brief Sets where contact constraints are enforced
+ *
+ * \param [in] cs_id coupling scheme id
+ * \param [in] location contact enforcement location
+ *
+ * \note Quadrature-point computation is only supported with penalty enforcement.
+ */
+void setEnforcementLocation( IndexT cs_id, EnforcementLocation location );
+
 /*!
  *
  * \brief sets the timestep scale factor applied to the timestep vote

src/tribol/mesh/CouplingScheme.cpp

@@ -1155,6 +1155,18 @@ int CouplingScheme::apply( int cycle, RealT t, RealT& dt )
 
 }  // end CouplingScheme::apply()
 
+//------------------------------------------------------------------------------
+void CouplingScheme::updateContactFormulation()
+{
+  if ( m_contactMethod == ENERGY_MORTAR ) {
+    this->setContactFormulation( createContactFormulation( this ) );
+  } else {
+    SLIC_WARNING_ROOT(
+        "tribol::CouplingScheme::updateContactFormulation(): rebuilding is only supported for ENERGY_MORTAR at this "
+        "time." );
+  }
+}
+
 //------------------------------------------------------------------------------
 bool CouplingScheme::init()
 {