diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt index 118ca0d9..271c9832 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_enzyme_element_mortar.cpp b/src/tests/tribol_enzyme_element_mortar.cpp index fed9b902..562b038a 100644 --- a/src/tests/tribol_enzyme_element_mortar.cpp +++ b/src/tests/tribol_enzyme_element_mortar.cpp @@ -974,6 +974,43 @@ TEST_F( EnzymeElementMortarTest, NoOverlap ) FDCheck( x1, x2, n1, p1, x1_stencil, x2_stencil, num_nodes, check_scale, len_collapse_ratio ); } +TEST_F( EnzymeElementMortarTest, ExactOverlapResidualGap ) +{ + // clang-format off + // Setup two overlapping elements with a gap of 0.1 + double x1[12] = { 0.0, 1.0, 1.0, 0.0, + 0.0, 0.0, 1.0, 1.0, + 0.0, 0.0, 0.0, 0.0 }; + // Define x2 in CW order so ElemReverse makes it CCW + double x2[12] = { 0.0, 0.0, 1.0, 1.0, + 0.0, 1.0, 1.0, 0.0, + 0.1, 0.1, 0.1, 0.1 }; + double n1[12] = { 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, + 1.0, 1.0, 1.0, 1.0 }; + double p1[4] = { 1.0, 1.0, 1.0, 1.0 }; + // clang-format on + + double residual_gap = 0.15; + + double f1[12], f2[12], g1[4]; + int num_nodes = 4; + + // Call with residual_gap + tribol::ComputeMortarForceEnzyme( x1, n1, p1, f1, g1, num_nodes, x2, f2, num_nodes, 1.0e-8, residual_gap ); + + // Kinematic gap is 0.1. Effective gap is 0.1 - 0.15 = -0.05. + // Integrated gap for unit square (overlap area 1.0) distributed to 4 nodes: + // Each node i has integrated gap G_i = \int \phi_i (g_kin - g_r) dA. + // For constant g_kin and g_r: G_i = (g_kin - g_r) \int \phi_i dA. + // \int \phi_i dA for bilinear quad node is 0.25 (1/4 of element area). + // So G_i = -0.05 * 0.25 = -0.0125. + + for ( int i = 0; i < 4; ++i ) { + EXPECT_NEAR( g1[i], -0.0125, 1e-12 ); + } +} + } // namespace tribol //------------------------------------------------------------------------------ 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/tests/tribol_mfem_mortar_lm.cpp b/src/tests/tribol_mfem_mortar_lm.cpp index 98db248b..46f6f613 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 } ); @@ -144,6 +146,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" 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 df31b3d2..bf9b532f 100644 --- a/src/tribol/common/Parameters.hpp +++ b/src/tribol/common/Parameters.hpp @@ -505,6 +505,8 @@ struct Parameters { // Note, auto-contact will require registration of element thicknesses. bool auto_contact_check = false; ///! True if auto-contact checks should be enabled + RealT residual_gap = 0.0; ///! User defined residual gap constraint + EnforcementLocation enforcement_location = EnforcementLocation::QuadraturePoint; ///! Defaults to quadrature-point enforcement }; diff --git a/src/tribol/geom/CompGeom.hpp b/src/tribol/geom/CompGeom.hpp index cffbb483..b50575fb 100644 --- a/src/tribol/geom/CompGeom.hpp +++ b/src/tribol/geom/CompGeom.hpp @@ -1034,7 +1034,7 @@ TRIBOL_HOST_DEVICE inline bool CommonPlanePair::exceedsMaxAutoInterpen( const Me RealT max_interpen = -1. * params.auto_contact_pen_frac * axom::utilities::min( mesh1.getElementData().m_thickness[faceId1], mesh2.getElementData().m_thickness[faceId2] ); - if ( gap < max_interpen ) { + if ( ( gap - params.residual_gap ) < max_interpen ) { return true; } } diff --git a/src/tribol/interface/mfem_tribol.cpp b/src/tribol/interface/mfem_tribol.cpp index d130084f..3fc7ed74 100644 --- a/src/tribol/interface/mfem_tribol.cpp +++ b/src/tribol/interface/mfem_tribol.cpp @@ -775,9 +775,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 aca5eb42..cc5fa28f 100644 --- a/src/tribol/interface/tribol.cpp +++ b/src/tribol/interface/tribol.cpp @@ -352,6 +352,22 @@ 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; + if ( cs->getContactFormulation() ) { + cs->getContactFormulation()->updateResidualGap( 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 306c876a..51aa214d 100644 --- a/src/tribol/interface/tribol.hpp +++ b/src/tribol/interface/tribol.hpp @@ -225,6 +225,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 7114e067..9fe8a7bf 100644 --- a/src/tribol/mesh/CouplingScheme.hpp +++ b/src/tribol/mesh/CouplingScheme.hpp @@ -1050,8 +1050,9 @@ TRIBOL_HOST_DEVICE inline RealT CouplingScheme::Viewer::getGapTol( IndexT fid1, 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 e3cda2d8..64f35cf5 100644 --- a/src/tribol/mesh/MfemData.cpp +++ b/src/tribol/mesh/MfemData.cpp @@ -721,7 +721,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; @@ -771,10 +772,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; } @@ -1046,19 +1047,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 0b0a4d14..f5ffd922 100644 --- a/src/tribol/mesh/MfemData.hpp +++ b/src/tribol/mesh/MfemData.hpp @@ -654,11 +654,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 @@ -1134,11 +1136,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 fb216298..f1fe4e67 100644 --- a/src/tribol/physics/AlignedMortar.cpp +++ b/src/tribol/physics/AlignedMortar.cpp @@ -72,7 +72,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; @@ -140,7 +140,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<>() @@ -173,6 +173,8 @@ void ComputeAlignedMortarGaps( CouplingScheme* cs ) Array2D mortarX( numNodesPerFace, dim ); Array2D nonmortarX( numNodesPerFace, dim ); + RealT residual_gap = cs->getParameters().residual_gap; + //////////////////////////// // compute nonmortar gaps // //////////////////////////// @@ -214,7 +216,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 5278e3b5..66ec2190 100644 --- a/src/tribol/physics/CommonPlane.cpp +++ b/src/tribol/physics/CommonPlane.cpp @@ -206,7 +206,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/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..5d0d472a 100644 --- a/src/tribol/physics/EnergyMortar.hpp +++ b/src/tribol/physics/EnergyMortar.hpp @@ -20,10 +20,11 @@ struct QuadPoints { /// Parameters controlling ENERGY_MORTAR contact evaluation. struct ContactParams { - double del; ///< Smoothing length used for integration bounds. - double k; ///< Penalty stiffness. - int N; ///< Number of quadrature points. - bool enzyme_quadrature; ///< Whether Enzyme differentiates the quadrature construction. + double del; ///< Smoothing length used for integration bounds. + 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. @@ -79,8 +80,9 @@ struct FiniteDiffResult { /// Stores fixed quadrature data passed to differentiated kernels. 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. + 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