diff --git a/docs/common_plane_contact_mass_timestep_plan.md b/docs/common_plane_contact_mass_timestep_plan.md new file mode 100644 index 00000000..aea1fc0d --- /dev/null +++ b/docs/common_plane_contact_mass_timestep_plan.md @@ -0,0 +1,75 @@ +# Plan: Common Plane Contact Mass-Based Timestep Limit + +## Summary + +Refactor the `COMMON_PLANE + PENALTY` timestep stability limit from a host-`dt` scaling into an absolute contact timestep estimate based on host-provided **lumped nodal masses**. Implement the core registration path in `tribol.hpp` and add an MFEM convenience wrapper in `mfem_tribol.hpp`. + +This work should be implemented **after** higher-order Common Plane quadrature is added, since both changes touch overlap integration, contact weighting, and how overlap contributions are accumulated. + +## Key Changes + +- Replace the current CFL-like stability vote with an absolute contact timestep for active Common Plane penalty contact: + \[ + \Delta t_{\text{contact}} = \alpha \, 2 \sqrt{\frac{m_{\text{pair}}}{k_{\text{pair}}}} + \] + where `alpha = timestep_scale`, `k_pair = k_c * overlap_area`, and `k_c` is the existing interface stiffness per area. +- Compute side effective mass from **lumped nodal masses** using the same contact weights already used to distribute force: + \[ + m_{i,\text{eff}} = \left(\sum_a \frac{\phi_{i,a}^2}{m_{i,a}}\right)^{-1} + \] + using the overlap-based weak-form weights for each side. +- Combine the two sides into a pair effective mass with harmonic reduction: + \[ + m_{\text{pair}} = \left(\frac{1}{m_{1,\text{eff}}} + \frac{1}{m_{2,\text{eff}}}\right)^{-1} + \] +- Do not add extra overlap-area fraction scaling for mass in v1. The overlap localization is handled through the contact weights. +- Keep the existing gap-based and velocity-projection timestep checks unchanged. Only the stability-limit bucket changes from a `dt_in` scaler to an absolute contact vote. + +## Public API / Interface Changes + +- Add a core registration function in `tribol.hpp` / `tribol.cpp`: + - `registerLumpedNodalMass( IndexT mesh_id, const RealT* lumped_nodal_mass )` +- Add storage and validation for lumped nodal mass in mesh nodal data. + - Scalar field, length = number of contact surface nodes on the registered mesh. + - All values must be strictly positive. +- Add an MFEM convenience wrapper in `mfem_tribol.hpp` / `mfem_tribol.cpp`: + - `registerMfemLumpedNodalMass( IndexT cs_id, const mfem::ParGridFunction& lumped_mass )` +- MFEM helper behavior: + - Accept a scalar nodal `ParGridFunction` on the parent mesh. + - Transfer/project it onto the Tribol contact surface nodes using the existing parent-to-submesh machinery. +- Validation rule: + - If `enableTimestepStabilityLimits()` is on for a Common Plane penalty scheme, lumped nodal mass must be registered and valid on both meshes. + - Missing or nonpositive nodal mass makes the coupling scheme invalid during `init()`, with a warning, matching existing invalid-enforcement-data behavior. + +## Implementation Notes + +- Reuse the existing overlap weights for the effective-mass calculation so the inertial reduction follows the same contact discretization as the force calculation. +- Use the existing stiffness path: + - `KINEMATIC_CONSTANT`: `k_i = pen_scale * const_penalty` + - `KINEMATIC_ELEMENT`: `k_i = pen_scale * mat_mod / thickness` + - `k_c = (k1 * k2) / (k1 + k2)` + - `k_pair = k_c * plane.m_area` +- Guard all effective-mass calculations against zero or invalid mass before division. +- Keep the existing invalid-thickness early return in Common Plane penalty. +- Do not support consistent mass in v1 and do not implement internal mass lumping. The host is responsible for providing already-lumped nodal masses. +- When higher-order quadrature lands, generalize the single-point weight usage to a quadrature-point accumulation rather than reintroducing a single-point assumption. + +## Test Plan + +- Add/update core tests for Common Plane penalty: + - Constant-penalty contact with registered nodal masses: verify absolute contact timestep matches `2*sqrt(m_pair/k_pair)` times `timestep_scale`. + - Element-penalty contact with registered nodal masses: same verification using element-based stiffness. + - Stationary contact with stability limits enabled: verify the absolute contact limit is computed from mass/stiffness and is not a function of incoming `dt`. + - Missing nodal mass with stability limits enabled: verify the coupling scheme is skipped as invalid. + - Nonpositive nodal mass with stability limits enabled: verify invalid scheme and no contact response. +- Add MFEM-path coverage: + - Register MFEM lumped nodal mass through the new helper and verify the resulting timestep vote matches the core formula on a simple contact case. +- Re-run existing `tribol_timestep_vote` and `tribol_common_plane_penalty` suites to ensure the penetration-based votes still behave the same. + +## Assumptions + +- V1 supports only **lumped nodal masses** as inertial input. +- V1 computes contact mass from contact weights and does not apply extra overlap-fraction scaling. +- V1 adds both the core Tribol API support and an MFEM convenience registration helper. +- `timestep_scale` remains the user-facing safety factor and multiplies the absolute contact critical timestep. +- The stability constant is `2`, corresponding to the lumped-mass single-DOF spring stability bound. diff --git a/src/tests/tribol_common_plane_penalty.cpp b/src/tests/tribol_common_plane_penalty.cpp index 5f264725..d847d0df 100644 --- a/src/tests/tribol_common_plane_penalty.cpp +++ b/src/tests/tribol_common_plane_penalty.cpp @@ -173,6 +173,26 @@ void checkForceSense( tribol::CouplingScheme const* cs, bool isTied = false ) } } // end checkForceSense() +void checkZeroResponses( tribol::CouplingScheme const* cs, const RealT tol = 1.e-14 ) +{ + const auto mesh1 = const_cast( cs )->getMesh1().getView(); + const auto mesh2 = const_cast( cs )->getMesh2().getView(); + + for ( int i = 0; i < 2; ++i ) { + auto& mesh = ( i == 0 ) ? mesh1 : mesh2; + for ( tribol::IndexT kf = 0; kf < mesh.numberOfElements(); ++kf ) { + for ( tribol::IndexT a = 0; a < mesh.numberOfNodesPerElement(); ++a ) { + int node_id = mesh.getGlobalNodeId( kf, a ); + EXPECT_NEAR( mesh.getResponse()[0][node_id], 0., tol ); + EXPECT_NEAR( mesh.getResponse()[1][node_id], 0., tol ); + if ( mesh.spatialDimension() == 3 ) { + EXPECT_NEAR( mesh.getResponse()[2][node_id], 0., tol ); + } + } + } + } +} // end checkZeroResponses() + /*! * Test fixture class with some setup necessary to test * the COMMON_PLANE + PENALTY implementation @@ -449,6 +469,76 @@ TEST_F( CommonPlaneTest, element_penalty_check ) tribol::finalize(); } +TEST_F( CommonPlaneTest, element_penalty_nonpositive_thickness_returns_error_without_forces ) +{ + this->m_mesh.mortarMeshId = 0; + this->m_mesh.nonmortarMeshId = 1; + + int nMortarElems = 4; + int nElemsXM = nMortarElems; + int nElemsYM = nMortarElems; + int nElemsZM = nMortarElems; + + int nNonmortarElems = 5; + int nElemsXS = nNonmortarElems; + int nElemsYS = nNonmortarElems; + int nElemsZS = nNonmortarElems; + + RealT x_min1 = 0.; + RealT y_min1 = 0.; + RealT z_min1 = 0.; + RealT x_max1 = 1.; + RealT y_max1 = 1.; + RealT z_max1 = 1.05; + + RealT x_min2 = 0.; + RealT y_min2 = 0.; + RealT z_min2 = 0.95; + RealT x_max2 = 1.; + RealT y_max2 = 1.; + RealT z_max2 = 2.; + + RealT element_thickness2 = ( z_max2 - z_min2 ) / nElemsZS; + + this->m_mesh.setupContactMeshHex( nElemsXM, nElemsYM, nElemsZM, x_min1, y_min1, z_min1, x_max1, y_max1, z_max1, + nElemsXS, nElemsYS, nElemsZS, x_min2, y_min2, z_min2, x_max2, y_max2, z_max2, 0., + 0. ); + + RealT bulk_mod1 = 1.0; + RealT bulk_mod2 = 1.0; + RealT velX1 = 0.; + RealT velY1 = 0.; + RealT velZ1 = 0.; + RealT velX2 = 0.; + RealT velY2 = 0.; + RealT velZ2 = 0.; + + this->m_mesh.allocateAndSetVelocities( m_mesh.mortarMeshId, velX1, velY1, velZ1 ); + this->m_mesh.allocateAndSetVelocities( m_mesh.nonmortarMeshId, velX2, velY2, -velZ2 ); + + this->m_mesh.allocateAndSetElementThickness( m_mesh.mortarMeshId, 0. ); + this->m_mesh.allocateAndSetBulkModulus( m_mesh.mortarMeshId, bulk_mod1 ); + this->m_mesh.allocateAndSetElementThickness( m_mesh.nonmortarMeshId, element_thickness2 ); + this->m_mesh.allocateAndSetBulkModulus( m_mesh.nonmortarMeshId, bulk_mod2 ); + + tribol::TestControlParameters parameters; + parameters.penalty_ratio = true; + parameters.const_penalty = 0.75; + parameters.dt = 1.e-3; + + int test_mesh_update_err = this->m_mesh.tribolSetupAndUpdate( + tribol::COMMON_PLANE, tribol::PENALTY, tribol::FRICTIONLESS, tribol::NO_CASE, false, parameters ); + + EXPECT_EQ( test_mesh_update_err, 0 ); + + tribol::CouplingSchemeManager& couplingSchemeManager = tribol::CouplingSchemeManager::getInstance(); + tribol::CouplingScheme* couplingScheme = &couplingSchemeManager.at( 0 ); + EXPECT_FALSE( couplingScheme->init() ); + checkZeroResponses( couplingScheme ); + + tribol::finalize(); +} + TEST_F( CommonPlaneTest, tied_contact_check ) { this->m_mesh.mortarMeshId = 0; diff --git a/src/tests/tribol_timestep_vote.cpp b/src/tests/tribol_timestep_vote.cpp index 4e2b0f82..0721ee61 100644 --- a/src/tests/tribol_timestep_vote.cpp +++ b/src/tests/tribol_timestep_vote.cpp @@ -305,7 +305,7 @@ TEST_F( CommonPlaneTest, numerically_zero_velocity_small_gap ) EXPECT_EQ( test_mesh_update_err, 0 ); - // expect no change in dt + // expect no change in dt because velocity is numerically zero and stability limits are disabled by default EXPECT_EQ( parameters.dt, dt ); tribol::finalize(); @@ -390,10 +390,7 @@ TEST_F( CommonPlaneTest, zero_velocity_large_gap ) tribol::COMMON_PLANE, tribol::PENALTY, tribol::FRICTIONLESS, tribol::NO_CASE, false, parameters ); EXPECT_EQ( test_mesh_update_err, 0 ); - // note that with very small velocity, the dt estimate will be - // very large, and won't change the timestep, even if the gap is large. This - // allows for a soft contact response in the presence of a small velocity that - // won't actually correct too much interpen with a contact dt vote + // expect no change in dt because velocity is zero and stability limits are disabled by default EXPECT_EQ( parameters.dt, dt ); tribol::finalize(); @@ -652,7 +649,7 @@ TEST_F( CommonPlaneTest, separation_velocity_small_gap ) EXPECT_EQ( test_mesh_update_err, 0 ); - // no change in dt because of separation velocities + // expect no change in dt because velocities are separating and stability limits are disabled by default EXPECT_EQ( parameters.dt, dt ); tribol::finalize(); @@ -917,6 +914,100 @@ TEST_F( CommonPlaneTest, large_velocity_small_separation_set_alpha ) tribol::finalize(); } +TEST_F( CommonPlaneTest, critical_timestep_cfl_limit ) +{ + // This test isolates and validates the Courant (CFL) stability limit. + // By setting a small timestep_scale (alpha = 0.1), we ensure that the CFL limit + // becomes the governing minimum timestep limit. + this->m_mesh.mortarMeshId = 0; + this->m_mesh.nonmortarMeshId = 1; + + constexpr int nMortarElems = 4; + constexpr int nElemsXM = nMortarElems; + constexpr int nElemsYM = nMortarElems; + constexpr int nElemsZM = nMortarElems; + + constexpr int nNonmortarElems = 5; + constexpr int nElemsXS = nNonmortarElems; + constexpr int nElemsYS = nNonmortarElems; + constexpr int nElemsZS = nNonmortarElems; + + constexpr RealT x_min1 = 0.; + constexpr RealT y_min1 = 0.; + constexpr RealT z_min1 = 0.; + constexpr RealT x_max1 = 1.; + constexpr RealT y_max1 = 1.; + constexpr RealT z_max1 = 1.005; + + constexpr RealT x_min2 = 0.; + constexpr RealT y_min2 = 0.; + constexpr RealT z_min2 = 0.95; + constexpr RealT x_max2 = 1.; + constexpr RealT y_max2 = 1.; + constexpr RealT z_max2 = 2.; + + constexpr RealT element_thickness1 = ( z_max1 - z_min1 ) / nElemsZM; + constexpr RealT element_thickness2 = ( z_max2 - z_min2 ) / nElemsZS; + + this->m_mesh.setupContactMeshHex( nElemsXM, nElemsYM, nElemsZM, x_min1, y_min1, z_min1, x_max1, y_max1, z_max1, + nElemsXS, nElemsYS, nElemsZS, x_min2, y_min2, z_min2, x_max2, y_max2, z_max2, 0., + 0. ); + + // Stationary contact to bypass impact limit + constexpr RealT dt = 1.0; + RealT bulk_mod1 = 1.0; + RealT bulk_mod2 = 1.0; + RealT velX1 = 0.; + RealT velY1 = 0.; + RealT velZ1 = 0.; + RealT velX2 = 0.; + RealT velY2 = 0.; + RealT velZ2 = 0.; + + this->m_mesh.allocateAndSetVelocities( m_mesh.mortarMeshId, velX1, velY1, velZ1 ); + this->m_mesh.allocateAndSetVelocities( m_mesh.nonmortarMeshId, velX2, velY2, -velZ2 ); + + this->m_mesh.allocateAndSetElementThickness( m_mesh.mortarMeshId, element_thickness1 ); + this->m_mesh.allocateAndSetBulkModulus( m_mesh.mortarMeshId, bulk_mod1 ); + this->m_mesh.allocateAndSetElementThickness( m_mesh.nonmortarMeshId, element_thickness2 ); + this->m_mesh.allocateAndSetBulkModulus( m_mesh.nonmortarMeshId, bulk_mod2 ); + + tribol::TestControlParameters parameters; + parameters.penalty_ratio = true; + parameters.const_penalty = 0.75; + parameters.dt = dt; + parameters.enable_timestep_vote = true; + parameters.enable_timestep_stability_limits = true; + parameters.timestep_pen_frac = 0.3; + parameters.timestep_scale = 0.1; // alpha = 0.1 + + int test_mesh_update_err = this->m_mesh.tribolSetupAndUpdate( + tribol::COMMON_PLANE, tribol::PENALTY, tribol::FRICTIONLESS, tribol::NO_CASE, false, parameters ); + + EXPECT_EQ( test_mesh_update_err, 0 ); + + // Expected CFL value: + // Stiffness under KINEMATIC_ELEMENT is computed as: bulk_modulus / thickness + // for each contacting face, where bulk_modulus is 1.0. + // The contact interface spring is modeled as mesh 1 and mesh 2's springs in series: + // penalty_stiff_per_area = (stiffness1 * stiffness2) / (stiffness1 + stiffness2) + constexpr RealT stiffness1 = 1.0 / element_thickness1; + constexpr RealT stiffness2 = 1.0 / element_thickness2; + constexpr RealT penalty_stiff_per_area = ( stiffness1 * stiffness2 ) / ( stiffness1 + stiffness2 ); + + // CFL scaling factors based on added interface stiffness: + // f_scale = 1 / sqrt( 1 + penalty_stiff_per_area / stiffness ) + const RealT f_scale1 = 1.0 / std::sqrt( 1.0 + ( penalty_stiff_per_area / stiffness1 ) ); + const RealT f_scale2 = 1.0 / std::sqrt( 1.0 + ( penalty_stiff_per_area / stiffness2 ) ); + + // dt_CFL = alpha * dt * min( f_scale1, f_scale2 ) + const RealT expected_dt = parameters.timestep_scale * dt * std::min( f_scale1, f_scale2 ); + RealT dt_tol = 1.e-6; + EXPECT_NEAR( parameters.dt, expected_dt, dt_tol ); + + tribol::finalize(); +} + int main( int argc, char* argv[] ) { int result = 0; diff --git a/src/tribol/common/Parameters.hpp b/src/tribol/common/Parameters.hpp index 6bb52366..aa661f31 100644 --- a/src/tribol/common/Parameters.hpp +++ b/src/tribol/common/Parameters.hpp @@ -483,6 +483,8 @@ struct Parameters { int vis_cycle_incr = 100; ///! Frequency for visualizations dumps VisType vis_type = VIS_NONE; ///! Type of interface physics visualization output bool enable_timestep_vote = false; ///! True if host-code desires the timestep vote to be calculated and returned + bool enable_timestep_stability_limits = + false; ///! True if host-code desires dynamic contact stability limits to be included in timestep vote bool auto_interpen_check = false; ///! True if the auto-contact interpenetration check is used for full-overlap pairs diff --git a/src/tribol/interface/tribol.cpp b/src/tribol/interface/tribol.cpp index 895381ab..46f39336 100644 --- a/src/tribol/interface/tribol.cpp +++ b/src/tribol/interface/tribol.cpp @@ -350,6 +350,19 @@ void enableTimestepVote( IndexT cs_id, const bool enable ) } // end enableTimestepVote() +//------------------------------------------------------------------------------ +void enableTimestepStabilityLimits( IndexT cs_id, const bool enable ) +{ + auto cs = CouplingSchemeManager::getInstance().findData( cs_id ); + + // check to see if coupling scheme exists + SLIC_ERROR_ROOT_IF( !cs, "tribol::enableTimestepStabilityLimits(): call tribol::registerCouplingScheme() " + << "prior to calling this routine." ); + + cs->getParameters().enable_timestep_stability_limits = enable; + +} // end enableTimestepStabilityLimits() + //------------------------------------------------------------------------------ void enableEnzyme( IndexT cs_id, [[maybe_unused]] bool use_enzyme ) { diff --git a/src/tribol/interface/tribol.hpp b/src/tribol/interface/tribol.hpp index 36e6fa03..9a7e55ff 100644 --- a/src/tribol/interface/tribol.hpp +++ b/src/tribol/interface/tribol.hpp @@ -226,6 +226,17 @@ void setBinningProximityScale( IndexT cs_id, RealT binning_proximity_scale ); */ void enableTimestepVote( IndexT cs_id, const bool enable ); +/*! + * \brief Enable CFL-like contact stiffness limits in the timestep vote + * + * \param [in] cs_id coupling scheme id + * \param [in] enable the contact stiffness limit will be included in the timestep vote if true + * + * \note default behavior is to not enable the contact stiffness limit + * + */ +void enableTimestepStabilityLimits( IndexT cs_id, const bool enable ); + /** * @brief Enable Enzyme AD for Jacobian calculations * diff --git a/src/tribol/mesh/CouplingScheme.cpp b/src/tribol/mesh/CouplingScheme.cpp index 99fca093..235729b0 100644 --- a/src/tribol/mesh/CouplingScheme.cpp +++ b/src/tribol/mesh/CouplingScheme.cpp @@ -27,6 +27,7 @@ #include "tribol/common/Parameters.hpp" #include "tribol/physics/Physics.hpp" #include "tribol/physics/ContactFormulationFactory.hpp" +#include "tribol/physics/CommonPlane.hpp" #include "tribol/integ/FE.hpp" namespace tribol { @@ -1349,365 +1350,377 @@ void CouplingScheme::computeCommonPlaneTimeStep( RealT& dt ) // all candidates, not necessarily ones that are deemed to be in contact per // the gap constraint. auto cs_view = getView(); - ArrayT dt_temp_data( { dt, dt }, getAllocatorId() ); + ArrayT dt_temp_data( { dt, dt, dt }, getAllocatorId() ); ArrayViewT dt_temp = dt_temp_data; // [0]: exceed_max_gap1, [1]: exceed_max_gap2, [2]: neg_dt_gap_msg, [3]: neg_dt_vel_proj_msg ArrayT msg_data( { static_cast( false ), static_cast( false ), static_cast( false ), static_cast( false ) }, getAllocatorId() ); ArrayViewT msg = msg_data; - forAllExec( getExecutionMode(), getNumActivePairs(), - [cs_view, dim, proj_ratio, msg, dt_temp, dt] TRIBOL_HOST_DEVICE( IndexT i ) { - auto& cg_view = cs_view.getCompGeomView(); - auto& plane = cg_view.getCommonPlane( i ); - - auto& mesh1 = cs_view.getMesh1View(); - auto& mesh2 = cs_view.getMesh2View(); - - // get pair indices - IndexT index1 = plane.getCpElementId1(); - IndexT index2 = plane.getCpElementId2(); - - constexpr int max_dim = 3; - constexpr int max_nodes_per_elem = 4; - StackArrayT x1; - StackArrayT v1; - mesh1.getFaceCoords( index1, x1 ); - mesh1.getFaceVelocities( index1, v1 ); - - StackArrayT x2; - StackArrayT v2; - mesh2.getFaceCoords( index2, x2 ); - mesh2.getFaceVelocities( index2, v2 ); - - ///////////////////////////////////////////////////////////// - // calculate face velocities at projected overlap centroid // - ///////////////////////////////////////////////////////////// - StackArrayT vel_f1; - StackArrayT vel_f2; - initRealArray( vel_f1, dim, 0.0 ); - initRealArray( vel_f2, dim, 0.0 ); - - // interpolate nodal velocity at overlap centroid as projected - // onto face 1 - RealT cXf1 = plane.m_cXf1; - RealT cYf1 = plane.m_cYf1; - RealT cZf1 = ( dim == 3 ) ? plane.m_cZf1 : 0.; - GalerkinEval( x1, cXf1, cYf1, cZf1, LINEAR, PHYSICAL, dim, dim, v1, vel_f1 ); - // interpolate nodal velocity at overlap centroid as projected - // onto face 2 - RealT cXf2 = plane.m_cXf2; - RealT cYf2 = plane.m_cYf2; - RealT cZf2 = ( dim == 3 ) ? plane.m_cZf2 : 0.; - GalerkinEval( x2, cXf2, cYf2, cZf2, LINEAR, PHYSICAL, dim, dim, v2, vel_f2 ); - - //////////////////////////////////////////////// - // // - // Compute Timestep Vote Based on a Few Cases // - // // - //////////////////////////////////////////////// - - /////////////////////////////////////////////// - // compute data common to all timestep votes // - /////////////////////////////////////////////// - - // compute velocity projections: - // compute the dot product between the face velocities - // at the overlap-centroid-to-face projected centroid and each - // face's outward unit normal AND the overlap normal. - RealT v1_dot_n, v2_dot_n, v1_dot_n1, v2_dot_n2; - RealT overlapNormal[max_dim]; - overlapNormal[0] = plane.m_nX; - overlapNormal[1] = plane.m_nY; - if ( dim == 3 ) { - overlapNormal[2] = plane.m_nZ; - } + forAllExec( + getExecutionMode(), getNumActivePairs(), + [cs_view, dim, proj_ratio, msg, dt_temp, dt] TRIBOL_HOST_DEVICE( IndexT i ) { + auto& cg_view = cs_view.getCompGeomView(); + auto& plane = cg_view.getCommonPlane( i ); + + auto& mesh1 = cs_view.getMesh1View(); + auto& mesh2 = cs_view.getMesh2View(); + + // get pair indices + IndexT index1 = plane.getCpElementId1(); + IndexT index2 = plane.getCpElementId2(); + + constexpr int max_dim = 3; + constexpr int max_nodes_per_elem = 4; + StackArrayT x1; + StackArrayT v1; + mesh1.getFaceCoords( index1, x1 ); + mesh1.getFaceVelocities( index1, v1 ); + + StackArrayT x2; + StackArrayT v2; + mesh2.getFaceCoords( index2, x2 ); + mesh2.getFaceVelocities( index2, v2 ); + + ///////////////////////////////////////////////////////////// + // calculate face velocities at projected overlap centroid // + ///////////////////////////////////////////////////////////// + StackArrayT vel_f1; + StackArrayT vel_f2; + initRealArray( vel_f1, dim, 0.0 ); + initRealArray( vel_f2, dim, 0.0 ); + + // interpolate nodal velocity at overlap centroid as projected + // onto face 1 + RealT cXf1 = plane.m_cXf1; + RealT cYf1 = plane.m_cYf1; + RealT cZf1 = ( dim == 3 ) ? plane.m_cZf1 : 0.; + GalerkinEval( x1, cXf1, cYf1, cZf1, LINEAR, PHYSICAL, dim, dim, v1, vel_f1 ); + // interpolate nodal velocity at overlap centroid as projected + // onto face 2 + RealT cXf2 = plane.m_cXf2; + RealT cYf2 = plane.m_cYf2; + RealT cZf2 = ( dim == 3 ) ? plane.m_cZf2 : 0.; + GalerkinEval( x2, cXf2, cYf2, cZf2, LINEAR, PHYSICAL, dim, dim, v2, vel_f2 ); + + //////////////////////////////////////////////// + // // + // Compute Timestep Vote Based on a Few Cases // + // // + //////////////////////////////////////////////// + + /////////////////////////////////////////////// + // compute data common to all timestep votes // + /////////////////////////////////////////////// + + // compute velocity projections: + // compute the dot product between the face velocities + // at the overlap-centroid-to-face projected centroid and each + // face's outward unit normal AND the overlap normal. + RealT v1_dot_n, v2_dot_n, v1_dot_n1, v2_dot_n2; + RealT overlapNormal[max_dim]; + overlapNormal[0] = plane.m_nX; + overlapNormal[1] = plane.m_nY; + if ( dim == 3 ) { + overlapNormal[2] = plane.m_nZ; + } - // get face normals - RealT fn1[max_dim], fn2[max_dim]; - mesh1.getFaceNormal( index1, fn1 ); - mesh2.getFaceNormal( index2, fn2 ); - - // compute projections - v1_dot_n = dotProd( vel_f1, overlapNormal, dim ); - v2_dot_n = dotProd( vel_f2, overlapNormal, dim ); - v1_dot_n1 = dotProd( vel_f1, fn1, dim ); - v2_dot_n2 = dotProd( vel_f2, fn2, dim ); - - // Keep debug print statements. This routine is still in the testing phase - // std::cout << "face 1 normal: " << fn1[0] << ", " << fn1[1] << ", " << fn1[2] << std::endl; - // std::cout << "face 2 normal: " << fn2[0] << ", " << fn2[1] << ", " << fn2[2] << std::endl; - // std::cout << " " << std::endl; - // std::cout << "face 1 vel: " << vel_f1[0] << ", " << vel_f1[1] << ", " << vel_f1[2] << std::endl; - // std::cout << "face 2 vel: " << vel_f2[0] << ", " << vel_f2[1] << ", " << vel_f2[2] << std::endl; - // std::cout << " " << std::endl; - // std::cout << "First v1_dot_n1 calc: " << v1_dot_n1 << std::endl; - // std::cout << "First v2_dot_n2 calc: " << v2_dot_n2 << std::endl; - // std::cout << "First v1_dot_n: " << v1_dot_n << std::endl; - // std::cout << "First v2_dot_n: " << v2_dot_n << std::endl; - - // add tiny amount to velocity projections to avoid division by zero. - // Note that if these projections are close to zero, there may be - // stationary interactions or tangential motion. In this case, any - // timestep estimate will be very large, and not control the simulation - RealT tiny = 1.e-12; - RealT tiny1 = ( v1_dot_n >= 0. ) ? tiny : -1. * tiny; - RealT tiny2 = ( v2_dot_n >= 0. ) ? tiny : -1. * tiny; - v1_dot_n += tiny1; - v2_dot_n += tiny2; - // reset tiny velocity based on face normal projections. - tiny1 = ( v1_dot_n1 >= 0. ) ? tiny : -1. * tiny; - tiny2 = ( v2_dot_n2 >= 0. ) ? tiny : -1. * tiny; - v1_dot_n1 += tiny1; - v2_dot_n2 += tiny2; - - // Keep debug print statements. This routine is still in the testing phase - // std::cout << "Second v1_dot_n1 calc: " << v1_dot_n1 << std::endl; - // std::cout << "Second v2_dot_n2 calc: " << v2_dot_n2 << std::endl; - // std::cout << "Second v1_dot_n: " << v1_dot_n << std::endl; - // std::cout << "Second v2_dot_n: " << v2_dot_n << std::endl; - - // get volume element thicknesses associated with each face in this pair - RealT t1 = mesh1.getElementData().m_thickness[index1]; - RealT t2 = mesh2.getElementData().m_thickness[index2]; - - // compute the existing gap vector (recall gap is x1-x2 by convention) - RealT gapVec[max_dim]; - gapVec[0] = plane.m_cXf1 - plane.m_cXf2; - gapVec[1] = plane.m_cYf1 - plane.m_cYf2; - if ( dim == 3 ) { - gapVec[2] = plane.m_cZf1 - plane.m_cZf2; - } + // get face normals + RealT fn1[max_dim], fn2[max_dim]; + mesh1.getFaceNormal( index1, fn1 ); + mesh2.getFaceNormal( index2, fn2 ); + + // compute projections + v1_dot_n = dotProd( vel_f1, overlapNormal, dim ); + v2_dot_n = dotProd( vel_f2, overlapNormal, dim ); + v1_dot_n1 = dotProd( vel_f1, fn1, dim ); + v2_dot_n2 = dotProd( vel_f2, fn2, dim ); + + // Velocity tolerance for divide-by-zero threshold guarding + constexpr RealT tiny = 1.e-12; + // Large fallback value used when the time-to-threshold is effectively infinite. + constexpr RealT fallback_dt = 1.e6; + + // get volume element thicknesses associated with each face in this pair + RealT t1 = mesh1.getElementData().m_thickness[index1]; + RealT t2 = mesh2.getElementData().m_thickness[index2]; + + // compute the existing gap vector (recall gap is x1-x2 by convention) + RealT gapVec[max_dim]; + gapVec[0] = plane.m_cXf1 - plane.m_cXf2; + gapVec[1] = plane.m_cYf1 - plane.m_cYf2; + if ( dim == 3 ) { + gapVec[2] = plane.m_cZf1 - plane.m_cZf2; + } - // compute the dot product between gap vector and the outward unit face normals. - RealT gap_f1_n1 = dotProd( gapVec, fn1, dim ); - RealT gap_f2_n2 = dotProd( gapVec, fn2, dim ); - - RealT dt1 = 1.e6; // initialize as large number - RealT dt2 = 1.e6; // initialize as large number - RealT alpha = cs_view.getTimestepScale(); // multiplier on timestep estimate - bool dt1_check1 = false; - bool dt2_check1 = false; - bool dt1_vel_check = false; - bool dt2_vel_check = false; - - // maximum allowable interpenetration in the normal direction of each element - RealT max_delta1 = proj_ratio * t1; - RealT max_delta2 = proj_ratio * t2; - - // Separation or interpenetration trigger for check 1 and 2: - // check if there is further interpen or separation based on the - // velocity projection in the direction of the common-plane normal, - // which is in the direction of face-2 normal. - // The two cases are: - // if v1*n < 0 there is interpen - // if v2*n > 0 there is interpen - // - // Note: we compare strictly to 0. here since a 'tiny' value was - // appropriately added to the velocity projections, which is akin - // to some tolerancing effect - dt1_vel_check = ( v1_dot_n < 0. ) ? true : false; - dt2_vel_check = ( v2_dot_n > 0. ) ? true : false; - - ////////////////////////////////////////////////////////////////////////// - // Check 1. Current interpenetration gap exceeds max allowable interpen // - ////////////////////////////////////////////////////////////////////////// - - // check if face-pair is in contact (i.e. gap < gap_tol), which is determined - // in Common Plane ApplyNormal<>() routine - if ( plane.m_inContact ) { - // compute the difference between the 'face-gaps' and the max allowable - // interpen as a function of element thickness. Note, we have to use the - // gap projected onto the outward unit face-normal to check against the - // max allowable gap as a factor of the thickness in the element normal - // direction - RealT delta1 = max_delta1 - gap_f1_n1; // >0 not exceeding max allowable - RealT delta2 = max_delta2 + gap_f2_n2; // >0 not exceeding max allowable - - auto exceed_max_gap1 = ( delta1 < 0. ) ? true : false; - auto exceed_max_gap2 = ( delta2 < 0. ) ? true : false; - - // if velocity projection indicates further interpenetration, and the gaps - // EXCEED max allowable, then compute time step estimates to reduce overlap - dt1_check1 = ( dt1_vel_check ) ? exceed_max_gap1 : false; - dt2_check1 = ( dt2_vel_check ) ? exceed_max_gap2 : false; + // compute the dot product between gap vector and the outward unit face normals. + RealT gap_f1_n1 = dotProd( gapVec, fn1, dim ); + RealT gap_f2_n2 = dotProd( gapVec, fn2, dim ); + + RealT dt1 = fallback_dt; + RealT dt2 = fallback_dt; + RealT alpha = cs_view.getTimestepScale(); // multiplier on timestep estimate + bool dt1_check1 = false; + bool dt2_check1 = false; + bool dt1_vel_check = false; + bool dt2_vel_check = false; + + // maximum allowable interpenetration in the normal direction of each element + RealT max_delta1 = proj_ratio * t1; + RealT max_delta2 = proj_ratio * t2; + + // Separation or interpenetration trigger for check 1 and 2: + // check if there is further interpen or separation based on the + // velocity projection in the direction of the common-plane normal, + // which is in the direction of face-2 normal. + // The two cases are: + // if v1*n < 0 there is interpen + // if v2*n > 0 there is interpen + dt1_vel_check = ( v1_dot_n < 0. ) ? true : false; + dt2_vel_check = ( v2_dot_n > 0. ) ? true : false; + + ////////////////////////////////////////////////////////////////////////// + // Check 1. Current interpenetration gap exceeds max allowable interpen // + ////////////////////////////////////////////////////////////////////////// + + // check if face-pair is in contact (i.e. gap < gap_tol), which is determined + // in Common Plane ApplyNormal<>() routine + if ( plane.m_inContact ) { + // compute the difference between the 'face-gaps' and the max allowable + // interpen as a function of element thickness. Note, we have to use the + // gap projected onto the outward unit face-normal to check against the + // max allowable gap as a factor of the thickness in the element normal + // direction + RealT delta1 = max_delta1 - gap_f1_n1; // >0 not exceeding max allowable + RealT delta2 = max_delta2 + gap_f2_n2; // >0 not exceeding max allowable + + auto exceed_max_gap1 = ( delta1 < 0. ) ? true : false; + auto exceed_max_gap2 = ( delta2 < 0. ) ? true : false; + + // if velocity projection indicates further interpenetration, and the gaps + // EXCEED max allowable, then compute time step estimates to reduce overlap + dt1_check1 = ( dt1_vel_check ) ? exceed_max_gap1 : false; + dt2_check1 = ( dt2_vel_check ) ? exceed_max_gap2 : false; #ifdef TRIBOL_USE_RAJA - RAJA::atomicMax( &msg[0], static_cast( exceed_max_gap1 ) ); - RAJA::atomicMax( &msg[1], static_cast( exceed_max_gap2 ) ); + RAJA::atomicMax( &msg[0], static_cast( exceed_max_gap1 ) ); + RAJA::atomicMax( &msg[1], static_cast( exceed_max_gap2 ) ); #else - msg[0] = exceed_max_gap1; - msg[1] = exceed_max_gap2; + msg[0] = exceed_max_gap1; + msg[1] = exceed_max_gap2; #endif - // compute dt for face 1 and 2 based on the velocity and gap projections onto - // the face-normals for faces where currect gap exceeds max allowable gap. - // - // NOTE: - // - // This calculation RESETS the current gap to be g = 0, and computes a timestep - // such that the velocity projection of the overlap-to-face projected overlap - // centroid does not exceed the max allowable gap. - // - // This avoid a timestep crash in the case that the current gap barely exceeds - // the max allowable and also allows a soft contact response with interpen - // in excess of the max allowable gap without causing timestep crashes. - // - // v1_dot_n1 > 0 and v2_dot_n2 > 0 for further interpen - dt1 = ( dt1_check1 ) ? alpha * max_delta1 / v1_dot_n1 : dt1; - dt2 = ( dt2_check1 ) ? alpha * max_delta2 / v2_dot_n2 : dt2; - - // Keep debug print statements. This routine is still in the testing phase - // std::cout << "dt1_check1, delta1 and v1_dot_n1: " << dt1_check1 << ", " << max_delta1 << ", " << - // v1_dot_n1 << std::endl; std::cout << "dt2_check1, delta2 and v2_dot_n2: " << dt2_check1 << ", " << - // max_delta2 << ", " << v2_dot_n2 << std::endl; std::cout << "dt1 and dt2: " << dt1 << ", " << dt2 << - // std::endl; - - // update dt_temp1 only for positive dt1 and/or dt2 - if ( dt1 > 0. ) { + // compute dt for face 1 and 2 based on the velocity and gap projections onto + // the face-normals for faces where currect gap exceeds max allowable gap. + // + // NOTE: + // + // This calculation RESETS the current gap to be g = 0, and computes a timestep + // such that the velocity projection of the overlap-to-face projected overlap + // centroid does not exceed the max allowable gap. + // + // This avoid a timestep crash in the case that the current gap barely exceeds + // the max allowable and also allows a soft contact response with interpen + // in excess of the max allowable gap without causing timestep crashes. + // + // v1_dot_n1 > 0 and v2_dot_n2 > 0 for further interpen + if ( dt1_check1 ) { + // If velocity is effectively zero, the time to exceed gap is effectively infinite. + dt1 = ( std::abs( v1_dot_n1 ) > tiny ) ? ( alpha * max_delta1 / v1_dot_n1 ) : fallback_dt; + } + if ( dt2_check1 ) { + // If velocity is effectively zero, the time to exceed gap is effectively infinite. + dt2 = ( std::abs( v2_dot_n2 ) > tiny ) ? ( alpha * max_delta2 / v2_dot_n2 ) : fallback_dt; + } + + // update dt_temp1 only for positive dt1 and/or dt2 + if ( dt1 > 0. ) { #ifdef TRIBOL_USE_RAJA - RAJA::atomicMin( &dt_temp[0], axom::utilities::min( dt1, 1.e6 ) ); + RAJA::atomicMin( &dt_temp[0], axom::utilities::min( dt1, fallback_dt ) ); #else - dt_temp[0] = axom::utilities::min(dt_temp[0], axom::utilities::min(dt1, 1.e6)); + dt_temp[0] = axom::utilities::min( dt_temp[0], axom::utilities::min( dt1, fallback_dt ) ); #endif - } - if ( dt2 > 0. ) { + } + if ( dt2 > 0. ) { #ifdef TRIBOL_USE_RAJA - RAJA::atomicMin( &dt_temp[0], axom::utilities::min( 1.e6, dt2 ) ); + RAJA::atomicMin( &dt_temp[0], axom::utilities::min( fallback_dt, dt2 ) ); #else - dt_temp[0] = axom::utilities::min(dt_temp[0], axom::utilities::min(1.e6, dt2)); + dt_temp[0] = axom::utilities::min( dt_temp[0], axom::utilities::min( fallback_dt, dt2 ) ); #endif - } + } - if ( dt1 < 0. || dt2 < 0. ) { + if ( dt1 < 0. || dt2 < 0. ) { #ifdef TRIBOL_USE_RAJA - RAJA::atomicMax( &msg[2], static_cast( true ) ); + RAJA::atomicMax( &msg[2], static_cast( true ) ); #else - msg[2] = true; + msg[2] = true; #endif - } - } // end case 1 - - //////////////////////////////////////////////////////////////////////// - // 2. Velocity projection exceeds max interpenetration // - // // - // Note: This is performed for all contact candidates even if they // - // are not 'in contact' per the common-plane method. Every // - // contact candidate has a contact plane // - //////////////////////////////////////////////////////////////////////// - - { - // compute the delta between the velocity projection of each face-projected - // common plane centroid location - // - // First project each face-projected common plane centroid using linear velocity - // projection as approximation of configuration next cycle - RealT proj_delta_x1 = plane.m_cXf1 + dt * vel_f1[0]; - RealT proj_delta_y1 = plane.m_cYf1 + dt * vel_f1[1]; - RealT proj_delta_z1 = 0.; - - RealT proj_delta_x2 = plane.m_cXf2 + dt * vel_f2[0]; - RealT proj_delta_y2 = plane.m_cYf2 + dt * vel_f2[1]; - RealT proj_delta_z2 = 0.; - - // Second compute the amount of interpenetration as the difference between the two - // velocity projected points - RealT proj_delta_x1_fixed = proj_delta_x1; - RealT proj_delta_y1_fixed = proj_delta_y1; - - proj_delta_x1 -= proj_delta_x2; - proj_delta_y1 -= proj_delta_y2; - - proj_delta_x2 -= proj_delta_x1_fixed; - proj_delta_y2 -= proj_delta_y1_fixed; - - // compute the dot product between each face's delta and the OTHER - // face's outward unit normal. This is the magnitude of interpenetration - // of one face's projected overlap-centroid in the 'thickness-direction' - // of the other face (with whom in may be in contact currently, or in - // a velocity projected sense). - RealT proj_delta_n_1 = proj_delta_x1 * fn2[0] + proj_delta_y1 * fn2[1]; - RealT proj_delta_n_2 = proj_delta_x2 * fn1[0] + proj_delta_y2 * fn1[1]; - - if ( dim == 3 ) { - // project the z-component - proj_delta_z1 = plane.m_cZf1 + dt * vel_f1[2]; - proj_delta_z2 = plane.m_cZf2 + dt * vel_f2[2]; - - RealT proj_delta_z1_fixed = proj_delta_z1; - - // compute difference between each projected point - proj_delta_z1 -= proj_delta_z2; - proj_delta_z2 -= proj_delta_z1_fixed; - - // add the z-component of the projection onto the face normal - proj_delta_n_1 += proj_delta_z1 * fn2[2]; - proj_delta_n_2 += proj_delta_z2 * fn1[2]; - } - - // Reset the dt velocity check only for faces with continued interpen that exceeds the - // max allowable gap AND where the current gap did NOT exceed that face's max allowable - // gap per check 1 (would result in same dt calc). - // - // Note: - // If proj_delta_n_i < 0, (i=1,2) there is interpen from the velocity projection. - // Check this interpen against the maximum allowable to determine if a velocity projection - // timestep estimate is still required. - if ( dt1_vel_check && !dt1_check1 ) // continued interpen - { - dt1_vel_check = ( proj_delta_n_1 < 0. ) - ? ( ( std::abs( proj_delta_n_1 ) > max_delta1 ) ? true : false ) - : false; - } - - if ( dt2_vel_check && !dt2_check1 ) // continued interpen - { - dt2_vel_check = ( proj_delta_n_2 < 0. ) - ? ( ( std::abs( proj_delta_n_2 ) > max_delta2 ) ? true : false ) - : false; - } - - // compute velocity projection based dt (check 2) using a RESET gap (g=0) such that - // the velocity projected gap does not exceed the max allowable gap. This avoid timestep - // crashes for velocity projected gaps slightly in excess of the max allowable and still - // allows for a soft contact response without a timestep crash. - // - // v1_dot_n1 > 0 and v2_dot_n2 > 0 for further interpen - dt1 = ( dt1_vel_check ) ? alpha * max_delta1 / v1_dot_n1 : dt1; - dt2 = ( dt2_vel_check ) ? alpha * max_delta2 / v2_dot_n2 : dt2; - - // Keep debug print statements. This routine is still in the testing phase - // std::cout << "dt1_vel_check, (proj_delta_n_1+max_delta1), v1_dot_n1: " << dt1_vel_check << ", " - // << proj_delta_n_1+max_delta1 << ", " << v1_dot_n1 << std::endl; - // std::cout << "dt2_vel_check, (proj_delta_n_2+max_delta2), v2_dot_n2: " << dt2_vel_check << ", " - // << proj_delta_n_2+max_delta2 << ", " << v2_dot_n2 << std::endl; - // std::cout << "dt1 and dt2: " << dt1 << ", " << dt2 << std::endl; - - // update dt_temp2 only for positive dt1 and/or dt2 - if ( dt1 > 0. ) { + } + } // end case 1 + + //////////////////////////////////////////////////////////////////////// + // 2. Velocity projection exceeds max interpenetration // + // // + // Note: This is performed for all contact candidates even if they // + // are not 'in contact' per the common-plane method. Every // + // contact candidate has a contact plane // + //////////////////////////////////////////////////////////////////////// + + { + // compute the delta between the velocity projection of each face-projected + // common plane centroid location + // + // First project each face-projected common plane centroid using linear velocity + // projection as approximation of configuration next cycle + RealT proj_delta_x1 = plane.m_cXf1 + dt * vel_f1[0]; + RealT proj_delta_y1 = plane.m_cYf1 + dt * vel_f1[1]; + RealT proj_delta_z1 = 0.; + + RealT proj_delta_x2 = plane.m_cXf2 + dt * vel_f2[0]; + RealT proj_delta_y2 = plane.m_cYf2 + dt * vel_f2[1]; + RealT proj_delta_z2 = 0.; + + // Second compute the amount of interpenetration as the difference between the two + // velocity projected points + RealT proj_delta_x1_fixed = proj_delta_x1; + RealT proj_delta_y1_fixed = proj_delta_y1; + + proj_delta_x1 -= proj_delta_x2; + proj_delta_y1 -= proj_delta_y2; + + proj_delta_x2 -= proj_delta_x1_fixed; + proj_delta_y2 -= proj_delta_y1_fixed; + + // compute the dot product between each face's delta and the OTHER + // face's outward unit normal. This is the magnitude of interpenetration + // of one face's projected overlap-centroid in the 'thickness-direction' + // of the other face (with whom in may be in contact currently, or in + // a velocity projected sense). + RealT proj_delta_n_1 = proj_delta_x1 * fn2[0] + proj_delta_y1 * fn2[1]; + RealT proj_delta_n_2 = proj_delta_x2 * fn1[0] + proj_delta_y2 * fn1[1]; + + if ( dim == 3 ) { + // project the z-component + proj_delta_z1 = plane.m_cZf1 + dt * vel_f1[2]; + proj_delta_z2 = plane.m_cZf2 + dt * vel_f2[2]; + + RealT proj_delta_z1_fixed = proj_delta_z1; + + // compute difference between each projected point + proj_delta_z1 -= proj_delta_z2; + proj_delta_z2 -= proj_delta_z1_fixed; + + // add the z-component of the projection onto the face normal + proj_delta_n_1 += proj_delta_z1 * fn2[2]; + proj_delta_n_2 += proj_delta_z2 * fn1[2]; + } + + // Reset the dt velocity check only for faces with continued interpen that exceeds the + // max allowable gap AND where the current gap did NOT exceed that face's max allowable + // gap per check 1 (would result in same dt calc). + // + // Note: + // If proj_delta_n_i < 0, (i=1,2) there is interpen from the velocity projection. + // Check this interpen against the maximum allowable to determine if a velocity projection + // timestep estimate is still required. + if ( dt1_vel_check && !dt1_check1 ) // continued interpen + { + dt1_vel_check = + ( proj_delta_n_1 < 0. ) ? ( ( std::abs( proj_delta_n_1 ) > max_delta1 ) ? true : false ) : false; + } + + if ( dt2_vel_check && !dt2_check1 ) // continued interpen + { + dt2_vel_check = + ( proj_delta_n_2 < 0. ) ? ( ( std::abs( proj_delta_n_2 ) > max_delta2 ) ? true : false ) : false; + } + + // compute velocity projection based dt (check 2) using a RESET gap (g=0) such that + // the velocity projected gap does not exceed the max allowable gap. This avoid timestep + // crashes for velocity projected gaps slightly in excess of the max allowable and still + // allows for a soft contact response without a timestep crash. + // + // v1_dot_n1 > 0 and v2_dot_n2 > 0 for further interpen + if ( dt1_vel_check ) { + // If velocity is effectively zero, the time to exceed gap is effectively infinite. + dt1 = ( std::abs( v1_dot_n1 ) > tiny ) ? ( alpha * max_delta1 / v1_dot_n1 ) : fallback_dt; + } + if ( dt2_vel_check ) { + // If velocity is effectively zero, the time to exceed gap is effectively infinite. + dt2 = ( std::abs( v2_dot_n2 ) > tiny ) ? ( alpha * max_delta2 / v2_dot_n2 ) : fallback_dt; + } + + // update dt_temp2 only for positive dt1 and/or dt2 + if ( dt1 > 0. ) { #ifdef TRIBOL_USE_RAJA - RAJA::atomicMin( &dt_temp[1], axom::utilities::min( dt1, 1.e6 ) ); + RAJA::atomicMin( &dt_temp[1], axom::utilities::min( dt1, fallback_dt ) ); #else - dt_temp[1] = axom::utilities::min(dt_temp[1], axom::utilities::min(dt1, 1.e6)); + dt_temp[1] = axom::utilities::min( dt_temp[1], axom::utilities::min( dt1, fallback_dt ) ); #endif - } - if ( dt2 > 0. ) { + } + if ( dt2 > 0. ) { #ifdef TRIBOL_USE_RAJA - RAJA::atomicMin( &dt_temp[1], axom::utilities::min( 1.e6, dt2 ) ); + RAJA::atomicMin( &dt_temp[1], axom::utilities::min( fallback_dt, dt2 ) ); #else - dt_temp[1] = axom::utilities::min(dt_temp[1], axom::utilities::min(1.e6, dt2)); + dt_temp[1] = axom::utilities::min( dt_temp[1], axom::utilities::min( fallback_dt, dt2 ) ); #endif - } - if ( dt1 < 0. || dt2 < 0. ) { + } + if ( dt1 < 0. || dt2 < 0. ) { #ifdef TRIBOL_USE_RAJA - RAJA::atomicMax( &msg[3], static_cast( true ) ); + RAJA::atomicMax( &msg[3], static_cast( true ) ); #else - msg[3] = true; + msg[3] = true; #endif - } + } - } // end check 2 - } ); + } // end check 2 + + //////////////////////////////////////////////////////////////////////// + // 3. Contact Dynamics and Stability Limits // + // // + // These limits constrain the explicit time integration based on // + // the added stiffness of active contact springs. // + //////////////////////////////////////////////////////////////////////// + + if ( plane.m_inContact && cs_view.enableTimestepStabilityLimits() ) { + RealT stiffness1 = 0.; + RealT stiffness2 = 0.; + auto& enforcement_options = cs_view.getEnforcementOptions(); + const PenaltyEnforcementOptions& pen_enfrc_options = enforcement_options.penalty_options; + RealT pen_scale1 = mesh1.getElementData().m_penalty_scale; + RealT pen_scale2 = mesh2.getElementData().m_penalty_scale; + RealT mat_mod1 = 0.; + RealT mat_mod2 = 0.; + + if ( pen_enfrc_options.kinematic_calculation == KINEMATIC_ELEMENT ) { + mat_mod1 = mesh1.getElementData().m_mat_mod[index1]; + mat_mod2 = mesh2.getElementData().m_mat_mod[index2]; + } + + ComputeUncoupledStiffness( pen_enfrc_options.kinematic_calculation, t1, pen_enfrc_options.tiny_length, + pen_scale1, mat_mod1, mesh1.getElementData().m_penalty_stiffness, stiffness1 ); + ComputeUncoupledStiffness( pen_enfrc_options.kinematic_calculation, t2, pen_enfrc_options.tiny_length, + pen_scale2, mat_mod2, mesh2.getElementData().m_penalty_stiffness, stiffness2 ); + + if ( stiffness1 > 0. && stiffness2 > 0. ) { + RealT penalty_stiff_per_area = ComputePenaltyStiffnessPerArea( stiffness1, stiffness2 ); + if ( penalty_stiff_per_area > 0. ) { + // Rescale the host timestep by the added interface stiffness. + RealT f_scale1 = 1.0 / std::sqrt( 1.0 + ( penalty_stiff_per_area / stiffness1 ) ); + RealT f_scale2 = 1.0 / std::sqrt( 1.0 + ( penalty_stiff_per_area / stiffness2 ) ); + RealT dt_crit = alpha * dt * axom::utilities::min( f_scale1, f_scale2 ); + + if ( dt_crit > 0. ) { +#ifdef TRIBOL_USE_RAJA + RAJA::atomicMin( &dt_temp[2], axom::utilities::min( dt_crit, fallback_dt ) ); +#else + dt_temp[2] = axom::utilities::min( dt_temp[2], axom::utilities::min( dt_crit, fallback_dt ) ); +#endif + } + } + } + } // end contact dynamics and stability limits + } ); // print general messages once // Can we output this message on root? SRW @@ -1726,7 +1739,7 @@ void CouplingScheme::computeCommonPlaneTimeStep( RealT& dt ) << "velocity projection calculation." ); ArrayT dt_temp_host( dt_temp_data ); - dt = axom::utilities::min( dt_temp_host[0], dt_temp_host[1] ); + dt = axom::utilities::min( dt_temp_host[0], axom::utilities::min( dt_temp_host[1], dt_temp_host[2] ) ); } //------------------------------------------------------------------------------ diff --git a/src/tribol/mesh/CouplingScheme.hpp b/src/tribol/mesh/CouplingScheme.hpp index df8ec888..bdbc8b11 100644 --- a/src/tribol/mesh/CouplingScheme.hpp +++ b/src/tribol/mesh/CouplingScheme.hpp @@ -187,6 +187,16 @@ class CouplingScheme { */ TRIBOL_HOST_DEVICE RealT getTimestepScale() const { return m_parameters.timestep_scale; } + /** + * @brief Return whether contact dynamics stability limits are enabled + * + * @return true if stability limits are enabled, false otherwise + */ + TRIBOL_HOST_DEVICE bool enableTimestepStabilityLimits() const + { + return m_parameters.enable_timestep_stability_limits; + } + /** * @brief Get the gap tolerance that determines in contact face-pairs * diff --git a/src/tribol/physics/CommonPlane.cpp b/src/tribol/physics/CommonPlane.cpp index 39eae36a..d976f04a 100644 --- a/src/tribol/physics/CommonPlane.cpp +++ b/src/tribol/physics/CommonPlane.cpp @@ -167,37 +167,31 @@ int ApplyNormal( CouplingScheme* cs ) const PenaltyEnforcementOptions& pen_enfrc_options = enforcement_options.penalty_options; RealT pen_scale1 = mesh1.getElementData().m_penalty_scale; RealT pen_scale2 = mesh2.getElementData().m_penalty_scale; - switch ( pen_enfrc_options.kinematic_calculation ) { - case KINEMATIC_CONSTANT: { - // pre-multiply each spring stiffness by each mesh's penalty scale - auto stiffness1 = pen_scale1 * mesh1.getElementData().m_penalty_stiffness; - auto stiffness2 = pen_scale2 * mesh2.getElementData().m_penalty_stiffness; - // compute the equivalent contact penalty spring stiffness per area - penalty_stiff_per_area = ComputePenaltyStiffnessPerArea( stiffness1, stiffness2 ); - break; - } - case KINEMATIC_ELEMENT: { - // add tiny_length to element thickness to avoid division by zero - auto t1 = mesh1.getElementData().m_thickness[index1] + pen_enfrc_options.tiny_length; - auto t2 = mesh2.getElementData().m_thickness[index2] + pen_enfrc_options.tiny_length; - - if ( t1 < 0. || t2 < 0. ) { - neg_thickness[0] = true; - err[0] = 1; - } - - // compute each element spring stiffness. Pre-multiply the material modulus - // (i.e. material stiffness) by each mesh's penalty scale - auto stiffness1 = pen_scale1 * mesh1.getElementData().m_mat_mod[index1] / t1; - auto stiffness2 = pen_scale2 * mesh2.getElementData().m_mat_mod[index2] / t2; - // compute the equivalent contact penalty spring stiffness per area - penalty_stiff_per_area = ComputePenaltyStiffnessPerArea( stiffness1, stiffness2 ); - break; + RealT stiffness1 = 0.; + RealT stiffness2 = 0.; + RealT t1 = 0.; + RealT t2 = 0.; + RealT mat_mod1 = 0.; + RealT mat_mod2 = 0.; + + if ( pen_enfrc_options.kinematic_calculation == KINEMATIC_ELEMENT ) { + t1 = mesh1.getElementData().m_thickness[index1]; + t2 = mesh2.getElementData().m_thickness[index2]; + mat_mod1 = mesh1.getElementData().m_mat_mod[index1]; + mat_mod2 = mesh2.getElementData().m_mat_mod[index2]; + if ( t1 <= 0. || t2 <= 0. ) { + neg_thickness[0] = true; + err[0] = 1; + return; } - default: - // no-op, quiet compiler - break; - } // end switch on kinematic penalty calculation option + } + + ComputeUncoupledStiffness( pen_enfrc_options.kinematic_calculation, t1, pen_enfrc_options.tiny_length, pen_scale1, + mat_mod1, mesh1.getElementData().m_penalty_stiffness, stiffness1 ); + ComputeUncoupledStiffness( pen_enfrc_options.kinematic_calculation, t2, pen_enfrc_options.tiny_length, pen_scale2, + mat_mod2, mesh2.getElementData().m_penalty_stiffness, stiffness2 ); + + penalty_stiff_per_area = ComputePenaltyStiffnessPerArea( stiffness1, stiffness2 ); //////////////////////////////////////////////////// // Compute contact pressure(s) on current overlap // diff --git a/src/tribol/physics/CommonPlane.hpp b/src/tribol/physics/CommonPlane.hpp index a7b0ef8c..11bd3d9b 100644 --- a/src/tribol/physics/CommonPlane.hpp +++ b/src/tribol/physics/CommonPlane.hpp @@ -7,8 +7,44 @@ #define SRC_TRIBOL_PHYSICS_COMMONPLANE_HPP_ #include "Physics.hpp" +#include "tribol/common/Parameters.hpp" namespace tribol { +/*! + * + * \brief computes the individual (uncoupled) kinematic penalty spring stiffness for a single face + * + * \details This routine calculates the stiffness of a single contact spring on one side of + * the contact pair (either mortar or non-mortar), completely independent of (uncoupled from) + * the other side. The resulting uncoupled stiffnesses from both sides are later combined + * in series (harmonic mean) via ComputePenaltyStiffnessPerArea() to obtain the coupled + * equivalent interface stiffness per unit area. + * + * \param [in] kinematic_calc calculation option (constant or element-based) + * \param [in] thickness element thickness + * \param [in] tiny_length numeric safeguard offset to avoid division by zero + * \param [in] pen_scale mesh penalty scale + * \param [in] mat_mod material modulus (bulk modulus) + * \param [in] const_penalty constant penalty value + * \param [out] stiffness computed uncoupled spring stiffness for this face/side + * + */ +TRIBOL_HOST_DEVICE inline void ComputeUncoupledStiffness( const KinematicPenaltyCalculation kinematic_calc, + const RealT thickness, const RealT tiny_length, + const RealT pen_scale, const RealT mat_mod, + const RealT const_penalty, RealT& stiffness ) +{ + if ( kinematic_calc == KINEMATIC_CONSTANT ) { + stiffness = pen_scale * const_penalty; + } else if ( kinematic_calc == KINEMATIC_ELEMENT ) { + RealT denom = thickness; + if ( denom < tiny_length ) { + denom = tiny_length; + } + stiffness = pen_scale * mat_mod / denom; + } +} + /*! * * \brief computes penalty stiffness for Common Plane + Penalty diff --git a/src/tribol/utils/TestUtils.cpp b/src/tribol/utils/TestUtils.cpp index 6da885f4..9bda6fee 100644 --- a/src/tribol/utils/TestUtils.cpp +++ b/src/tribol/utils/TestUtils.cpp @@ -1109,6 +1109,7 @@ int TestMesh::tribolSetupAndUpdate( ContactMethod method, EnforcementMethod enfo model, enforcement, BINNING_GRID, ExecutionMode::Sequential ); enableTimestepVote( csIndex, params.enable_timestep_vote ); + enableTimestepStabilityLimits( csIndex, params.enable_timestep_stability_limits ); setTimestepPenFrac( csIndex, params.timestep_pen_frac ); setTimestepScale( csIndex, params.timestep_scale ); @@ -1917,4 +1918,4 @@ mfem::Vector ExplicitMechanics::ComputeInvMass( mfem::ParFiniteElementSpace& fes #endif -} // namespace mfem_ext \ No newline at end of file +} // namespace mfem_ext diff --git a/src/tribol/utils/TestUtils.hpp b/src/tribol/utils/TestUtils.hpp index 52272275..d0fc109d 100644 --- a/src/tribol/utils/TestUtils.hpp +++ b/src/tribol/utils/TestUtils.hpp @@ -48,6 +48,7 @@ struct TestControlParameters { rate_penalty_ratio( 0.0 ), const_penalty( 1.0 ), enable_timestep_vote( false ), + enable_timestep_stability_limits( false ), timestep_pen_frac( 0.30 ), timestep_scale( 1.0 ) { @@ -69,6 +70,7 @@ struct TestControlParameters { RealT rate_penalty_ratio; RealT const_penalty; bool enable_timestep_vote; + bool enable_timestep_stability_limits; RealT timestep_pen_frac; RealT timestep_scale; };