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 ) { diff --git a/src/shared/math/ParSparseMat.cpp b/src/shared/math/ParSparseMat.cpp index e12d177b..c36822f7 100644 --- a/src/shared/math/ParSparseMat.cpp +++ b/src/shared/math/ParSparseMat.cpp @@ -97,10 +97,11 @@ void ParSparseMatView::eliminateRows( const mfem::Array& rows ) invokeHypreMethod( [&]() { mat_->EliminateRows( rows ); } ); } -void ParSparseMatView::eliminateRowsCols( const mfem::Array& rows_cols ) +ParSparseMat ParSparseMatView::eliminateRowsCols( const mfem::Array& rows_cols ) { ensureHostMemory( *this ); - invokeHypreMethod( [&]() { mat_->EliminateRowsCols( rows_cols ); } ); + return ParSparseMat( + createHypreParMatrix( [&]() { return mat_->EliminateRowsCols( rows_cols ); } ) ); } ParSparseMat ParSparseMatView::eliminateCols( const mfem::Array& cols ) diff --git a/src/shared/math/ParSparseMat.hpp b/src/shared/math/ParSparseMat.hpp index 4487960e..320f767c 100644 --- a/src/shared/math/ParSparseMat.hpp +++ b/src/shared/math/ParSparseMat.hpp @@ -137,8 +137,9 @@ class ParSparseMatView { * @brief Eliminates chosen rows and columns from the matrix * * @param rows_cols Array of rows/columns to eliminate + * @return Matrix containing the eliminated entries */ - void eliminateRowsCols( const mfem::Array& rows_cols ); + ParSparseMat eliminateRowsCols( const mfem::Array& rows_cols ); /** * @brief Eliminates chosen columns from the matrix diff --git a/src/tests/shared_par_sparse_mat.cpp b/src/tests/shared_par_sparse_mat.cpp index 720bc1d1..77361a31 100644 --- a/src/tests/shared_par_sparse_mat.cpp +++ b/src/tests/shared_par_sparse_mat.cpp @@ -349,12 +349,12 @@ TEST_F( ParSparseMatTest, Elimination ) // Eliminate row 0 (globally) // Determine if I own row 0 - mfem::Array rows_to_elim; + mfem::Array indices_to_elim; int row_starts_idx = HYPRE_AssumedPartitionCheck() ? 0 : rank; if ( row_starts[row_starts_idx] == 0 ) { - rows_to_elim.Append( 0 ); + indices_to_elim.Append( 0 ); } - A.eliminateRows( rows_to_elim ); + A.eliminateRows( indices_to_elim ); // Check if row 0 is identity (or zero with diagonal 1) // Diagonal matrix means we can just check multiplication @@ -377,6 +377,15 @@ TEST_F( ParSparseMatTest, Elimination ) } } + A = shared::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 3.0 ); + shared::ParSparseMat Ae = A.eliminateRowsCols( indices_to_elim ); + x.fill( 1.0 ); + y = A * x; + shared::ParVector ye = Ae * x; + for ( int i = 0; i < y.size(); ++i ) { + EXPECT_NEAR( y[i] + ye[i], 3.0, 1e-12 ); + } + A = shared::ParSparseMat::diagonalMatrix( MPI_COMM_WORLD, 10, row_starts, 3.0 ); int num_procs; MPI_Comm_size( MPI_COMM_WORLD, &num_procs ); @@ -385,7 +394,7 @@ TEST_F( ParSparseMatTest, Elimination ) auto last_local_col = A.width() - 1; mfem::Array cols_to_elim( { last_local_col } ); - shared::ParSparseMat Ae = A.eliminateCols( cols_to_elim ); + Ae = A.eliminateCols( cols_to_elim ); // Now check A * e_last = 0 // Create vector with 1 at last_local_col, 0 elsewhere @@ -396,7 +405,7 @@ TEST_F( ParSparseMatTest, Elimination ) EXPECT_NEAR( y_last[last_local_col], 0.0, 1e-12 ); // Check Ae * e_last = original value - shared::ParVector ye = Ae * x_last; + ye = Ae * x_last; double expected_val = 3.0; EXPECT_NEAR( ye[last_local_col], expected_val, 1e-12 ); diff --git a/src/tests/tribol_energy_mortar_patch.cpp b/src/tests/tribol_energy_mortar_patch.cpp index f1f17caa..45cee5da 100644 --- a/src/tests/tribol_energy_mortar_patch.cpp +++ b/src/tests/tribol_energy_mortar_patch.cpp @@ -39,7 +39,7 @@ * u_y(x,y) = eps_yy * y * u_x(x,y) = eps_xx * x */ -class MfemMortarEnergyPatchTest : public testing::TestWithParam> { +class MfemMortarEnergyPatchTest : public testing::TestWithParam> { protected: tribol::RealT max_disp_; double l2_err_vec_; @@ -56,6 +56,7 @@ class MfemMortarEnergyPatchTest : public testing::TestWithParam> void SetUp() override { int ref_levels = std::get<0>( GetParam() ); + tribol::EnforcementLocation enforcement_location = std::get<1>( GetParam() ); int order = 1; auto mortar_attrs = std::set( { 5 } ); @@ -175,7 +176,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::setEnforcementLocation( cs_id, enforcement_location ); + tribol::setMfemKinematicConstantPenalty( cs_id, 100.0, 100.0 ); mfem::Vector X( par_fe_space.GetTrueVSize() ); X = 0.0; @@ -356,7 +358,10 @@ TEST_P( MfemMortarEnergyPatchTest, check_patch_test ) MPI_Barrier( MPI_COMM_WORLD ); } -INSTANTIATE_TEST_SUITE_P( tribol, MfemMortarEnergyPatchTest, testing::Values( std::make_tuple( 2 ) ) ); +INSTANTIATE_TEST_SUITE_P( tribol, MfemMortarEnergyPatchTest, + testing::Combine( testing::Values( 2 ), + testing::Values( tribol::EnforcementLocation::Nodal, + tribol::EnforcementLocation::QuadraturePoint ) ) ); //------------------------------------------------------------------------------ #include "axom/slic/core/SimpleLogger.hpp" diff --git a/src/tests/tribol_finite_diff_energy_mortar.cpp b/src/tests/tribol_finite_diff_energy_mortar.cpp index 450a5596..fb565b6c 100644 --- a/src/tests/tribol_finite_diff_energy_mortar.cpp +++ b/src/tests/tribol_finite_diff_energy_mortar.cpp @@ -372,6 +372,122 @@ FiniteDiffResult EnergyMortarCalculator::validate_hessian( const InterfacePair& return result; } +TEST( QuadraturePointPenaltyCheck, OpenGapIsInactive ) +{ + RealT x1[2] = { 0.0, 1.0 }; + RealT y1[2] = { 0.0, 0.0 }; + IndexT conn1[2] = { 1, 0 }; + MeshData mesh1( 0, 1, 2, conn1, LINEAR_EDGE, x1, y1, nullptr, MemorySpace::Host ); + + RealT x2[2] = { 0.2, 0.8 }; + RealT y2[2] = { 0.1, 0.1 }; + IndexT conn2[2] = { 0, 1 }; + MeshData mesh2( 1, 1, 2, conn2, LINEAR_EDGE, x2, y2, nullptr, MemorySpace::Host ); + + ContactParams params; + params.del = 0.1; + params.k = 3.0; + params.N = 3; + params.enzyme_quadrature = true; + + EnergyMortarCalculator evaluator( params ); + const auto result = + evaluator.compute_quadrature_point_penalty_data( InterfacePair( 0, 0 ), mesh1.getView(), mesh2.getView() ); + EXPECT_EQ( result.energy, 0.0 ); +} + +TEST( QuadraturePointPenaltyCheck, DerivativesMatchFiniteDifference ) +{ + RealT x1[2] = { 0.0, 1.0 }; + RealT y1[2] = { 0.0, 0.0 }; + IndexT conn1[2] = { 1, 0 }; + MeshData mesh1( 0, 1, 2, conn1, LINEAR_EDGE, x1, y1, nullptr, MemorySpace::Host ); + + RealT x2[2] = { 0.2, 0.8 }; + RealT y2[2] = { -0.1, -0.1 }; + IndexT conn2[2] = { 0, 1 }; + MeshData mesh2( 1, 1, 2, conn2, LINEAR_EDGE, x2, y2, nullptr, MemorySpace::Host ); + + ContactParams params; + params.del = 0.1; + params.k = 3.0; + params.N = 3; + params.enzyme_quadrature = true; + + EnergyMortarCalculator evaluator( params ); + const InterfacePair pair( 0, 0 ); + const auto analytical = evaluator.compute_quadrature_point_penalty_data( pair, mesh1.getView(), mesh2.getView() ); + ASSERT_GT( analytical.energy, 0.0 ); + + const std::array 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 +602,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 5ec6403b..d2ad26a1 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::setEnforcementLocation( cs_id, tribol::EnforcementLocation::QuadraturePoint ); tribol::setPenaltyOptions( cs_id, tribol::KINEMATIC, tribol::KINEMATIC_CONSTANT ); tribol::setMfemKinematicConstantPenalty( cs_id, 1.0, 1.0 ); diff --git a/src/tribol/common/Parameters.hpp b/src/tribol/common/Parameters.hpp index 6bb52366..df31b3d2 100644 --- a/src/tribol/common/Parameters.hpp +++ b/src/tribol/common/Parameters.hpp @@ -155,6 +155,15 @@ enum EnforcementMethod NUM_ENFORCEMENT_METHODS }; +/*! + * \brief Enumerates where contact constraints are enforced + */ +enum class EnforcementLocation +{ + Nodal, ///! Enforce contact using assembled nodal gaps + QuadraturePoint ///! Enforce contact independently at quadrature points +}; + /*! * \brief Enumerates the available spatial binning methods */ @@ -495,6 +504,9 @@ struct Parameters { // constituent face elements, then we don't consider the face-pair a contact candidate. // Note, auto-contact will require registration of element thicknesses. bool auto_contact_check = false; ///! True if auto-contact checks should be enabled + + EnforcementLocation enforcement_location = + EnforcementLocation::QuadraturePoint; ///! Defaults to quadrature-point enforcement }; } // namespace tribol diff --git a/src/tribol/interface/mfem_tribol.cpp b/src/tribol/interface/mfem_tribol.cpp index 39f2fa70..d130084f 100644 --- a/src/tribol/interface/mfem_tribol.cpp +++ b/src/tribol/interface/mfem_tribol.cpp @@ -14,7 +14,6 @@ #include "tribol.hpp" #include "tribol/mesh/CouplingScheme.hpp" -#include "tribol/physics/ContactFormulationFactory.hpp" namespace tribol { @@ -293,6 +292,7 @@ void registerMfemCouplingScheme( IndexT cs_id, int mesh_id_1, int mesh_id_2, con auto& cs = CouplingSchemeManager::getInstance().at( cs_id ); cs.setMPIComm( mesh.GetComm() ); if ( contact_method == ENERGY_MORTAR && enforcement_method == LAGRANGE_MULTIPLIER ) { + cs.getParameters().enforcement_location = EnforcementLocation::Nodal; SLIC_WARNING_ROOT( "ENERGY_MORTAR with Lagrange multiplier enforcement is experimental, has no testing, and has " "no support from Tribol developers." ); @@ -338,7 +338,7 @@ void registerMfemCouplingScheme( IndexT cs_id, int mesh_id_1, int mesh_id_2, con } cs.setMfemMeshData( std::move( mfem_data ) ); if ( contact_method == ENERGY_MORTAR ) { - cs.setContactFormulation( createContactFormulation( &cs ) ); + cs.updateContactFormulation(); } } diff --git a/src/tribol/interface/tribol.cpp b/src/tribol/interface/tribol.cpp index 895381ab..aca5eb42 100644 --- a/src/tribol/interface/tribol.cpp +++ b/src/tribol/interface/tribol.cpp @@ -26,6 +26,7 @@ #include "axom/slic.hpp" // C/C++ includes + #include #include #include @@ -170,6 +171,20 @@ void setTimestepPenFrac( IndexT cs_id, RealT frac ) } // end setTimestepPenFrac() +//------------------------------------------------------------------------------ +void setEnforcementLocation( IndexT cs_id, EnforcementLocation location ) +{ + auto cs = CouplingSchemeManager::getInstance().findData( cs_id ); + + SLIC_ERROR_ROOT_IF( !cs, "tribol::setEnforcementLocation(): call tribol::registerCouplingScheme() " + << "prior to calling this routine." ); + + cs->getParameters().enforcement_location = location; + + // Automatically rebuild the formulation to reflect the new setting + cs->updateContactFormulation(); +} + //------------------------------------------------------------------------------ void setTimestepScale( IndexT cs_id, RealT scale ) { diff --git a/src/tribol/interface/tribol.hpp b/src/tribol/interface/tribol.hpp index 36e6fa03..306c876a 100644 --- a/src/tribol/interface/tribol.hpp +++ b/src/tribol/interface/tribol.hpp @@ -131,6 +131,16 @@ void setAutoContactPenScale( IndexT cs_id, RealT scale ); */ void setTimestepPenFrac( IndexT cs_id, RealT frac ); +/*! + * \brief Sets where contact constraints are enforced + * + * \param [in] cs_id coupling scheme id + * \param [in] location contact enforcement location + * + * \note Quadrature-point computation is only supported with penalty enforcement. + */ +void setEnforcementLocation( IndexT cs_id, EnforcementLocation location ); + /*! * * \brief sets the timestep scale factor applied to the timestep vote diff --git a/src/tribol/mesh/CouplingScheme.cpp b/src/tribol/mesh/CouplingScheme.cpp index 5c23c4e8..93db4c81 100644 --- a/src/tribol/mesh/CouplingScheme.cpp +++ b/src/tribol/mesh/CouplingScheme.cpp @@ -1155,6 +1155,18 @@ int CouplingScheme::apply( int cycle, RealT t, RealT& dt ) } // end CouplingScheme::apply() +//------------------------------------------------------------------------------ +void CouplingScheme::updateContactFormulation() +{ + if ( m_contactMethod == ENERGY_MORTAR ) { + this->setContactFormulation( createContactFormulation( this ) ); + } else { + SLIC_WARNING_ROOT( + "tribol::CouplingScheme::updateContactFormulation(): rebuilding is only supported for ENERGY_MORTAR at this " + "time." ); + } +} + //------------------------------------------------------------------------------ bool CouplingScheme::init() { diff --git a/src/tribol/mesh/CouplingScheme.hpp b/src/tribol/mesh/CouplingScheme.hpp index b053d126..7114e067 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/ContactFormulation.hpp b/src/tribol/physics/ContactFormulation.hpp index c48e37f3..b64a2332 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 + * @brief Returns a reference to the MFEM dual t-dof vector on the submesh * * @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/ContactFormulationFactory.cpp b/src/tribol/physics/ContactFormulationFactory.cpp index 22b682dc..54d857e8 100644 --- a/src/tribol/physics/ContactFormulationFactory.cpp +++ b/src/tribol/physics/ContactFormulationFactory.cpp @@ -43,9 +43,16 @@ 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 ); + const auto enforcement_location = cs->getParameters().enforcement_location; + if ( enforcement_location == EnforcementLocation::QuadraturePoint ) { + 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.cpp b/src/tribol/physics/EnergyMortar.cpp index 87df94e9..52b4a0cd 100644 --- a/src/tribol/physics/EnergyMortar.cpp +++ b/src/tribol/physics/EnergyMortar.cpp @@ -17,12 +17,20 @@ namespace tribol { namespace { // This MUST match what the ContactParams struct has in EnergyMortarAdapter -// Theese had to be saved locally in order for enzyme to work correctly +// These had to be saved locally in order for enzyme to work correctly struct KernelParams { - int N = 3; // No. of quadrature points - double del = 0.1; // Smoothing parameter + int N{ 3 }; // No. of quadrature points + double del{ 0.1 }; // Smoothing parameter + double k{ 1.0 }; // Penalty stiffness }; +// Return the line-element mapping Jacobian. Local edge coordinates span [-0.5, 0.5], so the Jacobian is the physical +// length of the segment from A0 to A1. +TRIBOL_ENZYME_INLINE double line_jacobian( const double* A0, const double* A1 ) +{ + 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 ) { @@ -228,9 +236,10 @@ TRIBOL_ENZYME_INLINE void gtilde_kernel( const double* x, Gparams* gp, double* g find_normal( A0, A1, nA ); // Only keep the contribution when the edge normals oppose each other. + // NOTE: geomFilter already rejects pairs with co-oriented normals (dot > 0), + // but the clamp is retained for defensive correctness in tests and direct calls. double dot = nB[0] * nA[0] + nB[1] * nA[1]; - // double eta = ( dot < 0 ) ? dot : 0.0; //Normal smoothing - double eta = dot; + double eta = ( dot < 0 ) ? dot : 0.0; double g1 = 0.0, g2 = 0.0; double AI_1 = 0.0, AI_2 = 0.0; @@ -294,9 +303,10 @@ TRIBOL_ENZYME_INLINE void gtilde_kernel_quad( const double* x, const Gparams* gp double nA[2]; find_normal( A0, A1, nA ); // Only keep the contribution when the edge normals oppose each other. + // NOTE: geomFilter already rejects pairs with co-oriented normals (dot > 0), + // but the clamp is retained for defensive correctness in tests and direct calls. double dot = nB[0] * nA[0] + nB[1] * nA[1]; - // double eta = ( dot < 0 ) ? dot : 0.0; - double eta = dot; + double eta = ( dot < 0 ) ? dot : 0.0; double g1 = 0.0, g2 = 0.0; double AI_1 = 0.0, AI_2 = 0.0; @@ -386,9 +396,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 +415,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 +443,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 +460,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 +470,86 @@ 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 ); + // Only keep the contribution when the edge normals oppose each other. + // NOTE: geomFilter already rejects pairs with co-oriented normals (dot > 0), + // but the clamp is retained for defensive correctness in tests and direct calls. + const double dot = nA[0] * nB[0] + nA[1] * nB[1]; + const double eta = ( dot < 0 ) ? dot : 0.0; + 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 +837,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 +870,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 +902,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 +939,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 +950,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..21316132 100644 --- a/src/tribol/physics/EnergyMortar.hpp +++ b/src/tribol/physics/EnergyMortar.hpp @@ -12,23 +12,37 @@ namespace tribol { #ifdef TRIBOL_USE_ENZYME -// EnergyMortar uses a 3-point Gauss-Legendre quad rule +/// Stores quadrature-point locations and weights for the supported Gauss-Legendre rule. struct QuadPoints { - std::array qp; // qp locations - std::array w; // weights + std::array qp; ///< Quadrature-point locations in the local coordinate of the integration edge. + std::array w; ///< Quadrature weights mapped to the local integration interval. }; +/// Parameters controlling ENERGY_MORTAR contact evaluation. struct ContactParams { - double del; // Smoothing Parameter - double k; // Penalty - int N; // Quadrature Points - bool enzyme_quadrature; // Determines how enzyming is performed (default = True) + double del; ///< Smoothing length used for integration bounds. + double k; ///< Penalty stiffness. + int N; ///< Number of quadrature points. + bool enzyme_quadrature; ///< Whether Enzyme differentiates the quadrature construction. }; -// Weighted gap and trib area +/// Stores quadrature-point penalty energy derivatives for one interface pair. +struct QuadraturePointPenaltyData { + static constexpr int dim = 2; ///< Spatial dimension. + static constexpr int max_nodes_per_elem = 2; ///< Maximum nodes on each line element. + static constexpr int pair_size = 2; ///< Number of elements in an interface pair. + static constexpr int num_force_dofs = dim * max_nodes_per_elem * pair_size; ///< Pair coordinate degrees of freedom. + static constexpr int num_stiffness_entries = num_force_dofs * num_force_dofs; ///< Flattened stiffness size. + + double energy{ 0.0 }; ///< Penalty energy for the interface pair. + std::array force{}; ///< Derivative with respect to pair coordinates. + std::array stiffness{}; ///< Flattened force derivative matrix. +}; + +/// Stores weighted nodal gaps and tributary areas for one interface pair. struct NodalContactData { - std::array AI; // Trib area - std::array g_tilde; // Weighted gap + std::array AI; ///< Tributary areas for the two integration-edge nodes. + std::array g_tilde; ///< Weighted gaps for the two integration-edge nodes. }; /// Stores finite-difference and analytical derivative data for validation tests. @@ -63,9 +77,10 @@ struct FiniteDiffResult { double g_tilde2_baseline{ 0.0 }; }; +/// Stores fixed quadrature data passed to differentiated kernels. struct Gparams { - std::array qp; - std::array w; + std::array qp; ///< Quadrature-point locations in the integration-edge local coordinate. + std::array w; ///< Quadrature weights mapped to the local integration interval. }; /// Provides smoothing operations for the Energy Mortar contact formulation. @@ -170,6 +185,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..b58fe61a 100644 --- a/src/tribol/physics/EnergyMortarAdapter.cpp +++ b/src/tribol/physics/EnergyMortarAdapter.cpp @@ -11,9 +11,10 @@ namespace tribol { #ifdef TRIBOL_USE_ENZYME -EnergyMortarAdapter::EnergyMortarAdapter( MfemMeshData& mesh_data, MfemSubmeshData& submesh_data, - MfemJacobianData& jac_data, double k, double delta, int N, - bool enzyme_quadrature, bool use_penalty ) +template