From 6f3b6cba37f21cf8d5e92dc065fbcad758d6d959 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Sun, 22 Feb 2026 22:02:09 -1000 Subject: [PATCH 1/7] initial implementation --- src/tests/CMakeLists.txt | 1 + src/tests/tribol_residual_gap.cpp | 179 ++++++++++++++++++++++ src/tribol/common/Parameters.hpp | 2 + src/tribol/geom/CompGeom.hpp | 2 +- src/tribol/interface/mfem_tribol.cpp | 4 +- src/tribol/interface/tribol.cpp | 13 ++ src/tribol/interface/tribol.hpp | 9 ++ src/tribol/mesh/CouplingScheme.hpp | 5 +- src/tribol/mesh/MfemData.cpp | 17 +- src/tribol/mesh/MfemData.hpp | 7 +- src/tribol/physics/AlignedMortar.cpp | 8 +- src/tribol/physics/AlignedMortar.hpp | 17 +- src/tribol/physics/CommonPlane.cpp | 3 +- src/tribol/physics/Mortar.cpp | 33 ++-- src/tribol/physics/Mortar.hpp | 12 +- src/tribol/search/InterfacePairFinder.cpp | 86 ++++------- src/tribol/search/InterfacePairFinder.hpp | 83 +++++++--- 17 files changed, 357 insertions(+), 124 deletions(-) create mode 100644 src/tests/tribol_residual_gap.cpp diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt index 150f1dfa..7ea037ff 100644 --- a/src/tests/CMakeLists.txt +++ b/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 diff --git a/src/tests/tribol_residual_gap.cpp b/src/tests/tribol_residual_gap.cpp new file mode 100644 index 00000000..5985b617 --- /dev/null +++ b/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( tribol::LINEAR_QUAD ), &x1[0], &y1[0], + &z1[0], tribol::MemorySpace::Host ); + tribol::registerMesh( 1, numElems, numVerts, &conn2[0], static_cast( 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; +} diff --git a/src/tribol/common/Parameters.hpp b/src/tribol/common/Parameters.hpp index d5dd2e6b..1dcbba71 100644 --- a/src/tribol/common/Parameters.hpp +++ b/src/tribol/common/Parameters.hpp @@ -496,6 +496,8 @@ 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 + + RealT residual_gap = 0.0; ///! User defined residual gap constraint }; } // namespace tribol diff --git a/src/tribol/geom/CompGeom.hpp b/src/tribol/geom/CompGeom.hpp index 80a2f460..d5657db1 100644 --- a/src/tribol/geom/CompGeom.hpp +++ b/src/tribol/geom/CompGeom.hpp @@ -1039,7 +1039,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; } } diff --git a/src/tribol/interface/mfem_tribol.cpp b/src/tribol/interface/mfem_tribol.cpp index 68874175..60ed363a 100644 --- a/src/tribol/interface/mfem_tribol.cpp +++ b/src/tribol/interface/mfem_tribol.cpp @@ -402,9 +402,11 @@ void updateMfemParallelDecomposition( int n_ranks, bool force_new_redecomp ) if ( mfem_data->GetLORFactor() > 1 ) { effective_binning_proximity *= static_cast( mfem_data->GetLORFactor() ); } + auto residual_gap = cs.getParameters().residual_gap; // 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 new_redecomp = + mfem_data->UpdateMfemMeshData( effective_binning_proximity, n_ranks, force_new_redecomp, residual_gap ); auto coord_ptrs = mfem_data->GetRedecompCoordsPtrs(); registerMesh( mesh_ids[0], mfem_data->GetMesh1NE(), mfem_data->GetNV(), mfem_data->GetMesh1Conn(), diff --git a/src/tribol/interface/tribol.cpp b/src/tribol/interface/tribol.cpp index bad3cbe3..134a5aa2 100644 --- a/src/tribol/interface/tribol.cpp +++ b/src/tribol/interface/tribol.cpp @@ -337,6 +337,19 @@ void setBinningProximityScale( IndexT cs_id, RealT binning_proximity_scale ) } // end setBinningProximityScale() +//------------------------------------------------------------------------------ +void setResidualGap( IndexT cs_id, RealT residual_gap ) +{ + auto cs = CouplingSchemeManager::getInstance().findData( cs_id ); + + // check to see if coupling scheme exists + SLIC_ERROR_ROOT_IF( + !cs, "tribol::setResidualGap(): call tribol::registerCouplingScheme() " << "prior to calling this routine." ); + + cs->getParameters().residual_gap = residual_gap; + +} // end setResidualGap() + //------------------------------------------------------------------------------ void enableTimestepVote( IndexT cs_id, const bool enable ) { diff --git a/src/tribol/interface/tribol.hpp b/src/tribol/interface/tribol.hpp index 4ac9e3be..1240a3f4 100644 --- a/src/tribol/interface/tribol.hpp +++ b/src/tribol/interface/tribol.hpp @@ -211,6 +211,15 @@ void setLoggingLevel( IndexT cs_id, LoggingLevel log_level ); */ void setBinningProximityScale( IndexT cs_id, RealT binning_proximity_scale ); +/*! + * @brief Sets the residual gap for a coupling scheme + * + * @param [in] cs_id coupling scheme id + * @param [in] residual_gap the gap offset. Positive values shift the contact surface away from the mesh + * surface, making contact occur earlier (with a gap). Effective gap = kinematic gap - residual gap. + */ +void setResidualGap( IndexT cs_id, RealT residual_gap ); + /*! * \brief Enable the contact timestep vote * diff --git a/src/tribol/mesh/CouplingScheme.hpp b/src/tribol/mesh/CouplingScheme.hpp index 23e7911f..1fcc2723 100644 --- a/src/tribol/mesh/CouplingScheme.hpp +++ b/src/tribol/mesh/CouplingScheme.hpp @@ -1013,8 +1013,9 @@ TRIBOL_HOST_DEVICE inline RealT CouplingScheme::Viewer::getGapTol( int fid1, int break; default: - gap_tol = -1. * m_parameters.gap_tol_ratio * - axom::utilities::max( m_mesh1.getFaceRadius()[fid1], m_mesh2.getFaceRadius()[fid2] ); + gap_tol = m_parameters.residual_gap - + m_parameters.gap_tol_ratio * + axom::utilities::max( m_mesh1.getFaceRadius()[fid1], m_mesh2.getFaceRadius()[fid2] ); break; } // end switch over m_contactModel diff --git a/src/tribol/mesh/MfemData.cpp b/src/tribol/mesh/MfemData.cpp index 3e552c3d..6c2e4d3b 100644 --- a/src/tribol/mesh/MfemData.cpp +++ b/src/tribol/mesh/MfemData.cpp @@ -348,7 +348,8 @@ void MfemMeshData::SetParentReferenceCoords( const mfem::ParGridFunction& refere } } -bool MfemMeshData::UpdateMfemMeshData( RealT binning_proximity_scale, int n_ranks, bool force_new_redecomp ) +bool MfemMeshData::UpdateMfemMeshData( RealT binning_proximity_scale, int n_ranks, bool force_new_redecomp, + RealT residual_gap ) { TRIBOL_MARK_FUNCTION; @@ -399,10 +400,10 @@ bool MfemMeshData::UpdateMfemMeshData( RealT binning_proximity_scale, int n_rank submesh_lor_xfer_->SubmeshToLOR( *submesh_nodes, *lor_nodes ); TRIBOL_MARK_END( "Update LOR coords" ); } - update_data_ = - std::make_unique( submesh_, lor_mesh_.get(), *coords_.GetParentGridFn().ParFESpace(), - submesh_xfer_gridfn_, submesh_lor_xfer_.get(), attributes_1_, attributes_2_, - binning_proximity_scale, n_ranks, allocator_id_, redecomp_trigger_displacement_ ); + update_data_ = std::make_unique( submesh_, lor_mesh_.get(), *coords_.GetParentGridFn().ParFESpace(), + submesh_xfer_gridfn_, submesh_lor_xfer_.get(), attributes_1_, + attributes_2_, binning_proximity_scale, n_ranks, allocator_id_, + redecomp_trigger_displacement_, residual_gap ); rebuilt = true; } @@ -644,19 +645,19 @@ MfemMeshData::UpdateData::UpdateData( mfem::ParSubMesh& submesh, mfem::ParMesh* mfem::ParGridFunction& submesh_gridfn, SubmeshLORTransfer* submesh_lor_xfer, const std::set& attributes_1, const std::set& attributes_2, RealT binning_proximity_scale, int n_ranks, int allocator_id, - RealT redecomp_trigger_displacement ) + RealT redecomp_trigger_displacement, RealT residual_gap ) : redecomp_mesh_{ lor_mesh ? redecomp::RedecompMesh( *lor_mesh, binning_proximity_scale * redecomp::RedecompMesh::MaxElementSize( *lor_mesh, redecomp::MPIUtility( lor_mesh->GetComm() ) ) + - redecomp_trigger_displacement, + redecomp_trigger_displacement + residual_gap, redecomp::RedecompMesh::RCB, n_ranks ) : redecomp::RedecompMesh( submesh, binning_proximity_scale * redecomp::RedecompMesh::MaxElementSize( submesh, redecomp::MPIUtility( submesh.GetComm() ) ) + - redecomp_trigger_displacement, + redecomp_trigger_displacement + residual_gap, redecomp::RedecompMesh::RCB, n_ranks ) }, vector_xfer_{ parent_fes, submesh_gridfn, submesh_lor_xfer, redecomp_mesh_ }, allocator_id_{ allocator_id } diff --git a/src/tribol/mesh/MfemData.hpp b/src/tribol/mesh/MfemData.hpp index de2846d8..88c011bc 100644 --- a/src/tribol/mesh/MfemData.hpp +++ b/src/tribol/mesh/MfemData.hpp @@ -639,11 +639,13 @@ class MfemMeshData { * @param n_ranks Number of ranks in the parallel decomposition * @param force_new_redecomp If true, construct a new RedecompMesh even if displacement threshold is not met (default * = false) + * @param residual_gap Residual gap to add to ghost layer thickness (default = 0.0) * @return True if a new RedecompMesh is created by this method * * @note This method should be called after the coordinate grid function is updated. */ - bool UpdateMfemMeshData( RealT binning_proximity_scale, int n_ranks, bool force_new_redecomp = false ); + bool UpdateMfemMeshData( RealT binning_proximity_scale, int n_ranks, bool force_new_redecomp = false, + RealT residual_gap = 0.0 ); /** * @brief Get the integer identifier for the first Tribol registered mesh @@ -1119,11 +1121,12 @@ class MfemMeshData { * @param allocator_id Allocation space ID for Tribol memory * @param redecomp_trigger_displacement Additional length to add to redecomp ghost length equal to the * displacement required to trigger a new RedecompMesh to be built. + * @param residual_gap user-defined residual gap */ UpdateData( mfem::ParSubMesh& submesh, mfem::ParMesh* lor_mesh, const mfem::ParFiniteElementSpace& parent_fes, mfem::ParGridFunction& submesh_gridfn, SubmeshLORTransfer* submesh_lor_xfer, const std::set& attributes_1, const std::set& attributes_2, RealT binning_proximity_scale, - int n_ranks, int allocator_id, RealT redecomp_trigger_displacement ); + int n_ranks, int allocator_id, RealT redecomp_trigger_displacement, RealT residual_gap ); /** * @brief Redecomposed boundary element mesh diff --git a/src/tribol/physics/AlignedMortar.cpp b/src/tribol/physics/AlignedMortar.cpp index a9aac9ec..5e5f471c 100644 --- a/src/tribol/physics/AlignedMortar.cpp +++ b/src/tribol/physics/AlignedMortar.cpp @@ -70,7 +70,7 @@ void ComputeAlignedMortarWeights( SurfaceContactElem& elem ) //------------------------------------------------------------------------------ template <> -void ComputeNodalGap( SurfaceContactElem& elem ) +void ComputeNodalGap( SurfaceContactElem& elem, RealT residual_gap ) { // get pointer to mesh view to store gaps on mesh data object auto& nonmortarMesh = *elem.m_mesh2; @@ -135,7 +135,7 @@ void ComputeNodalGap( SurfaceContactElem& elem ) v[1] = elem.faceCoords1[elem.dim * mortarNodeId + 1] - elem.faceCoords2[elem.dim * a + 1]; v[2] = elem.faceCoords1[elem.dim * mortarNodeId + 2] - elem.faceCoords2[elem.dim * a + 2]; - nonmortarMesh.getNodalFields().m_node_gap[glbId] += dotProd( &v[0], &nrml_a[0], elem.dim ); + nonmortarMesh.getNodalFields().m_node_gap[glbId] += dotProd( &v[0], &nrml_a[0], elem.dim ) - residual_gap; } } // end of ComputeNodalGap<>() @@ -168,6 +168,8 @@ void ComputeAlignedMortarGaps( CouplingScheme* cs ) Array2D mortarX( numNodesPerFace, dim ); Array2D nonmortarX( numNodesPerFace, dim ); + RealT residual_gap = cs->getParameters().residual_gap; + //////////////////////////// // compute nonmortar gaps // //////////////////////////// @@ -209,7 +211,7 @@ void ComputeAlignedMortarGaps( CouplingScheme* cs ) ///////////////////////// // compute mortar gaps // ///////////////////////// - ComputeNodalGap( elem_for_gap ); + ComputeNodalGap( elem_for_gap, residual_gap ); // HAVE TO set the number of active constraints. For now set to // all nonmortar face nodes. diff --git a/src/tribol/physics/AlignedMortar.hpp b/src/tribol/physics/AlignedMortar.hpp index c9998149..12ec72fd 100644 --- a/src/tribol/physics/AlignedMortar.hpp +++ b/src/tribol/physics/AlignedMortar.hpp @@ -43,23 +43,14 @@ void ComputeAlignedMortarWeights( SurfaceContactElem& elem ); * * \brief compute a contact element's contribution to nodal gaps * - * \param [in] cs pointer to the coupling scheme - * - */ -template -void ComputeNodalGap( SurfaceContactElem& elem ); - -/*! + * \note explicit specialization for aligned mortar method * - * \brief compute a contact element's contribution to nodal gaps - * - * \note explicit specialization for single mortar method - * - * \param [in] cs pointer to the coupling scheme + * \param [in] elem surface contact element object for contact face-pair + * \param [in] residual_gap user-defined residual gap * */ template <> -void ComputeNodalGap( SurfaceContactElem& elem ); +void ComputeNodalGap( SurfaceContactElem& elem, RealT residual_gap ); /*! * diff --git a/src/tribol/physics/CommonPlane.cpp b/src/tribol/physics/CommonPlane.cpp index 39eae36a..f5973718 100644 --- a/src/tribol/physics/CommonPlane.cpp +++ b/src/tribol/physics/CommonPlane.cpp @@ -205,7 +205,8 @@ int ApplyNormal( CouplingScheme* cs ) // compute total pressure based on constraint type RealT totalPressure = 0.; - plane.m_pressure = gap * penalty_stiff_per_area; // kinematic contribution + RealT residual_gap = cs_view.getParameters().residual_gap; + plane.m_pressure = ( gap - residual_gap ) * penalty_stiff_per_area; // kinematic contribution switch ( pen_enfrc_options.constraint_type ) { case KINEMATIC_AND_RATE: { // kinematic contribution diff --git a/src/tribol/physics/Mortar.cpp b/src/tribol/physics/Mortar.cpp index 91e85f1d..044a4d64 100644 --- a/src/tribol/physics/Mortar.cpp +++ b/src/tribol/physics/Mortar.cpp @@ -113,7 +113,7 @@ void ComputeMortarWeights( SurfaceContactElem& elem ) //------------------------------------------------------------------------------ template <> -void ComputeNodalGap( SurfaceContactElem& elem ) +void ComputeNodalGap( SurfaceContactElem& elem, RealT residual_gap ) { // check to make sure mortar weights have been computed locally // for the SurfaceContactElem object @@ -161,6 +161,9 @@ void ComputeNodalGap( SurfaceContactElem& elem ) g1 += dotProd( nrml_a.data(), &elem.faceCoords1[elem.dim * b], elem.dim ) * nab_1; g2 += dotProd( nrml_a.data(), &elem.faceCoords2[elem.dim * b], elem.dim ) * nab_2; + + // also integrate the residual gap contribution + g2 += residual_gap * nab_2; } // store local gap @@ -200,6 +203,8 @@ void ComputeSingleMortarGaps( CouplingScheme* cs ) Array2D mortarX_bar( numNodesPerFace, dim ); Array2D nonmortarX_bar( numNodesPerFace, dim ); + RealT residual_gap = cs->getParameters().residual_gap; + //////////////////////////////////////////////////////////////////// // compute nonmortar gaps to determine active set of contact dofs // //////////////////////////////////////////////////////////////////// @@ -247,7 +252,7 @@ void ComputeSingleMortarGaps( CouplingScheme* cs ) elem.faceCoords1 = mortarX.data(); elem.faceCoords2 = nonmortarX.data(); - ComputeNodalGap( elem ); + ComputeNodalGap( elem, residual_gap ); // TODO: fix this to register the actual number of active nonmortar gaps. // This is not the appropriate data structure to put this information in @@ -693,7 +698,8 @@ int ApplyNormalEnzyme( CouplingScheme* cs ) ComputeMortarJacobianEnzyme( x1, n1, p1, f1, blockJ( 0, 0 ).data(), blockJ( 0, 1 ).data(), blockJ_n( 0, 0 ).data(), blockJ( 0, 2 ).data(), g1, blockJ( 2, 0 ).data(), blockJ( 2, 1 ).data(), blockJ_n( 2, 0 ).data(), size1, x2, f2, blockJ( 1, 0 ).data(), - blockJ( 1, 1 ).data(), blockJ_n( 1, 0 ).data(), blockJ( 1, 2 ).data(), size2 ); + blockJ( 1, 1 ).data(), blockJ_n( 1, 0 ).data(), blockJ( 1, 2 ).data(), size2, + cs->getParameters().residual_gap ); if ( lm_opts.sparse_mode == SparseMode::MFEM_ELEMENT_DENSE ) { cs->getMethodData()->storeElemBlockJ( { elem1, elem2, elem1 }, blockJ ); @@ -704,7 +710,7 @@ int ApplyNormalEnzyme( CouplingScheme* cs ) } } else if ( lm_opts.eval_mode == ImplicitEvalMode::MORTAR_GAP || lm_opts.eval_mode == ImplicitEvalMode::MORTAR_RESIDUAL ) { - ComputeMortarForceEnzyme( x1, n1, p1, f1, g1, size1, x2, f2, size2 ); + ComputeMortarForceEnzyme( x1, n1, p1, f1, g1, size1, x2, f2, size2, cs->getParameters().residual_gap ); } for ( int i{ 0 }; i < size1; ++i ) { int node_id = mesh1.getGlobalNodeId( elem1, i ); @@ -726,7 +732,7 @@ int ApplyNormalEnzyme( CouplingScheme* cs ) //------------------------------------------------------------------------------ void ComputeMortarForceEnzyme( const RealT* x1, const RealT* n1, const RealT* p1, RealT* f1, RealT* g1, int size1, - const RealT* x2, RealT* f2, int size2 ) + const RealT* x2, RealT* f2, int size2, RealT residual_gap ) { // convention: elem1 = nonmortar element // elem2 = mortar element @@ -973,6 +979,9 @@ void ComputeMortarForceEnzyme( const RealT* x1, const RealT* n1, const RealT* p1 for ( int d{ 0 }; d < 3; ++d ) { g1[i] += n1[d * size1 + i] * gap_v[d]; } + for ( int j{ 0 }; j < size1; ++j ) { + g1[i] -= residual_gap * mortar_mat1[i * size1 + j]; + } } // compute nonmortar force contributions @@ -1004,7 +1013,7 @@ void ComputeMortarForceEnzyme( const RealT* x1, const RealT* n1, const RealT* p1 void ComputeMortarJacobianEnzyme( const RealT* x1, const RealT* n1, const RealT* p1, RealT* f1, RealT* df1dx1, RealT* df1dx2, RealT* df1dn1, RealT* df1dp1, RealT* g1, RealT* dg1dx1, RealT* dg1dx2, RealT* dg1dn1, int size1, const RealT* x2, RealT* f2, RealT* df2dx1, RealT* df2dx2, - RealT* df2dn1, RealT* df2dp1, int size2 ) + RealT* df2dn1, RealT* df2dp1, int size2, RealT residual_gap ) { RealT x1_dot[12] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }; for ( int i{ 0 }; i < size1 * 3; ++i ) { @@ -1019,7 +1028,8 @@ void ComputeMortarJacobianEnzyme( const RealT* x1, const RealT* n1, const RealT* TRIBOL_ENZYME_CONST, size1, TRIBOL_ENZYME_CONST, x2, TRIBOL_ENZYME_DUP, f2, &df2dx1[size1*3*i], - TRIBOL_ENZYME_CONST, size2); + TRIBOL_ENZYME_CONST, size2, + TRIBOL_ENZYME_CONST, residual_gap); // clang-format on x1_dot[i] = 0.0; } @@ -1036,7 +1046,8 @@ void ComputeMortarJacobianEnzyme( const RealT* x1, const RealT* n1, const RealT* TRIBOL_ENZYME_CONST, size1, TRIBOL_ENZYME_CONST, x2, TRIBOL_ENZYME_DUP, f2, &df2dn1[size1*3*i], - TRIBOL_ENZYME_CONST, size2); + TRIBOL_ENZYME_CONST, size2, + TRIBOL_ENZYME_CONST, residual_gap); // clang-format on n1_dot[i] = 0.0; } @@ -1053,7 +1064,8 @@ void ComputeMortarJacobianEnzyme( const RealT* x1, const RealT* n1, const RealT* TRIBOL_ENZYME_CONST, size1, TRIBOL_ENZYME_CONST, x2, TRIBOL_ENZYME_DUP, f2, &df2dp1[size1*3*i], - TRIBOL_ENZYME_CONST, size2); + TRIBOL_ENZYME_CONST, size2, + TRIBOL_ENZYME_CONST, residual_gap); // clang-format on p1_dot[i] = 0.0; } @@ -1070,7 +1082,8 @@ void ComputeMortarJacobianEnzyme( const RealT* x1, const RealT* n1, const RealT* TRIBOL_ENZYME_CONST, size1, TRIBOL_ENZYME_DUP, x2, x2_dot, TRIBOL_ENZYME_DUP, f2, &df2dx2[size2*3*i], - TRIBOL_ENZYME_CONST, size2); + TRIBOL_ENZYME_CONST, size2, + TRIBOL_ENZYME_CONST, residual_gap); // clang-format on x2_dot[i] = 0.0; } diff --git a/src/tribol/physics/Mortar.hpp b/src/tribol/physics/Mortar.hpp index efec3a9f..559eb5c5 100644 --- a/src/tribol/physics/Mortar.hpp +++ b/src/tribol/physics/Mortar.hpp @@ -49,10 +49,11 @@ void ComputeSingleMortarGaps( CouplingScheme* cs ); * \brief compute a contact element's contribution to nodal gaps * * \param [in] elem surface contact element object for contact face-pair + * \param [in] residual_gap user-defined residual gap * */ template -void ComputeNodalGap( SurfaceContactElem& elem ); +void ComputeNodalGap( SurfaceContactElem& elem, RealT residual_gap = 0. ); /*! * @@ -61,10 +62,11 @@ void ComputeNodalGap( SurfaceContactElem& elem ); * \note explicit specialization for single mortar method * * \param [in] elem surface contact element object for contact face-pair + * \param [in] residual_gap user-defined residual gap * */ template <> -void ComputeNodalGap( SurfaceContactElem& elem ); +void ComputeNodalGap( SurfaceContactElem& elem, RealT residual_gap ); /*! * @@ -185,9 +187,10 @@ int ApplyNormalEnzyme( CouplingScheme* cs ); * z3]) * @param [out] f2 Nodal forces for element 2 (stored by node, e.g. [x0, x1, x2, x3, y0, y1, y2, y3, z0, z1, z2, z3]) * @param [in] size2 Number of nodes on element 2 + * @param [in] residual_gap user-defined residual gap */ void ComputeMortarForceEnzyme( const RealT* x1, const RealT* n1, const RealT* p1, RealT* f1, RealT* g1, int size1, - const RealT* x2, RealT* f2, int size2 ); + const RealT* x2, RealT* f2, int size2, RealT residual_gap = 0.0 ); /** * @brief Computes the frictionless mortar forces and Jacobian for a 3D quad element (following Puso and Laursen (2004)) @@ -226,11 +229,12 @@ void ComputeMortarForceEnzyme( const RealT* x1, const RealT* n1, const RealT* p1 * @param [out] df2dp1 Derivative of element 2 nodal forces with respect to element 1 nodal pressures (size = * num_nodes_per_elem^2 x spatial dim) * @param [in] size2 Number of nodes on element 2 + * @param [in] residual_gap user-defined residual gap */ void ComputeMortarJacobianEnzyme( const RealT* x1, const RealT* n1, const RealT* p1, RealT* f1, RealT* df1dx1, RealT* df1dx2, RealT* df1dn1, RealT* df1dp1, RealT* g1, RealT* dg1dx1, RealT* dg1dx2, RealT* dg1dn1, int size1, const RealT* x2, RealT* f2, RealT* df2dx1, RealT* df2dx2, - RealT* df2dn1, RealT* df2dp1, int size2 ); + RealT* df2dn1, RealT* df2dp1, int size2, RealT residual_gap = 0.0 ); #endif /*! diff --git a/src/tribol/search/InterfacePairFinder.cpp b/src/tribol/search/InterfacePairFinder.cpp index 469457c2..50b1af21 100644 --- a/src/tribol/search/InterfacePairFinder.cpp +++ b/src/tribol/search/InterfacePairFinder.cpp @@ -23,27 +23,6 @@ namespace spin = axom::spin; namespace tribol { -/*! - * \brief Base class to compute the candidate pairs for a coupling scheme - * - * \a initialize() must be called prior to \a findInterfacePairs() - * - */ -class SearchBase { - public: - SearchBase() {}; - virtual ~SearchBase() {}; - /*! - * Prepares the object for spatial searches - */ - virtual void initialize() = 0; - - /*! - * Find candidates in first mesh for each element in second mesh of coupling scheme. - */ - virtual void findInterfacePairs() = 0; -}; - /////////////////////////////////////////////////////////////////////////////// /*! @@ -64,7 +43,7 @@ class CartesianProduct : public SearchBase { * Constructs a CartesianProduct instance over CouplingScheme \a couplingScheme * \pre couplingScheme is not null */ - CartesianProduct( CouplingScheme* couplingScheme ) : m_coupling_scheme( couplingScheme ) {} + CartesianProduct( CouplingScheme* couplingScheme ) : SearchBase( couplingScheme ) {} void initialize() override {} @@ -92,10 +71,9 @@ class CartesianProduct : public SearchBase { ArrayT countArray( 1, 1, m_coupling_scheme->getAllocatorId() ); int* pCount = countArray.data(); - const auto cs_view = m_coupling_scheme->getView(); // count how many pairs are proximate forAllExec( m_coupling_scheme->getExecutionMode(), maxNumPairs, - [cs_view, mesh1NumElems, mesh2NumElems, is_symm, isProximate, pCount] TRIBOL_HOST_DEVICE( IndexT i ) { + [this, mesh1NumElems, mesh2NumElems, is_symm, isProximate, pCount] TRIBOL_HOST_DEVICE( IndexT i ) { IndexT fromIdx = i / mesh2NumElems; IndexT toIdx = i % mesh2NumElems; if ( is_symm ) { @@ -104,7 +82,7 @@ class CartesianProduct : public SearchBase { fromIdx = row; toIdx = i - offset; } - isProximate[i] = geomFilter( cs_view, fromIdx, toIdx ); + isProximate[i] = geomFilter( fromIdx, toIdx ); #ifdef TRIBOL_USE_RAJA RAJA::atomicAdd( pCount, static_cast( isProximate[i] ) ); #else @@ -156,8 +134,6 @@ class CartesianProduct : public SearchBase { << "." ); } - private: - CouplingScheme* m_coupling_scheme; }; // End of CartesianProduct definition /////////////////////////////////////////////////////////////////////////////// @@ -191,7 +167,7 @@ class GridSearch : public SearchBase { * \pre couplingScheme is not null */ GridSearch( CouplingScheme* couplingScheme ) - : m_coupling_scheme( couplingScheme ), + : SearchBase( couplingScheme ), m_mesh1( m_coupling_scheme->getMesh1().getView() ), m_mesh2( m_coupling_scheme->getMesh2().getView() ) { @@ -231,9 +207,10 @@ class GridSearch : public SearchBase { // * Find the average extents (range) of the boxes Assumption is that elements are roughly the same size // * Grid resolution for each dimension is overall box width divided by half the average element width SpaceVec ranges; + RealT residual_gap = m_coupling_scheme->getParameters().residual_gap; for ( int i = 0; i < m_mesh1.numberOfElements(); ++i ) { auto& bbox = m_meshBBoxes1[i]; - inflateBBox( bbox, e_binning_proximity_scale ); + inflateBBox( bbox, e_binning_proximity_scale, residual_gap ); ranges += bbox.range(); @@ -294,9 +271,10 @@ class GridSearch : public SearchBase { // Find matches in first mesh (with index 'fromIdx') // with candidate elements in second mesh (with index 'toIdx') // int k = 0; // Debug only + RealT residual_gap = m_coupling_scheme->getParameters().residual_gap; for ( int toIdx = 0; toIdx < m_mesh2.numberOfElements(); ++toIdx ) { SpatialBoundingBox bbox = elementBoundingBox( m_mesh2, toIdx ); - inflateBBox( bbox, e_binning_proximity_scale ); + inflateBBox( bbox, e_binning_proximity_scale, residual_gap ); // Query the mesh auto candidateBits = m_grid.getCandidates( bbox ); @@ -313,7 +291,7 @@ class GridSearch : public SearchBase { // TODO: Add extra filter by bbox // Preliminary geometry/proximity checks, SRW - bool contact = geomFilter( m_coupling_scheme->getView(), fromIdx, toIdx ); + bool contact = geomFilter( fromIdx, toIdx ); if ( contact ) { contactPairs.emplace_back( fromIdx, toIdx, true ); @@ -347,16 +325,15 @@ class GridSearch : public SearchBase { return box; } /*! - * Expands bounding box by range_multiplier * the longest dimension's range + * Expands bounding box by range_multiplier * the longest dimension's range + residual gap */ - void inflateBBox( SpatialBoundingBox& bbox, RealT range_multiplier ) + void inflateBBox( SpatialBoundingBox& bbox, RealT range_multiplier, RealT residual_gap = 0.0 ) { int d = bbox.getLongestDimension(); - const RealT expansionFac = range_multiplier * bbox.range()[d]; + const RealT expansionFac = range_multiplier * bbox.range()[d] + residual_gap; bbox.expand( expansionFac ); } - CouplingScheme* m_coupling_scheme; const MeshData::Viewer m_mesh1; const MeshData::Viewer m_mesh2; @@ -395,7 +372,7 @@ class BvhSearch : public SearchBase { * \pre couplingScheme is not null */ BvhSearch( CouplingScheme* coupling_scheme ) - : m_coupling_scheme( coupling_scheme ), + : SearchBase( coupling_scheme ), m_mesh1( m_coupling_scheme->getMesh1().getView() ), m_mesh2( m_coupling_scheme->getMesh2().getView() ), m_boxes1( axom::ArrayOptions::Uninitialized{}, m_mesh1.numberOfElements(), m_mesh1.numberOfElements(), @@ -416,11 +393,11 @@ class BvhSearch : public SearchBase { void initialize() override { // we want binning proximity scaled by LOR factor on HO meshes, i.e. the effective binning proximity - buildMeshBBoxes( m_boxes1, m_coupling_scheme->getMesh1().getView(), - m_coupling_scheme->getEffectiveBinningProximityScale() ); + auto e_binning_proximity_scale = m_coupling_scheme->getEffectiveBinningProximityScale(); + auto residual_gap = m_coupling_scheme->getParameters().residual_gap; + buildMeshBBoxes( m_boxes1, m_coupling_scheme->getMesh1().getView(), e_binning_proximity_scale, residual_gap ); // we want binning proximity scaled by LOR factor on HO meshes, i.e. the effective binning proximity - buildMeshBBoxes( m_boxes2, m_coupling_scheme->getMesh2().getView(), - m_coupling_scheme->getEffectiveBinningProximityScale() ); + buildMeshBBoxes( m_boxes2, m_coupling_scheme->getMesh2().getView(), e_binning_proximity_scale, residual_gap ); } // end initialize() /*! @@ -446,23 +423,21 @@ class BvhSearch : public SearchBase { // with device kernels ArrayT filtered_candidates_data( 1, 1, m_coupling_scheme->getAllocatorId() ); auto filtered_candidates = filtered_candidates_data.view(); - const auto cs_view = m_coupling_scheme->getView(); // count the number of filtered proximate pairs - forAllExec( - m_coupling_scheme->getExecutionMode(), m_candidates.size(), - [cs_view, offsets_view, counts_view, candidates_view, filtered_candidates] TRIBOL_HOST_DEVICE( IndexT i ) { - auto mesh1_elem = candidates_view[i]; - auto mesh2_elem = algorithm::binarySearch( offsets_view, counts_view, i ); - if ( geomFilter( cs_view, mesh1_elem, mesh2_elem ) ) { + forAllExec( m_coupling_scheme->getExecutionMode(), m_candidates.size(), + [this, offsets_view, counts_view, candidates_view, filtered_candidates] TRIBOL_HOST_DEVICE( IndexT i ) { + auto mesh1_elem = candidates_view[i]; + auto mesh2_elem = algorithm::binarySearch( offsets_view, counts_view, i ); + if ( geomFilter( mesh1_elem, mesh2_elem ) ) { #ifdef TRIBOL_USE_RAJA - RAJA::atomicInc( filtered_candidates.data() ); + RAJA::atomicInc( filtered_candidates.data() ); #else ++filtered_candidates[0]; #endif - } else { - candidates_view[i] = -1; - } - } ); + } else { + candidates_view[i] = -1; + } + } ); ArrayT filtered_candidates_host( filtered_candidates_data ); m_coupling_scheme->getInterfacePairs().resize( filtered_candidates_host[0] ); @@ -493,11 +468,11 @@ class BvhSearch : public SearchBase { } ); } // end findInterfacePairs() - void buildMeshBBoxes( ArrayT& boxes, const MeshData::Viewer& mesh, RealT binning_proximity ) + void buildMeshBBoxes( ArrayT& boxes, const MeshData::Viewer& mesh, RealT binning_proximity, RealT residual_gap ) { auto boxes_view = boxes.view(); forAllExec( m_coupling_scheme->getExecutionMode(), mesh.numberOfElements(), - [this, mesh, boxes_view, binning_proximity] TRIBOL_HOST_DEVICE( IndexT i ) { + [this, mesh, boxes_view, binning_proximity, residual_gap] TRIBOL_HOST_DEVICE( IndexT i ) { BoxT box; auto num_nodes_per_elem = mesh.numberOfNodesPerElement(); for ( IndexT j{ 0 }; j < num_nodes_per_elem; ++j ) { @@ -513,7 +488,7 @@ class BvhSearch : public SearchBase { mesh.getFaceNormal( i, vnorm ); VectorT faceNormal( vnorm ); RealT faceRadius = mesh.getFaceRadius()[i]; - expandBBoxNormal( box, faceNormal, binning_proximity * faceRadius ); + expandBBoxNormal( box, faceNormal, binning_proximity * faceRadius + residual_gap ); boxes_view[i] = std::move( box ); } ); } @@ -541,7 +516,6 @@ class BvhSearch : public SearchBase { */ TRIBOL_HOST_DEVICE void inflateBBox( BoxT& bbox, const RealT faceRadius ) { bbox.expand( faceRadius ); } - CouplingScheme* m_coupling_scheme; const MeshData::Viewer m_mesh1; const MeshData::Viewer m_mesh2; ArrayT m_boxes1; diff --git a/src/tribol/search/InterfacePairFinder.hpp b/src/tribol/search/InterfacePairFinder.hpp index 8a028b94..c8b5798a 100644 --- a/src/tribol/search/InterfacePairFinder.hpp +++ b/src/tribol/search/InterfacePairFinder.hpp @@ -6,34 +6,35 @@ #ifndef SRC_TRIBOL_SEARCH_INTERFACE_PAIR_FINDER_HPP_ #define SRC_TRIBOL_SEARCH_INTERFACE_PAIR_FINDER_HPP_ -#include "tribol/common/Parameters.hpp" -#include "tribol/mesh/MeshData.hpp" +#include "tribol/common/BasicTypes.hpp" #include "tribol/mesh/CouplingScheme.hpp" +#include "tribol/utils/Math.hpp" namespace tribol { -// Forward Declarations class SearchBase; -/// Free functions - -/*! - * \brief Basic geometry/proximity checks for face pairs - * - * \param [in] cs_view View of the coupling scheme - * \param [in] element_id1 id of 1st element in pair - * \param [in] element_id2 id of 2nd element in pair +/** + * @brief Performs a geometric filter on two elements to determine if they are a contact candidate * + * @param cs_view View of the coupling scheme + * @param fid1 ID of element in the first mesh + * @param fid2 ID of element in the second mesh + * @return true Elements are a contact candidate + * @return false Elements are NOT a contact candidate */ -TRIBOL_HOST_DEVICE inline bool geomFilter( const CouplingScheme::Viewer& cs_view, IndexT element_id1, - IndexT element_id2 ) +TRIBOL_HOST_DEVICE inline bool geomFilter( const CouplingScheme::Viewer& cs_view, const IndexT element_id1, + const IndexT element_id2 ) { auto& mesh1 = cs_view.getMesh1View(); auto& mesh2 = cs_view.getMesh2View(); bool auto_contact_check = cs_view.getParameters().auto_contact_check; // we want binning proximity scaled by LOR factor on HO meshes, i.e. the effective binning proximity - auto element_radius_multiplier = cs_view.getEffectiveBinningProximityScale(); + auto binning_proximity_scale = cs_view.getEffectiveBinningProximityScale(); + auto dim = cs_view.spatialDimension(); auto mode = cs_view.getContactMode(); + auto& params = cs_view.getParameters(); + RealT residual_gap = params.residual_gap; /// CHECK #1: Check to make sure the two face ids are not the same /// and the two mesh ids are not the same. @@ -41,8 +42,6 @@ TRIBOL_HOST_DEVICE inline bool geomFilter( const CouplingScheme::Viewer& cs_view return false; } - int dim = mesh1.spatialDimension(); - /// CHECK #2: Auto-contact precludes faces that share a common /// node(s). We want to preclude two adjacent faces from interacting /// due to problematic configurations, such as corners where the @@ -85,20 +84,21 @@ TRIBOL_HOST_DEVICE inline bool geomFilter( const CouplingScheme::Viewer& cs_view /// The face radii are taken to be the magnitude of the /// longest vector from that face's vertex averaged /// centroid to one its nodes. - RealT offset_tol = 0.05; + constexpr RealT offset_tol = 0.05; + if ( dim == 3 ) { + // get face radius off the mesh data RealT r1 = mesh1.getFaceRadius()[element_id1]; RealT r2 = mesh2.getFaceRadius()[element_id2]; - // set maximum offset of face centroids for inclusion - RealT distMax = element_radius_multiplier * ( r1 + r2 ); // default is sum of face radii + RealT distMax = binning_proximity_scale * ( r1 + r2 ) + residual_gap; // check if the contact mode is conforming, in which case the // faces are supposed to be aligned if ( mode == SURFACE_TO_SURFACE_CONFORMING ) { // use 5% of max face radius for conforming case as // tolerance on face offsets - distMax *= offset_tol; + distMax = offset_tol * binning_proximity_scale * ( r1 + r2 ) + residual_gap; } // compute the distance between the two face centroids @@ -117,14 +117,14 @@ TRIBOL_HOST_DEVICE inline bool geomFilter( const CouplingScheme::Viewer& cs_view RealT e1 = 0.5 * mesh1.getElementAreas()[element_id1]; RealT e2 = 0.5 * mesh2.getElementAreas()[element_id2]; - RealT distMax = element_radius_multiplier * ( e1 + e2 ); + RealT distMax = binning_proximity_scale * ( e1 + e2 ) + residual_gap; // check if the contact mode is conforming, in which case the // edges are supposed to be aligned if ( mode == SURFACE_TO_SURFACE_CONFORMING ) { // use 5% of max face radius for conforming case as // tolerance on face offsets - distMax *= offset_tol; + distMax = offset_tol * binning_proximity_scale * ( e1 + e2 ) + residual_gap; } // compute the distance between the two edge centroids @@ -147,8 +147,45 @@ TRIBOL_HOST_DEVICE inline bool geomFilter( const CouplingScheme::Viewer& cs_view // if we made it here we passed all checks return true; +} + +/*! + * \class SearchBase + * + * \brief This is the base class for the candidate search + */ +class SearchBase { + public: + SearchBase( CouplingScheme* cs ) : m_coupling_scheme( cs ) {} + + virtual ~SearchBase() {} -} // end geomFilter() + /*! + * Initializes the search strategy + */ + virtual void initialize() = 0; + + /*! + * Performs the candidate search + */ + virtual void findInterfacePairs() = 0; + + protected: + CouplingScheme* m_coupling_scheme; + + /** + * @brief Performs a geometric filter on two elements to determine if they are a contact candidate + * + * @param fid1 ID of element in the first mesh + * @param fid2 ID of element in the second mesh + * @return true Elements are a contact candidate + * @return false Elements are NOT a contact candidate + */ + TRIBOL_HOST_DEVICE inline bool geomFilter( const IndexT element_id1, const IndexT element_id2 ) const + { + return tribol::geomFilter( m_coupling_scheme->getView(), element_id1, element_id2 ); + } +}; /*! * \class InterfacePairFinder From 6f3a9eddfd66478b783fe414bcd7f745815739ce Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Sun, 22 Feb 2026 22:02:24 -1000 Subject: [PATCH 2/7] add mortar tests --- src/tests/tribol_enzyme_element_mortar.cpp | 37 ++++++++++++++++++++++ src/tests/tribol_mfem_mortar_lm.cpp | 14 +++++--- 2 files changed, 47 insertions(+), 4 deletions(-) diff --git a/src/tests/tribol_enzyme_element_mortar.cpp b/src/tests/tribol_enzyme_element_mortar.cpp index 34b05912..58fb108b 100644 --- a/src/tests/tribol_enzyme_element_mortar.cpp +++ b/src/tests/tribol_enzyme_element_mortar.cpp @@ -969,6 +969,43 @@ TEST_F( EnzymeElementMortarTest, NoOverlap ) FDCheck( x1, x2, n1, p1, x1_stencil, x2_stencil ); } +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, 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 //------------------------------------------------------------------------------ diff --git a/src/tests/tribol_mfem_mortar_lm.cpp b/src/tests/tribol_mfem_mortar_lm.cpp index 3c193e64..5dc33dfe 100644 --- a/src/tests/tribol_mfem_mortar_lm.cpp +++ b/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> { +class MfemMortarTest : public testing::TestWithParam> { protected: tribol::RealT max_disp_; void SetUp() override @@ -46,6 +46,8 @@ class MfemMortarTest : public testing::TestWithParam( GetParam() ); + // fixed options // boundary element attributes of mortar surface auto mortar_attrs = std::set( { 4 } ); @@ -145,6 +147,7 @@ class MfemMortarTest : public testing::TestWithParam( 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" From 61244dead72028cd0301034bab8b97f23247b277 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Sun, 22 Feb 2026 22:27:43 -1000 Subject: [PATCH 3/7] naming fixes --- src/tribol/search/InterfacePairFinder.hpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/tribol/search/InterfacePairFinder.hpp b/src/tribol/search/InterfacePairFinder.hpp index c8b5798a..5c52e1c0 100644 --- a/src/tribol/search/InterfacePairFinder.hpp +++ b/src/tribol/search/InterfacePairFinder.hpp @@ -18,13 +18,13 @@ class SearchBase; * @brief Performs a geometric filter on two elements to determine if they are a contact candidate * * @param cs_view View of the coupling scheme - * @param fid1 ID of element in the first mesh - * @param fid2 ID of element in the second mesh + * @param element_id1 ID of element in the first mesh + * @param element_id2 ID of element in the second mesh * @return true Elements are a contact candidate * @return false Elements are NOT a contact candidate */ -TRIBOL_HOST_DEVICE inline bool geomFilter( const CouplingScheme::Viewer& cs_view, const IndexT element_id1, - const IndexT element_id2 ) +TRIBOL_HOST_DEVICE inline bool geomFilter( const CouplingScheme::Viewer& cs_view, IndexT element_id1, + IndexT element_id2 ) { auto& mesh1 = cs_view.getMesh1View(); auto& mesh2 = cs_view.getMesh2View(); @@ -176,8 +176,8 @@ class SearchBase { /** * @brief Performs a geometric filter on two elements to determine if they are a contact candidate * - * @param fid1 ID of element in the first mesh - * @param fid2 ID of element in the second mesh + * @param element_id1 ID of element in the first mesh + * @param element_id2 ID of element in the second mesh * @return true Elements are a contact candidate * @return false Elements are NOT a contact candidate */ From 16f4ffdc01dfd10432c76b95129b6c41294024c4 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Sun, 22 Feb 2026 22:38:35 -1000 Subject: [PATCH 4/7] remove unneeded changes --- src/tribol/search/InterfacePairFinder.cpp | 61 ++++++++++++----- src/tribol/search/InterfacePairFinder.hpp | 80 +++++++---------------- 2 files changed, 66 insertions(+), 75 deletions(-) diff --git a/src/tribol/search/InterfacePairFinder.cpp b/src/tribol/search/InterfacePairFinder.cpp index 50b1af21..d21d7b54 100644 --- a/src/tribol/search/InterfacePairFinder.cpp +++ b/src/tribol/search/InterfacePairFinder.cpp @@ -23,6 +23,27 @@ namespace spin = axom::spin; namespace tribol { +/*! + * \brief Base class to compute the candidate pairs for a coupling scheme + * + * \a initialize() must be called prior to \a findInterfacePairs() + * + */ +class SearchBase { + public: + SearchBase() {}; + virtual ~SearchBase() {}; + /*! + * Prepares the object for spatial searches + */ + virtual void initialize() = 0; + + /*! + * Find candidates in first mesh for each element in second mesh of coupling scheme. + */ + virtual void findInterfacePairs() = 0; +}; + /////////////////////////////////////////////////////////////////////////////// /*! @@ -43,7 +64,7 @@ class CartesianProduct : public SearchBase { * Constructs a CartesianProduct instance over CouplingScheme \a couplingScheme * \pre couplingScheme is not null */ - CartesianProduct( CouplingScheme* couplingScheme ) : SearchBase( couplingScheme ) {} + CartesianProduct( CouplingScheme* couplingScheme ) : m_coupling_scheme( couplingScheme ) {} void initialize() override {} @@ -61,7 +82,6 @@ class CartesianProduct : public SearchBase { if ( is_symm ) { // account for symmetry: the max number of pairs when the meshes are the // same is the upper triangular portion of the cartesian product pair - // matrix maxNumPairs = mesh1NumElems * ( mesh1NumElems + 1 ) / 2; } ArrayT proximityArray( maxNumPairs, maxNumPairs, m_coupling_scheme->getAllocatorId() ); @@ -71,9 +91,10 @@ class CartesianProduct : public SearchBase { ArrayT countArray( 1, 1, m_coupling_scheme->getAllocatorId() ); int* pCount = countArray.data(); + const auto cs_view = m_coupling_scheme->getView(); // count how many pairs are proximate forAllExec( m_coupling_scheme->getExecutionMode(), maxNumPairs, - [this, mesh1NumElems, mesh2NumElems, is_symm, isProximate, pCount] TRIBOL_HOST_DEVICE( IndexT i ) { + [cs_view, mesh1NumElems, mesh2NumElems, is_symm, isProximate, pCount] TRIBOL_HOST_DEVICE( IndexT i ) { IndexT fromIdx = i / mesh2NumElems; IndexT toIdx = i % mesh2NumElems; if ( is_symm ) { @@ -82,7 +103,7 @@ class CartesianProduct : public SearchBase { fromIdx = row; toIdx = i - offset; } - isProximate[i] = geomFilter( fromIdx, toIdx ); + isProximate[i] = geomFilter( cs_view, fromIdx, toIdx ); #ifdef TRIBOL_USE_RAJA RAJA::atomicAdd( pCount, static_cast( isProximate[i] ) ); #else @@ -134,6 +155,8 @@ class CartesianProduct : public SearchBase { << "." ); } + private: + CouplingScheme* m_coupling_scheme; }; // End of CartesianProduct definition /////////////////////////////////////////////////////////////////////////////// @@ -167,7 +190,7 @@ class GridSearch : public SearchBase { * \pre couplingScheme is not null */ GridSearch( CouplingScheme* couplingScheme ) - : SearchBase( couplingScheme ), + : m_coupling_scheme( couplingScheme ), m_mesh1( m_coupling_scheme->getMesh1().getView() ), m_mesh2( m_coupling_scheme->getMesh2().getView() ) { @@ -291,7 +314,7 @@ class GridSearch : public SearchBase { // TODO: Add extra filter by bbox // Preliminary geometry/proximity checks, SRW - bool contact = geomFilter( fromIdx, toIdx ); + bool contact = geomFilter( m_coupling_scheme->getView(), fromIdx, toIdx ); if ( contact ) { contactPairs.emplace_back( fromIdx, toIdx, true ); @@ -334,6 +357,7 @@ class GridSearch : public SearchBase { bbox.expand( expansionFac ); } + CouplingScheme* m_coupling_scheme; const MeshData::Viewer m_mesh1; const MeshData::Viewer m_mesh2; @@ -372,7 +396,7 @@ class BvhSearch : public SearchBase { * \pre couplingScheme is not null */ BvhSearch( CouplingScheme* coupling_scheme ) - : SearchBase( coupling_scheme ), + : m_coupling_scheme( coupling_scheme ), m_mesh1( m_coupling_scheme->getMesh1().getView() ), m_mesh2( m_coupling_scheme->getMesh2().getView() ), m_boxes1( axom::ArrayOptions::Uninitialized{}, m_mesh1.numberOfElements(), m_mesh1.numberOfElements(), @@ -423,21 +447,23 @@ class BvhSearch : public SearchBase { // with device kernels ArrayT filtered_candidates_data( 1, 1, m_coupling_scheme->getAllocatorId() ); auto filtered_candidates = filtered_candidates_data.view(); + const auto cs_view = m_coupling_scheme->getView(); // count the number of filtered proximate pairs - forAllExec( m_coupling_scheme->getExecutionMode(), m_candidates.size(), - [this, offsets_view, counts_view, candidates_view, filtered_candidates] TRIBOL_HOST_DEVICE( IndexT i ) { - auto mesh1_elem = candidates_view[i]; - auto mesh2_elem = algorithm::binarySearch( offsets_view, counts_view, i ); - if ( geomFilter( mesh1_elem, mesh2_elem ) ) { + forAllExec( + m_coupling_scheme->getExecutionMode(), m_candidates.size(), + [cs_view, offsets_view, counts_view, candidates_view, filtered_candidates] TRIBOL_HOST_DEVICE( IndexT i ) { + auto mesh1_elem = candidates_view[i]; + auto mesh2_elem = algorithm::binarySearch( offsets_view, counts_view, i ); + if ( geomFilter( cs_view, mesh1_elem, mesh2_elem ) ) { #ifdef TRIBOL_USE_RAJA - RAJA::atomicInc( filtered_candidates.data() ); + RAJA::atomicInc( filtered_candidates.data() ); #else ++filtered_candidates[0]; #endif - } else { - candidates_view[i] = -1; - } - } ); + } else { + candidates_view[i] = -1; + } + } ); ArrayT filtered_candidates_host( filtered_candidates_data ); m_coupling_scheme->getInterfacePairs().resize( filtered_candidates_host[0] ); @@ -516,6 +542,7 @@ class BvhSearch : public SearchBase { */ TRIBOL_HOST_DEVICE void inflateBBox( BoxT& bbox, const RealT faceRadius ) { bbox.expand( faceRadius ); } + CouplingScheme* m_coupling_scheme; const MeshData::Viewer m_mesh1; const MeshData::Viewer m_mesh2; ArrayT m_boxes1; diff --git a/src/tribol/search/InterfacePairFinder.hpp b/src/tribol/search/InterfacePairFinder.hpp index 5c52e1c0..bbb625ec 100644 --- a/src/tribol/search/InterfacePairFinder.hpp +++ b/src/tribol/search/InterfacePairFinder.hpp @@ -6,22 +6,24 @@ #ifndef SRC_TRIBOL_SEARCH_INTERFACE_PAIR_FINDER_HPP_ #define SRC_TRIBOL_SEARCH_INTERFACE_PAIR_FINDER_HPP_ -#include "tribol/common/BasicTypes.hpp" +#include "tribol/common/Parameters.hpp" +#include "tribol/mesh/MeshData.hpp" #include "tribol/mesh/CouplingScheme.hpp" -#include "tribol/utils/Math.hpp" namespace tribol { +// Forward Declarations class SearchBase; -/** - * @brief Performs a geometric filter on two elements to determine if they are a contact candidate +/// Free functions + +/*! + * \brief Basic geometry/proximity checks for face pairs + * + * \param [in] cs_view View of the coupling scheme + * \param [in] element_id1 id of 1st element in pair + * \param [in] element_id2 id of 2nd element in pair * - * @param cs_view View of the coupling scheme - * @param element_id1 ID of element in the first mesh - * @param element_id2 ID of element in the second mesh - * @return true Elements are a contact candidate - * @return false Elements are NOT a contact candidate */ TRIBOL_HOST_DEVICE inline bool geomFilter( const CouplingScheme::Viewer& cs_view, IndexT element_id1, IndexT element_id2 ) @@ -30,11 +32,9 @@ TRIBOL_HOST_DEVICE inline bool geomFilter( const CouplingScheme::Viewer& cs_view auto& mesh2 = cs_view.getMesh2View(); bool auto_contact_check = cs_view.getParameters().auto_contact_check; // we want binning proximity scaled by LOR factor on HO meshes, i.e. the effective binning proximity - auto binning_proximity_scale = cs_view.getEffectiveBinningProximityScale(); - auto dim = cs_view.spatialDimension(); + auto element_radius_multiplier = cs_view.getEffectiveBinningProximityScale(); auto mode = cs_view.getContactMode(); - auto& params = cs_view.getParameters(); - RealT residual_gap = params.residual_gap; + RealT residual_gap = cs_view.getParameters().residual_gap; /// CHECK #1: Check to make sure the two face ids are not the same /// and the two mesh ids are not the same. @@ -42,6 +42,8 @@ TRIBOL_HOST_DEVICE inline bool geomFilter( const CouplingScheme::Viewer& cs_view return false; } + int dim = mesh1.spatialDimension(); + /// CHECK #2: Auto-contact precludes faces that share a common /// node(s). We want to preclude two adjacent faces from interacting /// due to problematic configurations, such as corners where the @@ -84,21 +86,20 @@ TRIBOL_HOST_DEVICE inline bool geomFilter( const CouplingScheme::Viewer& cs_view /// The face radii are taken to be the magnitude of the /// longest vector from that face's vertex averaged /// centroid to one its nodes. - constexpr RealT offset_tol = 0.05; - + RealT offset_tol = 0.05; if ( dim == 3 ) { - // get face radius off the mesh data RealT r1 = mesh1.getFaceRadius()[element_id1]; RealT r2 = mesh2.getFaceRadius()[element_id2]; - RealT distMax = binning_proximity_scale * ( r1 + r2 ) + residual_gap; + // set maximum offset of face centroids for inclusion + RealT distMax = element_radius_multiplier * ( r1 + r2 ) + residual_gap; // default is sum of face radii // check if the contact mode is conforming, in which case the // faces are supposed to be aligned if ( mode == SURFACE_TO_SURFACE_CONFORMING ) { // use 5% of max face radius for conforming case as // tolerance on face offsets - distMax = offset_tol * binning_proximity_scale * ( r1 + r2 ) + residual_gap; + distMax = offset_tol * element_radius_multiplier * ( r1 + r2 ) + residual_gap; } // compute the distance between the two face centroids @@ -117,14 +118,14 @@ TRIBOL_HOST_DEVICE inline bool geomFilter( const CouplingScheme::Viewer& cs_view RealT e1 = 0.5 * mesh1.getElementAreas()[element_id1]; RealT e2 = 0.5 * mesh2.getElementAreas()[element_id2]; - RealT distMax = binning_proximity_scale * ( e1 + e2 ) + residual_gap; + RealT distMax = element_radius_multiplier * ( e1 + e2 ) + residual_gap; // check if the contact mode is conforming, in which case the // edges are supposed to be aligned if ( mode == SURFACE_TO_SURFACE_CONFORMING ) { // use 5% of max face radius for conforming case as // tolerance on face offsets - distMax = offset_tol * binning_proximity_scale * ( e1 + e2 ) + residual_gap; + distMax = offset_tol * element_radius_multiplier * ( e1 + e2 ) + residual_gap; } // compute the distance between the two edge centroids @@ -147,45 +148,8 @@ TRIBOL_HOST_DEVICE inline bool geomFilter( const CouplingScheme::Viewer& cs_view // if we made it here we passed all checks return true; -} - -/*! - * \class SearchBase - * - * \brief This is the base class for the candidate search - */ -class SearchBase { - public: - SearchBase( CouplingScheme* cs ) : m_coupling_scheme( cs ) {} - - virtual ~SearchBase() {} - /*! - * Initializes the search strategy - */ - virtual void initialize() = 0; - - /*! - * Performs the candidate search - */ - virtual void findInterfacePairs() = 0; - - protected: - CouplingScheme* m_coupling_scheme; - - /** - * @brief Performs a geometric filter on two elements to determine if they are a contact candidate - * - * @param element_id1 ID of element in the first mesh - * @param element_id2 ID of element in the second mesh - * @return true Elements are a contact candidate - * @return false Elements are NOT a contact candidate - */ - TRIBOL_HOST_DEVICE inline bool geomFilter( const IndexT element_id1, const IndexT element_id2 ) const - { - return tribol::geomFilter( m_coupling_scheme->getView(), element_id1, element_id2 ); - } -}; +} // end geomFilter() /*! * \class InterfacePairFinder From ef79c317fc47cab1431a65d2575104babb253b18 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Sun, 22 Feb 2026 22:59:18 -1000 Subject: [PATCH 5/7] formatting --- src/tests/tribol_enzyme_element_mortar.cpp | 10 +++++----- src/tribol/physics/Mortar.cpp | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/tests/tribol_enzyme_element_mortar.cpp b/src/tests/tribol_enzyme_element_mortar.cpp index 58fb108b..7b000706 100644 --- a/src/tests/tribol_enzyme_element_mortar.cpp +++ b/src/tests/tribol_enzyme_element_mortar.cpp @@ -987,21 +987,21 @@ TEST_F( EnzymeElementMortarTest, ExactOverlapResidualGap ) // 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, 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) { + + for ( int i = 0; i < 4; ++i ) { EXPECT_NEAR( g1[i], -0.0125, 1e-12 ); } } diff --git a/src/tribol/physics/Mortar.cpp b/src/tribol/physics/Mortar.cpp index 044a4d64..ee1e1640 100644 --- a/src/tribol/physics/Mortar.cpp +++ b/src/tribol/physics/Mortar.cpp @@ -161,7 +161,7 @@ void ComputeNodalGap( SurfaceContactElem& elem, RealT residual_ga g1 += dotProd( nrml_a.data(), &elem.faceCoords1[elem.dim * b], elem.dim ) * nab_1; g2 += dotProd( nrml_a.data(), &elem.faceCoords2[elem.dim * b], elem.dim ) * nab_2; - + // also integrate the residual gap contribution g2 += residual_gap * nab_2; } From 66d598d3bd0e03242381ad13d2963b66c5c09025 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Tue, 1 Sep 2026 12:34:05 -0700 Subject: [PATCH 6/7] propagate residual gap calculation to energy mortar --- .../tribol_finite_diff_energy_mortar.cpp | 71 +++++++++++++++++++ src/tribol/interface/tribol.cpp | 3 + src/tribol/physics/ContactFormulation.hpp | 5 ++ .../physics/ContactFormulationFactory.cpp | 5 +- src/tribol/physics/EnergyMortar.cpp | 38 ++++++---- src/tribol/physics/EnergyMortar.hpp | 2 + src/tribol/physics/EnergyMortarAdapter.cpp | 11 ++- src/tribol/physics/EnergyMortarAdapter.hpp | 10 ++- 8 files changed, 126 insertions(+), 19 deletions(-) diff --git a/src/tests/tribol_finite_diff_energy_mortar.cpp b/src/tests/tribol_finite_diff_energy_mortar.cpp index 3dc059e8..4588b272 100644 --- a/src/tests/tribol_finite_diff_energy_mortar.cpp +++ b/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_ ); diff --git a/src/tribol/interface/tribol.cpp b/src/tribol/interface/tribol.cpp index e44f8e19..cc5fa28f 100644 --- a/src/tribol/interface/tribol.cpp +++ b/src/tribol/interface/tribol.cpp @@ -362,6 +362,9 @@ void setResidualGap( IndexT cs_id, RealT residual_gap ) !cs, "tribol::setResidualGap(): call tribol::registerCouplingScheme() " << "prior to calling this routine." ); cs->getParameters().residual_gap = residual_gap; + if ( cs->getContactFormulation() ) { + cs->getContactFormulation()->updateResidualGap( residual_gap ); + } } // end setResidualGap() diff --git a/src/tribol/physics/ContactFormulation.hpp b/src/tribol/physics/ContactFormulation.hpp index b64a2332..f0717140 100644 --- a/src/tribol/physics/ContactFormulation.hpp +++ b/src/tribol/physics/ContactFormulation.hpp @@ -119,6 +119,11 @@ class ContactFormulation { */ virtual void updateConstantPenaltyStiffness( double /*mesh1_penalty*/, double /*mesh2_penalty*/ ) {} + /** + * @brief Update the residual-gap offset on formulations that cache contact parameters + */ + virtual void updateResidualGap( RealT /*residual_gap*/ ) {} + #ifdef BUILD_REDECOMP /** * @brief Returns t-dof vector of forces on parent mesh diff --git a/src/tribol/physics/ContactFormulationFactory.cpp b/src/tribol/physics/ContactFormulationFactory.cpp index 54d857e8..373ac5e0 100644 --- a/src/tribol/physics/ContactFormulationFactory.cpp +++ b/src/tribol/physics/ContactFormulationFactory.cpp @@ -44,14 +44,15 @@ std::unique_ptr createContactFormulation( CouplingScheme* cs SLIC_ERROR_ROOT_IF( !cs->hasMfemJacobianData(), "ENERGY_MORTAR requires MFEM Jacobian data." ); const auto enforcement_location = cs->getParameters().enforcement_location; + const auto residual_gap = cs->getParameters().residual_gap; if ( enforcement_location == EnforcementLocation::QuadraturePoint ) { return std::make_unique>( *cs->getMfemMeshData(), *cs->getMfemSubmeshData(), *cs->getMfemJacobianData(), k, delta, N, - enzyme_quadrature, use_penalty ); + enzyme_quadrature, use_penalty, residual_gap ); } else { return std::make_unique>( *cs->getMfemMeshData(), *cs->getMfemSubmeshData(), *cs->getMfemJacobianData(), k, delta, N, enzyme_quadrature, - use_penalty ); + use_penalty, residual_gap ); } #else SLIC_ERROR_ROOT( "ENERGY_MORTAR requires Enzyme and redecomp to be built." ); diff --git a/src/tribol/physics/EnergyMortar.cpp b/src/tribol/physics/EnergyMortar.cpp index fa8d2424..41df0ea9 100644 --- a/src/tribol/physics/EnergyMortar.cpp +++ b/src/tribol/physics/EnergyMortar.cpp @@ -19,11 +19,17 @@ namespace { // This MUST match what the ContactParams struct has in EnergyMortarAdapter // These had to be saved locally in order for enzyme to work correctly struct KernelParams { - int N{ 3 }; // No. of quadrature points - double del{ 0.1 }; // Smoothing parameter - double k{ 1.0 }; // Penalty stiffness + int N{ 3 }; // No. of quadrature points + double del{ 0.1 }; // Smoothing parameter + double k{ 1.0 }; // Penalty stiffness + double residual_gap{ 0.0 }; // User-defined gap offset }; +TRIBOL_ENZYME_INLINE double effective_gap( double normal_gap, double normal_alignment, double residual_gap ) +{ + return normal_gap * normal_alignment - residual_gap; +} + // Return the line-element mapping Jacobian. Local edge coordinates span [-0.5, 0.5], so the Jacobian is the physical // length of the segment from A0 to A1. TRIBOL_ENZYME_INLINE double line_jacobian( const double* A0, const double* A1 ) @@ -286,7 +292,7 @@ TRIBOL_ENZYME_INLINE void gtilde_kernel( const double* x, Gparams* gp, double* g // lagged normal on B const double gn = -( dx * nB[0] + dy * nB[1] ); - const double g = gn * eta; + const double g = effective_gap( gn, eta, gp->residual_gap ); g1 += w * N1 * g * J; g2 += w * N2 * g * J; @@ -353,7 +359,7 @@ TRIBOL_ENZYME_INLINE void gtilde_kernel_quad( const double* x, const Gparams* gp // lagged normal on B const double gn = -( dx * nB[0] + dy * nB[1] ); - const double g = gn * eta; + const double g = effective_gap( gn, eta, gp->residual_gap ); g1 += w * N1 * g * J; g2 += w * N2 * g * J; @@ -444,6 +450,7 @@ static void kernel_out_enzyme( const double* x, const void* kp_void, double* out EnergyMortarCalculator::compute_quadrature( xi_bounds, kp->N, &qp ); Gparams gp; + gp.residual_gap = kp->residual_gap; for ( std::size_t i = 0; i < qp.qp.size(); ++i ) { gp.qp[i] = qp.qp[i]; gp.w[i] = qp.w[i]; @@ -503,7 +510,7 @@ void d2_kernel( const double* x, const KernelParams* kp, double* H ) // Isolate loop-local arrays to avoid a leak in Enzyme's reverse-mode tape. TRIBOL_ENZYME_INLINE double qp_penalty_kernel_qp_energy( double xiA, double w, const double* A0, const double* A1, const double* B0, const double* B1, const double* nB, - double eta, double penalty, double J ) + double eta, double residual_gap, double penalty, double J ) { double x1[2]; iso_map( A0, A1, xiA, x1 ); @@ -514,7 +521,7 @@ TRIBOL_ENZYME_INLINE double qp_penalty_kernel_qp_energy( double xiA, double w, c const double dx = x1[0] - x2[0]; const double dy = x1[1] - x2[1]; const double gn = -( dx * nB[0] + dy * nB[1] ); - const double gap = gn * eta; + const double gap = effective_gap( gn, eta, residual_gap ); return gap < 0.0 ? 0.5 * penalty * gap * gap * w * J : 0.0; } @@ -548,7 +555,7 @@ TRIBOL_ENZYME_INLINE void qp_penalty_kernel( const double* x, const KernelParams double value = 0.0; for ( int i = 0; i < kp->N; ++i ) { - value += qp_penalty_kernel_qp_energy( qp.qp[i], qp.w[i], A0, A1, B0, B1, nB, eta, kp->k, J ); + value += qp_penalty_kernel_qp_energy( qp.qp[i], qp.w[i], A0, A1, B0, B1, nB, eta, kp->residual_gap, kp->k, J ); } *energy = value; @@ -638,6 +645,7 @@ Gparams EnergyMortarCalculator::construct_gparams( const InterfacePair& pair, co } Gparams gp; + gp.residual_gap = p_.residual_gap; // int N = eval.get_N(); for ( std::size_t i = 0; i < qp.qp.size(); ++i ) { @@ -750,7 +758,7 @@ double EnergyMortarCalculator::compute_weighted_normal_gap( const InterfacePair& double dot = nB[0] * nA[0] + nB[1] * nA[1]; double eta = ( dot < 0 ) ? dot : 0.0; - return gn * eta; + return effective_gap( gn, eta, p_.residual_gap ); } // Assemble nodal gap and tributary area data for the current interface pair. @@ -844,7 +852,7 @@ void EnergyMortarCalculator::grad_gtilde( const InterfacePair& pair, const MeshD } else { // Differentiate through the geometry-dependent quadrature construction. - const KernelParams kp{ p_.N, p_.del, p_.k }; + const KernelParams kp{ p_.N, p_.del, p_.k, p_.residual_gap }; grad_kernel_enzyme( x, &kp, dg1_du ); grad_kernel_enzyme( x, &kp, dg2_du ); } @@ -877,7 +885,7 @@ void EnergyMortarCalculator::grad_trib_area( const InterfacePair& pair, const Me grad_kernel( x, &gp, dA2_dx ); } else { // Differentiate through the geometry-dependent quadrature construction. - const KernelParams kp{ p_.N, p_.del, p_.k }; + const KernelParams kp{ p_.N, p_.del, p_.k, p_.residual_gap }; grad_kernel_enzyme( x, &kp, dA1_dx ); grad_kernel_enzyme( x, &kp, dA2_dx ); } @@ -909,7 +917,7 @@ void EnergyMortarCalculator::d2_g2tilde( const InterfacePair& pair, const MeshDa } else { // Differentiate through the geometry-dependent quadrature construction. - const KernelParams kp{ p_.N, p_.del, p_.k }; + const KernelParams kp{ p_.N, p_.del, p_.k, p_.residual_gap }; d2_kernel( x, &kp, d2g1_d2u ); d2_kernel( x, &kp, d2g2_d2u ); } @@ -946,7 +954,7 @@ void EnergyMortarCalculator::compute_d2A_d2u( const InterfacePair& pair, const M d2_kernel_quad( x, &gp, d2A2_d2u ); } else { // Differentiate through the geometry-dependent quadrature construction. - const KernelParams kp{ p_.N, p_.del, p_.k }; + const KernelParams kp{ p_.N, p_.del, p_.k, p_.residual_gap }; d2_kernel( x, &kp, d2A1_d2u ); d2_kernel( x, &kp, d2A2_d2u ); } @@ -967,7 +975,7 @@ double EnergyMortarCalculator::compute_quadrature_point_penalty_energy( const In endpoints( mesh2, pair.m_element_id2, B0, B1 ); const double x[8] = { A0[0], A0[1], A1[0], A1[1], B0[0], B0[1], B1[0], B1[1] }; - const KernelParams kp{ p_.N, p_.del, p_.k }; + const KernelParams kp{ p_.N, p_.del, p_.k, p_.residual_gap }; double energy = 0.0; qp_penalty_kernel( x, &kp, &energy ); return energy; @@ -982,7 +990,7 @@ QuadraturePointPenaltyData EnergyMortarCalculator::compute_quadrature_point_pena endpoints( mesh2, pair.m_element_id2, B0, B1 ); const double x[8] = { A0[0], A0[1], A1[0], A1[1], B0[0], B0[1], B1[0], B1[1] }; - const KernelParams kp{ p_.N, p_.del, p_.k }; + const KernelParams kp{ p_.N, p_.del, p_.k, p_.residual_gap }; QuadraturePointPenaltyData result; qp_penalty_kernel( x, &kp, &result.energy ); diff --git a/src/tribol/physics/EnergyMortar.hpp b/src/tribol/physics/EnergyMortar.hpp index cf0f3217..044b22fb 100644 --- a/src/tribol/physics/EnergyMortar.hpp +++ b/src/tribol/physics/EnergyMortar.hpp @@ -24,6 +24,7 @@ struct ContactParams { double k; ///< Penalty stiffness. int N; ///< Number of quadrature points. bool enzyme_quadrature; ///< Whether Enzyme differentiates the quadrature construction. + double residual_gap{ 0.0 }; ///< User-defined gap offset subtracted from the kinematic gap. }; /// Stores quadrature-point penalty energy derivatives for one interface pair. @@ -81,6 +82,7 @@ struct FiniteDiffResult { struct Gparams { std::array qp; ///< Quadrature-point locations in the integration-edge local coordinate. std::array w; ///< Quadrature weights mapped to the local integration interval. + double residual_gap{ 0.0 }; ///< User-defined gap offset subtracted from the kinematic gap. }; /// Provides smoothing operations for the Energy Mortar contact formulation. diff --git a/src/tribol/physics/EnergyMortarAdapter.cpp b/src/tribol/physics/EnergyMortarAdapter.cpp index b58fe61a..fef91de5 100644 --- a/src/tribol/physics/EnergyMortarAdapter.cpp +++ b/src/tribol/physics/EnergyMortarAdapter.cpp @@ -14,7 +14,8 @@ namespace tribol { template