From 363e7e8fc8243b3ac69391a10b5b283ec96f58fa Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Mon, 27 Jul 2026 15:31:53 -0700 Subject: [PATCH 01/14] bug fix --- src/redecomp/transfer/MatrixTransfer.cpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/redecomp/transfer/MatrixTransfer.cpp b/src/redecomp/transfer/MatrixTransfer.cpp index b4c7218c..e9d22f86 100644 --- a/src/redecomp/transfer/MatrixTransfer.cpp +++ b/src/redecomp/transfer/MatrixTransfer.cpp @@ -13,6 +13,8 @@ #include "shared/math/ParSparseMat.hpp" #include "shared/infrastructure/Profiling.hpp" +#include + namespace redecomp { void MatrixTransfer::validateTransferInputs( const axom::Array& test_elem_idx, @@ -302,9 +304,9 @@ shared::ParSparseMat MatrixTransfer::TransferToParallel( const axom::Array& } auto num_rows = parent_test_fes_.GetVSize(); - int* I_ptr = new int[num_rows + 1]; - HYPRE_BigInt* J_ptr = new HYPRE_BigInt[num_unique_nonzeros]; - double* data_ptr = new double[num_unique_nonzeros]; + std::vector I_ptr( num_rows + 1, 0 ); + std::vector J_ptr( num_unique_nonzeros ); + std::vector data_ptr( num_unique_nonzeros ); // Initialize I_ptr with zeros for ( int i = 0; i <= num_rows; ++i ) { @@ -337,7 +339,7 @@ shared::ParSparseMat MatrixTransfer::TransferToParallel( const axom::Array& // Construct rectangular HypreParMatrix shared::ParSparseMat J_full( getMPIUtility().MPIComm(), num_rows, parent_test_fes_.GlobalVSize(), - parent_trial_fes_.GlobalVSize(), I_ptr, J_ptr, data_ptr, + parent_trial_fes_.GlobalVSize(), I_ptr.data(), J_ptr.data(), data_ptr.data(), parent_test_fes_.GetDofOffsets(), parent_trial_fes_.GetDofOffsets() ); if ( !parallel_assemble ) { From 0aaece2d14a0f0aecd40a292210001b5ce32fa19 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Mon, 27 Jul 2026 15:33:54 -0700 Subject: [PATCH 02/14] add quadrature point gap option --- src/tribol/common/Parameters.hpp | 12 ++ src/tribol/interface/tribol.cpp | 14 ++ src/tribol/interface/tribol.hpp | 8 + src/tribol/physics/ContactFormulation.hpp | 5 + .../physics/ContactFormulationFactory.cpp | 2 +- src/tribol/physics/EnergyMortar.cpp | 153 +++++++++++++++--- src/tribol/physics/EnergyMortar.hpp | 17 ++ src/tribol/physics/EnergyMortarAdapter.cpp | 107 +++++++++++- src/tribol/physics/EnergyMortarAdapter.hpp | 15 +- 9 files changed, 311 insertions(+), 22 deletions(-) diff --git a/src/tribol/common/Parameters.hpp b/src/tribol/common/Parameters.hpp index 6bb52366..ef034898 100644 --- a/src/tribol/common/Parameters.hpp +++ b/src/tribol/common/Parameters.hpp @@ -122,6 +122,15 @@ enum ContactMethod // all mortar methods go first NUM_CONTACT_METHODS }; +/*! + * \brief Enumerates penalty enforcement options for ENERGY_MORTAR + */ +enum class EnergyMortarPenaltyMode +{ + QUADRATURE_POINT_GAP, ///! Penalize gaps independently at quadrature points + NODAL_GAP ///! Penalize assembled nodal gaps +}; + /*! * \brief Enumerates the available contact model options. * @@ -495,6 +504,9 @@ struct Parameters { // constituent face elements, then we don't consider the face-pair a contact candidate. // Note, auto-contact will require registration of element thicknesses. bool auto_contact_check = false; ///! True if auto-contact checks should be enabled + + EnergyMortarPenaltyMode energy_mortar_penalty_mode = + EnergyMortarPenaltyMode::QUADRATURE_POINT_GAP; ///! Penalty enforcement mode used by ENERGY_MORTAR }; } // namespace tribol diff --git a/src/tribol/interface/tribol.cpp b/src/tribol/interface/tribol.cpp index 895381ab..8f057afe 100644 --- a/src/tribol/interface/tribol.cpp +++ b/src/tribol/interface/tribol.cpp @@ -170,6 +170,20 @@ void setTimestepPenFrac( IndexT cs_id, RealT frac ) } // end setTimestepPenFrac() +//------------------------------------------------------------------------------ +void setEnergyMortarPenaltyMode( IndexT cs_id, EnergyMortarPenaltyMode mode ) +{ + auto cs = CouplingSchemeManager::getInstance().findData( cs_id ); + + SLIC_ERROR_ROOT_IF( !cs, "tribol::setEnergyMortarPenaltyMode(): call tribol::registerCouplingScheme() " + << "prior to calling this routine." ); + + cs->getParameters().energy_mortar_penalty_mode = mode; + if ( cs->hasContactFormulation() ) { + cs->getContactFormulation()->updateEnergyMortarPenaltyMode( mode ); + } +} + //------------------------------------------------------------------------------ void setTimestepScale( IndexT cs_id, RealT scale ) { diff --git a/src/tribol/interface/tribol.hpp b/src/tribol/interface/tribol.hpp index 36e6fa03..6e006218 100644 --- a/src/tribol/interface/tribol.hpp +++ b/src/tribol/interface/tribol.hpp @@ -131,6 +131,14 @@ void setAutoContactPenScale( IndexT cs_id, RealT scale ); */ void setTimestepPenFrac( IndexT cs_id, RealT frac ); +/*! + * \brief Sets the penalty enforcement mode for ENERGY_MORTAR + * + * \param [in] cs_id coupling scheme id + * \param [in] mode penalty mode + */ +void setEnergyMortarPenaltyMode( IndexT cs_id, EnergyMortarPenaltyMode mode ); + /*! * * \brief sets the timestep scale factor applied to the timestep vote diff --git a/src/tribol/physics/ContactFormulation.hpp b/src/tribol/physics/ContactFormulation.hpp index c48e37f3..40392127 100644 --- a/src/tribol/physics/ContactFormulation.hpp +++ b/src/tribol/physics/ContactFormulation.hpp @@ -118,6 +118,11 @@ class ContactFormulation { */ virtual void updateConstantPenaltyStiffness( double /*mesh1_penalty*/, double /*mesh2_penalty*/ ) {} + /** + * @brief Update ENERGY_MORTAR penalty mode on formulations that support it + */ + virtual void updateEnergyMortarPenaltyMode( EnergyMortarPenaltyMode /*mode*/ ) {} + #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 beaed413..7f695292 100644 --- a/src/tribol/physics/ContactFormulationFactory.cpp +++ b/src/tribol/physics/ContactFormulationFactory.cpp @@ -45,7 +45,7 @@ std::unique_ptr createContactFormulation( CouplingScheme* cs return std::make_unique( *cs->getMfemMeshData(), *cs->getMfemSubmeshData(), *cs->getMfemJacobianData(), k, delta, N, enzyme_quadrature, - use_penalty_ ); + use_penalty_, cs->getParameters().energy_mortar_penalty_mode ); #else SLIC_ERROR_ROOT( "ENERGY_MORTAR requires Enzyme and redecomp to be built." ); return nullptr; diff --git a/src/tribol/physics/EnergyMortar.cpp b/src/tribol/physics/EnergyMortar.cpp index 87df94e9..a362c329 100644 --- a/src/tribol/physics/EnergyMortar.cpp +++ b/src/tribol/physics/EnergyMortar.cpp @@ -19,10 +19,16 @@ namespace { // This MUST match what the ContactParams struct has in EnergyMortarAdapter // Theese 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 + int N{ 3 }; // No. of quadrature points + double del{ 0.1 }; // Smoothing parameter + double k{ 1.0 }; // Penalty stiffness }; +TRIBOL_ENZYME_INLINE double line_jacobian( const double* A0, const double* A1 ) +{ + return std::sqrt( ( A1[0] - A0[0] ) * ( A1[0] - A0[0] ) + ( A1[1] - A0[1] ) * ( A1[1] - A0[1] ) ); +} + // Compute a unit normal vector for the line segment from coord1 to coord2 TRIBOL_ENZYME_INLINE void find_normal( const double* coord1, const double* coord2, double* normal ) { @@ -386,9 +392,9 @@ void grad_kernel( const double* x, const Gparams* gp, double* dout_du ) // Wrap the varying-quadrature kernel as a scalar-valued function for Enzyme. template -static void kernel_out_enzyme( const double* x, double* out ) +static void kernel_out_enzyme( const double* x, const void* kp_void, double* out ) { - KernelParams kp; + const KernelParams* kp = static_cast( kp_void ); // x stores the two endpoints of edge A followed by the two endpoints of edge B. double A0[2], A1[2], B0[2], B1[2]; A0[0] = x[0]; @@ -405,10 +411,10 @@ static void kernel_out_enzyme( const double* x, double* out ) std::array projections = { projs[0], projs[1] }; // Recompute the integration bounds and quadrature from the current geometry. - auto bounds = ContactSmoothing::bounds_from_projections( projections, kp.del ); - auto xi_bounds = ContactSmoothing::smooth_bounds( bounds, kp.del ); + auto bounds = ContactSmoothing::bounds_from_projections( projections, kp->del ); + auto xi_bounds = ContactSmoothing::smooth_bounds( bounds, kp->del ); - auto qp = EnergyMortarCalculator::compute_quadrature( xi_bounds, kp.N ); + auto qp = EnergyMortarCalculator::compute_quadrature( xi_bounds, kp->N ); Gparams gp; for ( std::size_t i = 0; i < qp.qp.size(); ++i ) { @@ -433,14 +439,15 @@ static void kernel_out_enzyme( const double* x, double* out ) // Differentiate the selected varying-quadrature scalar kernel with respect to the 8 endpoint coordinates. template -void grad_kernel_enzyme( const double* x, double* dout_du ) +void grad_kernel_enzyme( const double* x, const KernelParams* kp, double* dout_du ) { double dx[8] = { 0.0 }; double out = 0.0; double dout = 1.0; // Seed the scalar output with 1.0 so Enzyme accumulates dOutput/dx into dx. - __enzyme_autodiff( (void*)kernel_out_enzyme, enzyme_dup, x, dx, enzyme_dup, &out, &dout ); + __enzyme_autodiff( (void*)kernel_out_enzyme, enzyme_dup, x, dx, enzyme_const, (const void*)kp, + enzyme_dup, &out, &dout ); for ( int i = 0; i < 8; ++i ) { dout_du[i] = dx[i]; @@ -449,7 +456,7 @@ void grad_kernel_enzyme( const double* x, double* dout_du ) // Compute the Hessian of the selected varying-quadrature scalar kernel. template -void d2_kernel( const double* x, double* H ) +void d2_kernel( const double* x, const KernelParams* kp, double* H ) { for ( int col = 0; col < 8; ++col ) { double dx[8] = { 0.0 }; @@ -459,12 +466,82 @@ void d2_kernel( const double* x, double* H ) double dgrad[8] = { 0.0 }; // Differentiate the gradient in coordinate direction col to form one Hessian column. - __enzyme_fwddiff( (void*)grad_kernel_enzyme, enzyme_dup, x, dx, enzyme_dup, grad, dgrad ); + __enzyme_fwddiff( (void*)grad_kernel_enzyme, enzyme_dup, x, dx, enzyme_const, (const void*)kp, + enzyme_dup, grad, dgrad ); for ( int row = 0; row < 8; ++row ) H[row * 8 + col] = dgrad[row]; } } +TRIBOL_ENZYME_INLINE void qp_penalty_kernel( const double* x, const KernelParams* kp, double* energy ) +{ + double A0[2] = { x[0], x[1] }; + double A1[2] = { x[2], x[3] }; + double B0[2] = { x[4], x[5] }; + double B1[2] = { x[6], x[7] }; + + double projs[2] = { 0.0, 0.0 }; + get_projections( A0, A1, B0, B1, projs ); + const std::array projections{ projs[0], projs[1] }; + const auto bounds = ContactSmoothing::bounds_from_projections( projections, kp->del ); + const auto xi_bounds = ContactSmoothing::smooth_bounds( bounds, kp->del ); + const auto qp = EnergyMortarCalculator::compute_quadrature( xi_bounds, kp->N ); + + double nB[2]; + find_normal( B0, B1, nB ); + double nA[2]; + find_normal( A0, A1, nA ); + const double eta = nA[0] * nB[0] + nA[1] * nB[1]; + const double J = line_jacobian( A0, A1 ); + + double value = 0.0; + for ( int i = 0; i < kp->N; ++i ) { + double x1[2]; + iso_map( A0, A1, qp.qp[i], x1 ); + + double x2[2]; + find_intersection( B0, B1, x1, nB, x2 ); + + 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; + if ( gap < 0.0 ) { + value += 0.5 * kp->k * gap * gap * qp.w[i] * J; + } + } + + *energy = value; +} + +void grad_qp_penalty_kernel( const double* x, const KernelParams* kp, double* dout_du ) +{ + double dx[8] = { 0.0 }; + double out = 0.0; + double dout = 1.0; + __enzyme_autodiff( (void*)qp_penalty_kernel, enzyme_dup, x, dx, enzyme_const, (const void*)kp, enzyme_dup, + &out, &dout ); + + for ( int i = 0; i < 8; ++i ) { + dout_du[i] = dx[i]; + } +} + +void d2_qp_penalty_kernel( const double* x, const KernelParams* kp, double* H ) +{ + for ( int col = 0; col < 8; ++col ) { + double dx[8] = { 0.0 }; + dx[col] = 1.0; + double grad[8] = { 0.0 }; + double dgrad[8] = { 0.0 }; + __enzyme_fwddiff( (void*)grad_qp_penalty_kernel, enzyme_dup, x, dx, enzyme_const, (const void*)kp, + enzyme_dup, grad, dgrad ); + for ( int row = 0; row < 8; ++row ) { + H[row * 8 + col] = dgrad[row]; + } + } +} + // Compute the Hessian of the selected fixed-quadrature scalar kernel. template void d2_kernel_quad( const double* x, const Gparams* gp, double* H ) @@ -752,8 +829,9 @@ void EnergyMortarCalculator::grad_gtilde( const InterfacePair& pair, const MeshD } else { // Differentiate through the geometry-dependent quadrature construction. - grad_kernel_enzyme( x, dg1_du ); - grad_kernel_enzyme( x, dg2_du ); + const KernelParams kp{ p_.N, p_.del, p_.k }; + grad_kernel_enzyme( x, &kp, dg1_du ); + grad_kernel_enzyme( x, &kp, dg2_du ); } for ( int i = 0; i < 8; ++i ) { @@ -784,8 +862,9 @@ void EnergyMortarCalculator::grad_trib_area( const InterfacePair& pair, const Me grad_kernel( x, &gp, dA2_dx ); } else { // Differentiate through the geometry-dependent quadrature construction. - grad_kernel_enzyme( x, dA1_dx ); - grad_kernel_enzyme( x, dA2_dx ); + const KernelParams kp{ p_.N, p_.del, p_.k }; + grad_kernel_enzyme( x, &kp, dA1_dx ); + grad_kernel_enzyme( x, &kp, dA2_dx ); } } @@ -815,8 +894,9 @@ void EnergyMortarCalculator::d2_g2tilde( const InterfacePair& pair, const MeshDa } else { // Differentiate through the geometry-dependent quadrature construction. - d2_kernel( x, d2g1_d2u ); - d2_kernel( x, d2g2_d2u ); + const KernelParams kp{ p_.N, p_.del, p_.k }; + d2_kernel( x, &kp, d2g1_d2u ); + d2_kernel( x, &kp, d2g2_d2u ); } for ( int i = 0; i < 64; ++i ) { @@ -851,8 +931,9 @@ 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. - d2_kernel( x, d2A1_d2u ); - d2_kernel( x, d2A2_d2u ); + const KernelParams kp{ p_.N, p_.del, p_.k }; + d2_kernel( x, &kp, d2A1_d2u ); + d2_kernel( x, &kp, d2A2_d2u ); } for ( int i = 0; i < 64; ++i ) { @@ -861,6 +942,40 @@ void EnergyMortarCalculator::compute_d2A_d2u( const InterfacePair& pair, const M } } +double EnergyMortarCalculator::compute_quadrature_point_penalty_energy( const InterfacePair& pair, + const MeshData::Viewer& mesh1, + const MeshData::Viewer& mesh2 ) const +{ + double A0[2], A1[2], B0[2], B1[2]; + + endpoints( mesh1, pair.m_element_id1, A0, A1 ); + 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 }; + double energy = 0.0; + qp_penalty_kernel( x, &kp, &energy ); + return energy; +} + +QuadraturePointPenaltyData EnergyMortarCalculator::compute_quadrature_point_penalty_data( + const InterfacePair& pair, const MeshData::Viewer& mesh1, const MeshData::Viewer& mesh2 ) const +{ + double A0[2], A1[2], B0[2], B1[2]; + + endpoints( mesh1, pair.m_element_id1, A0, A1 ); + 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 }; + + QuadraturePointPenaltyData result; + qp_penalty_kernel( x, &kp, &result.energy ); + grad_qp_penalty_kernel( x, &kp, result.force.data() ); + d2_qp_penalty_kernel( x, &kp, result.stiffness.data() ); + return result; +} + #endif // TRIBOL_USE_ENZYME } // namespace tribol diff --git a/src/tribol/physics/EnergyMortar.hpp b/src/tribol/physics/EnergyMortar.hpp index 66961fe7..08a5335c 100644 --- a/src/tribol/physics/EnergyMortar.hpp +++ b/src/tribol/physics/EnergyMortar.hpp @@ -23,6 +23,14 @@ struct ContactParams { double k; // Penalty int N; // Quadrature Points bool enzyme_quadrature; // Determines how enzyming is performed (default = True) + EnergyMortarPenaltyMode penalty_mode{ EnergyMortarPenaltyMode::QUADRATURE_POINT_GAP }; // Penalty enforcement mode +}; + +/// Stores quadrature-point penalty energy derivatives for one interface pair. +struct QuadraturePointPenaltyData { + double energy{ 0.0 }; + std::array force{}; + std::array stiffness{}; }; // Weighted gap and trib area @@ -170,6 +178,15 @@ class EnergyMortarCalculator { void compute_d2A_d2u( const InterfacePair& pair, const MeshData::Viewer& mesh1, const MeshData::Viewer& mesh2, double dgt1_dx[64], double dgt2_dx[64] ) const; + /// Compute local energy, force, and stiffness for quadrature-point penalty enforcement. + QuadraturePointPenaltyData compute_quadrature_point_penalty_data( const InterfacePair& pair, + const MeshData::Viewer& mesh1, + const MeshData::Viewer& mesh2 ) const; + + /// Evaluate only the local quadrature-point penalty energy. + double compute_quadrature_point_penalty_energy( const InterfacePair& pair, const MeshData::Viewer& mesh1, + const MeshData::Viewer& mesh2 ) const; + /// Evaluate and return the two nodal smoothed gap integrals. /// /// This is a convenience wrapper for obtaining only the gap integral diff --git a/src/tribol/physics/EnergyMortarAdapter.cpp b/src/tribol/physics/EnergyMortarAdapter.cpp index 3be37931..28aad634 100644 --- a/src/tribol/physics/EnergyMortarAdapter.cpp +++ b/src/tribol/physics/EnergyMortarAdapter.cpp @@ -13,7 +13,8 @@ namespace tribol { EnergyMortarAdapter::EnergyMortarAdapter( MfemMeshData& mesh_data, MfemSubmeshData& submesh_data, MfemJacobianData& jac_data, double k, double delta, int N, - bool enzyme_quadrature, bool use_penalty ) + bool enzyme_quadrature, bool use_penalty, + EnergyMortarPenaltyMode penalty_mode ) // NOTE: mesh1 maps to mesh2_ and mesh2 maps to mesh1_. This is to keep consistent with mesh1_ being non-mortar and // mesh2_ being mortar as is typical in the literature, but different from Tribol convention. : use_penalty_( use_penalty ), mesh_data_( mesh_data ), submesh_data_( submesh_data ), jac_data_( jac_data ) @@ -22,6 +23,7 @@ EnergyMortarAdapter::EnergyMortarAdapter( MfemMeshData& mesh_data, MfemSubmeshDa params_.del = delta; params_.N = N; params_.enzyme_quadrature = enzyme_quadrature; + params_.penalty_mode = penalty_mode; evaluator_ = std::make_unique( params_ ); @@ -43,6 +45,13 @@ void EnergyMortarAdapter::updateConstantPenaltyStiffness( double mesh1_penalty, { use_penalty_ = true; params_.k = 0.5 * ( mesh1_penalty + mesh2_penalty ); + evaluator_ = std::make_unique( params_ ); +} + +void EnergyMortarAdapter::updateEnergyMortarPenaltyMode( EnergyMortarPenaltyMode mode ) +{ + params_.penalty_mode = mode; + evaluator_ = std::make_unique( params_ ); } const mfem::HypreParVector& EnergyMortarAdapter::getMfemGap() const @@ -211,6 +220,11 @@ void EnergyMortarAdapter::updateNodalForces() // NOTE: user should have called updateNodalGaps() with updated coords before calling this if ( use_penalty_ ) { + if ( params_.penalty_mode == EnergyMortarPenaltyMode::QUADRATURE_POINT_GAP ) { + updateQuadraturePointPenaltyForces(); + return; + } + // Penalty mode: p = k * (g_tilde / A) pressure_vec_ = params_.k * gap_vec_; } else { @@ -452,6 +466,97 @@ shared::ParSparseMat EnergyMortarAdapter::computeDfDxSecondDerivativesPenalty( mesh_data_.GetParentCoords().ParFESpace(), df_contribs ); } +void EnergyMortarAdapter::updateQuadraturePointPenaltyForces() +{ + const bool use_lor = ( mesh_data_.GetLORMesh() != nullptr ); + const auto& displacement_surface_fes = use_lor ? *mesh_data_.GetLORMeshFESpace() : mesh_data_.GetSubmeshFESpace(); + const auto& displacement_redecomp_fes = *mesh_data_.GetRedecompResponse().FESpace(); + const auto& mortar_elem_map = mesh_data_.GetElemMap1(); + const auto& nonmortar_elem_map = mesh_data_.GetElemMap2(); + + PackedPairJacobianContribs df_nm_nm( displacement_surface_fes, displacement_surface_fes, displacement_redecomp_fes, + displacement_redecomp_fes, nonmortar_elem_map, nonmortar_elem_map ); + PackedPairJacobianContribs df_nm_m( displacement_surface_fes, displacement_surface_fes, displacement_redecomp_fes, + displacement_redecomp_fes, nonmortar_elem_map, mortar_elem_map ); + PackedPairJacobianContribs df_m_nm( displacement_surface_fes, displacement_surface_fes, displacement_redecomp_fes, + displacement_redecomp_fes, mortar_elem_map, nonmortar_elem_map ); + PackedPairJacobianContribs df_m_m( displacement_surface_fes, displacement_surface_fes, displacement_redecomp_fes, + displacement_redecomp_fes, mortar_elem_map, mortar_elem_map ); + + df_nm_nm.reserve( pairs_.size(), 16 ); + df_nm_m.reserve( pairs_.size(), 16 ); + df_m_nm.reserve( pairs_.size(), 16 ); + df_m_m.reserve( pairs_.size(), 16 ); + + mfem::GridFunction redecomp_force( const_cast( &displacement_redecomp_fes ) ); + redecomp_force = 0.0; + const int scalar_size = redecomp_force.FESpace()->GetVSize() / redecomp_force.FESpace()->GetVDim(); + energy_ = 0.0; + + const int node_idx[8] = { 0, 2, 1, 3, 4, 6, 5, 7 }; + + SLIC_ERROR_ROOT_IF( mesh1_ == nullptr || mesh2_ == nullptr, "ENERGY_MORTAR meshes not set." ); + auto mesh1_view = mesh1_->getView(); + auto mesh2_view = mesh2_->getView(); + + for ( const auto& pair : pairs_ ) { + InterfacePair flipped_pair( pair.m_element_id2, pair.m_element_id1 ); + const auto elem1 = static_cast( flipped_pair.m_element_id1 ); + const auto elem2 = static_cast( flipped_pair.m_element_id2 ); + const auto qp_data = evaluator_->compute_quadrature_point_penalty_data( flipped_pair, mesh1_view, mesh2_view ); + + if ( qp_data.energy == 0.0 ) { + continue; + } + + energy_ += qp_data.energy; + + auto A_conn = mesh1_view.getConnectivity()( elem1 ); + auto B_conn = mesh2_view.getConnectivity()( elem2 ); + redecomp_force( A_conn[0] ) += qp_data.force[0]; + redecomp_force( scalar_size + A_conn[0] ) += qp_data.force[1]; + redecomp_force( A_conn[1] ) += qp_data.force[2]; + redecomp_force( scalar_size + A_conn[1] ) += qp_data.force[3]; + redecomp_force( B_conn[0] ) += qp_data.force[4]; + redecomp_force( scalar_size + B_conn[0] ) += qp_data.force[5]; + redecomp_force( B_conn[1] ) += qp_data.force[6]; + redecomp_force( scalar_size + B_conn[1] ) += qp_data.force[7]; + + double df_dx_blocks[2][2][16]; + for ( int i{ 0 }; i < 2; ++i ) { + for ( int j{ 0 }; j < 2; ++j ) { + for ( int k{ 0 }; k < 4; ++k ) { + for ( int l{ 0 }; l < 4; ++l ) { + const auto idx = node_idx[l + i * 4] + node_idx[k + j * 4] * 8; + df_dx_blocks[i][j][l + k * 4] = qp_data.stiffness[idx]; + } + } + } + } + + df_nm_nm.append( elem1, elem1, df_dx_blocks[0][0], 16 ); + df_nm_m.append( elem1, elem2, df_dx_blocks[0][1], 16 ); + df_m_nm.append( elem2, elem1, df_dx_blocks[1][0], 16 ); + df_m_m.append( elem2, elem2, df_dx_blocks[1][1], 16 ); + } + + auto* parent_fes = mesh_data_.GetParentCoords().ParFESpace(); + force_vec_ = shared::ParVector( const_cast( parent_fes ) ); + force_vec_.fill( 0.0 ); + mfem::Vector parent_force( parent_fes->GetVSize() ); + parent_force = 0.0; + mesh_data_.GetParentRedecompTransfer().RedecompToParent( redecomp_force, parent_force ); + parent_fes->GetProlongationMatrix()->MultTranspose( parent_force, force_vec_.get() ); + + std::vector df_contribs; + df_contribs.reserve( 4 ); + df_contribs.push_back( std::move( df_nm_nm ) ); + df_contribs.push_back( std::move( df_nm_m ) ); + df_contribs.push_back( std::move( df_m_nm ) ); + df_contribs.push_back( std::move( df_m_m ) ); + df_dx_ = jac_data_.GetMfemJacobian( parent_fes, parent_fes, df_contribs ); +} + std::unique_ptr EnergyMortarAdapter::getMfemDfDx() const { return std::unique_ptr( df_dx_.release() ); diff --git a/src/tribol/physics/EnergyMortarAdapter.hpp b/src/tribol/physics/EnergyMortarAdapter.hpp index 41e7d599..d6f57e3a 100644 --- a/src/tribol/physics/EnergyMortarAdapter.hpp +++ b/src/tribol/physics/EnergyMortarAdapter.hpp @@ -50,7 +50,8 @@ class EnergyMortarAdapter : public ContactFormulation { * relative to the order of the meshes provided here. */ EnergyMortarAdapter( MfemMeshData& mesh_data, MfemSubmeshData& submesh_data, MfemJacobianData& jac_data, double k, - double delta, int N, bool enzyme_quadrature, bool use_penalty = true ); + double delta, int N, bool enzyme_quadrature, bool use_penalty = true, + EnergyMortarPenaltyMode penalty_mode = EnergyMortarPenaltyMode::QUADRATURE_POINT_GAP ); /** * @brief Default destructor @@ -132,6 +133,13 @@ class EnergyMortarAdapter : public ContactFormulation { */ void updateConstantPenaltyStiffness( double mesh1_penalty, double mesh2_penalty ) override; + /** + * @brief Update the EnergyMortar penalty mode + * + * @param mode Penalty mode + */ + void updateEnergyMortarPenaltyMode( EnergyMortarPenaltyMode mode ) override; + #ifdef BUILD_REDECOMP /** * @brief Return the parent true-dof force vector @@ -333,6 +341,11 @@ class EnergyMortarAdapter : public ContactFormulation { shared::ParSparseMat computeDfDxSecondDerivativesPenalty( const mfem::GridFunction& redecomp_pressure, const mfem::GridFunction& redecomp_g_tilde, const mfem::GridFunction& redecomp_A ); + + /** + * @brief Assemble QP-gap penalty force and Jacobian contributions + */ + void updateQuadraturePointPenaltyForces(); }; #endif // TRIBOL_USE_ENZYME From a829c3a57024ed1330788f0b435d326a4e492dc1 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Mon, 27 Jul 2026 23:30:57 -0700 Subject: [PATCH 03/14] update tests for quad pt gap --- src/tests/tribol_energy_mortar_patch.cpp | 3 +- .../tribol_finite_diff_energy_mortar.cpp | 117 +++++++++++++++++- src/tests/tribol_mfem_jacobian.cpp | 1 + 3 files changed, 119 insertions(+), 2 deletions(-) diff --git a/src/tests/tribol_energy_mortar_patch.cpp b/src/tests/tribol_energy_mortar_patch.cpp index f1f17caa..6f44ac93 100644 --- a/src/tests/tribol_energy_mortar_patch.cpp +++ b/src/tests/tribol_energy_mortar_patch.cpp @@ -175,7 +175,8 @@ class MfemMortarEnergyPatchTest : public testing::TestWithParam> tribol::registerMfemCouplingScheme( cs_id, mesh1_id, mesh2_id, mesh, coords, mortar_attrs, nonmortar_attrs, tribol::SURFACE_TO_SURFACE, tribol::NO_SLIDING, tribol::ENERGY_MORTAR, tribol::FRICTIONLESS, tribol::PENALTY, tribol::BINNING_GRID ); - tribol::setMfemKinematicConstantPenalty( cs_id, 10000.0, 10000.0 ); + tribol::setEnergyMortarPenaltyMode( cs_id, tribol::EnergyMortarPenaltyMode::QUADRATURE_POINT_GAP ); + tribol::setMfemKinematicConstantPenalty( cs_id, 100.0, 100.0 ); mfem::Vector X( par_fe_space.GetTrueVSize() ); X = 0.0; diff --git a/src/tests/tribol_finite_diff_energy_mortar.cpp b/src/tests/tribol_finite_diff_energy_mortar.cpp index 450a5596..3a75e59b 100644 --- a/src/tests/tribol_finite_diff_energy_mortar.cpp +++ b/src/tests/tribol_finite_diff_energy_mortar.cpp @@ -372,6 +372,121 @@ FiniteDiffResult EnergyMortarCalculator::validate_hessian( const InterfacePair& return result; } + +TEST( QuadraturePointPenaltyCheck, OpenGapIsInactive ) +{ + RealT x1[2] = { 0.0, 1.0 }; + RealT y1[2] = { 0.0, 0.0 }; + IndexT conn1[2] = { 1, 0 }; + MeshData mesh1( 0, 1, 2, conn1, LINEAR_EDGE, x1, y1, nullptr, MemorySpace::Host ); + + RealT x2[2] = { 0.2, 0.8 }; + RealT y2[2] = { 0.1, 0.1 }; + IndexT conn2[2] = { 0, 1 }; + MeshData mesh2( 1, 1, 2, conn2, LINEAR_EDGE, x2, y2, nullptr, MemorySpace::Host ); + + ContactParams params; + params.del = 0.1; + params.k = 3.0; + params.N = 3; + params.enzyme_quadrature = true; + params.penalty_mode = EnergyMortarPenaltyMode::QUADRATURE_POINT_GAP; + + EnergyMortarCalculator evaluator( params ); + const auto result = evaluator.compute_quadrature_point_penalty_data( InterfacePair( 0, 0 ), mesh1.getView(), + mesh2.getView() ); + EXPECT_EQ( result.energy, 0.0 ); +} + +TEST( QuadraturePointPenaltyCheck, DerivativesMatchFiniteDifference ) +{ + RealT x1[2] = { 0.0, 1.0 }; + RealT y1[2] = { 0.0, 0.0 }; + IndexT conn1[2] = { 1, 0 }; + MeshData mesh1( 0, 1, 2, conn1, LINEAR_EDGE, x1, y1, nullptr, MemorySpace::Host ); + + RealT x2[2] = { 0.2, 0.8 }; + RealT y2[2] = { -0.1, -0.1 }; + IndexT conn2[2] = { 0, 1 }; + MeshData mesh2( 1, 1, 2, conn2, LINEAR_EDGE, x2, y2, nullptr, MemorySpace::Host ); + + ContactParams params; + params.del = 0.1; + params.k = 3.0; + params.N = 3; + params.enzyme_quadrature = true; + params.penalty_mode = EnergyMortarPenaltyMode::QUADRATURE_POINT_GAP; + + EnergyMortarCalculator evaluator( params ); + const InterfacePair pair( 0, 0 ); + const auto analytical = evaluator.compute_quadrature_point_penalty_data( pair, mesh1.getView(), mesh2.getView() ); + ASSERT_GT( analytical.energy, 0.0 ); + + const std::array x1_orig{ x1[0], x1[1] }; + const std::array y1_orig{ y1[0], y1[1] }; + const std::array x2_orig{ x2[0], x2[1] }; + const std::array y2_orig{ y2[0], y2[1] }; + + auto restore = [&]() { + x1[0] = x1_orig[0]; + x1[1] = x1_orig[1]; + y1[0] = y1_orig[0]; + y1[1] = y1_orig[1]; + x2[0] = x2_orig[0]; + x2[1] = x2_orig[1]; + y2[0] = y2_orig[0]; + y2[1] = y2_orig[1]; + mesh1.setPosition( x1, y1, nullptr ); + mesh2.setPosition( x2, y2, nullptr ); + }; + + auto perturb = [&]( int dof, double delta ) { + if ( dof < 4 ) { + const int endpoint = dof / 2; + const int component = dof % 2; + const int node = conn1[endpoint]; + ( component == 0 ? x1[node] : y1[node] ) += delta; + mesh1.setPosition( x1, y1, nullptr ); + } else { + const int endpoint = ( dof - 4 ) / 2; + const int component = ( dof - 4 ) % 2; + const int node = conn2[endpoint]; + ( component == 0 ? x2[node] : y2[node] ) += delta; + mesh2.setPosition( x2, y2, nullptr ); + } + }; + + const double gradient_eps = 1.0e-7; + const double gradient_tol = 1.0e-6; + for ( int dof = 0; dof < 8; ++dof ) { + restore(); + perturb( dof, gradient_eps ); + const double energy_plus = evaluator.compute_quadrature_point_penalty_energy( pair, mesh1.getView(), mesh2.getView() ); + restore(); + perturb( dof, -gradient_eps ); + const double energy_minus = evaluator.compute_quadrature_point_penalty_energy( pair, mesh1.getView(), mesh2.getView() ); + const double fd_force = ( energy_plus - energy_minus ) / ( 2.0 * gradient_eps ); + EXPECT_NEAR( fd_force, analytical.force[dof], gradient_tol ) << "force mismatch at dof " << dof; + } + + const double hessian_eps = 1.0e-6; + const double hessian_tol = 1.0e-4; + for ( int col = 0; col < 8; ++col ) { + restore(); + perturb( col, hessian_eps ); + const auto force_plus = evaluator.compute_quadrature_point_penalty_data( pair, mesh1.getView(), mesh2.getView() ).force; + restore(); + perturb( col, -hessian_eps ); + const auto force_minus = evaluator.compute_quadrature_point_penalty_data( pair, mesh1.getView(), mesh2.getView() ).force; + for ( int row = 0; row < 8; ++row ) { + const double fd_stiffness = ( force_plus[row] - force_minus[row] ) / ( 2.0 * hessian_eps ); + EXPECT_NEAR( fd_stiffness, analytical.stiffness[row * 8 + col], hessian_tol ) + << "stiffness mismatch at row " << row << ", col " << col; + } + } + restore(); +} + TEST( GradientCheck, GtildeFDvsAD ) { // ── Geometry: two facing LINEAR_EDGE segments ──────────────────────────── @@ -486,4 +601,4 @@ TEST( HessianCheck, GtildeFDvsAD ) } } -} // namespace tribol \ No newline at end of file +} // namespace tribol diff --git a/src/tests/tribol_mfem_jacobian.cpp b/src/tests/tribol_mfem_jacobian.cpp index f181799c..78ae5c0d 100644 --- a/src/tests/tribol_mfem_jacobian.cpp +++ b/src/tests/tribol_mfem_jacobian.cpp @@ -1053,6 +1053,7 @@ TEST_F( MfemJacobianTest, mfem_penalty_jacobian_retrieval ) tribol::registerMfemCouplingScheme( cs_id, mesh1_id, mesh2_id, mesh, coords, mortar_attrs, nonmortar_attrs, tribol::SURFACE_TO_SURFACE, tribol::NO_CASE, tribol::ENERGY_MORTAR, tribol::FRICTIONLESS, tribol::PENALTY, tribol::BINNING_GRID ); + tribol::setEnergyMortarPenaltyMode( cs_id, tribol::EnergyMortarPenaltyMode::QUADRATURE_POINT_GAP ); tribol::setPenaltyOptions( cs_id, tribol::KINEMATIC, tribol::KINEMATIC_CONSTANT ); tribol::setMfemKinematicConstantPenalty( cs_id, 1.0, 1.0 ); From 66a7f775645b03a3c041d1dd7dfaf30aa504774d Mon Sep 17 00:00:00 2001 From: "E. B. Chin" Date: Wed, 5 Aug 2026 15:58:36 -0700 Subject: [PATCH 04/14] make enforcement a template class --- src/tests/tribol_energy_mortar_patch.cpp | 2 +- src/tests/tribol_energy_mortar_patch_LM.cpp | 1 + .../tribol_finite_diff_energy_mortar.cpp | 19 +- src/tests/tribol_mfem_jacobian.cpp | 2 +- src/tribol/common/Parameters.hpp | 10 +- src/tribol/interface/tribol.cpp | 32 +- src/tribol/interface/tribol.hpp | 17 +- src/tribol/physics/ContactFormulation.hpp | 5 - .../physics/ContactFormulationFactory.cpp | 12 +- src/tribol/physics/EnergyMortar.hpp | 1 - src/tribol/physics/EnergyMortarAdapter.cpp | 334 +++++++++--------- src/tribol/physics/EnergyMortarAdapter.hpp | 67 ++-- 12 files changed, 262 insertions(+), 240 deletions(-) diff --git a/src/tests/tribol_energy_mortar_patch.cpp b/src/tests/tribol_energy_mortar_patch.cpp index 6f44ac93..755d1b9b 100644 --- a/src/tests/tribol_energy_mortar_patch.cpp +++ b/src/tests/tribol_energy_mortar_patch.cpp @@ -175,7 +175,7 @@ class MfemMortarEnergyPatchTest : public testing::TestWithParam> tribol::registerMfemCouplingScheme( cs_id, mesh1_id, mesh2_id, mesh, coords, mortar_attrs, nonmortar_attrs, tribol::SURFACE_TO_SURFACE, tribol::NO_SLIDING, tribol::ENERGY_MORTAR, tribol::FRICTIONLESS, tribol::PENALTY, tribol::BINNING_GRID ); - tribol::setEnergyMortarPenaltyMode( cs_id, tribol::EnergyMortarPenaltyMode::QUADRATURE_POINT_GAP ); + tribol::setEnergyMortarEnforcementOption( cs_id, tribol::EnergyMortarEnforcementOption::QUADRATURE_POINT_GAP ); tribol::setMfemKinematicConstantPenalty( cs_id, 100.0, 100.0 ); mfem::Vector X( par_fe_space.GetTrueVSize() ); diff --git a/src/tests/tribol_energy_mortar_patch_LM.cpp b/src/tests/tribol_energy_mortar_patch_LM.cpp index a56c7dce..d1f85ded 100644 --- a/src/tests/tribol_energy_mortar_patch_LM.cpp +++ b/src/tests/tribol_energy_mortar_patch_LM.cpp @@ -195,6 +195,7 @@ class MfemMortarEnergyLagrangePatchTest : public testing::TestWithParam +#include + +#include "tribol/physics/ContactFormulationFactory.hpp" #include #include #include @@ -171,16 +175,34 @@ void setTimestepPenFrac( IndexT cs_id, RealT frac ) } // end setTimestepPenFrac() //------------------------------------------------------------------------------ -void setEnergyMortarPenaltyMode( IndexT cs_id, EnergyMortarPenaltyMode mode ) +void setEnergyMortarEnforcementOption( IndexT cs_id, EnergyMortarEnforcementOption mode ) { auto cs = CouplingSchemeManager::getInstance().findData( cs_id ); - SLIC_ERROR_ROOT_IF( !cs, "tribol::setEnergyMortarPenaltyMode(): call tribol::registerCouplingScheme() " + SLIC_ERROR_ROOT_IF( !cs, "tribol::setEnergyMortarEnforcementOption(): call tribol::registerCouplingScheme() " << "prior to calling this routine." ); - cs->getParameters().energy_mortar_penalty_mode = mode; - if ( cs->hasContactFormulation() ) { - cs->getContactFormulation()->updateEnergyMortarPenaltyMode( mode ); + cs->getParameters().energy_mortar_enforcement_option = mode; + + // Automatically rebuild the formulation to reflect the new setting + if ( cs->getContactMethod() == ENERGY_MORTAR ) { + rebuildContactFormulation( cs_id ); + } +} + +//------------------------------------------------------------------------------ +void rebuildContactFormulation( IndexT cs_id ) +{ + auto cs = CouplingSchemeManager::getInstance().findData( cs_id ); + + SLIC_ERROR_ROOT_IF( !cs, "tribol::rebuildContactFormulation(): call tribol::registerCouplingScheme() " + << "prior to calling this routine." ); + + if ( cs->getContactMethod() == ENERGY_MORTAR ) { + cs->setContactFormulation( createContactFormulation( cs ) ); + } else { + SLIC_WARNING_ROOT( + "tribol::rebuildContactFormulation(): rebuilding is only supported for ENERGY_MORTAR at this time." ); } } diff --git a/src/tribol/interface/tribol.hpp b/src/tribol/interface/tribol.hpp index 6e006218..691dfddc 100644 --- a/src/tribol/interface/tribol.hpp +++ b/src/tribol/interface/tribol.hpp @@ -132,12 +132,23 @@ void setAutoContactPenScale( IndexT cs_id, RealT scale ); void setTimestepPenFrac( IndexT cs_id, RealT frac ); /*! - * \brief Sets the penalty enforcement mode for ENERGY_MORTAR + * \brief Sets the enforcement mode for ENERGY_MORTAR * * \param [in] cs_id coupling scheme id - * \param [in] mode penalty mode + * \param [in] mode enforcement option */ -void setEnergyMortarPenaltyMode( IndexT cs_id, EnergyMortarPenaltyMode mode ); +void setEnergyMortarEnforcementOption( IndexT cs_id, EnergyMortarEnforcementOption mode ); + +/*! + * \brief Rebuilds the contact formulation for a coupling scheme. + * + * This function can be called after modifying coupling scheme settings + * (e.g. setEnergyMortarEnforcementOption) to recreate the underlying formulation + * instance with the new settings. + * + * \param [in] cs_id coupling scheme id + */ +void rebuildContactFormulation( IndexT cs_id ); /*! * diff --git a/src/tribol/physics/ContactFormulation.hpp b/src/tribol/physics/ContactFormulation.hpp index 40392127..c48e37f3 100644 --- a/src/tribol/physics/ContactFormulation.hpp +++ b/src/tribol/physics/ContactFormulation.hpp @@ -118,11 +118,6 @@ class ContactFormulation { */ virtual void updateConstantPenaltyStiffness( double /*mesh1_penalty*/, double /*mesh2_penalty*/ ) {} - /** - * @brief Update ENERGY_MORTAR penalty mode on formulations that support it - */ - virtual void updateEnergyMortarPenaltyMode( EnergyMortarPenaltyMode /*mode*/ ) {} - #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 7f695292..c0212092 100644 --- a/src/tribol/physics/ContactFormulationFactory.cpp +++ b/src/tribol/physics/ContactFormulationFactory.cpp @@ -43,9 +43,15 @@ std::unique_ptr createContactFormulation( CouplingScheme* cs SLIC_ERROR_ROOT_IF( !cs->hasMfemSubmeshData(), "ENERGY_MORTAR requires MFEM submesh data." ); SLIC_ERROR_ROOT_IF( !cs->hasMfemJacobianData(), "ENERGY_MORTAR requires MFEM Jacobian data." ); - return std::make_unique( *cs->getMfemMeshData(), *cs->getMfemSubmeshData(), - *cs->getMfemJacobianData(), k, delta, N, enzyme_quadrature, - use_penalty_, cs->getParameters().energy_mortar_penalty_mode ); + if ( cs->getParameters().energy_mortar_enforcement_option == EnergyMortarEnforcementOption::QUADRATURE_POINT_GAP ) { + return std::make_unique>( + *cs->getMfemMeshData(), *cs->getMfemSubmeshData(), *cs->getMfemJacobianData(), k, delta, N, enzyme_quadrature, + use_penalty_ ); + } else { + return std::make_unique>( + *cs->getMfemMeshData(), *cs->getMfemSubmeshData(), *cs->getMfemJacobianData(), k, delta, N, enzyme_quadrature, + use_penalty_ ); + } #else SLIC_ERROR_ROOT( "ENERGY_MORTAR requires Enzyme and redecomp to be built." ); return nullptr; diff --git a/src/tribol/physics/EnergyMortar.hpp b/src/tribol/physics/EnergyMortar.hpp index 08a5335c..9a98bef5 100644 --- a/src/tribol/physics/EnergyMortar.hpp +++ b/src/tribol/physics/EnergyMortar.hpp @@ -23,7 +23,6 @@ struct ContactParams { double k; // Penalty int N; // Quadrature Points bool enzyme_quadrature; // Determines how enzyming is performed (default = True) - EnergyMortarPenaltyMode penalty_mode{ EnergyMortarPenaltyMode::QUADRATURE_POINT_GAP }; // Penalty enforcement mode }; /// Stores quadrature-point penalty energy derivatives for one interface pair. diff --git a/src/tribol/physics/EnergyMortarAdapter.cpp b/src/tribol/physics/EnergyMortarAdapter.cpp index 28aad634..7885bf6d 100644 --- a/src/tribol/physics/EnergyMortarAdapter.cpp +++ b/src/tribol/physics/EnergyMortarAdapter.cpp @@ -11,10 +11,10 @@ namespace tribol { #ifdef TRIBOL_USE_ENZYME -EnergyMortarAdapter::EnergyMortarAdapter( MfemMeshData& mesh_data, MfemSubmeshData& submesh_data, +template +EnergyMortarAdapter::EnergyMortarAdapter( MfemMeshData& mesh_data, MfemSubmeshData& submesh_data, MfemJacobianData& jac_data, double k, double delta, int N, - bool enzyme_quadrature, bool use_penalty, - EnergyMortarPenaltyMode penalty_mode ) + bool enzyme_quadrature, bool use_penalty ) // NOTE: mesh1 maps to mesh2_ and mesh2 maps to mesh1_. This is to keep consistent with mesh1_ being non-mortar and // mesh2_ being mortar as is typical in the literature, but different from Tribol convention. : use_penalty_( use_penalty ), mesh_data_( mesh_data ), submesh_data_( submesh_data ), jac_data_( jac_data ) @@ -23,8 +23,7 @@ EnergyMortarAdapter::EnergyMortarAdapter( MfemMeshData& mesh_data, MfemSubmeshDa params_.del = delta; params_.N = N; params_.enzyme_quadrature = enzyme_quadrature; - params_.penalty_mode = penalty_mode; - + evaluator_ = std::make_unique( params_ ); // Allocate the (pressure) true-dof vector early so host code can set it via tribol::getMfemContactPressure() after @@ -34,61 +33,61 @@ EnergyMortarAdapter::EnergyMortarAdapter( MfemMeshData& mesh_data, MfemSubmeshDa pressure_vec_.fill( 0.0 ); } -void EnergyMortarAdapter::updateMeshes( MeshData& mesh1, MeshData& mesh2 ) +template +void EnergyMortarAdapter::updateMeshes( MeshData& mesh1, MeshData& mesh2 ) { // Maintain the same "flipped" convention as the constructor. mesh1_ = &mesh2; mesh2_ = &mesh1; } -void EnergyMortarAdapter::updateConstantPenaltyStiffness( double mesh1_penalty, double mesh2_penalty ) +template +void EnergyMortarAdapter::updateConstantPenaltyStiffness( double mesh1_penalty, double mesh2_penalty ) { use_penalty_ = true; params_.k = 0.5 * ( mesh1_penalty + mesh2_penalty ); evaluator_ = std::make_unique( params_ ); } -void EnergyMortarAdapter::updateEnergyMortarPenaltyMode( EnergyMortarPenaltyMode mode ) -{ - params_.penalty_mode = mode; - evaluator_ = std::make_unique( params_ ); -} -const mfem::HypreParVector& EnergyMortarAdapter::getMfemGap() const +template +const mfem::HypreParVector& EnergyMortarAdapter::getMfemGap() const { // Penalty mode uses the normalized gap g = g_tilde / A. LM mode enforces the unnormalized constraint g_tilde = 0, // consistent with dg/dx returned by getMfemDgDx(). return use_penalty_ ? gap_vec_.get() : g_tilde_vec_.get(); } -void EnergyMortarAdapter::setInterfacePairs( ArrayT&& pairs, int /*check_level*/ ) +template +void EnergyMortarAdapter::setInterfacePairs( ArrayT&& pairs, int /*check_level*/ ) { // TODO: Consider design and how this interacts with binning and CG pairs_ = std::move( pairs ); } -void EnergyMortarAdapter::updateIntegrationRule() +template +void EnergyMortarAdapter::updateIntegrationRule() { SLIC_WARNING_ROOT( "Update integration rule not implemmented for any method" ); // TODO: break out integration rule as a separate method } -void EnergyMortarAdapter::updateNodalGaps() +void NodalGapEnforcement::updateNodalGaps(EnergyMortarAdapter* adapter) { // NOTE: user should have called updateMfemParallelDecomposition() with updated coords before calling this // Tribol level data structures for storing gap, area, and derivatives - auto& redecomp_gap = submesh_data_.GetRedecompGap(); + auto& redecomp_gap = adapter->submesh_data_.GetRedecompGap(); mfem::GridFunction redecomp_area( redecomp_gap.FESpace() ); redecomp_area = 0.0; - const bool use_lor = ( mesh_data_.GetLORMesh() != nullptr ); - const auto& displacement_surface_fes = use_lor ? *mesh_data_.GetLORMeshFESpace() : mesh_data_.GetSubmeshFESpace(); - const auto& pressure_surface_fes = use_lor ? *submesh_data_.GetLORMeshFESpace() : submesh_data_.GetSubmeshFESpace(); - const auto& displacement_redecomp_fes = *mesh_data_.GetRedecompResponse().FESpace(); - const auto& pressure_redecomp_fes = *submesh_data_.GetRedecompGap().FESpace(); - const auto& mortar_elem_map = mesh_data_.GetElemMap1(); - const auto& nonmortar_elem_map = mesh_data_.GetElemMap2(); + const bool use_lor = ( adapter->mesh_data_.GetLORMesh() != nullptr ); + const auto& displacement_surface_fes = use_lor ? *adapter->mesh_data_.GetLORMeshFESpace() : adapter->mesh_data_.GetSubmeshFESpace(); + const auto& pressure_surface_fes = use_lor ? *adapter->submesh_data_.GetLORMeshFESpace() : adapter->submesh_data_.GetSubmeshFESpace(); + const auto& displacement_redecomp_fes = *adapter->mesh_data_.GetRedecompResponse().FESpace(); + const auto& pressure_redecomp_fes = *adapter->submesh_data_.GetRedecompGap().FESpace(); + const auto& mortar_elem_map = adapter->mesh_data_.GetElemMap1(); + const auto& nonmortar_elem_map = adapter->mesh_data_.GetElemMap2(); PackedPairJacobianContribs dg_lm_nm( pressure_surface_fes, displacement_surface_fes, pressure_redecomp_fes, displacement_redecomp_fes, nonmortar_elem_map, nonmortar_elem_map ); @@ -99,19 +98,19 @@ void EnergyMortarAdapter::updateNodalGaps() PackedPairJacobianContribs dA_lm_m( pressure_surface_fes, displacement_surface_fes, pressure_redecomp_fes, displacement_redecomp_fes, nonmortar_elem_map, mortar_elem_map ); - dg_lm_nm.reserve( pairs_.size(), 8 ); - dg_lm_m.reserve( pairs_.size(), 8 ); - dA_lm_nm.reserve( pairs_.size(), 8 ); - dA_lm_m.reserve( pairs_.size(), 8 ); + dg_lm_nm.reserve( adapter->pairs_.size(), 8 ); + dg_lm_m.reserve( adapter->pairs_.size(), 8 ); + dA_lm_nm.reserve( adapter->pairs_.size(), 8 ); + dA_lm_m.reserve( adapter->pairs_.size(), 8 ); const int node_idx[8] = { 0, 2, 1, 3, 4, 6, 5, 7 }; - SLIC_ERROR_ROOT_IF( mesh1_ == nullptr || mesh2_ == nullptr, "ENERGY_MORTAR meshes not set." ); - auto mesh1_view = mesh1_->getView(); - auto mesh2_view = mesh2_->getView(); + SLIC_ERROR_ROOT_IF( adapter->mesh1_ == nullptr || adapter->mesh2_ == nullptr, "ENERGY_MORTAR meshes not set." ); + auto mesh1_view = adapter->mesh1_->getView(); + auto mesh2_view = adapter->mesh2_->getView(); // Compute local contributions - for ( const auto& pair : pairs_ ) { + for ( const auto& pair : adapter->pairs_ ) { // These need to be flipped, since the pairs are determined with element 1 associated with mesh 1, and we flipped // the mesh numbers to be consistent with the literature and since the underlying method integrates on element 1 InterfacePair flipped_pair( pair.m_element_id2, pair.m_element_id1 ); @@ -121,7 +120,7 @@ void EnergyMortarAdapter::updateNodalGaps() double g_tilde_elem[2]; double A_elem[2]; - evaluator_->compute_gtilde_and_area( flipped_pair, mesh1_view, mesh2_view, g_tilde_elem, A_elem ); + adapter->evaluator_->compute_gtilde_and_area( flipped_pair, mesh1_view, mesh2_view, g_tilde_elem, A_elem ); if ( A_elem[0] <= 0.0 && A_elem[1] <= 0.0 ) { continue; @@ -140,7 +139,7 @@ void EnergyMortarAdapter::updateNodalGaps() double dg_dx_node1[8]; double dg_dx_node2[8]; // TODO: make grad_gtilde return directly in dg_tilde_dx_blocks format - evaluator_->grad_gtilde( flipped_pair, mesh1_view, mesh2_view, dg_dx_node1, dg_dx_node2 ); + adapter->evaluator_->grad_gtilde( flipped_pair, mesh1_view, mesh2_view, dg_dx_node1, dg_dx_node2 ); double dg_tilde_dx_blocks[2][8]; for ( int i{ 0 }; i < 4; ++i ) { dg_tilde_dx_blocks[0][i * 2] = dg_dx_node1[node_idx[i]]; @@ -154,7 +153,7 @@ void EnergyMortarAdapter::updateNodalGaps() double dA_dx_node1[8]; double dA_dx_node2[8]; // TODO: make grad_trib_area return directly in dA_dx_blocks format - evaluator_->grad_trib_area( flipped_pair, mesh1_view, mesh2_view, dA_dx_node1, dA_dx_node2 ); + adapter->evaluator_->grad_trib_area( flipped_pair, mesh1_view, mesh2_view, dA_dx_node1, dA_dx_node2 ); double dA_dx_blocks[2][8]; for ( int i{ 0 }; i < 4; ++i ) { dA_dx_blocks[0][i * 2] = dA_dx_node1[node_idx[i]]; @@ -168,157 +167,153 @@ void EnergyMortarAdapter::updateNodalGaps() // Move gap and area to submesh level vectors mfem::ParLinearForm g_tilde_linear_form( - const_cast( &submesh_data_.GetSubmeshFESpace() ) ); - submesh_data_.GetSubmeshGap( g_tilde_linear_form ); - auto& P_submesh = *submesh_data_.GetSubmeshFESpace().GetProlongationMatrix(); - g_tilde_vec_ = shared::ParVector( const_cast( &submesh_data_.GetSubmeshFESpace() ) ); - g_tilde_vec_.fill( 0.0 ); - P_submesh.MultTranspose( g_tilde_linear_form, g_tilde_vec_.get() ); + const_cast( &adapter->submesh_data_.GetSubmeshFESpace() ) ); + adapter->submesh_data_.GetSubmeshGap( g_tilde_linear_form ); + auto& P_submesh = *adapter->submesh_data_.GetSubmeshFESpace().GetProlongationMatrix(); + adapter->g_tilde_vec_ = shared::ParVector( const_cast( &adapter->submesh_data_.GetSubmeshFESpace() ) ); + adapter->g_tilde_vec_.fill( 0.0 ); + P_submesh.MultTranspose( g_tilde_linear_form, adapter->g_tilde_vec_.get() ); mfem::Array rows_to_elim; - if ( !tied_contact_ && use_penalty_ ) { - rows_to_elim.Reserve( g_tilde_vec_.size() ); - for ( int i{ 0 }; i < g_tilde_vec_.size(); ++i ) { - if ( g_tilde_vec_[i] > 0.0 ) { - g_tilde_vec_[i] = 0.0; + if ( !adapter->tied_contact_ && adapter->use_penalty_ ) { + rows_to_elim.Reserve( adapter->g_tilde_vec_.size() ); + for ( int i{ 0 }; i < adapter->g_tilde_vec_.size(); ++i ) { + if ( adapter->g_tilde_vec_[i] > 0.0 ) { + adapter->g_tilde_vec_[i] = 0.0; rows_to_elim.push_back( i ); } } } - mfem::ParLinearForm A_linear_form( const_cast( &submesh_data_.GetSubmeshFESpace() ) ); - submesh_data_.GetPressureTransfer().RedecompToSubmesh( redecomp_area, A_linear_form ); - A_vec_ = shared::ParVector( const_cast( &submesh_data_.GetSubmeshFESpace() ) ); - A_vec_.fill( 0.0 ); - P_submesh.MultTranspose( A_linear_form, A_vec_.get() ); + mfem::ParLinearForm A_linear_form( const_cast( &adapter->submesh_data_.GetSubmeshFESpace() ) ); + adapter->submesh_data_.GetPressureTransfer().RedecompToSubmesh( redecomp_area, A_linear_form ); + adapter->A_vec_ = shared::ParVector( const_cast( &adapter->submesh_data_.GetSubmeshFESpace() ) ); + adapter->A_vec_.fill( 0.0 ); + P_submesh.MultTranspose( A_linear_form, adapter->A_vec_.get() ); - gap_vec_ = g_tilde_vec_.divide( A_vec_, area_tol_ ); + adapter->gap_vec_ = adapter->g_tilde_vec_.divide( adapter->A_vec_, adapter->area_tol_ ); // Move gap and area derivatives to (pressure true-dof rows, displacement true-dof cols) std::vector dg_contribs; dg_contribs.reserve( 2 ); dg_contribs.push_back( std::move( dg_lm_nm ) ); dg_contribs.push_back( std::move( dg_lm_m ) ); - dg_tilde_dx_ = jac_data_.GetMfemJacobian( &submesh_data_.GetSubmeshFESpace(), - mesh_data_.GetParentCoords().ParFESpace(), dg_contribs ); - if ( !tied_contact_ && use_penalty_ ) { + adapter->dg_tilde_dx_ = adapter->jac_data_.GetMfemJacobian( &adapter->submesh_data_.GetSubmeshFESpace(), + adapter->mesh_data_.GetParentCoords().ParFESpace(), dg_contribs ); + if ( !adapter->tied_contact_ && adapter->use_penalty_ ) { // technically, we should do this on all the vectors/matrices below, but it looks like the mutliplication operators // below will zero them out anyway - dg_tilde_dx_.eliminateRows( rows_to_elim ); + adapter->dg_tilde_dx_.eliminateRows( rows_to_elim ); } std::vector dA_contribs; dA_contribs.reserve( 2 ); dA_contribs.push_back( std::move( dA_lm_nm ) ); dA_contribs.push_back( std::move( dA_lm_m ) ); - dA_dx_ = jac_data_.GetMfemJacobian( &submesh_data_.GetSubmeshFESpace(), mesh_data_.GetParentCoords().ParFESpace(), + adapter->dA_dx_ = adapter->jac_data_.GetMfemJacobian( &adapter->submesh_data_.GetSubmeshFESpace(), adapter->mesh_data_.GetParentCoords().ParFESpace(), dA_contribs ); } -void EnergyMortarAdapter::updateNodalForces() +void NodalGapEnforcement::updateNodalForces(EnergyMortarAdapter* adapter) { // NOTE: user should have called updateNodalGaps() with updated coords before calling this - if ( use_penalty_ ) { - if ( params_.penalty_mode == EnergyMortarPenaltyMode::QUADRATURE_POINT_GAP ) { - updateQuadraturePointPenaltyForces(); - return; - } - + if ( adapter->use_penalty_ ) { // Penalty mode: p = k * (g_tilde / A) - pressure_vec_ = params_.k * gap_vec_; + adapter->pressure_vec_ = adapter->params_.k * adapter->gap_vec_; } else { - // LM mode: pressure_vec_ is treated as the Lagrange multiplier vector (lambda) - SLIC_ERROR_ROOT_IF( submesh_data_.GetSubmeshFESpace().GetTrueVSize() != pressure_vec_.size(), + // LM mode: adapter->pressure_vec_ is treated as the Lagrange multiplier vector (lambda) + SLIC_ERROR_ROOT_IF( adapter->submesh_data_.GetSubmeshFESpace().GetTrueVSize() != adapter->pressure_vec_.size(), "LM vector is not initialized. Call tribol::update() once to initialize the formulation." ); - SLIC_ERROR_ROOT_IF( pressure_vec_.size() != g_tilde_vec_.size(), + SLIC_ERROR_ROOT_IF( adapter->pressure_vec_.size() != adapter->g_tilde_vec_.size(), "LM vector size mismatch with contact dofs (g_tilde)." ); } - energy_ = pressure_vec_.dot( g_tilde_vec_ ); + adapter->energy_ = adapter->pressure_vec_.dot( adapter->g_tilde_vec_ ); - if ( !use_penalty_ ) { + if ( !adapter->use_penalty_ ) { // ------------------------------------------------------------------------- // LM mode: force = G^T * lambda and df/dx = lambda · d^2(g_tilde)/dx^2 // ------------------------------------------------------------------------- - force_vec_ = pressure_vec_ * dg_tilde_dx_; + adapter->force_vec_ = adapter->pressure_vec_ * adapter->dg_tilde_dx_; - mfem::GridFunction redecomp_lambda( submesh_data_.GetRedecompGap() ); + mfem::GridFunction redecomp_lambda( adapter->submesh_data_.GetRedecompGap() ); mfem::ParGridFunction submesh_lambda( - const_cast( &submesh_data_.GetSubmeshFESpace() ) ); - submesh_lambda.SetFromTrueDofs( pressure_vec_.get() ); - submesh_data_.GetPressureTransfer().SubmeshToRedecomp( submesh_lambda, redecomp_lambda ); + const_cast( &adapter->submesh_data_.GetSubmeshFESpace() ) ); + submesh_lambda.SetFromTrueDofs( adapter->pressure_vec_.get() ); + adapter->submesh_data_.GetPressureTransfer().SubmeshToRedecomp( submesh_lambda, redecomp_lambda ); - df_dx_ = computeDfDxSecondDerivativesLM( redecomp_lambda ); + adapter->df_dx_ = computeDfDxSecondDerivativesLM(adapter, redecomp_lambda ); return; } // --------------------------------------------------------------------------- // Penalty mode: force and Jacobian include pressure/area coupling terms // --------------------------------------------------------------------------- - auto k_over_a = params_.k * A_vec_.inverse( area_tol_ ); - auto p_over_a = pressure_vec_.divide( A_vec_, area_tol_ ); + auto k_over_a = adapter->params_.k * adapter->A_vec_.inverse( adapter->area_tol_ ); + auto p_over_a = adapter->pressure_vec_.divide( adapter->A_vec_, adapter->area_tol_ ); - shared::ParSparseMat dp_dx( dg_tilde_dx_.get() ); + shared::ParSparseMat dp_dx( adapter->dg_tilde_dx_.get() ); dp_dx->ScaleRows( k_over_a.get() ); - shared::ParSparseMat dp_dx_temp( dA_dx_.get() ); + shared::ParSparseMat dp_dx_temp( adapter->dA_dx_.get() ); dp_dx_temp->ScaleRows( p_over_a.get() ); dp_dx -= dp_dx_temp; - force_vec_ = ( pressure_vec_ * dg_tilde_dx_ ) + ( g_tilde_vec_ * dp_dx ); + adapter->force_vec_ = ( adapter->pressure_vec_ * adapter->dg_tilde_dx_ ) + ( adapter->g_tilde_vec_ * dp_dx ); // TODO (EBC): Move transfer path-specific logic out of this file - mfem::GridFunction redecomp_pressure( submesh_data_.GetRedecompGap() ); + mfem::GridFunction redecomp_pressure( adapter->submesh_data_.GetRedecompGap() ); mfem::ParGridFunction submesh_pressure( - const_cast( &submesh_data_.GetSubmeshFESpace() ) ); - submesh_pressure.SetFromTrueDofs( pressure_vec_.get() ); - submesh_data_.GetPressureTransfer().SubmeshToRedecomp( submesh_pressure, redecomp_pressure ); + const_cast( &adapter->submesh_data_.GetSubmeshFESpace() ) ); + submesh_pressure.SetFromTrueDofs( adapter->pressure_vec_.get() ); + adapter->submesh_data_.GetPressureTransfer().SubmeshToRedecomp( submesh_pressure, redecomp_pressure ); - mfem::GridFunction redecomp_g_tilde( submesh_data_.GetRedecompGap() ); + mfem::GridFunction redecomp_g_tilde( adapter->submesh_data_.GetRedecompGap() ); mfem::ParGridFunction submesh_g_tilde( - const_cast( &submesh_data_.GetSubmeshFESpace() ) ); - submesh_g_tilde.SetFromTrueDofs( g_tilde_vec_.get() ); - submesh_data_.GetPressureTransfer().SubmeshToRedecomp( submesh_g_tilde, redecomp_g_tilde ); + const_cast( &adapter->submesh_data_.GetSubmeshFESpace() ) ); + submesh_g_tilde.SetFromTrueDofs( adapter->g_tilde_vec_.get() ); + adapter->submesh_data_.GetPressureTransfer().SubmeshToRedecomp( submesh_g_tilde, redecomp_g_tilde ); - mfem::GridFunction redecomp_A( submesh_data_.GetRedecompGap() ); - mfem::ParGridFunction submesh_A( const_cast( &submesh_data_.GetSubmeshFESpace() ) ); - submesh_A.SetFromTrueDofs( A_vec_.get() ); - submesh_data_.GetPressureTransfer().SubmeshToRedecomp( submesh_A, redecomp_A ); + mfem::GridFunction redecomp_A( adapter->submesh_data_.GetRedecompGap() ); + mfem::ParGridFunction submesh_A( const_cast( &adapter->submesh_data_.GetSubmeshFESpace() ) ); + submesh_A.SetFromTrueDofs( adapter->A_vec_.get() ); + adapter->submesh_data_.GetPressureTransfer().SubmeshToRedecomp( submesh_A, redecomp_A ); - df_dx_ = computeDfDxSecondDerivativesPenalty( redecomp_pressure, redecomp_g_tilde, redecomp_A ); + adapter->df_dx_ = computeDfDxSecondDerivativesPenalty(adapter, redecomp_pressure, redecomp_g_tilde, redecomp_A ); - auto pg2_over_asq = ( 2.0 * pressure_vec_ ) - .multiplyInPlace( g_tilde_vec_ ) - .divideInPlace( A_vec_, area_tol_ ) - .divideInPlace( A_vec_, area_tol_ ); + auto pg2_over_asq = ( 2.0 * adapter->pressure_vec_ ) + .multiplyInPlace( adapter->g_tilde_vec_ ) + .divideInPlace( adapter->A_vec_, adapter->area_tol_ ) + .divideInPlace( adapter->A_vec_, adapter->area_tol_ ); - auto& submesh_fes = submesh_data_.GetSubmeshFESpace(); + auto& submesh_fes = adapter->submesh_data_.GetSubmeshFESpace(); auto p_over_a_diag = shared::ParSparseMat::diagonalMatrix( submesh_fes.GetComm(), submesh_fes.GlobalTrueVSize(), submesh_fes.GetTrueDofOffsets(), p_over_a.get() ); auto pg2_over_asq_diag = shared::ParSparseMat::diagonalMatrix( submesh_fes.GetComm(), submesh_fes.GlobalTrueVSize(), submesh_fes.GetTrueDofOffsets(), pg2_over_asq.get() ); - df_dx_ -= shared::ParSparseMat::rap( dg_tilde_dx_, p_over_a_diag, dA_dx_ ); - df_dx_ -= shared::ParSparseMat::rap( dA_dx_, p_over_a_diag, dg_tilde_dx_ ); - df_dx_ += shared::ParSparseMat::rap( dA_dx_, pg2_over_asq_diag, dg_tilde_dx_ ); - df_dx_ += dp_dx.transpose() * dg_tilde_dx_; - df_dx_ += dg_tilde_dx_.transpose() * dp_dx; + adapter->df_dx_ -= shared::ParSparseMat::rap( adapter->dg_tilde_dx_, p_over_a_diag, adapter->dA_dx_ ); + adapter->df_dx_ -= shared::ParSparseMat::rap( adapter->dA_dx_, p_over_a_diag, adapter->dg_tilde_dx_ ); + adapter->df_dx_ += shared::ParSparseMat::rap( adapter->dA_dx_, pg2_over_asq_diag, adapter->dg_tilde_dx_ ); + adapter->df_dx_ += dp_dx.transpose() * adapter->dg_tilde_dx_; + adapter->df_dx_ += adapter->dg_tilde_dx_.transpose() * dp_dx; } -RealT EnergyMortarAdapter::computeTimeStep() +template +RealT EnergyMortarAdapter::computeTimeStep() { SLIC_INFO_ROOT( "computeTimestep() not implemented for EnergyMortar" ); // TODO: implement timestep calculation return 1.0; } -shared::ParSparseMat EnergyMortarAdapter::computeDfDxSecondDerivativesLM( const mfem::GridFunction& redecomp_lambda ) +shared::ParSparseMat NodalGapEnforcement::computeDfDxSecondDerivativesLM(EnergyMortarAdapter* adapter, const mfem::GridFunction& redecomp_lambda ) { - const bool use_lor = ( mesh_data_.GetLORMesh() != nullptr ); - const auto& displacement_surface_fes = use_lor ? *mesh_data_.GetLORMeshFESpace() : mesh_data_.GetSubmeshFESpace(); - const auto& displacement_redecomp_fes = *mesh_data_.GetRedecompResponse().FESpace(); - const auto& mortar_elem_map = mesh_data_.GetElemMap1(); - const auto& nonmortar_elem_map = mesh_data_.GetElemMap2(); + const bool use_lor = ( adapter->mesh_data_.GetLORMesh() != nullptr ); + const auto& displacement_surface_fes = use_lor ? *adapter->mesh_data_.GetLORMeshFESpace() : adapter->mesh_data_.GetSubmeshFESpace(); + const auto& displacement_redecomp_fes = *adapter->mesh_data_.GetRedecompResponse().FESpace(); + const auto& mortar_elem_map = adapter->mesh_data_.GetElemMap1(); + const auto& nonmortar_elem_map = adapter->mesh_data_.GetElemMap2(); PackedPairJacobianContribs df_nm_nm( displacement_surface_fes, displacement_surface_fes, displacement_redecomp_fes, displacement_redecomp_fes, nonmortar_elem_map, nonmortar_elem_map ); @@ -329,18 +324,18 @@ shared::ParSparseMat EnergyMortarAdapter::computeDfDxSecondDerivativesLM( const PackedPairJacobianContribs df_m_m( displacement_surface_fes, displacement_surface_fes, displacement_redecomp_fes, displacement_redecomp_fes, mortar_elem_map, mortar_elem_map ); - df_nm_nm.reserve( pairs_.size(), 16 ); - df_nm_m.reserve( pairs_.size(), 16 ); - df_m_nm.reserve( pairs_.size(), 16 ); - df_m_m.reserve( pairs_.size(), 16 ); + df_nm_nm.reserve( adapter->pairs_.size(), 16 ); + df_nm_m.reserve( adapter->pairs_.size(), 16 ); + df_m_nm.reserve( adapter->pairs_.size(), 16 ); + df_m_m.reserve( adapter->pairs_.size(), 16 ); const int node_idx[8] = { 0, 2, 1, 3, 4, 6, 5, 7 }; - SLIC_ERROR_ROOT_IF( mesh1_ == nullptr || mesh2_ == nullptr, "ENERGY_MORTAR meshes not set." ); - auto mesh1_view = mesh1_->getView(); - auto mesh2_view = mesh2_->getView(); + SLIC_ERROR_ROOT_IF( adapter->mesh1_ == nullptr || adapter->mesh2_ == nullptr, "ENERGY_MORTAR meshes not set." ); + auto mesh1_view = adapter->mesh1_->getView(); + auto mesh2_view = adapter->mesh2_->getView(); - for ( auto& pair : pairs_ ) { + for ( auto& pair : adapter->pairs_ ) { InterfacePair flipped_pair( pair.m_element_id2, pair.m_element_id1 ); const auto elem1 = static_cast( flipped_pair.m_element_id1 ); const auto node11 = mesh1_view.getConnectivity()( elem1, 0 ); @@ -352,7 +347,7 @@ shared::ParSparseMat EnergyMortarAdapter::computeDfDxSecondDerivativesLM( const double d2g_dx2_node1[64]; double d2g_dx2_node2[64]; - evaluator_->d2_g2tilde( flipped_pair, mesh1_view, mesh2_view, d2g_dx2_node1, d2g_dx2_node2 ); + adapter->evaluator_->d2_g2tilde( flipped_pair, mesh1_view, mesh2_view, d2g_dx2_node1, d2g_dx2_node2 ); double df_dx_blocks[2][2][16]; for ( int i{ 0 }; i < 2; ++i ) { @@ -378,19 +373,19 @@ shared::ParSparseMat EnergyMortarAdapter::computeDfDxSecondDerivativesLM( const df_contribs.push_back( std::move( df_nm_m ) ); df_contribs.push_back( std::move( df_m_nm ) ); df_contribs.push_back( std::move( df_m_m ) ); - return jac_data_.GetMfemJacobian( mesh_data_.GetParentCoords().ParFESpace(), - mesh_data_.GetParentCoords().ParFESpace(), df_contribs ); + return adapter->jac_data_.GetMfemJacobian( adapter->mesh_data_.GetParentCoords().ParFESpace(), + adapter->mesh_data_.GetParentCoords().ParFESpace(), df_contribs ); } -shared::ParSparseMat EnergyMortarAdapter::computeDfDxSecondDerivativesPenalty( +shared::ParSparseMat NodalGapEnforcement::computeDfDxSecondDerivativesPenalty(EnergyMortarAdapter* adapter, const mfem::GridFunction& redecomp_pressure, const mfem::GridFunction& redecomp_g_tilde, const mfem::GridFunction& redecomp_A ) { - const bool use_lor = ( mesh_data_.GetLORMesh() != nullptr ); - const auto& displacement_surface_fes = use_lor ? *mesh_data_.GetLORMeshFESpace() : mesh_data_.GetSubmeshFESpace(); - const auto& displacement_redecomp_fes = *mesh_data_.GetRedecompResponse().FESpace(); - const auto& mortar_elem_map = mesh_data_.GetElemMap1(); - const auto& nonmortar_elem_map = mesh_data_.GetElemMap2(); + const bool use_lor = ( adapter->mesh_data_.GetLORMesh() != nullptr ); + const auto& displacement_surface_fes = use_lor ? *adapter->mesh_data_.GetLORMeshFESpace() : adapter->mesh_data_.GetSubmeshFESpace(); + const auto& displacement_redecomp_fes = *adapter->mesh_data_.GetRedecompResponse().FESpace(); + const auto& mortar_elem_map = adapter->mesh_data_.GetElemMap1(); + const auto& nonmortar_elem_map = adapter->mesh_data_.GetElemMap2(); PackedPairJacobianContribs df_nm_nm( displacement_surface_fes, displacement_surface_fes, displacement_redecomp_fes, displacement_redecomp_fes, nonmortar_elem_map, nonmortar_elem_map ); @@ -401,18 +396,18 @@ shared::ParSparseMat EnergyMortarAdapter::computeDfDxSecondDerivativesPenalty( PackedPairJacobianContribs df_m_m( displacement_surface_fes, displacement_surface_fes, displacement_redecomp_fes, displacement_redecomp_fes, mortar_elem_map, mortar_elem_map ); - df_nm_nm.reserve( pairs_.size(), 16 ); - df_nm_m.reserve( pairs_.size(), 16 ); - df_m_nm.reserve( pairs_.size(), 16 ); - df_m_m.reserve( pairs_.size(), 16 ); + df_nm_nm.reserve( adapter->pairs_.size(), 16 ); + df_nm_m.reserve( adapter->pairs_.size(), 16 ); + df_m_nm.reserve( adapter->pairs_.size(), 16 ); + df_m_m.reserve( adapter->pairs_.size(), 16 ); const int node_idx[8] = { 0, 2, 1, 3, 4, 6, 5, 7 }; - SLIC_ERROR_ROOT_IF( mesh1_ == nullptr || mesh2_ == nullptr, "ENERGY_MORTAR meshes not set." ); - auto mesh1_view = mesh1_->getView(); - auto mesh2_view = mesh2_->getView(); + SLIC_ERROR_ROOT_IF( adapter->mesh1_ == nullptr || adapter->mesh2_ == nullptr, "ENERGY_MORTAR meshes not set." ); + auto mesh1_view = adapter->mesh1_->getView(); + auto mesh2_view = adapter->mesh2_->getView(); - for ( auto& pair : pairs_ ) { + for ( auto& pair : adapter->pairs_ ) { InterfacePair flipped_pair( pair.m_element_id2, pair.m_element_id1 ); const auto elem1 = static_cast( flipped_pair.m_element_id1 ); const auto node11 = mesh1_view.getConnectivity()( elem1, 0 ); @@ -431,11 +426,11 @@ shared::ParSparseMat EnergyMortarAdapter::computeDfDxSecondDerivativesPenalty( double d2g_dx2_node1[64]; double d2g_dx2_node2[64]; - evaluator_->d2_g2tilde( flipped_pair, mesh1_view, mesh2_view, d2g_dx2_node1, d2g_dx2_node2 ); + adapter->evaluator_->d2_g2tilde( flipped_pair, mesh1_view, mesh2_view, d2g_dx2_node1, d2g_dx2_node2 ); double d2A_dx2_node1[64]; double d2A_dx2_node2[64]; - evaluator_->compute_d2A_d2u( flipped_pair, mesh1_view, mesh2_view, d2A_dx2_node1, d2A_dx2_node2 ); + adapter->evaluator_->compute_d2A_d2u( flipped_pair, mesh1_view, mesh2_view, d2A_dx2_node1, d2A_dx2_node2 ); double df_dx_blocks[2][2][16]; for ( int i{ 0 }; i < 2; ++i ) { @@ -462,17 +457,17 @@ shared::ParSparseMat EnergyMortarAdapter::computeDfDxSecondDerivativesPenalty( df_contribs.push_back( std::move( df_nm_m ) ); df_contribs.push_back( std::move( df_m_nm ) ); df_contribs.push_back( std::move( df_m_m ) ); - return jac_data_.GetMfemJacobian( mesh_data_.GetParentCoords().ParFESpace(), - mesh_data_.GetParentCoords().ParFESpace(), df_contribs ); + return adapter->jac_data_.GetMfemJacobian( adapter->mesh_data_.GetParentCoords().ParFESpace(), + adapter->mesh_data_.GetParentCoords().ParFESpace(), df_contribs ); } -void EnergyMortarAdapter::updateQuadraturePointPenaltyForces() +void QuadraturePointEnforcement::updateNodalForces(EnergyMortarAdapter* adapter) { - const bool use_lor = ( mesh_data_.GetLORMesh() != nullptr ); - const auto& displacement_surface_fes = use_lor ? *mesh_data_.GetLORMeshFESpace() : mesh_data_.GetSubmeshFESpace(); - const auto& displacement_redecomp_fes = *mesh_data_.GetRedecompResponse().FESpace(); - const auto& mortar_elem_map = mesh_data_.GetElemMap1(); - const auto& nonmortar_elem_map = mesh_data_.GetElemMap2(); + const bool use_lor = ( adapter->mesh_data_.GetLORMesh() != nullptr ); + const auto& displacement_surface_fes = use_lor ? *adapter->mesh_data_.GetLORMeshFESpace() : adapter->mesh_data_.GetSubmeshFESpace(); + const auto& displacement_redecomp_fes = *adapter->mesh_data_.GetRedecompResponse().FESpace(); + const auto& mortar_elem_map = adapter->mesh_data_.GetElemMap1(); + const auto& nonmortar_elem_map = adapter->mesh_data_.GetElemMap2(); PackedPairJacobianContribs df_nm_nm( displacement_surface_fes, displacement_surface_fes, displacement_redecomp_fes, displacement_redecomp_fes, nonmortar_elem_map, nonmortar_elem_map ); @@ -483,33 +478,33 @@ void EnergyMortarAdapter::updateQuadraturePointPenaltyForces() PackedPairJacobianContribs df_m_m( displacement_surface_fes, displacement_surface_fes, displacement_redecomp_fes, displacement_redecomp_fes, mortar_elem_map, mortar_elem_map ); - df_nm_nm.reserve( pairs_.size(), 16 ); - df_nm_m.reserve( pairs_.size(), 16 ); - df_m_nm.reserve( pairs_.size(), 16 ); - df_m_m.reserve( pairs_.size(), 16 ); + df_nm_nm.reserve( adapter->pairs_.size(), 16 ); + df_nm_m.reserve( adapter->pairs_.size(), 16 ); + df_m_nm.reserve( adapter->pairs_.size(), 16 ); + df_m_m.reserve( adapter->pairs_.size(), 16 ); mfem::GridFunction redecomp_force( const_cast( &displacement_redecomp_fes ) ); redecomp_force = 0.0; const int scalar_size = redecomp_force.FESpace()->GetVSize() / redecomp_force.FESpace()->GetVDim(); - energy_ = 0.0; + adapter->energy_ = 0.0; const int node_idx[8] = { 0, 2, 1, 3, 4, 6, 5, 7 }; - SLIC_ERROR_ROOT_IF( mesh1_ == nullptr || mesh2_ == nullptr, "ENERGY_MORTAR meshes not set." ); - auto mesh1_view = mesh1_->getView(); - auto mesh2_view = mesh2_->getView(); + SLIC_ERROR_ROOT_IF( adapter->mesh1_ == nullptr || adapter->mesh2_ == nullptr, "ENERGY_MORTAR meshes not set." ); + auto mesh1_view = adapter->mesh1_->getView(); + auto mesh2_view = adapter->mesh2_->getView(); - for ( const auto& pair : pairs_ ) { + for ( const auto& pair : adapter->pairs_ ) { InterfacePair flipped_pair( pair.m_element_id2, pair.m_element_id1 ); const auto elem1 = static_cast( flipped_pair.m_element_id1 ); const auto elem2 = static_cast( flipped_pair.m_element_id2 ); - const auto qp_data = evaluator_->compute_quadrature_point_penalty_data( flipped_pair, mesh1_view, mesh2_view ); + const auto qp_data = adapter->evaluator_->compute_quadrature_point_penalty_data( flipped_pair, mesh1_view, mesh2_view ); if ( qp_data.energy == 0.0 ) { continue; } - energy_ += qp_data.energy; + adapter->energy_ += qp_data.energy; auto A_conn = mesh1_view.getConnectivity()( elem1 ); auto B_conn = mesh2_view.getConnectivity()( elem2 ); @@ -540,13 +535,13 @@ void EnergyMortarAdapter::updateQuadraturePointPenaltyForces() df_m_m.append( elem2, elem2, df_dx_blocks[1][1], 16 ); } - auto* parent_fes = mesh_data_.GetParentCoords().ParFESpace(); - force_vec_ = shared::ParVector( const_cast( parent_fes ) ); - force_vec_.fill( 0.0 ); + auto* parent_fes = adapter->mesh_data_.GetParentCoords().ParFESpace(); + adapter->force_vec_ = shared::ParVector( const_cast( parent_fes ) ); + adapter->force_vec_.fill( 0.0 ); mfem::Vector parent_force( parent_fes->GetVSize() ); parent_force = 0.0; - mesh_data_.GetParentRedecompTransfer().RedecompToParent( redecomp_force, parent_force ); - parent_fes->GetProlongationMatrix()->MultTranspose( parent_force, force_vec_.get() ); + adapter->mesh_data_.GetParentRedecompTransfer().RedecompToParent( redecomp_force, parent_force ); + parent_fes->GetProlongationMatrix()->MultTranspose( parent_force, adapter->force_vec_.get() ); std::vector df_contribs; df_contribs.reserve( 4 ); @@ -554,20 +549,23 @@ void EnergyMortarAdapter::updateQuadraturePointPenaltyForces() df_contribs.push_back( std::move( df_nm_m ) ); df_contribs.push_back( std::move( df_m_nm ) ); df_contribs.push_back( std::move( df_m_m ) ); - df_dx_ = jac_data_.GetMfemJacobian( parent_fes, parent_fes, df_contribs ); + adapter->df_dx_ = adapter->jac_data_.GetMfemJacobian( parent_fes, parent_fes, df_contribs ); } -std::unique_ptr EnergyMortarAdapter::getMfemDfDx() const +template +std::unique_ptr EnergyMortarAdapter::getMfemDfDx() const { return std::unique_ptr( df_dx_.release() ); } -std::unique_ptr EnergyMortarAdapter::getMfemDgDx() const +template +std::unique_ptr EnergyMortarAdapter::getMfemDgDx() const { return std::unique_ptr( dg_tilde_dx_.release() ); } -std::unique_ptr EnergyMortarAdapter::getMfemDfDp() const +template +std::unique_ptr EnergyMortarAdapter::getMfemDfDp() const { if ( use_penalty_ ) { return nullptr; @@ -579,6 +577,10 @@ std::unique_ptr EnergyMortarAdapter::getMfemDfDp() const return std::unique_ptr( df_dlambda.release() ); } + +template class EnergyMortarAdapter; +template class EnergyMortarAdapter; + #endif // TRIBOL_USE_ENZYME } // namespace tribol diff --git a/src/tribol/physics/EnergyMortarAdapter.hpp b/src/tribol/physics/EnergyMortarAdapter.hpp index d6f57e3a..08032efb 100644 --- a/src/tribol/physics/EnergyMortarAdapter.hpp +++ b/src/tribol/physics/EnergyMortarAdapter.hpp @@ -30,6 +30,27 @@ namespace tribol { * the dual field is interpreted as the multiplier vector `lambda` and the formulation computes `f = G^T * lambda` where * `G = d(g_tilde)/dx`. */ + +template +class EnergyMortarAdapter; + +struct NodalGapEnforcement { + static void updateNodalGaps(EnergyMortarAdapter* adapter); + static void updateNodalForces(EnergyMortarAdapter* adapter); + + static shared::ParSparseMat computeDfDxSecondDerivativesLM(EnergyMortarAdapter* adapter, const mfem::GridFunction& redecomp_lambda); + static shared::ParSparseMat computeDfDxSecondDerivativesPenalty(EnergyMortarAdapter* adapter, + const mfem::GridFunction& redecomp_pressure, + const mfem::GridFunction& redecomp_g_tilde, + const mfem::GridFunction& redecomp_A); +}; + +struct QuadraturePointEnforcement { + static void updateNodalGaps(EnergyMortarAdapter* /*adapter*/) {} + static void updateNodalForces(EnergyMortarAdapter* adapter); +}; + +template class EnergyMortarAdapter : public ContactFormulation { public: /** @@ -50,8 +71,7 @@ class EnergyMortarAdapter : public ContactFormulation { * relative to the order of the meshes provided here. */ EnergyMortarAdapter( MfemMeshData& mesh_data, MfemSubmeshData& submesh_data, MfemJacobianData& jac_data, double k, - double delta, int N, bool enzyme_quadrature, bool use_penalty = true, - EnergyMortarPenaltyMode penalty_mode = EnergyMortarPenaltyMode::QUADRATURE_POINT_GAP ); + double delta, int N, bool enzyme_quadrature, bool use_penalty = true ); /** * @brief Default destructor @@ -80,7 +100,7 @@ class EnergyMortarAdapter : public ContactFormulation { * In penalty mode this assembles both `g_tilde` and `g = g_tilde / A`. In LM mode the constraint is `g_tilde = 0` and * derivatives are taken with respect to `g_tilde`. */ - void updateNodalGaps() override; + void updateNodalGaps() override { EnforcementPolicy::updateNodalGaps( this ); } /** * @brief Assemble nodal forces/residual and Jacobian contributions @@ -88,7 +108,7 @@ class EnergyMortarAdapter : public ContactFormulation { * In penalty mode this computes pressure from the current gap and penalty stiffness. In LM mode this interprets the * stored pressure vector as the Lagrange multiplier vector `lambda`. */ - void updateNodalForces() override; + void updateNodalForces() override { EnforcementPolicy::updateNodalForces( this ); } /** * @brief Reports if formulation has a maximum allowable timestep calculation @@ -133,13 +153,6 @@ class EnergyMortarAdapter : public ContactFormulation { */ void updateConstantPenaltyStiffness( double mesh1_penalty, double mesh2_penalty ) override; - /** - * @brief Update the EnergyMortar penalty mode - * - * @param mode Penalty mode - */ - void updateEnergyMortarPenaltyMode( EnergyMortarPenaltyMode mode ) override; - #ifdef BUILD_REDECOMP /** * @brief Return the parent true-dof force vector @@ -201,6 +214,8 @@ class EnergyMortarAdapter : public ContactFormulation { #endif private: + friend EnforcementPolicy; + /** * @brief Controls penalty vs. Lagrange multiplier (LM) mode * @@ -316,36 +331,6 @@ class EnergyMortarAdapter : public ContactFormulation { * @brief Derivative df/dx assembled on parent displacement true-dofs */ mutable shared::ParSparseMat df_dx_; - - /** - * @brief Assemble the LM second-derivative df/dx contribution - * - * Assembles `df/dx = lambda · d²(g_tilde)/dx²`. - * - * @param redecomp_lambda Lagrange multiplier on the redecomp mesh - * @return Assembled df/dx contribution on parent true-dofs - */ - shared::ParSparseMat computeDfDxSecondDerivativesLM( const mfem::GridFunction& redecomp_lambda ); - - /** - * @brief Assemble the penalty second-derivative df/dx contribution - * - * Assembles the second-derivative penalty contribution: - * `p · d²(g_tilde)/dx² - (g_tilde p / A) · d²A/dx²`. - * - * @param redecomp_pressure Pressure field on the redecomp mesh - * @param redecomp_g_tilde g_tilde on the redecomp mesh - * @param redecomp_A Area weighting A on the redecomp mesh - * @return Assembled df/dx contribution on parent true-dofs - */ - shared::ParSparseMat computeDfDxSecondDerivativesPenalty( const mfem::GridFunction& redecomp_pressure, - const mfem::GridFunction& redecomp_g_tilde, - const mfem::GridFunction& redecomp_A ); - - /** - * @brief Assemble QP-gap penalty force and Jacobian contributions - */ - void updateQuadraturePointPenaltyForces(); }; #endif // TRIBOL_USE_ENZYME From 976efdcd33205234ca140d38bba9a7b5df869515 Mon Sep 17 00:00:00 2001 From: "E. B. Chin" Date: Wed, 5 Aug 2026 16:22:19 -0700 Subject: [PATCH 05/14] update naming and clean up methods --- src/tests/tribol_energy_mortar_patch.cpp | 2 +- src/tests/tribol_energy_mortar_patch_LM.cpp | 2 +- src/tests/tribol_mfem_jacobian.cpp | 2 +- src/tribol/common/Parameters.hpp | 6 +++--- src/tribol/interface/mfem_tribol.cpp | 4 ++-- src/tribol/interface/tribol.cpp | 21 +------------------ src/tribol/interface/tribol.hpp | 11 +--------- src/tribol/mesh/CouplingScheme.cpp | 10 +++++++++ src/tribol/mesh/CouplingScheme.hpp | 5 +++++ .../physics/ContactFormulationFactory.cpp | 2 +- 10 files changed, 26 insertions(+), 39 deletions(-) diff --git a/src/tests/tribol_energy_mortar_patch.cpp b/src/tests/tribol_energy_mortar_patch.cpp index 755d1b9b..f5061200 100644 --- a/src/tests/tribol_energy_mortar_patch.cpp +++ b/src/tests/tribol_energy_mortar_patch.cpp @@ -175,7 +175,7 @@ class MfemMortarEnergyPatchTest : public testing::TestWithParam> tribol::registerMfemCouplingScheme( cs_id, mesh1_id, mesh2_id, mesh, coords, mortar_attrs, nonmortar_attrs, tribol::SURFACE_TO_SURFACE, tribol::NO_SLIDING, tribol::ENERGY_MORTAR, tribol::FRICTIONLESS, tribol::PENALTY, tribol::BINNING_GRID ); - tribol::setEnergyMortarEnforcementOption( cs_id, tribol::EnergyMortarEnforcementOption::QUADRATURE_POINT_GAP ); + tribol::setEnergyMortarEnforcementOption( cs_id, tribol::EnergyMortarEnforcementOption::QuadraturePointGap ); tribol::setMfemKinematicConstantPenalty( cs_id, 100.0, 100.0 ); mfem::Vector X( par_fe_space.GetTrueVSize() ); diff --git a/src/tests/tribol_energy_mortar_patch_LM.cpp b/src/tests/tribol_energy_mortar_patch_LM.cpp index d1f85ded..caf3ae8c 100644 --- a/src/tests/tribol_energy_mortar_patch_LM.cpp +++ b/src/tests/tribol_energy_mortar_patch_LM.cpp @@ -195,7 +195,7 @@ class MfemMortarEnergyLagrangePatchTest : public testing::TestWithParam #include -#include "tribol/physics/ContactFormulationFactory.hpp" #include #include #include @@ -185,25 +184,7 @@ void setEnergyMortarEnforcementOption( IndexT cs_id, EnergyMortarEnforcementOpti cs->getParameters().energy_mortar_enforcement_option = mode; // Automatically rebuild the formulation to reflect the new setting - if ( cs->getContactMethod() == ENERGY_MORTAR ) { - rebuildContactFormulation( cs_id ); - } -} - -//------------------------------------------------------------------------------ -void rebuildContactFormulation( IndexT cs_id ) -{ - auto cs = CouplingSchemeManager::getInstance().findData( cs_id ); - - SLIC_ERROR_ROOT_IF( !cs, "tribol::rebuildContactFormulation(): call tribol::registerCouplingScheme() " - << "prior to calling this routine." ); - - if ( cs->getContactMethod() == ENERGY_MORTAR ) { - cs->setContactFormulation( createContactFormulation( cs ) ); - } else { - SLIC_WARNING_ROOT( - "tribol::rebuildContactFormulation(): rebuilding is only supported for ENERGY_MORTAR at this time." ); - } + cs->updateContactFormulation(); } //------------------------------------------------------------------------------ diff --git a/src/tribol/interface/tribol.hpp b/src/tribol/interface/tribol.hpp index 691dfddc..5cbe6c1b 100644 --- a/src/tribol/interface/tribol.hpp +++ b/src/tribol/interface/tribol.hpp @@ -139,16 +139,7 @@ void setTimestepPenFrac( IndexT cs_id, RealT frac ); */ void setEnergyMortarEnforcementOption( IndexT cs_id, EnergyMortarEnforcementOption mode ); -/*! - * \brief Rebuilds the contact formulation for a coupling scheme. - * - * This function can be called after modifying coupling scheme settings - * (e.g. setEnergyMortarEnforcementOption) to recreate the underlying formulation - * instance with the new settings. - * - * \param [in] cs_id coupling scheme id - */ -void rebuildContactFormulation( IndexT cs_id ); + /*! * diff --git a/src/tribol/mesh/CouplingScheme.cpp b/src/tribol/mesh/CouplingScheme.cpp index 5c23c4e8..74374db8 100644 --- a/src/tribol/mesh/CouplingScheme.cpp +++ b/src/tribol/mesh/CouplingScheme.cpp @@ -1155,6 +1155,16 @@ int CouplingScheme::apply( int cycle, RealT t, RealT& dt ) } // end CouplingScheme::apply() +//------------------------------------------------------------------------------ +void CouplingScheme::updateContactFormulation() +{ + if ( m_contactMethod == ENERGY_MORTAR ) { + this->setContactFormulation( createContactFormulation( this ) ); + } else { + SLIC_WARNING_ROOT( "tribol::CouplingScheme::updateContactFormulation(): rebuilding is only supported for ENERGY_MORTAR at this time." ); + } +} + //------------------------------------------------------------------------------ bool CouplingScheme::init() { diff --git a/src/tribol/mesh/CouplingScheme.hpp b/src/tribol/mesh/CouplingScheme.hpp index e7dbe61a..1d6bb7da 100644 --- a/src/tribol/mesh/CouplingScheme.hpp +++ b/src/tribol/mesh/CouplingScheme.hpp @@ -743,6 +743,11 @@ class CouplingScheme { m_formulation = std::move( formulation ); } + /** + * @brief Rebuilds the contact formulation based on the current scheme parameters + */ + void updateContactFormulation(); + /** * @brief Check if a ContactFormulation implementation is set * diff --git a/src/tribol/physics/ContactFormulationFactory.cpp b/src/tribol/physics/ContactFormulationFactory.cpp index c0212092..2bb006ec 100644 --- a/src/tribol/physics/ContactFormulationFactory.cpp +++ b/src/tribol/physics/ContactFormulationFactory.cpp @@ -43,7 +43,7 @@ std::unique_ptr createContactFormulation( CouplingScheme* cs SLIC_ERROR_ROOT_IF( !cs->hasMfemSubmeshData(), "ENERGY_MORTAR requires MFEM submesh data." ); SLIC_ERROR_ROOT_IF( !cs->hasMfemJacobianData(), "ENERGY_MORTAR requires MFEM Jacobian data." ); - if ( cs->getParameters().energy_mortar_enforcement_option == EnergyMortarEnforcementOption::QUADRATURE_POINT_GAP ) { + if ( cs->getParameters().energy_mortar_enforcement_option == EnergyMortarEnforcementOption::QuadraturePointGap ) { return std::make_unique>( *cs->getMfemMeshData(), *cs->getMfemSubmeshData(), *cs->getMfemJacobianData(), k, delta, N, enzyme_quadrature, use_penalty_ ); From bf2f119ed26582282479fe85405b060f62c481c9 Mon Sep 17 00:00:00 2001 From: "E. B. Chin" Date: Wed, 5 Aug 2026 16:51:22 -0700 Subject: [PATCH 06/14] design classes with parameterized inheritance --- src/tribol/physics/ContactFormulation.hpp | 39 ++++-- src/tribol/physics/EnergyMortarAdapter.cpp | 120 ++++++++++--------- src/tribol/physics/EnergyMortarAdapter.hpp | 132 +++++++-------------- 3 files changed, 138 insertions(+), 153 deletions(-) diff --git a/src/tribol/physics/ContactFormulation.hpp b/src/tribol/physics/ContactFormulation.hpp index c48e37f3..f699e7af 100644 --- a/src/tribol/physics/ContactFormulation.hpp +++ b/src/tribol/physics/ContactFormulation.hpp @@ -7,6 +7,7 @@ #define SRC_TRIBOL_PHYSICS_CONTACTFORMULATION_HPP_ #include "tribol/config.hpp" +#include #include "tribol/common/Parameters.hpp" #include "tribol/common/ArrayTypes.hpp" @@ -124,23 +125,33 @@ class ContactFormulation { * * @note Requires updateNodalForces() to be called first. */ - virtual const mfem::HypreParVector& getMfemForce() const = 0; + virtual const mfem::HypreParVector& getMfemForce() const + { + SLIC_ERROR_ROOT( "getMfemForce() is not supported by this formulation." ); + throw std::runtime_error( "Not supported" ); + } /** * @brief Returns t-dof vector of gaps on submesh * * @note Requires updateNodalGaps() to be called first. */ - virtual const mfem::HypreParVector& getMfemGap() const = 0; + virtual const mfem::HypreParVector& getMfemGap() const + { + SLIC_ERROR_ROOT( "getMfemGap() is not supported by this formulation." ); + throw std::runtime_error( "Not supported" ); + } /** * @brief Returns a reference to the MFEM dual t-dof vector * * @return Reference to the dual t-dof vector (e.g. pressure in penalty mode, or Lagrange multiplier in LM mode) - * - * TODO: specify what mesh object this is define on. */ - virtual mfem::HypreParVector& getMfemPressure() = 0; + virtual mfem::HypreParVector& getMfemPressure() + { + SLIC_ERROR_ROOT( "getMfemPressure() is not supported by this formulation." ); + throw std::runtime_error( "Not supported" ); + } /** * @brief Get the derivative of force with respect to displacement @@ -149,7 +160,11 @@ class ContactFormulation { * * @note Requires updateNodalForces() to be called first. */ - virtual std::unique_ptr getMfemDfDx() const = 0; + virtual std::unique_ptr getMfemDfDx() const + { + SLIC_ERROR_ROOT( "getMfemDfDx() is not supported by this formulation." ); + return nullptr; + } /** * @brief Get the derivative of the gap constraint with respect to displacement @@ -158,7 +173,11 @@ class ContactFormulation { * * @note Requires updateNodalGaps() to be called first. */ - virtual std::unique_ptr getMfemDgDx() const = 0; + virtual std::unique_ptr getMfemDgDx() const + { + SLIC_ERROR_ROOT( "getMfemDgDx() is not supported by this formulation." ); + return nullptr; + } /** * @brief Get the derivative of force with respect to the dual variable @@ -167,7 +186,11 @@ class ContactFormulation { * * @note Requires updateNodalForces() to be called first. */ - virtual std::unique_ptr getMfemDfDp() const = 0; + virtual std::unique_ptr getMfemDfDp() const + { + SLIC_ERROR_ROOT( "getMfemDfDp() is not supported by this formulation." ); + return nullptr; + } #endif }; diff --git a/src/tribol/physics/EnergyMortarAdapter.cpp b/src/tribol/physics/EnergyMortarAdapter.cpp index 7885bf6d..7d179c19 100644 --- a/src/tribol/physics/EnergyMortarAdapter.cpp +++ b/src/tribol/physics/EnergyMortarAdapter.cpp @@ -11,10 +11,10 @@ namespace tribol { #ifdef TRIBOL_USE_ENZYME -template +template