From 916978a3992f66fd385950f19684089929c9b1f4 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Wed, 20 May 2026 14:41:51 -0700 Subject: [PATCH 1/6] implement polygon triangulation and triangle integration for common plane --- src/tests/tribol_common_plane_penalty.cpp | 175 +++++++- src/tests/tribol_enforcement_options.cpp | 21 + src/tests/tribol_mfem_common_plane.cpp | 19 +- src/tribol/common/Parameters.hpp | 2 + src/tribol/integ/FE.hpp | 114 ++++- src/tribol/integ/Integration.hpp | 218 ++++++---- src/tribol/interface/tribol.cpp | 22 + src/tribol/interface/tribol.hpp | 10 + src/tribol/physics/CommonPlane.cpp | 489 +++++++++++++++------- src/tribol/utils/TestUtils.cpp | 3 +- src/tribol/utils/TestUtils.hpp | 4 + 11 files changed, 836 insertions(+), 241 deletions(-) diff --git a/src/tests/tribol_common_plane_penalty.cpp b/src/tests/tribol_common_plane_penalty.cpp index 5f264725..b0be687b 100644 --- a/src/tests/tribol_common_plane_penalty.cpp +++ b/src/tests/tribol_common_plane_penalty.cpp @@ -186,9 +186,75 @@ class CommonPlaneTest : public ::testing::Test { void TearDown() override { this->m_mesh.clear(); } - protected: +protected: +}; + +struct WarpedQuadForceResult { + int err{ -1 }; + RealT gap{ 0. }; + RealT total_abs_force{ 0. }; + RealT total_force_z{ 0. }; }; +WarpedQuadForceResult runWarpedQuadForceCase( tribol::PolyInteg rule ) +{ + constexpr int numVerts = 4; + + RealT x1[numVerts] = { 0.0, 1.0, 1.0, 0.0 }; + RealT y1[numVerts] = { 0.0, 0.0, 1.0, 1.0 }; + RealT z1[numVerts] = { 0.0, 0.0, 0.0, 0.0 }; + + RealT x2[numVerts] = { 0.0, 0.0, 1.0, 1.0 }; + RealT y2[numVerts] = { 0.0, 1.0, 1.0, 0.0 }; + RealT z2[numVerts] = { -0.20, 1.00, 1.00, 1.00 }; + + tribol::IndexT conn1[numVerts] = { 0, 1, 2, 3 }; + tribol::IndexT conn2[numVerts] = { 0, 1, 2, 3 }; + + tribol::registerMesh( 0, 1, numVerts, &conn1[0], (int)( tribol::LINEAR_QUAD ), &x1[0], &y1[0], &z1[0], + tribol::MemorySpace::Host ); + tribol::registerMesh( 1, 1, numVerts, &conn2[0], (int)( tribol::LINEAR_QUAD ), &x2[0], &y2[0], &z2[0], + tribol::MemorySpace::Host ); + + RealT fx1[numVerts] = { 0., 0., 0., 0. }; + RealT fy1[numVerts] = { 0., 0., 0., 0. }; + RealT fz1[numVerts] = { 0., 0., 0., 0. }; + RealT fx2[numVerts] = { 0., 0., 0., 0. }; + RealT fy2[numVerts] = { 0., 0., 0., 0. }; + RealT fz2[numVerts] = { 0., 0., 0., 0. }; + + tribol::registerNodalResponse( 0, &fx1[0], &fy1[0], &fz1[0] ); + tribol::registerNodalResponse( 1, &fx2[0], &fy2[0], &fz2[0] ); + + tribol::setKinematicConstantPenalty( 0, 1.0 ); + tribol::setKinematicConstantPenalty( 1, 1.0 ); + + tribol::registerCouplingScheme( 0, 0, 1, tribol::SURFACE_TO_SURFACE, tribol::NO_CASE, tribol::COMMON_PLANE, + tribol::FRICTIONLESS, tribol::PENALTY, tribol::BINNING_GRID, + tribol::ExecutionMode::Sequential ); + + tribol::setPenaltyOptions( 0, tribol::KINEMATIC, tribol::KINEMATIC_CONSTANT ); + tribol::setCommonPlaneIntegrationOptions( 0, rule, 3 ); + tribol::setContactAreaFrac( 0, 1.e-12 ); + + WarpedQuadForceResult result; + RealT dt = 1.; + result.err = tribol::update( 1, 1., dt ); + + auto& cs = tribol::CouplingSchemeManager::getInstance().at( 0 ); + EXPECT_EQ( 1, cs.getNumActivePairs() ); + result.gap = cs.getCompGeom().getCommonPlane( 0 ).m_gap; + + for ( int i = 0; i < numVerts; ++i ) { + result.total_abs_force += std::abs( fx1[i] ) + std::abs( fy1[i] ) + std::abs( fz1[i] ); + result.total_abs_force += std::abs( fx2[i] ) + std::abs( fy2[i] ) + std::abs( fz2[i] ); + result.total_force_z += fz1[i] + fz2[i]; + } + + tribol::finalize(); + return result; +} + TEST_F( CommonPlaneTest, penetration_gap_check ) { this->m_mesh.mortarMeshId = 0; @@ -245,6 +311,104 @@ TEST_F( CommonPlaneTest, penetration_gap_check ) tribol::finalize(); } +TEST_F( CommonPlaneTest, full_triangle_decomp_quad_execution ) +{ + this->m_mesh.mortarMeshId = 0; + this->m_mesh.nonmortarMeshId = 1; + + const int nElems = 2; + const RealT x_min1 = 0.; + const RealT y_min1 = 0.; + const RealT z_min1 = 0.; + const RealT x_max1 = 1.; + const RealT y_max1 = 1.; + const RealT z_max1 = 1.05; + + const RealT x_min2 = 0.; + const RealT y_min2 = 0.; + const RealT z_min2 = 0.95; + const RealT x_max2 = 1.; + const RealT y_max2 = 1.; + const RealT z_max2 = 2.; + + this->m_mesh.setupContactMeshHex( nElems, nElems, nElems, x_min1, y_min1, z_min1, x_max1, y_max1, z_max1, nElems, + nElems, nElems, x_min2, y_min2, z_min2, x_max2, y_max2, z_max2, 0., 0. ); + + tribol::TestControlParameters parameters; + parameters.dt = 1.e-3; + parameters.const_penalty = 1.0; + parameters.common_plane_rule = tribol::FULL_TRI_DECOMP; + parameters.common_plane_triangle_order = 3; + + int err = this->m_mesh.tribolSetupAndUpdate( tribol::COMMON_PLANE, tribol::PENALTY, tribol::FRICTIONLESS, + tribol::NO_CASE, false, parameters ); + + EXPECT_EQ( err, 0 ); + + tribol::CouplingScheme* couplingScheme = &tribol::CouplingSchemeManager::getInstance().at( 0 ); + compareGaps( couplingScheme, z_min2 - z_max1, 1.E-8, "kinematic_penetration" ); + checkForceSense( couplingScheme ); + + tribol::finalize(); +} + +TEST_F( CommonPlaneTest, full_triangle_decomp_triangle_execution ) +{ + this->m_mesh.mortarMeshId = 0; + this->m_mesh.nonmortarMeshId = 1; + + const int nElems = 2; + const RealT x_min1 = 0.; + const RealT y_min1 = 0.; + const RealT z_min1 = 0.; + const RealT x_max1 = 1.; + const RealT y_max1 = 1.; + const RealT z_max1 = 1.05; + + const RealT x_min2 = 0.; + const RealT y_min2 = 0.; + const RealT z_min2 = 0.95; + const RealT x_max2 = 1.; + const RealT y_max2 = 1.; + const RealT z_max2 = 2.; + + this->m_mesh.setupContactMeshTet( nElems, nElems, nElems, x_min1, y_min1, z_min1, x_max1, y_max1, z_max1, nElems, + nElems, nElems, x_min2, y_min2, z_min2, x_max2, y_max2, z_max2, 0., 0. ); + + tribol::TestControlParameters parameters; + parameters.dt = 1.e-3; + parameters.const_penalty = 1.0; + parameters.common_plane_rule = tribol::FULL_TRI_DECOMP; + parameters.common_plane_triangle_order = 3; + + int err = this->m_mesh.tribolSetupAndUpdate( tribol::COMMON_PLANE, tribol::PENALTY, tribol::FRICTIONLESS, + tribol::NO_CASE, false, parameters ); + + EXPECT_EQ( err, 0 ); + + tribol::CouplingScheme* couplingScheme = &tribol::CouplingSchemeManager::getInstance().at( 0 ); + compareGaps( couplingScheme, z_min2 - z_max1, 1.E-8, "kinematic_penetration" ); + checkForceSense( couplingScheme ); + + tribol::finalize(); +} + +TEST_F( CommonPlaneTest, full_triangle_decomp_warped_quad_local_contact ) +{ + const auto single_point = runWarpedQuadForceCase( tribol::SINGLE_POINT ); + const auto full_tri = runWarpedQuadForceCase( tribol::FULL_TRI_DECOMP ); + + EXPECT_EQ( single_point.err, 0 ); + EXPECT_EQ( full_tri.err, 0 ); + + EXPECT_GT( single_point.gap, 0. ); + EXPECT_GT( full_tri.gap, 0. ); + + EXPECT_NEAR( single_point.total_abs_force, 0., 1.e-12 ); + EXPECT_GT( full_tri.total_abs_force, 1.e-6 ); + EXPECT_NEAR( full_tri.total_force_z, 0., 1.e-12 ); +} + TEST_F( CommonPlaneTest, separation_gap_check ) { this->m_mesh.mortarMeshId = 0; @@ -671,11 +835,12 @@ TEST_F( CommonPlaneTest, common_plane_viscous_tangential_2d ) // by the gap and then divided amongst the edge nodes RealT force_y = 0.5 * gap / numVerts; RealT force_x = visc_coeff * ( vx1[0] - vx2[0] ) / numVerts; + constexpr RealT tol = 5.e-5; for ( int i = 0; i < numVerts; ++i ) { - EXPECT_NEAR( fx1[i], -force_x, 1.e-10 ); - EXPECT_NEAR( fy1[i], -force_y, 1.e-10 ); - EXPECT_NEAR( fx2[i], force_x, 1.e-10 ); - EXPECT_NEAR( fy2[i], force_y, 1.e-10 ); + EXPECT_NEAR( fx1[i], -force_x, tol ); + EXPECT_NEAR( fy1[i], -force_y, tol ); + EXPECT_NEAR( fx2[i], force_x, tol ); + EXPECT_NEAR( fy2[i], force_y, tol ); } } diff --git a/src/tests/tribol_enforcement_options.cpp b/src/tests/tribol_enforcement_options.cpp index 12113bc3..0676df8e 100644 --- a/src/tests/tribol_enforcement_options.cpp +++ b/src/tests/tribol_enforcement_options.cpp @@ -162,6 +162,27 @@ TEST_F( EnforcementOptionsTest, penalty_kinematic_constant_error ) delete mesh; } +TEST_F( EnforcementOptionsTest, common_plane_integration_options_are_stored ) +{ + tribol::TestMesh* mesh = new tribol::TestMesh(); + SetupTest( mesh ); + + RealT penalty = 1.0; + tribol::setKinematicConstantPenalty( 0, penalty ); + tribol::setKinematicConstantPenalty( 1, penalty ); + tribol::setPenaltyOptions( 0, tribol::KINEMATIC, tribol::KINEMATIC_CONSTANT ); + tribol::setCommonPlaneIntegrationOptions( 0, tribol::FULL_TRI_DECOMP, 4 ); + + tribol::CouplingSchemeManager& csManager = tribol::CouplingSchemeManager::getInstance(); + tribol::CouplingScheme* scheme = &csManager.at( 0 ); + + EXPECT_EQ( scheme->getEnforcementOptions().penalty_options.common_plane_rule, tribol::FULL_TRI_DECOMP ); + EXPECT_EQ( scheme->getEnforcementOptions().penalty_options.common_plane_triangle_order, 4 ); + + tribol::finalize(); + delete mesh; +} + TEST_F( EnforcementOptionsTest, penalty_kinematic_element_error ) { // Setup boiler plate test data etc. diff --git a/src/tests/tribol_mfem_common_plane.cpp b/src/tests/tribol_mfem_common_plane.cpp index d6855d71..cdabf297 100644 --- a/src/tests/tribol_mfem_common_plane.cpp +++ b/src/tests/tribol_mfem_common_plane.cpp @@ -37,10 +37,12 @@ * difference explicit time integration scheme. * * Both the element penalty and a constant penalty are tested, with the constant penalty tuned to match the element - * penalty for this case. As a result, the test comparisons are the same for both penalty types. + * penalty for this case. The contact solve is also exercised with both the legacy single-point CommonPlane integration + * rule and the new full triangle-decomposition rule. * */ -class MfemCommonPlaneTest : public testing::TestWithParam> { +class MfemCommonPlaneTest + : public testing::TestWithParam> { protected: tribol::RealT max_disp_; void SetUp() override @@ -166,6 +168,7 @@ class MfemCommonPlaneTest : public testing::TestWithParam( GetParam() ) ); if ( std::get<1>( GetParam() ) == tribol::KINEMATIC_CONSTANT ) { tribol::setMfemKinematicConstantPenalty( coupling_scheme_id, p_kine, p_kine ); } else { @@ -213,10 +216,14 @@ TEST_P( MfemCommonPlaneTest, common_plane ) } INSTANTIATE_TEST_SUITE_P( tribol, MfemCommonPlaneTest, - testing::Values( std::make_tuple( 1, tribol::KINEMATIC_CONSTANT ), - std::make_tuple( 1, tribol::KINEMATIC_ELEMENT ), - std::make_tuple( 2, tribol::KINEMATIC_CONSTANT ), - std::make_tuple( 2, tribol::KINEMATIC_ELEMENT ) ) ); + testing::Values( std::make_tuple( 1, tribol::KINEMATIC_CONSTANT, tribol::SINGLE_POINT ), + std::make_tuple( 1, tribol::KINEMATIC_CONSTANT, tribol::FULL_TRI_DECOMP ), + std::make_tuple( 1, tribol::KINEMATIC_ELEMENT, tribol::SINGLE_POINT ), + std::make_tuple( 1, tribol::KINEMATIC_ELEMENT, tribol::FULL_TRI_DECOMP ), + std::make_tuple( 2, tribol::KINEMATIC_CONSTANT, tribol::SINGLE_POINT ), + std::make_tuple( 2, tribol::KINEMATIC_CONSTANT, tribol::FULL_TRI_DECOMP ), + std::make_tuple( 2, tribol::KINEMATIC_ELEMENT, tribol::SINGLE_POINT ), + std::make_tuple( 2, tribol::KINEMATIC_ELEMENT, tribol::FULL_TRI_DECOMP ) ) ); //------------------------------------------------------------------------------ int main( int argc, char* argv[] ) diff --git a/src/tribol/common/Parameters.hpp b/src/tribol/common/Parameters.hpp index 6bb52366..9aadc965 100644 --- a/src/tribol/common/Parameters.hpp +++ b/src/tribol/common/Parameters.hpp @@ -444,6 +444,8 @@ struct PenaltyEnforcementOptions { PenaltyConstraintType constraint_type; KinematicPenaltyCalculation kinematic_calculation; RatePenaltyCalculation rate_calculation; + PolyInteg common_plane_rule{ SINGLE_POINT }; + int common_plane_triangle_order{ 3 }; bool constraint_type_set{ false }; bool kinematic_calc_set{ false }; diff --git a/src/tribol/integ/FE.hpp b/src/tribol/integ/FE.hpp index 9281dc1e..494eed29 100644 --- a/src/tribol/integ/FE.hpp +++ b/src/tribol/integ/FE.hpp @@ -111,8 +111,8 @@ TRIBOL_HOST_DEVICE inline void SegmentBasis( const RealT* const x, const RealT p * x[2] is equal to 0. * */ -inline void InvIso( const RealT x[3], const RealT* xA, const RealT* yA, const RealT* zA, const int numNodes, - RealT xi[2] ) +TRIBOL_HOST_DEVICE inline void InvIso( const RealT x[3], const RealT* xA, const RealT* yA, const RealT* zA, + const int numNodes, RealT xi[2] ) { if ( numNodes == 4 ) { constexpr int kmax = 15; @@ -309,7 +309,7 @@ void FwdMapLinQuad( const RealT xi[2], RealT xa[4], RealT ya[4], RealT za[4], Re * of each node are as follows (-1,-1), (1,-1), (0,1). * */ -inline void LinIsoTriShapeFunc( const RealT xi, const RealT eta, const int a, RealT& phi ) +TRIBOL_HOST_DEVICE inline void LinIsoTriShapeFunc( const RealT xi, const RealT eta, const int a, RealT& phi ) { switch ( a ) { case 0: @@ -338,7 +338,7 @@ inline void LinIsoTriShapeFunc( const RealT xi, const RealT eta, const int a, Re * \param [in] xi array of length 2 holding parent coordinates * \param [in,out] phi shape function evaluation (array of length 3) */ -inline void LinIsoTriShapeFunc( const RealT* xi, RealT* phi ) +TRIBOL_HOST_DEVICE inline void LinIsoTriShapeFunc( const RealT* xi, RealT* phi ) { phi[0] = 1.0 - xi[0] - xi[1]; phi[1] = xi[0]; @@ -359,7 +359,7 @@ inline void LinIsoTriShapeFunc( const RealT* xi, RealT* phi ) * * */ -inline void LinIsoQuadShapeFunc( const RealT xi, const RealT eta, const int a, RealT& phi ) +TRIBOL_HOST_DEVICE inline void LinIsoQuadShapeFunc( const RealT xi, const RealT eta, const int a, RealT& phi ) { RealT xi_node, eta_node; switch ( a ) { @@ -419,6 +419,17 @@ void DetJQuad( const RealT xi, const RealT eta, const RealT* x, const int dim, R // Implementations //----------------------------------------------------------------------------- +/*! + * + * \brief Returns the number of nodes on a linear contact face or edge for the + * given problem dimension and face order. + * + * \param [in] dim the dimension of the contact problem + * \param [in] order_type the finite element face order/type + * + * \return number of nodes on the corresponding contact face or edge + * + */ TRIBOL_HOST_DEVICE inline int GetNumFaceNodes( int dim, FaceOrderType order_type ) { // SRW consider consolidating this to take a tribol topology and @@ -437,6 +448,99 @@ TRIBOL_HOST_DEVICE inline int GetNumFaceNodes( int dim, FaceOrderType order_type return numNodes; } +//----------------------------------------------------------------------------- +/*! + * + * \brief Evaluates a linear basis function directly on a physical edge, + * triangle, or quadrilateral face. + * + * \param [in] x pointer to stacked physical coordinates of the face vertices + * \param [in] pX x-coordinate of the evaluation point + * \param [in] pY y-coordinate of the evaluation point + * \param [in] pZ z-coordinate of the evaluation point + * \param [in] numNodes number of nodes on the face + * \param [in] vertexId local node id whose basis function is to be evaluated + * \param [in,out] phi evaluated basis function value + * + * \note For triangles and quadrilaterals this routine inverse-maps the physical + * point to parent space and then evaluates the corresponding linear + * isoparametric basis function. + * + */ +TRIBOL_HOST_DEVICE inline void EvalBasisOnPhysicalFace( const RealT* const x, const RealT pX, const RealT pY, + const RealT pZ, const int numNodes, const int vertexId, + RealT& phi ) +{ + if ( numNodes == 2 ) { + SegmentBasis( x, pX, pY, vertexId, phi ); + return; + } + +#ifdef TRIBOL_USE_HOST + SLIC_ERROR_IF( numNodes != 3 && numNodes != 4, + "EvalBasisOnPhysicalFace(): only linear triangle and quadrilateral faces are supported." ); +#endif + + constexpr int max_nodes_per_face = 4; + RealT xA[max_nodes_per_face] = { 0., 0., 0., 0. }; + RealT yA[max_nodes_per_face] = { 0., 0., 0., 0. }; + RealT zA[max_nodes_per_face] = { 0., 0., 0., 0. }; + for ( int i = 0; i < numNodes; ++i ) { + xA[i] = x[3 * i]; + yA[i] = x[3 * i + 1]; + zA[i] = x[3 * i + 2]; + } + + RealT xp[3] = { pX, pY, pZ }; + RealT xi[2] = { 0., 0. }; + InvIso( xp, xA, yA, zA, numNodes, xi ); + + if ( numNodes == 4 ) { + LinIsoQuadShapeFunc( xi[0], xi[1], vertexId, phi ); + } else { + LinIsoTriShapeFunc( xi[0], xi[1], vertexId, phi ); + } +} + +//----------------------------------------------------------------------------- +/*! + * + * \brief Evaluates a vector-valued Galerkin approximation directly on a physical + * edge, triangle, or quadrilateral face. + * + * \param [in] x pointer to stacked physical coordinates of the face vertices + * \param [in] pX x-coordinate of the evaluation point + * \param [in] pY y-coordinate of the evaluation point + * \param [in] pZ z-coordinate of the evaluation point + * \param [in] numNodes number of nodes on the face + * \param [in] galerkinDim vector dimension of the nodal coefficients + * \param [in] nodeVals stacked nodal coefficients for the Galerkin approximation + * \param [in,out] galerkinVal evaluated Galerkin approximation + * + * \note This helper is topology-aware and supports the linear segment, + * triangle, and quadrilateral basis evaluations used by CommonPlane. + * + */ +TRIBOL_HOST_DEVICE inline void GalerkinEvalOnPhysicalFace( const RealT* const x, const RealT pX, const RealT pY, + const RealT pZ, const int numNodes, const int galerkinDim, + RealT* nodeVals, RealT* galerkinVal ) +{ +#ifdef TRIBOL_USE_HOST + SLIC_ERROR_IF( x == nullptr, "GalerkinEvalOnPhysicalFace(): input pointer, x, is NULL." ); + SLIC_ERROR_IF( nodeVals == nullptr, "GalerkinEvalOnPhysicalFace(): input pointer, nodeVals, is NULL." ); + SLIC_ERROR_IF( galerkinVal == nullptr, "GalerkinEvalOnPhysicalFace(): galerkinVal pointer is NULL." ); + SLIC_ERROR_IF( galerkinDim < 1, "GalerkinEvalOnPhysicalFace(): scalar approximations not yet supported." ); +#endif + + for ( int nd = 0; nd < numNodes; ++nd ) { + RealT phi = 0.; + EvalBasisOnPhysicalFace( x, pX, pY, pZ, numNodes, nd, phi ); + for ( int i = 0; i < galerkinDim; ++i ) { + galerkinVal[i] += nodeVals[i + nd * galerkinDim] * phi; + } + } +} + //----------------------------------------------------------------------------- TRIBOL_HOST_DEVICE inline void GalerkinEval( const RealT* const x, const RealT pX, const RealT pY, const RealT pZ, FaceOrderType order_type, BasisEvalType basis_type, int dim, diff --git a/src/tribol/integ/Integration.hpp b/src/tribol/integ/Integration.hpp index 8f9a57eb..9f2bad8b 100644 --- a/src/tribol/integ/Integration.hpp +++ b/src/tribol/integ/Integration.hpp @@ -174,92 +174,168 @@ int NumTWBPointsPerTri( int order ); // Implementations //----------------------------------------------------------------------------- -template <> -TRIBOL_HOST_DEVICE inline void EvalWeakFormIntegral( SurfaceContactElem const& elem, - RealT* const integ1, - RealT* const integ2 ) +TRIBOL_HOST_DEVICE inline void GetCommonPlaneOverlapCentroid( SurfaceContactElem const& elem, RealT cx[3] ) { - // compute the area centroid of the overlap polygon, - // or vertex avg. centroid of the overlap segment, which - // serves as the single integration point - RealT cx[3] = { 0., 0., 0. }; + cx[0] = 0.; + cx[1] = 0.; + cx[2] = 0.; + if ( elem.dim == 2 ) { VertexAvgCentroid( elem.overlapCoords, elem.dim, elem.numPolyVert, cx[0], cx[1], cx[2] ); } else { PolyAreaCentroid( elem.overlapCoords, elem.dim, elem.numPolyVert, cx[0], cx[1], cx[2] ); } +} - // debug: leave commented out so we don't enter loop - { - // SLIC_DEBUG("Integration point: " << cx[0] << ", " << cx[1]); - - // SLIC_DEBUG("Overlap area: " << elem.overlapArea); - // SLIC_DEBUG("Overlap coords: "); - // for (int i=0; inumberOfNodesPerElement(); ++i ) { - ProjectPointToPlane( elem.faceCoords1[elem.dim * i], elem.faceCoords1[elem.dim * i + 1], - elem.faceCoords1[elem.dim * i + 2], elem.overlapNormal[0], elem.overlapNormal[1], - elem.overlapNormal[2], cx[0], cx[1], cx[2], projX1[elem.dim * i], projX1[elem.dim * i + 1], - projX1[elem.dim * i + 2] ); - - ProjectPointToPlane( elem.faceCoords2[elem.dim * i], elem.faceCoords2[elem.dim * i + 1], - elem.faceCoords2[elem.dim * i + 2], elem.overlapNormal[0], elem.overlapNormal[1], - elem.overlapNormal[2], cx[0], cx[1], cx[2], projX2[elem.dim * i], projX2[elem.dim * i + 1], - projX2[elem.dim * i + 2] ); + coords[0] = 0.1666666667; + coords[1] = 0.1666666667; + coords[2] = 0.6666666667; + coords[3] = 0.1666666667; + coords[4] = 0.1666666667; + coords[5] = 0.6666666667; + return 3; + case 3: + case 4: { + constexpr RealT wt1 = 0.109951743655322; + constexpr RealT wt2 = 0.223381589678011; + wts[0] = wt1; + wts[1] = wt1; + wts[2] = wt1; + wts[3] = wt2; + wts[4] = wt2; + wts[5] = wt2; + + constexpr RealT x1 = 0.091576213509771; + constexpr RealT x2 = 0.816847572980459; + constexpr RealT x3 = 0.108103018168070; + constexpr RealT x4 = 0.445948490915965; + coords[0] = x1; + coords[1] = x1; + coords[2] = x2; + coords[3] = x1; + coords[4] = x1; + coords[5] = x2; + coords[6] = x3; + coords[7] = x4; + coords[8] = x4; + coords[9] = x3; + coords[10] = x4; + coords[11] = x4; + return 6; } - } else { // dim == 2 - // loop over number of nodes per edge (same for each mesh) and project nodes to common plane. - // Can use the integration point as the point in the point-normal data. - for ( int i = 0; i < elem.m_mesh1->numberOfNodesPerElement(); ++i ) { - ProjectPointToSegment( elem.faceCoords1[elem.dim * i], elem.faceCoords1[elem.dim * i + 1], elem.overlapNormal[0], - elem.overlapNormal[1], cx[0], cx[1], projX1[elem.dim * i], projX1[elem.dim * i + 1] ); - - ProjectPointToSegment( elem.faceCoords2[elem.dim * i], elem.faceCoords2[elem.dim * i + 1], elem.overlapNormal[0], - elem.overlapNormal[1], cx[0], cx[1], projX2[elem.dim * i], projX2[elem.dim * i + 1] ); + default: +#ifdef TRIBOL_USE_HOST + SLIC_ERROR( "GetCommonPlaneTriangleRule(): only Gauss integration of order 2-4 is implemented." ); +#endif + return 0; + } +} + +TRIBOL_HOST_DEVICE inline void EvalWeakFormIntegralCommonPlaneFullTri( SurfaceContactElem const& elem, + const int tri_order, RealT* const integ1, + RealT* const integ2 ) +{ + if ( elem.dim != 3 ) { + RealT cx[3] = { 0., 0., 0. }; + GetCommonPlaneOverlapCentroid( elem, cx ); + AccumulateCommonPlaneIntegralAtPoint( elem, cx, 1.0, integ1, integ2 ); + return; + } + + constexpr int max_qpts = 6; + RealT rule_wts[max_qpts] = { 0., 0., 0., 0., 0., 0. }; + RealT rule_coords[2 * max_qpts] = { 0. }; + const int num_qpts = GetCommonPlaneTriangleRule( tri_order, rule_wts, rule_coords ); + + RealT centroid[3]; + GetCommonPlaneOverlapCentroid( elem, centroid ); + + RealT xTri[3]; + RealT yTri[3]; + RealT zTri[3]; + + for ( int j = 0; j < elem.numPolyVert; ++j ) { + const int next = ( j == elem.numPolyVert - 1 ) ? 0 : j + 1; + xTri[0] = elem.overlapCoords[elem.dim * j]; + yTri[0] = elem.overlapCoords[elem.dim * j + 1]; + zTri[0] = elem.overlapCoords[elem.dim * j + 2]; + xTri[1] = elem.overlapCoords[elem.dim * next]; + yTri[1] = elem.overlapCoords[elem.dim * next + 1]; + zTri[1] = elem.overlapCoords[elem.dim * next + 2]; + xTri[2] = centroid[0]; + yTri[2] = centroid[1]; + zTri[2] = centroid[2]; + + const RealT area = Area3DTri( xTri, yTri, zTri ); + if ( area <= 0. ) { + continue; + } + + for ( int qp = 0; qp < num_qpts; ++qp ) { + const RealT xi = rule_coords[2 * qp]; + const RealT eta = rule_coords[2 * qp + 1]; + const RealT n0 = 1. - xi - eta; + RealT x[3]; + x[0] = n0 * xTri[0] + xi * xTri[1] + eta * xTri[2]; + x[1] = n0 * yTri[0] + xi * yTri[1] + eta * yTri[2]; + x[2] = n0 * zTri[0] + xi * zTri[1] + eta * zTri[2]; + AccumulateCommonPlaneIntegralAtPoint( elem, x, area * rule_wts[qp], integ1, integ2 ); } } +} - // loop over nodes and compute nodal force integral - // contributions - for ( int a = 0; a < elem.numFaceVert; ++a ) { - EvalBasis( &projX1[0], cx[0], cx[1], cx[2], elem.numFaceVert, a, integ1[a] ); - EvalBasis( &projX2[0], cx[0], cx[1], cx[2], elem.numFaceVert, a, integ2[a] ); +TRIBOL_HOST_DEVICE inline void EvalWeakFormIntegralCommonPlane( SurfaceContactElem const& elem, const PolyInteg rule, + const int tri_order, RealT* const integ1, + RealT* const integ2 ) +{ + switch ( rule ) { + case SINGLE_POINT: { + RealT cx[3] = { 0., 0., 0. }; + GetCommonPlaneOverlapCentroid( elem, cx ); + AccumulateCommonPlaneIntegralAtPoint( elem, cx, 1.0, integ1, integ2 ); + break; + } + case FULL_TRI_DECOMP: + EvalWeakFormIntegralCommonPlaneFullTri( elem, tri_order, integ1, integ2 ); + break; + default: +#ifdef TRIBOL_USE_HOST + SLIC_ERROR( "EvalWeakFormIntegralCommonPlane(): unsupported polygon integration rule." ); +#endif + break; } +} - return; +template <> +TRIBOL_HOST_DEVICE inline void EvalWeakFormIntegral( SurfaceContactElem const& elem, + RealT* const integ1, + RealT* const integ2 ) +{ + RealT cx[3] = { 0., 0., 0. }; + GetCommonPlaneOverlapCentroid( elem, cx ); + AccumulateCommonPlaneIntegralAtPoint( elem, cx, 1.0, integ1, integ2 ); } } // end namespace tribol diff --git a/src/tribol/interface/tribol.cpp b/src/tribol/interface/tribol.cpp index 895381ab..760ebafb 100644 --- a/src/tribol/interface/tribol.cpp +++ b/src/tribol/interface/tribol.cpp @@ -96,6 +96,28 @@ void setPenaltyOptions( IndexT cs_id, PenaltyConstraintType pen_enfrc_option, } // end setPenaltyOptions() +//------------------------------------------------------------------------------ +void setCommonPlaneIntegrationOptions( IndexT cs_id, PolyInteg rule, int triangle_order ) +{ + auto cs = CouplingSchemeManager::getInstance().findData( cs_id ); + + SLIC_ERROR_ROOT_IF( !cs, "tribol::setCommonPlaneIntegrationOptions(): call tribol::registerCouplingScheme() prior " + << "to calling this routine." ); + + auto& penalty_options = cs->getEnforcementOptions().penalty_options; + + if ( !in_range( rule, NUM_INTEG_RULES ) ) { + SLIC_WARNING_ROOT( "tribol::setCommonPlaneIntegrationOptions(): polygon integration rule not available." ); + return; + } + + SLIC_ERROR_ROOT_IF( triangle_order < 2 || triangle_order > 4, + "tribol::setCommonPlaneIntegrationOptions(): triangle quadrature order must be in [2,4]." ); + + penalty_options.common_plane_rule = rule; + penalty_options.common_plane_triangle_order = triangle_order; +} + //------------------------------------------------------------------------------ void setKinematicConstantPenalty( IndexT mesh_id, RealT k ) { diff --git a/src/tribol/interface/tribol.hpp b/src/tribol/interface/tribol.hpp index 36e6fa03..f8c9b26d 100644 --- a/src/tribol/interface/tribol.hpp +++ b/src/tribol/interface/tribol.hpp @@ -65,6 +65,16 @@ void setPenaltyOptions( IndexT cs_id, PenaltyConstraintType pen_enfrc_option, KinematicPenaltyCalculation kinematic_calc, RatePenaltyCalculation rate_calc = NO_RATE_PENALTY ); +/*! + * \brief Sets CommonPlane overlap integration options + * + * \param [in] cs_id coupling scheme id + * \param [in] rule polygon integration rule for CommonPlane force integration + * \param [in] triangle_order order of the triangle quadrature used by FULL_TRI_DECOMP + * \pre user must register coupling scheme prior to setting CommonPlane integration options + */ +void setCommonPlaneIntegrationOptions( IndexT cs_id, PolyInteg rule, int triangle_order = 3 ); + /*! * \brief Sets the constant kinematic penalty stiffness * \param [in] mesh_id mesh id for penalty stiffness diff --git a/src/tribol/physics/CommonPlane.cpp b/src/tribol/physics/CommonPlane.cpp index 39eae36a..51bd794f 100644 --- a/src/tribol/physics/CommonPlane.cpp +++ b/src/tribol/physics/CommonPlane.cpp @@ -16,48 +16,141 @@ namespace tribol { -TRIBOL_HOST_DEVICE RealT ComputeGapRatePressure( CommonPlanePair& plane, const MeshData::Viewer& m1, - const MeshData::Viewer& m2, RealT element_penalty, - RatePenaltyCalculation rate_calc ) -{ - auto fId1 = plane.getCpElementId1(); - auto fId2 = plane.getCpElementId2(); +namespace { - const auto dim = plane.m_dim; +constexpr int max_dim = 3; +constexpr int max_nodes_per_face = 4; +constexpr int max_nodes_per_overlap = 10; - // compute the correct rate_penalty - RealT rate_penalty = 0.; +TRIBOL_HOST_DEVICE inline RealT ComputeRatePenalty( const MeshData::Viewer& m1, const MeshData::Viewer& m2, + RealT element_penalty, RatePenaltyCalculation rate_calc ) +{ switch ( rate_calc ) { case NO_RATE_PENALTY: { return 0.; } case RATE_CONSTANT: { - rate_penalty = - 0.5 * ( m1.getElementData().m_rate_penalty_stiffness + m2.getElementData().m_rate_penalty_stiffness ); - break; + return 0.5 * ( m1.getElementData().m_rate_penalty_stiffness + m2.getElementData().m_rate_penalty_stiffness ); } case RATE_PERCENT: { - rate_penalty = element_penalty * 0.5 * - ( m1.getElementData().m_rate_percent_stiffness + m2.getElementData().m_rate_percent_stiffness ); - break; + return element_penalty * 0.5 * + ( m1.getElementData().m_rate_percent_stiffness + m2.getElementData().m_rate_percent_stiffness ); } default: - // no-op, quiet compiler - break; + return 0.; + } +} + +TRIBOL_HOST_DEVICE inline void EvalLinearFaceAtPoint( const RealT* face_coords, const int num_nodes, + const RealT x_query[3], RealT x_face[3], RealT* phi, + const int value_dim = 0, const RealT* nodal_vals = nullptr, + RealT* values = nullptr ) +{ + RealT xA[max_nodes_per_face] = { 0., 0., 0., 0. }; + RealT yA[max_nodes_per_face] = { 0., 0., 0., 0. }; + RealT zA[max_nodes_per_face] = { 0., 0., 0., 0. }; + for ( int i = 0; i < num_nodes; ++i ) { + xA[i] = face_coords[3 * i]; + yA[i] = face_coords[3 * i + 1]; + zA[i] = face_coords[3 * i + 2]; + } + + RealT xi[2] = { 0., 0. }; + InvIso( x_query, xA, yA, zA, num_nodes, xi ); + + initRealArray( x_face, max_dim, 0. ); + initRealArray( phi, num_nodes, 0. ); + if ( values != nullptr ) { + initRealArray( values, value_dim, 0. ); + } + + for ( int a = 0; a < num_nodes; ++a ) { + if ( num_nodes == 4 ) { + LinIsoQuadShapeFunc( xi[0], xi[1], a, phi[a] ); + } else { + LinIsoTriShapeFunc( xi[0], xi[1], a, phi[a] ); + } + + x_face[0] += xA[a] * phi[a]; + x_face[1] += yA[a] * phi[a]; + x_face[2] += zA[a] * phi[a]; + + if ( values != nullptr ) { + for ( int i = 0; i < value_dim; ++i ) { + values[i] += nodal_vals[i + a * value_dim] * phi[a]; + } + } + } +} + +TRIBOL_HOST_DEVICE inline void AccumulateContactForce( const MeshData::Viewer& mesh1, const MeshData::Viewer& mesh2, + const IndexT index1, const IndexT index2, const int dim, + const int num_nodes_per_face, const RealT force_x, + const RealT force_y, const RealT force_z, const RealT* phi1, + const RealT* phi2 ) +{ + for ( IndexT a = 0; a < num_nodes_per_face; ++a ) { + IndexT node0 = mesh1.getGlobalNodeId( index1, a ); + IndexT node1 = mesh2.getGlobalNodeId( index2, a ); + + const RealT nodal_force_x1 = force_x * phi1[a]; + const RealT nodal_force_y1 = force_y * phi1[a]; + const RealT nodal_force_z1 = force_z * phi1[a]; + + const RealT nodal_force_x2 = force_x * phi2[a]; + const RealT nodal_force_y2 = force_y * phi2[a]; + const RealT nodal_force_z2 = force_z * phi2[a]; + +#ifdef TRIBOL_USE_RAJA + RAJA::atomicAdd( &mesh1.getResponse()[0][node0], -nodal_force_x1 ); + RAJA::atomicAdd( &mesh2.getResponse()[0][node1], nodal_force_x2 ); + + RAJA::atomicAdd( &mesh1.getResponse()[1][node0], -nodal_force_y1 ); + RAJA::atomicAdd( &mesh2.getResponse()[1][node1], nodal_force_y2 ); + + if ( dim == 3 ) { + RAJA::atomicAdd( &mesh1.getResponse()[2][node0], -nodal_force_z1 ); + RAJA::atomicAdd( &mesh2.getResponse()[2][node1], nodal_force_z2 ); + } +#else + mesh1.getResponse()[0][node0] -= nodal_force_x1; + mesh2.getResponse()[0][node1] += nodal_force_x2; + + mesh1.getResponse()[1][node0] -= nodal_force_y1; + mesh2.getResponse()[1][node1] += nodal_force_y2; + + if ( dim == 3 ) { + mesh1.getResponse()[2][node0] -= nodal_force_z1; + mesh2.getResponse()[2][node1] += nodal_force_z2; + } +#endif } // end switch on rate_calc +} - // compute the velocity gap and pressure contribution - constexpr int max_dim = 3; - constexpr int max_nodes_per_elem = 4; +} // namespace + +TRIBOL_HOST_DEVICE RealT ComputeGapRatePressure( CommonPlanePair& plane, const MeshData::Viewer& m1, + const MeshData::Viewer& m2, RealT element_penalty, + RatePenaltyCalculation rate_calc ) +{ + auto fId1 = plane.getCpElementId1(); + auto fId2 = plane.getCpElementId2(); + + const auto dim = plane.m_dim; + const RealT rate_penalty = ComputeRatePenalty( m1, m2, element_penalty, rate_calc ); + if ( rate_penalty == 0. ) { + return 0.; + } - StackArrayT x1; - StackArrayT v1; + // compute the velocity gap and pressure contribution + StackArrayT x1; + StackArrayT v1; auto numNodesPerFace1 = m1.numberOfNodesPerElement(); plane.getFace1Coords( x1, numNodesPerFace1 ); // get avg face coords off the contact plane m1.getFaceVelocities( fId1, v1 ); - StackArrayT x2; - StackArrayT v2; + StackArrayT x2; + StackArrayT v2; auto numNodesPerFace2 = m2.numberOfNodesPerElement(); plane.getFace2Coords( x2, numNodesPerFace2 ); // get avg face coords off the contact plane m2.getFaceVelocities( fId2, v2 ); @@ -76,14 +169,14 @@ TRIBOL_HOST_DEVICE RealT ComputeGapRatePressure( CommonPlanePair& plane, const M RealT cXf1 = plane.m_cXf1; RealT cYf1 = plane.m_cYf1; RealT cZf1 = ( dim == 3 ) ? plane.m_cZf1 : 0.; - GalerkinEval( x1, cXf1, cYf1, cZf1, LINEAR, PHYSICAL, dim, dim, v1, vel_f1 ); + GalerkinEvalOnPhysicalFace( x1, cXf1, cYf1, cZf1, numNodesPerFace1, dim, v1, vel_f1 ); // interpolate nodal velocity at overlap centroid as projected // onto face 2 RealT cXf2 = plane.m_cXf2; RealT cYf2 = plane.m_cYf2; RealT cZf2 = ( dim == 3 ) ? plane.m_cZf2 : 0.; - GalerkinEval( x2, cXf2, cYf2, cZf2, LINEAR, PHYSICAL, dim, dim, v2, vel_f2 ); + GalerkinEvalOnPhysicalFace( x2, cXf2, cYf2, cZf2, numNodesPerFace2, dim, v2, vel_f2 ); // compute velocity gap vector RealT velGap[max_dim]; @@ -148,14 +241,6 @@ int ApplyNormal( CouplingScheme* cs ) // allows for numerically zero interpenetration. RealT gap_tol = cs_view.getGapTol( index1, index2 ); - if ( gap > gap_tol ) { - // We are here if we have a pair that passes ALL geometric - // filter checks, BUT does not actually violate this method's - // gap constraint. - plane.m_inContact = false; - return; - } - // debug force sums // RealT dbg_sum_force1 {0.}; // RealT dbg_sum_force2 {0.}; @@ -165,6 +250,15 @@ int ApplyNormal( CouplingScheme* cs ) RealT penalty_stiff_per_area{ 0. }; auto& enforcement_options = cs_view.getEnforcementOptions(); const PenaltyEnforcementOptions& pen_enfrc_options = enforcement_options.penalty_options; + const bool use_full_tri = pen_enfrc_options.common_plane_rule == FULL_TRI_DECOMP && cs_view.spatialDimension() == 3; + if ( !use_full_tri && gap > gap_tol ) { + // We are here if we have a pair that passes ALL geometric + // filter checks, BUT does not actually violate this method's + // gap constraint. + plane.m_inContact = false; + return; + } + RealT pen_scale1 = mesh1.getElementData().m_penalty_scale; RealT pen_scale2 = mesh2.getElementData().m_penalty_scale; switch ( pen_enfrc_options.kinematic_calculation ) { @@ -236,9 +330,6 @@ int ApplyNormal( CouplingScheme* cs ) /////////////////////////////////////////// // construct array of nodal coordinates - constexpr int max_dim = 3; - constexpr int max_nodes_per_face = 4; - constexpr int max_nodes_per_overlap = 10; RealT xf1[max_dim * max_nodes_per_face]; RealT xf2[max_dim * max_nodes_per_face]; RealT xVert[max_dim * max_nodes_per_overlap]; @@ -298,7 +389,108 @@ int ApplyNormal( CouplingScheme* cs ) // Integration of contact integrals: integral of shape functions over // // contact overlap patch // //////////////////////////////////////////////////////////////////////// - EvalWeakFormIntegral( cntctElem, phi1, phi2 ); + EvalWeakFormIntegralCommonPlane( cntctElem, pen_enfrc_options.common_plane_rule, + pen_enfrc_options.common_plane_triangle_order, phi1, phi2 ); + + if ( use_full_tri ) { + StackArrayT actual_xf1; + StackArrayT actual_xf2; + mesh1.getFaceCoords( index1, actual_xf1 ); + mesh2.getFaceCoords( index2, actual_xf2 ); + + const bool use_rate = pen_enfrc_options.constraint_type == KINEMATIC_AND_RATE; + const RealT rate_penalty = + use_rate ? ComputeRatePenalty( mesh1, mesh2, penalty_stiff_per_area, pen_enfrc_options.rate_calculation ) : 0.; + + StackArrayT actual_vf1; + StackArrayT actual_vf2; + if ( use_rate ) { + mesh1.getFaceVelocities( index1, actual_vf1 ); + mesh2.getFaceVelocities( index2, actual_vf2 ); + } + + constexpr int max_qpts = 6; + RealT rule_wts[max_qpts] = { 0., 0., 0., 0., 0., 0. }; + RealT rule_coords[2 * max_qpts] = { 0. }; + const int num_qpts = GetCommonPlaneTriangleRule( pen_enfrc_options.common_plane_triangle_order, rule_wts, rule_coords ); + + RealT centroid[3]; + GetCommonPlaneOverlapCentroid( cntctElem, centroid ); + + RealT xTri[3]; + RealT yTri[3]; + RealT zTri[3]; + bool has_contact = false; + + for ( int j = 0; j < numPolyVert; ++j ) { + const int next = ( j == numPolyVert - 1 ) ? 0 : j + 1; + xTri[0] = xVert[dim * j]; + yTri[0] = xVert[dim * j + 1]; + zTri[0] = xVert[dim * j + 2]; + xTri[1] = xVert[dim * next]; + yTri[1] = xVert[dim * next + 1]; + zTri[1] = xVert[dim * next + 2]; + xTri[2] = centroid[0]; + yTri[2] = centroid[1]; + zTri[2] = centroid[2]; + + const RealT area = Area3DTri( xTri, yTri, zTri ); + if ( area <= 0. ) { + continue; + } + + for ( int qp = 0; qp < num_qpts; ++qp ) { + const RealT xi = rule_coords[2 * qp]; + const RealT eta = rule_coords[2 * qp + 1]; + const RealT n0 = 1. - xi - eta; + RealT x_q[3]; + x_q[0] = n0 * xTri[0] + xi * xTri[1] + eta * xTri[2]; + x_q[1] = n0 * yTri[0] + xi * yTri[1] + eta * yTri[2]; + x_q[2] = n0 * zTri[0] + xi * zTri[1] + eta * zTri[2]; + + RealT phi_q1[max_nodes_per_face] = { 0., 0., 0., 0. }; + RealT phi_q2[max_nodes_per_face] = { 0., 0., 0., 0. }; + RealT x_qf1[max_dim]; + RealT x_qf2[max_dim]; + RealT vel_q1[max_dim] = { 0., 0., 0. }; + RealT vel_q2[max_dim] = { 0., 0., 0. }; + + EvalLinearFaceAtPoint( actual_xf1, num_nodes_per_face, x_q, x_qf1, phi_q1, use_rate ? dim : 0, + use_rate ? &actual_vf1[0] : nullptr, use_rate ? vel_q1 : nullptr ); + EvalLinearFaceAtPoint( actual_xf2, num_nodes_per_face, x_q, x_qf2, phi_q2, use_rate ? dim : 0, + use_rate ? &actual_vf2[0] : nullptr, use_rate ? vel_q2 : nullptr ); + + RealT local_gap = ( x_qf1[0] - x_qf2[0] ) * overlapNormal[0] + ( x_qf1[1] - x_qf2[1] ) * overlapNormal[1] + + ( x_qf1[2] - x_qf2[2] ) * overlapNormal[2]; + if ( local_gap > gap_tol ) { + continue; + } + + has_contact = true; + + RealT local_pressure = local_gap * penalty_stiff_per_area; + if ( use_rate && rate_penalty > 0. ) { + RealT local_vel_gap = + ( vel_q1[0] - vel_q2[0] ) * overlapNormal[0] + ( vel_q1[1] - vel_q2[1] ) * overlapNormal[1] + + ( vel_q1[2] - vel_q2[2] ) * overlapNormal[2]; + if ( local_vel_gap <= 0. ) { + local_pressure += local_vel_gap * rate_penalty; + } + } + + const RealT weighted_force = area * rule_wts[qp] * local_pressure; + const RealT force_x = overlapNormal[0] * weighted_force; + const RealT force_y = overlapNormal[1] * weighted_force; + const RealT force_z = overlapNormal[2] * weighted_force; + + AccumulateContactForce( mesh1, mesh2, index1, index2, dim, num_nodes_per_face, force_x, force_y, force_z, + phi_q1, phi_q2 ); + } + } + + plane.m_inContact = has_contact; + return; + } /////////////////////////////////////////////////////////////////////// // Computation of full contact nodal force contributions // @@ -319,65 +511,7 @@ int ApplyNormal( CouplingScheme* cs ) force_z = overlapNormal[2] * contact_force; } - ////////////////////////////////////////////////////// - // loop over nodes and compute contact nodal forces // - ////////////////////////////////////////////////////// - for ( IndexT a = 0; a < num_nodes_per_face; ++a ) { - IndexT node0 = mesh1.getGlobalNodeId( index1, a ); - IndexT node1 = mesh2.getGlobalNodeId( index2, a ); - - // if (logLevel == TRIBOL_DEBUG) - // { - // phi_sum_1 += phi1[a]; - // phi_sum_2 += phi2[a]; - // } - - const RealT nodal_force_x1 = force_x * phi1[a]; - const RealT nodal_force_y1 = force_y * phi1[a]; - const RealT nodal_force_z1 = force_z * phi1[a]; - - const RealT nodal_force_x2 = force_x * phi2[a]; - const RealT nodal_force_y2 = force_y * phi2[a]; - const RealT nodal_force_z2 = force_z * phi2[a]; - - // if (logLevel == TRIBOL_DEBUG) - // { - // dbg_sum_force1 += magnitude( nodal_force_x1, - // nodal_force_y1, - // nodal_force_z1 ); - // dbg_sum_force2 += magnitude( nodal_force_x2, - // nodal_force_y2, - // nodal_force_z2 ); - // } - - // accumulate contributions in host code's registered nodal force arrays -#ifdef TRIBOL_USE_RAJA - RAJA::atomicAdd( &mesh1.getResponse()[0][node0], -nodal_force_x1 ); - RAJA::atomicAdd( &mesh2.getResponse()[0][node1], nodal_force_x2 ); - - RAJA::atomicAdd( &mesh1.getResponse()[1][node0], -nodal_force_y1 ); - RAJA::atomicAdd( &mesh2.getResponse()[1][node1], nodal_force_y2 ); - - // there is no z component for 2D - if ( dim == 3 ) { - RAJA::atomicAdd( &mesh1.getResponse()[2][node0], -nodal_force_z1 ); - RAJA::atomicAdd( &mesh2.getResponse()[2][node1], nodal_force_z2 ); - } -#else - mesh1.getResponse()[0][node0] -= nodal_force_x1; - mesh2.getResponse()[0][node1] += nodal_force_x2; - - mesh1.getResponse()[1][node0] -= nodal_force_y1; - mesh2.getResponse()[1][node1] += nodal_force_y2; - - // there is no z component for 2D - if (dim == 3) - { - mesh1.getResponse()[2][node0] -= nodal_force_z1; - mesh2.getResponse()[2][node1] += nodal_force_z2; - } -#endif - } // end for loop over face nodes + AccumulateContactForce( mesh1, mesh2, index1, index2, dim, num_nodes_per_face, force_x, force_y, force_z, phi1, phi2 ); // comment out debug logs; too much output during tests. Keep for easy // debugging if needed @@ -416,24 +550,23 @@ int ApplyTangential( CouplingScheme* const auto dim = plane.m_dim; auto& mesh1 = cs_view.getMesh1View(); auto& mesh2 = cs_view.getMesh2View(); + const PenaltyEnforcementOptions& pen_enfrc_options = cs_view.getEnforcementOptions().penalty_options; // get pair indices IndexT index1 = plane.getCpElementId1(); IndexT index2 = plane.getCpElementId2(); - // compute the velocity gap and pressure contribution - constexpr int max_dim = 3; - constexpr int max_nodes_per_elem = 4; - constexpr int max_nodes_per_overlap = 10; + const bool use_full_tri = pen_enfrc_options.common_plane_rule == FULL_TRI_DECOMP && dim == 3; - StackArrayT x1; - StackArrayT v1; + // compute the velocity gap and pressure contribution + StackArrayT x1; + StackArrayT v1; auto numNodesPerFace1 = mesh1.numberOfNodesPerElement(); plane.getFace1Coords( x1, numNodesPerFace1 ); // get avg face coords off the contact plane mesh1.getFaceVelocities( index1, v1 ); - StackArrayT x2; - StackArrayT v2; + StackArrayT x2; + StackArrayT v2; auto numNodesPerFace2 = mesh2.numberOfNodesPerElement(); plane.getFace2Coords( x2, numNodesPerFace2 ); // get avg face coords off the contact plane mesh2.getFaceVelocities( index2, v2 ); @@ -452,14 +585,14 @@ int ApplyTangential( CouplingScheme* RealT cXf1 = plane.m_cXf1; RealT cYf1 = plane.m_cYf1; RealT cZf1 = ( dim == 3 ) ? plane.m_cZf1 : 0.; - GalerkinEval( x1, cXf1, cYf1, cZf1, LINEAR, PHYSICAL, dim, dim, v1, vel_f1 ); + GalerkinEvalOnPhysicalFace( x1, cXf1, cYf1, cZf1, numNodesPerFace1, dim, v1, vel_f1 ); // interpolate nodal velocity at overlap centroid as projected // onto face 2 RealT cXf2 = plane.m_cXf2; RealT cYf2 = plane.m_cYf2; RealT cZf2 = ( dim == 3 ) ? plane.m_cZf2 : 0.; - GalerkinEval( x2, cXf2, cYf2, cZf2, LINEAR, PHYSICAL, dim, dim, v2, vel_f2 ); + GalerkinEvalOnPhysicalFace( x2, cXf2, cYf2, cZf2, numNodesPerFace2, dim, v2, vel_f2 ); // compute velocity gap vector RealT velGap[max_dim]; @@ -521,8 +654,8 @@ int ApplyTangential( CouplingScheme* cntctElem.overlapArea = plane.m_area; // create arrays to hold nodal residual weak form integral evaluations - RealT phi1[max_nodes_per_elem]; - RealT phi2[max_nodes_per_elem]; + RealT phi1[max_nodes_per_face]; + RealT phi2[max_nodes_per_face]; initRealArray( phi1, numNodesPerFace1, 0. ); initRealArray( phi2, numNodesPerFace2, 0. ); @@ -530,7 +663,99 @@ int ApplyTangential( CouplingScheme* // Integration of contact integrals: integral of shape functions over // // contact overlap patch // //////////////////////////////////////////////////////////////////////// - EvalWeakFormIntegral( cntctElem, phi1, phi2 ); + EvalWeakFormIntegralCommonPlane( cntctElem, pen_enfrc_options.common_plane_rule, + pen_enfrc_options.common_plane_triangle_order, phi1, phi2 ); + + if ( use_full_tri ) { + StackArrayT actual_xf1; + StackArrayT actual_xf2; + StackArrayT actual_vf1; + StackArrayT actual_vf2; + mesh1.getFaceCoords( index1, actual_xf1 ); + mesh2.getFaceCoords( index2, actual_xf2 ); + mesh1.getFaceVelocities( index1, actual_vf1 ); + mesh2.getFaceVelocities( index2, actual_vf2 ); + + constexpr int max_qpts = 6; + RealT rule_wts[max_qpts] = { 0., 0., 0., 0., 0., 0. }; + RealT rule_coords[2 * max_qpts] = { 0. }; + const int num_qpts = GetCommonPlaneTriangleRule( pen_enfrc_options.common_plane_triangle_order, rule_wts, rule_coords ); + + RealT centroid[3]; + GetCommonPlaneOverlapCentroid( cntctElem, centroid ); + + RealT xTri[3]; + RealT yTri[3]; + RealT zTri[3]; + const RealT visc = 0.5 * + ( mesh1.getElementData().m_viscous_damping_coeff + mesh2.getElementData().m_viscous_damping_coeff ); + + for ( int j = 0; j < numPolyVert; ++j ) { + const int next = ( j == numPolyVert - 1 ) ? 0 : j + 1; + xTri[0] = xVert[dim * j]; + yTri[0] = xVert[dim * j + 1]; + zTri[0] = xVert[dim * j + 2]; + xTri[1] = xVert[dim * next]; + yTri[1] = xVert[dim * next + 1]; + zTri[1] = xVert[dim * next + 2]; + xTri[2] = centroid[0]; + yTri[2] = centroid[1]; + zTri[2] = centroid[2]; + + const RealT area = Area3DTri( xTri, yTri, zTri ); + if ( area <= 0. ) { + continue; + } + + for ( int qp = 0; qp < num_qpts; ++qp ) { + const RealT xi = rule_coords[2 * qp]; + const RealT eta = rule_coords[2 * qp + 1]; + const RealT n0 = 1. - xi - eta; + RealT x_q[3]; + x_q[0] = n0 * xTri[0] + xi * xTri[1] + eta * xTri[2]; + x_q[1] = n0 * yTri[0] + xi * yTri[1] + eta * yTri[2]; + x_q[2] = n0 * zTri[0] + xi * zTri[1] + eta * zTri[2]; + + RealT phi_q1[max_nodes_per_face] = { 0., 0., 0., 0. }; + RealT phi_q2[max_nodes_per_face] = { 0., 0., 0., 0. }; + RealT x_qf1[max_dim]; + RealT x_qf2[max_dim]; + RealT vel_q1[max_dim]; + RealT vel_q2[max_dim]; + + EvalLinearFaceAtPoint( actual_xf1, numNodesPerFace1, x_q, x_qf1, phi_q1, dim, &actual_vf1[0], vel_q1 ); + EvalLinearFaceAtPoint( actual_xf2, numNodesPerFace2, x_q, x_qf2, phi_q2, dim, &actual_vf2[0], vel_q2 ); + + RealT local_gap = ( x_qf1[0] - x_qf2[0] ) * overlapNormal[0] + ( x_qf1[1] - x_qf2[1] ) * overlapNormal[1] + + ( x_qf1[2] - x_qf2[2] ) * overlapNormal[2]; + RealT gap_tol = cs_view.getGapTol( index1, index2 ); + if ( local_gap > gap_tol ) { + continue; + } + + RealT velGap_q[max_dim]; + velGap_q[0] = vel_q1[0] - vel_q2[0]; + velGap_q[1] = vel_q1[1] - vel_q2[1]; + velGap_q[2] = vel_q1[2] - vel_q2[2]; + + RealT velGap_dot_n = velGap_q[0] * overlapNormal[0] + velGap_q[1] * overlapNormal[1] + velGap_q[2] * overlapNormal[2]; + RealT velGapTan[max_dim]; + velGapTan[0] = velGap_q[0] - velGap_dot_n * overlapNormal[0]; + velGapTan[1] = velGap_q[1] - velGap_dot_n * overlapNormal[1]; + velGapTan[2] = velGap_q[2] - velGap_dot_n * overlapNormal[2]; + + const RealT weighted_force = area * rule_wts[qp] * visc; + const RealT force_x = weighted_force * velGapTan[0]; + const RealT force_y = weighted_force * velGapTan[1]; + const RealT force_z = weighted_force * velGapTan[2]; + + AccumulateContactForce( mesh1, mesh2, index1, index2, dim, numNodesPerFace1, force_x, force_y, force_z, + phi_q1, phi_q2 ); + } + } + + return; + } ///////////////////////////////////////////////////// // Computation of tangential viscous damping force // @@ -544,49 +769,7 @@ int ApplyTangential( CouplingScheme* force_z = visc * velGapTan[2]; } - ////////////////////////////////////////////////////// - // loop over nodes and compute contact nodal forces // - ////////////////////////////////////////////////////// - for ( IndexT a = 0; a < numNodesPerFace1; ++a ) { - IndexT node0 = mesh1.getGlobalNodeId( index1, a ); - IndexT node1 = mesh2.getGlobalNodeId( index2, a ); - - const RealT nodal_force_x1 = force_x * phi1[a]; - const RealT nodal_force_y1 = force_y * phi1[a]; - const RealT nodal_force_z1 = force_z * phi1[a]; - - const RealT nodal_force_x2 = force_x * phi2[a]; - const RealT nodal_force_y2 = force_y * phi2[a]; - const RealT nodal_force_z2 = force_z * phi2[a]; - - // accumulate contributions in host code's registered nodal force arrays -#ifdef TRIBOL_USE_RAJA - RAJA::atomicAdd( &mesh1.getResponse()[0][node0], -nodal_force_x1 ); - RAJA::atomicAdd( &mesh2.getResponse()[0][node1], nodal_force_x2 ); - - RAJA::atomicAdd( &mesh1.getResponse()[1][node0], -nodal_force_y1 ); - RAJA::atomicAdd( &mesh2.getResponse()[1][node1], nodal_force_y2 ); - - // there is no z component for 2D - if ( dim == 3 ) { - RAJA::atomicAdd( &mesh1.getResponse()[2][node0], -nodal_force_z1 ); - RAJA::atomicAdd( &mesh2.getResponse()[2][node1], nodal_force_z2 ); - } -#else - mesh1.getResponse()[0][node0] -= nodal_force_x1; - mesh2.getResponse()[0][node1] += nodal_force_x2; - - mesh1.getResponse()[1][node0] -= nodal_force_y1; - mesh2.getResponse()[1][node1] += nodal_force_y2; - - // there is no z component for 2D - if (dim == 3) - { - mesh1.getResponse()[2][node0] -= nodal_force_z1; - mesh2.getResponse()[2][node1] += nodal_force_z2; - } -#endif - } // end for loop over face nodes + AccumulateContactForce( mesh1, mesh2, index1, index2, dim, numNodesPerFace1, force_x, force_y, force_z, phi1, phi2 ); } ); return 0; diff --git a/src/tribol/utils/TestUtils.cpp b/src/tribol/utils/TestUtils.cpp index 6da885f4..22058f8d 100644 --- a/src/tribol/utils/TestUtils.cpp +++ b/src/tribol/utils/TestUtils.cpp @@ -1142,6 +1142,7 @@ int TestMesh::tribolSetupAndUpdate( ContactMethod method, EnforcementMethod enfo // set penalty options after registering coupling scheme setPenaltyOptions( csIndex, constraint_type, pen_calc, rate_calc ); + setCommonPlaneIntegrationOptions( csIndex, params.common_plane_rule, params.common_plane_triangle_order ); } else if ( ( method == SINGLE_MORTAR || method == ALIGNED_MORTAR ) && enforcement == LAGRANGE_MULTIPLIER ) { // note, eval modes and sparse modes not exposed in the interface to this class @@ -1917,4 +1918,4 @@ mfem::Vector ExplicitMechanics::ComputeInvMass( mfem::ParFiniteElementSpace& fes #endif -} // namespace mfem_ext \ No newline at end of file +} // namespace mfem_ext diff --git a/src/tribol/utils/TestUtils.hpp b/src/tribol/utils/TestUtils.hpp index 52272275..3dfe07f9 100644 --- a/src/tribol/utils/TestUtils.hpp +++ b/src/tribol/utils/TestUtils.hpp @@ -44,6 +44,8 @@ struct TestControlParameters { : penalty_ratio( false ), constant_rate_penalty( false ), percent_rate_penalty( false ), + common_plane_rule( SINGLE_POINT ), + common_plane_triangle_order( 3 ), rate_penalty( 1.0 ), rate_penalty_ratio( 0.0 ), const_penalty( 1.0 ), @@ -65,6 +67,8 @@ struct TestControlParameters { bool penalty_ratio; bool constant_rate_penalty; bool percent_rate_penalty; + PolyInteg common_plane_rule; + int common_plane_triangle_order; RealT rate_penalty; RealT rate_penalty_ratio; RealT const_penalty; From 89106e020786113e3ea1b0fcf5905a0796d3a685 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Thu, 21 May 2026 14:56:32 -0700 Subject: [PATCH 2/6] fix projection direction for common plane quad pts --- src/tribol/physics/CommonPlane.cpp | 237 +++++++++++++++++++++++++---- 1 file changed, 209 insertions(+), 28 deletions(-) diff --git a/src/tribol/physics/CommonPlane.cpp b/src/tribol/physics/CommonPlane.cpp index 51bd794f..a3b74a99 100644 --- a/src/tribol/physics/CommonPlane.cpp +++ b/src/tribol/physics/CommonPlane.cpp @@ -41,46 +41,213 @@ TRIBOL_HOST_DEVICE inline RealT ComputeRatePenalty( const MeshData::Viewer& m1, } } -TRIBOL_HOST_DEVICE inline void EvalLinearFaceAtPoint( const RealT* face_coords, const int num_nodes, - const RealT x_query[3], RealT x_face[3], RealT* phi, - const int value_dim = 0, const RealT* nodal_vals = nullptr, - RealT* values = nullptr ) +TRIBOL_HOST_DEVICE inline bool Solve3x3( const RealT A[3][3], const RealT b[3], RealT x[3] ) { - RealT xA[max_nodes_per_face] = { 0., 0., 0., 0. }; - RealT yA[max_nodes_per_face] = { 0., 0., 0., 0. }; - RealT zA[max_nodes_per_face] = { 0., 0., 0., 0. }; - for ( int i = 0; i < num_nodes; ++i ) { - xA[i] = face_coords[3 * i]; - yA[i] = face_coords[3 * i + 1]; - zA[i] = face_coords[3 * i + 2]; + const RealT detA = A[0][0] * ( A[1][1] * A[2][2] - A[1][2] * A[2][1] ) - + A[0][1] * ( A[1][0] * A[2][2] - A[1][2] * A[2][0] ) + + A[0][2] * ( A[1][0] * A[2][1] - A[1][1] * A[2][0] ); + constexpr RealT det_tol = 1.e-15; + if ( std::abs( detA ) <= det_tol ) { + return false; } - RealT xi[2] = { 0., 0. }; - InvIso( x_query, xA, yA, zA, num_nodes, xi ); + const RealT inv_detA = 1. / detA; + x[0] = inv_detA * ( b[0] * ( A[1][1] * A[2][2] - A[1][2] * A[2][1] ) - + A[0][1] * ( b[1] * A[2][2] - A[1][2] * b[2] ) + + A[0][2] * ( b[1] * A[2][1] - A[1][1] * b[2] ) ); + x[1] = inv_detA * ( A[0][0] * ( b[1] * A[2][2] - A[1][2] * b[2] ) - + b[0] * ( A[1][0] * A[2][2] - A[1][2] * A[2][0] ) + + A[0][2] * ( A[1][0] * b[2] - b[1] * A[2][0] ) ); + x[2] = inv_detA * ( A[0][0] * ( A[1][1] * b[2] - b[1] * A[2][1] ) - + A[0][1] * ( A[1][0] * b[2] - b[1] * A[2][0] ) + + b[0] * ( A[1][0] * A[2][1] - A[1][1] * A[2][0] ) ); + return true; +} +TRIBOL_HOST_DEVICE inline void AccumulateFaceInterpolation( const RealT* face_coords, const int num_nodes, const RealT* phi, + RealT x_face[3], const int value_dim = 0, + const RealT* nodal_vals = nullptr, RealT* values = nullptr ) +{ initRealArray( x_face, max_dim, 0. ); - initRealArray( phi, num_nodes, 0. ); if ( values != nullptr ) { initRealArray( values, value_dim, 0. ); } for ( int a = 0; a < num_nodes; ++a ) { - if ( num_nodes == 4 ) { - LinIsoQuadShapeFunc( xi[0], xi[1], a, phi[a] ); - } else { - LinIsoTriShapeFunc( xi[0], xi[1], a, phi[a] ); + x_face[0] += face_coords[3 * a] * phi[a]; + x_face[1] += face_coords[3 * a + 1] * phi[a]; + x_face[2] += face_coords[3 * a + 2] * phi[a]; + + if ( values != nullptr ) { + for ( int i = 0; i < value_dim; ++i ) { + values[i] += nodal_vals[i + a * value_dim] * phi[a]; + } + } + } +} + +TRIBOL_HOST_DEVICE inline bool EvalLinearFaceAtProjectedPoint( const RealT* face_coords, const int num_nodes, + const RealT x_query[3], const RealT projection_dir[3], + RealT x_face[3], RealT* phi, const int value_dim = 0, + const RealT* nodal_vals = nullptr, RealT* values = nullptr ) +{ + // CommonPlane quadrature points lie on the overlap polygon. Evaluate the face + // fields at the corresponding on-face point found by projecting along the + // common-plane normal rather than by an off-surface closest-point inverse map. + initRealArray( phi, num_nodes, 0. ); + + if ( num_nodes == 3 ) { + const RealT x0[3] = { face_coords[0], face_coords[1], face_coords[2] }; + const RealT e1[3] = { face_coords[3] - x0[0], face_coords[4] - x0[1], face_coords[5] - x0[2] }; + const RealT e2[3] = { face_coords[6] - x0[0], face_coords[7] - x0[1], face_coords[8] - x0[2] }; + + RealT n_face[3]; + crossProd( e1[0], e1[1], e1[2], e2[0], e2[1], e2[2], n_face[0], n_face[1], n_face[2] ); + + const RealT denom = dotProd( n_face[0], n_face[1], n_face[2], projection_dir[0], projection_dir[1], projection_dir[2] ); + constexpr RealT parallel_tol = 1.e-14; + if ( std::abs( denom ) <= parallel_tol ) { + return false; } - x_face[0] += xA[a] * phi[a]; - x_face[1] += yA[a] * phi[a]; - x_face[2] += zA[a] * phi[a]; + const RealT dx0 = x0[0] - x_query[0]; + const RealT dy0 = x0[1] - x_query[1]; + const RealT dz0 = x0[2] - x_query[2]; + const RealT step = + dotProd( n_face[0], n_face[1], n_face[2], dx0, dy0, dz0 ) / denom; + + x_face[0] = x_query[0] + step * projection_dir[0]; + x_face[1] = x_query[1] + step * projection_dir[1]; + x_face[2] = x_query[2] + step * projection_dir[2]; + + const RealT r[3] = { x_face[0] - x0[0], x_face[1] - x0[1], x_face[2] - x0[2] }; + const RealT gram11 = dotProd( e1[0], e1[1], e1[2], e1[0], e1[1], e1[2] ); + const RealT gram12 = dotProd( e1[0], e1[1], e1[2], e2[0], e2[1], e2[2] ); + const RealT gram22 = dotProd( e2[0], e2[1], e2[2], e2[0], e2[1], e2[2] ); + const RealT rhs1 = dotProd( e1[0], e1[1], e1[2], r[0], r[1], r[2] ); + const RealT rhs2 = dotProd( e2[0], e2[1], e2[2], r[0], r[1], r[2] ); + + const RealT gram_det = gram11 * gram22 - gram12 * gram12; + constexpr RealT gram_tol = 1.e-15; + if ( std::abs( gram_det ) <= gram_tol ) { + return false; + } - if ( values != nullptr ) { + const RealT xi = ( gram22 * rhs1 - gram12 * rhs2 ) / gram_det; + const RealT eta = ( gram11 * rhs2 - gram12 * rhs1 ) / gram_det; + phi[0] = 1. - xi - eta; + phi[1] = xi; + phi[2] = eta; + } else if ( num_nodes == 4 ) { + constexpr int max_iter = 25; + constexpr RealT step_tol = 1.e-12; + constexpr RealT residual_tol = 1.e-12; + constexpr RealT xi_tol = 1.e-8; + + RealT xi = 0.; + RealT eta = 0.; + RealT phi0[max_nodes_per_face] = { 0.25, 0.25, 0.25, 0.25 }; + RealT x_center[3]; + AccumulateFaceInterpolation( face_coords, num_nodes, phi0, x_center ); + RealT s = ( x_center[0] - x_query[0] ) * projection_dir[0] + ( x_center[1] - x_query[1] ) * projection_dir[1] + + ( x_center[2] - x_query[2] ) * projection_dir[2]; + bool converged = false; + + for ( int iter = 0; iter < max_iter; ++iter ) { + const RealT xi_node[4] = { 1., -1., -1., 1. }; + const RealT eta_node[4] = { 1., 1., -1., -1. }; + RealT dxdxi[3] = { 0., 0., 0. }; + RealT dxdeta[3] = { 0., 0., 0. }; + initRealArray( phi, num_nodes, 0. ); + initRealArray( x_face, max_dim, 0. ); + + for ( int a = 0; a < num_nodes; ++a ) { + phi[a] = 0.25 * ( 1. + xi_node[a] * xi ) * ( 1. + eta_node[a] * eta ); + const RealT dphi_dxi = 0.25 * xi_node[a] * ( 1. + eta_node[a] * eta ); + const RealT dphi_deta = 0.25 * eta_node[a] * ( 1. + xi_node[a] * xi ); + + const RealT xa = face_coords[3 * a]; + const RealT ya = face_coords[3 * a + 1]; + const RealT za = face_coords[3 * a + 2]; + + x_face[0] += xa * phi[a]; + x_face[1] += ya * phi[a]; + x_face[2] += za * phi[a]; + + dxdxi[0] += xa * dphi_dxi; + dxdxi[1] += ya * dphi_dxi; + dxdxi[2] += za * dphi_dxi; + + dxdeta[0] += xa * dphi_deta; + dxdeta[1] += ya * dphi_deta; + dxdeta[2] += za * dphi_deta; + } + + RealT residual[3] = { x_face[0] - x_query[0] - s * projection_dir[0], x_face[1] - x_query[1] - s * projection_dir[1], + x_face[2] - x_query[2] - s * projection_dir[2] }; + const RealT residual_norm = magnitude( residual[0], residual[1], residual[2] ); + if ( residual_norm <= residual_tol ) { + converged = true; + break; + } + + RealT J[3][3] = { { dxdxi[0], dxdeta[0], -projection_dir[0] }, + { dxdxi[1], dxdeta[1], -projection_dir[1] }, + { dxdxi[2], dxdeta[2], -projection_dir[2] } }; + RealT rhs[3] = { -residual[0], -residual[1], -residual[2] }; + RealT delta[3]; + if ( !Solve3x3( J, rhs, delta ) ) { + return false; + } + + xi += delta[0]; + eta += delta[1]; + s += delta[2]; + + const RealT step_norm = magnitude( delta[0], delta[1], delta[2] ); + if ( step_norm <= step_tol ) { + converged = true; + break; + } + } + + if ( !converged || xi < -1. - xi_tol || xi > 1. + xi_tol || eta < -1. - xi_tol || eta > 1. + xi_tol ) { + return false; + } + + initRealArray( phi, num_nodes, 0. ); + phi[0] = 0.25 * ( 1. + xi ) * ( 1. + eta ); + phi[1] = 0.25 * ( 1. - xi ) * ( 1. + eta ); + phi[2] = 0.25 * ( 1. - xi ) * ( 1. - eta ); + phi[3] = 0.25 * ( 1. + xi ) * ( 1. - eta ); + AccumulateFaceInterpolation( face_coords, num_nodes, phi, x_face ); + + const RealT line_residual = + magnitude( x_face[0] - x_query[0], x_face[1] - x_query[1], x_face[2] - x_query[2] ); + const RealT normal_step = + ( x_face[0] - x_query[0] ) * projection_dir[0] + ( x_face[1] - x_query[1] ) * projection_dir[1] + + ( x_face[2] - x_query[2] ) * projection_dir[2]; + const RealT projection_residual = + magnitude( x_face[0] - x_query[0] - normal_step * projection_dir[0], + x_face[1] - x_query[1] - normal_step * projection_dir[1], + x_face[2] - x_query[2] - normal_step * projection_dir[2] ); + if ( line_residual > 0. && projection_residual > residual_tol * line_residual ) { + return false; + } + } else { + return false; + } + + if ( values != nullptr ) { + initRealArray( values, value_dim, 0. ); + for ( int a = 0; a < num_nodes; ++a ) { for ( int i = 0; i < value_dim; ++i ) { values[i] += nodal_vals[i + a * value_dim] * phi[a]; } } } + + return true; } TRIBOL_HOST_DEVICE inline void AccumulateContactForce( const MeshData::Viewer& mesh1, const MeshData::Viewer& mesh2, @@ -455,10 +622,17 @@ int ApplyNormal( CouplingScheme* cs ) RealT vel_q1[max_dim] = { 0., 0., 0. }; RealT vel_q2[max_dim] = { 0., 0., 0. }; - EvalLinearFaceAtPoint( actual_xf1, num_nodes_per_face, x_q, x_qf1, phi_q1, use_rate ? dim : 0, - use_rate ? &actual_vf1[0] : nullptr, use_rate ? vel_q1 : nullptr ); - EvalLinearFaceAtPoint( actual_xf2, num_nodes_per_face, x_q, x_qf2, phi_q2, use_rate ? dim : 0, - use_rate ? &actual_vf2[0] : nullptr, use_rate ? vel_q2 : nullptr ); + const bool mapped_face1 = + EvalLinearFaceAtProjectedPoint( actual_xf1, num_nodes_per_face, x_q, overlapNormal, x_qf1, phi_q1, + use_rate ? dim : 0, use_rate ? &actual_vf1[0] : nullptr, + use_rate ? vel_q1 : nullptr ); + const bool mapped_face2 = + EvalLinearFaceAtProjectedPoint( actual_xf2, num_nodes_per_face, x_q, overlapNormal, x_qf2, phi_q2, + use_rate ? dim : 0, use_rate ? &actual_vf2[0] : nullptr, + use_rate ? vel_q2 : nullptr ); + if ( !mapped_face1 || !mapped_face2 ) { + continue; + } RealT local_gap = ( x_qf1[0] - x_qf2[0] ) * overlapNormal[0] + ( x_qf1[1] - x_qf2[1] ) * overlapNormal[1] + ( x_qf1[2] - x_qf2[2] ) * overlapNormal[2]; @@ -723,8 +897,15 @@ int ApplyTangential( CouplingScheme* RealT vel_q1[max_dim]; RealT vel_q2[max_dim]; - EvalLinearFaceAtPoint( actual_xf1, numNodesPerFace1, x_q, x_qf1, phi_q1, dim, &actual_vf1[0], vel_q1 ); - EvalLinearFaceAtPoint( actual_xf2, numNodesPerFace2, x_q, x_qf2, phi_q2, dim, &actual_vf2[0], vel_q2 ); + const bool mapped_face1 = + EvalLinearFaceAtProjectedPoint( actual_xf1, numNodesPerFace1, x_q, overlapNormal, x_qf1, phi_q1, dim, + &actual_vf1[0], vel_q1 ); + const bool mapped_face2 = + EvalLinearFaceAtProjectedPoint( actual_xf2, numNodesPerFace2, x_q, overlapNormal, x_qf2, phi_q2, dim, + &actual_vf2[0], vel_q2 ); + if ( !mapped_face1 || !mapped_face2 ) { + continue; + } RealT local_gap = ( x_qf1[0] - x_qf2[0] ) * overlapNormal[0] + ( x_qf1[1] - x_qf2[1] ) * overlapNormal[1] + ( x_qf1[2] - x_qf2[2] ) * overlapNormal[2]; From 339b92ffb9b0cf0a8c68a6cb3b3ff64bf435524e Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Thu, 21 May 2026 15:12:37 -0700 Subject: [PATCH 3/6] add higher order triangle rule --- src/tests/tribol_iso_integ.cpp | 71 ++++++++ src/tribol/integ/Integration.cpp | 94 ++--------- src/tribol/integ/Integration.hpp | 254 ++++++++++++++++++++++++++++- src/tribol/interface/tribol.cpp | 4 +- src/tribol/interface/tribol.hpp | 2 +- src/tribol/physics/CommonPlane.cpp | 8 +- src/tribol/physics/Mortar.cpp | 4 +- 7 files changed, 344 insertions(+), 93 deletions(-) diff --git a/src/tests/tribol_iso_integ.cpp b/src/tests/tribol_iso_integ.cpp index d34fdb23..62fa9d4e 100644 --- a/src/tests/tribol_iso_integ.cpp +++ b/src/tests/tribol_iso_integ.cpp @@ -18,6 +18,43 @@ using RealT = tribol::RealT; +namespace { + +RealT factorial( int n ) +{ + RealT result = 1.; + for ( int i = 2; i <= n; ++i ) { + result *= i; + } + return result; +} + +RealT referenceTriangleMoment( int px, int py ) +{ + return factorial( px ) * factorial( py ) / factorial( px + py + 2 ); +} + +RealT normalizedReferenceTriangleMoment( int px, int py ) +{ + return 2. * referenceTriangleMoment( px, py ); +} + +RealT evalTriangleRuleMoment( bool use_legacy, int order, int px, int py ) +{ + RealT wts[tribol::max_symmetric_triangle_qpts] = { 0. }; + RealT coords[2 * tribol::max_symmetric_triangle_qpts] = { 0. }; + const int num_qpts = + use_legacy ? tribol::GetLegacyTriangleRule( order, wts, coords ) : tribol::GetCommonPlaneTriangleRule( order, wts, coords ); + + RealT value = 0.; + for ( int qp = 0; qp < num_qpts; ++qp ) { + value += wts[qp] * std::pow( coords[2 * qp], px ) * std::pow( coords[2 * qp + 1], py ); + } + return value; +} + +} // namespace + /*! * Test fixture class with some setup necessary to use the * triangular decomposition of a quadrilateral with integration @@ -242,6 +279,40 @@ TEST_F( IsoIntegTest, nonaffine ) EXPECT_EQ( convrg, true ); } +TEST( TriangleRuleTest, legacy_and_symmetric_match_on_shared_orders ) +{ + for ( int order : { 2, 4 } ) { + EXPECT_NEAR( evalTriangleRuleMoment( true, order, 0, 0 ), evalTriangleRuleMoment( false, order, 0, 0 ), 2.e-10 ); + EXPECT_NEAR( evalTriangleRuleMoment( true, order, 2, 0 ), evalTriangleRuleMoment( false, order, 2, 0 ), 2.e-10 ); + EXPECT_NEAR( evalTriangleRuleMoment( true, order, 1, 1 ), evalTriangleRuleMoment( false, order, 1, 1 ), 2.e-10 ); + } +} + +TEST( TriangleRuleTest, gauss_poly_int_tri_supports_order_10 ) +{ + constexpr int dim = 3; + constexpr int num_nodes = 3; + RealT xyz[dim * num_nodes] = { 0., 0., 0., 1., 0., 0., 0., 1., 0. }; + + tribol::SurfaceContactElem elem( dim, xyz, xyz, xyz, num_nodes, num_nodes, nullptr, nullptr, 0, 0 ); + tribol::IntegPts integ; + tribol::GaussPolyIntTri( elem, integ, 10 ); + + RealT area = 0.; + RealT moment73 = 0.; + for ( int ip = 0; ip < integ.numIPs; ++ip ) { + const RealT x = integ.xy[dim * ip]; + const RealT y = integ.xy[dim * ip + 1]; + area += integ.wts[ip]; + moment73 += integ.wts[ip] * std::pow( x, 7 ) * std::pow( y, 3 ); + } + + EXPECT_EQ( integ.numIPs, 75 ); + EXPECT_NEAR( area, 0.5, 1.e-14 ); + EXPECT_NEAR( moment73, referenceTriangleMoment( 7, 3 ), 1.e-14 ); + EXPECT_NEAR( evalTriangleRuleMoment( false, 10, 7, 3 ), normalizedReferenceTriangleMoment( 7, 3 ), 1.e-14 ); +} + int main( int argc, char* argv[] ) { int result = 0; diff --git a/src/tribol/integ/Integration.cpp b/src/tribol/integ/Integration.cpp index a7bf5b20..63f89221 100644 --- a/src/tribol/integ/Integration.cpp +++ b/src/tribol/integ/Integration.cpp @@ -217,79 +217,19 @@ int NumTWBPointsPerTri( int order ) } //------------------------------------------------------------------------------ -void GaussPolyIntTri( SurfaceContactElem const& elem, IntegPts& integ, int k ) +void GaussPolyIntTri( SurfaceContactElem const& elem, IntegPts& integ, int k, TriangleQuadratureRuleFamily family ) { - // determine the number of integration points per triangle in the decomposed - // polygon and the total number of integration points on the polygon - int numTriPoints, numTotalPoints; - switch ( k ) { - case 2: - numTriPoints = 3; - numTotalPoints = numTriPoints * elem.numPolyVert; - break; - case 3: - // don't do anything, default to case 4 - case 4: - numTriPoints = 6; - numTotalPoints = numTriPoints * elem.numPolyVert; - break; - default: - SLIC_ERROR( "GaussPolyIntTri: only Gauss integration of order 2-4 is implemented." ); - return; + constexpr int parentDim = 2; + RealT rule_wts[max_symmetric_triangle_qpts] = { 0. }; + RealT rule_coords[parentDim * max_symmetric_triangle_qpts] = { 0. }; + const int numTriPoints = GetTriangleRule( k, family, rule_wts, rule_coords ); + if ( numTriPoints == 0 ) { + SLIC_ERROR( "GaussPolyIntTri: requested triangle integration rule is not available." ); + return; } - - int parentDim = 2; - + const int numTotalPoints = numTriPoints * elem.numPolyVert; integ.initialize( 3, numTotalPoints ); - // populate wts array and set parent space coordinates of - // integration points on triangle - RealT* coords; - switch ( k ) { - case 2: - for ( int i = 0; i < numTotalPoints; ++i ) { - integ.wts[i] = 0.3333333333; - } - coords = new RealT[6]; - coords[0] = 0.1666666667; - coords[1] = 0.1666666667; - coords[2] = 0.6666666667; - coords[3] = 0.1666666667; - coords[4] = 0.1666666667; - coords[5] = 0.6666666667; - break; - case 3: - case 4: - RealT wt1 = 0.109951743655322; - RealT wt2 = 0.223381589678011; - for ( int i = 0; i < elem.numPolyVert; ++i ) { - integ.wts[numTriPoints * i] = wt1; - integ.wts[numTriPoints * i + 1] = wt1; - integ.wts[numTriPoints * i + 2] = wt1; - integ.wts[numTriPoints * i + 3] = wt2; - integ.wts[numTriPoints * i + 4] = wt2; - integ.wts[numTriPoints * i + 5] = wt2; - } - RealT x1 = 0.091576213509771; - RealT x2 = 0.816847572980459; - RealT x3 = 0.108103018168070; - RealT x4 = 0.445948490915965; - coords = new RealT[12]; - coords[0] = x1; - coords[1] = x1; - coords[2] = x2; - coords[3] = x1; - coords[4] = x1; - coords[5] = x2; - coords[6] = x3; - coords[7] = x4; - coords[8] = x4; - coords[9] = x3; - coords[10] = x4; - coords[11] = x4; - break; - } - // compute area centroid of polygon RealT xTri[3] = { 0., 0., 0. }; RealT yTri[3] = { 0., 0., 0. }; @@ -311,31 +251,29 @@ void GaussPolyIntTri( SurfaceContactElem const& elem, IntegPts& integ, int k ) // compute area of triangle RealT area = Area3DTri( xTri, yTri, zTri ); - for ( int k = 0; k < numTriPoints; ++k ) { + for ( int ip = 0; ip < numTriPoints; ++ip ) { // NOTE: Per Puso 2004, the sum over integration point // evaluations per pallet are multiplied by the pallet area. // // multiply the integration point weights by the // triangle area (note: this is specific to how integrals // are computed on polygonal overlaps for Contact) - integ.wts[numTriPoints * j + k] *= area; + integ.wts[numTriPoints * j + ip] = area * rule_wts[ip]; // group parent space ip coordinates RealT xi[2]; - xi[0] = coords[parentDim * k]; - xi[1] = coords[parentDim * k + 1]; + xi[0] = rule_coords[parentDim * ip]; + xi[1] = rule_coords[parentDim * ip + 1]; // forward map parent space ip coords to physical space RealT x[3]; FwdMapLinTri( xi, xTri, yTri, zTri, x ); - integ.xy[( ( integ.ipDim ) * numTriPoints ) * j + ( integ.ipDim * k )] = x[0]; - integ.xy[( ( integ.ipDim ) * numTriPoints ) * j + ( integ.ipDim * k ) + 1] = x[1]; - integ.xy[( ( integ.ipDim ) * numTriPoints ) * j + ( integ.ipDim * k ) + 2] = x[2]; + integ.xy[( ( integ.ipDim ) * numTriPoints ) * j + ( integ.ipDim * ip )] = x[0]; + integ.xy[( ( integ.ipDim ) * numTriPoints ) * j + ( integ.ipDim * ip ) + 1] = x[1]; + integ.xy[( ( integ.ipDim ) * numTriPoints ) * j + ( integ.ipDim * ip ) + 2] = x[2]; } // end loop over number of ips per triangle } // end loop over triangles - - delete[] coords; } //------------------------------------------------------------------------------ diff --git a/src/tribol/integ/Integration.hpp b/src/tribol/integ/Integration.hpp index 9f2bad8b..30b7f890 100644 --- a/src/tribol/integ/Integration.hpp +++ b/src/tribol/integ/Integration.hpp @@ -90,6 +90,16 @@ template TRIBOL_HOST_DEVICE inline void EvalWeakFormIntegral( SurfaceContactElem const& elem, RealT* const integ1, RealT* const integ2 ); +/// Selector for the triangle quadrature family used by GaussPolyIntTri(). +enum TriangleQuadratureRuleFamily +{ + TRI_RULE_LEGACY, + TRI_RULE_SYMMETRIC +}; + +/// Maximum number of quadrature points in the built-in symmetric triangle rules. +constexpr int max_symmetric_triangle_qpts = 25; + /*! * * \brief Populates the integration points and weights on the IntegPts object @@ -120,13 +130,15 @@ void TWBPolyInt( SurfaceContactElem const& elem, IntegPts& integ, int k ); * \param [in] elem SurfaceContactElem object containing dimension and overlap vertices * \param [in,out] integ IntegPts object holding integration points and weights * \param [in] k order of integration + * \param [in] family selector for the triangle quadrature family * - * \pre order 2 <= k <= 3 + * \pre order 2 <= k <= 10 * \pre integ IntegPts object can be instantiated with no-op constructor. This routine * will allocate and populate necessary data. * */ -void GaussPolyIntTri( SurfaceContactElem const& elem, IntegPts& integ, int k ); +void GaussPolyIntTri( SurfaceContactElem const& elem, IntegPts& integ, int k, + TriangleQuadratureRuleFamily family = TRI_RULE_SYMMETRIC ); /*! * @@ -201,7 +213,18 @@ TRIBOL_HOST_DEVICE inline void AccumulateCommonPlaneIntegralAtPoint( SurfaceCont } } -TRIBOL_HOST_DEVICE inline int GetCommonPlaneTriangleRule( int order, RealT* wts, RealT* coords ) +/*! + * \brief Returns the legacy triangle quadrature rule historically used by + * CommonPlane and GaussPolyIntTri(). + * + * \param [in] order requested rule order + * \param [out] wts quadrature weights normalized so they sum to 1 on a triangle + * \param [out] coords quadrature coordinates stored as stacked (xi, eta) pairs + * + * \note This legacy rule is only available for the previously supported + * orders 2 and 3/4 and is kept for regression comparison tests. + */ +TRIBOL_HOST_DEVICE inline int GetLegacyTriangleRule( int order, RealT* wts, RealT* coords ) { switch ( order ) { case 2: @@ -247,7 +270,226 @@ TRIBOL_HOST_DEVICE inline int GetCommonPlaneTriangleRule( int order, RealT* wts, } default: #ifdef TRIBOL_USE_HOST - SLIC_ERROR( "GetCommonPlaneTriangleRule(): only Gauss integration of order 2-4 is implemented." ); + SLIC_ERROR( "GetLegacyTriangleRule(): only legacy Gauss integration of order 2-4 is implemented." ); +#endif + return 0; + } +} + +namespace detail { + +/*! + * \brief Minimal symmetric triangle quadrature orbit data. + * + * \note The compact orbit tables below are adapted from the symmetric triangle + * rules distributed in PETSc, which cite + * F.D. Witherden and P.E. Vincent, + * "On the identification of symmetric quadrature rules for finite element methods", + * Computers & Mathematics with Applications 69(10), 2015, + * doi:10.1016/j.camwa.2015.03.017. + * + * PETSc stores weights for a reference triangle of area 2. Tribol uses + * weights normalized so the weights sum to 1 and the physical triangle + * area is applied separately, so the imported weights are scaled by 1/2 + * during expansion. + */ +struct SymmetricTriangleRuleData +{ + int num_centroid_orbits; + int num_edge_orbits; + int num_general_orbits; + const RealT* weights; + const RealT* orbits; +}; + +constexpr RealT symmetric_triangle_weight_scale = 0.5; + +constexpr RealT tri_deg2_weights[] = { 6.66666666666666666666666666666666635e-01 }; +constexpr RealT tri_deg2_orbits[] = { 1.66666666666666666666666666666666659e-01, 6.66666666666666666666666666666666635e-01 }; + +constexpr RealT tri_deg4_weights[] = { 4.46763179356022931390014016866245598e-01, 2.19903487310643735276652649800421061e-01 }; +constexpr RealT tri_deg4_orbits[] = { 4.45948490915964886318329253883051984e-01, 1.08103018168070227363341492233896033e-01, + 9.15762135097707434595714634022014804e-02, 8.16847572980458513080857073195597039e-01 }; + +constexpr RealT tri_deg5_weights[] = { 4.50000000000000000000000000000000010e-01, 2.51878361089654305191367891000362687e-01, + 2.64788305577012361475298775666303977e-01 }; +constexpr RealT tri_deg5_orbits[] = { 3.33333333333333333333333333333333317e-01, 1.01286507323456338800987361915123836e-01, + 7.97426985353087322398025276169752328e-01, 4.70142064105115089770441209513447613e-01, + 5.97158717897698204591175809731048219e-02 }; + +constexpr RealT tri_deg6_weights[] = { 1.01689812740413633841873618213737963e-01, 2.33572551452758732050579222771158894e-01, + 1.65702151236747150387106912840884901e-01 }; +constexpr RealT tri_deg6_orbits[] = { 6.30890144915022283403316028708191300e-02, 8.73821971016995543319336794258361644e-01, + 2.49286745170910421291638553107019076e-01, 5.01426509658179157416722893785961848e-01, + 5.31450498448169473532496716313981651e-02, 6.36502499121398647230142594412049640e-01, + 3.10352451033784405416607733956552146e-01 }; + +constexpr RealT tri_deg7_weights[] = { 3.30901002215842620719558969458348911e-02, 2.55888342460311145565802470369292636e-01, + 1.54173292371972135669643041667482776e-01, 1.11757465806399561679632628842028190e-01 }; +constexpr RealT tri_deg7_orbits[] = { 3.37306485545878487149717263008162317e-02, 9.32538702890824302570056547398367537e-01, + 2.41577382595403558950186769837781999e-01, 5.16845234809192882099626460324436002e-01, + 4.74309692504718234209580735949185780e-01, 5.13806149905635315808385281016284391e-02, + 4.70366446525952333414099753568849895e-02, 7.54280040550053177356239324628119970e-01, + 1.98683314797351589302350700014995040e-01 }; + +constexpr RealT tri_deg8_weights[] = { 2.88631215355574336502182220978129237e-01, 1.90183268534569249587792208777168633e-01, + 2.06434741069436500563583100584258068e-01, 6.49169952463961606218518566835611904e-02, + 5.44606283488699885296893801478178481e-02 }; +constexpr RealT tri_deg8_orbits[] = { 3.33333333333333333333333333333333317e-01, 4.59292588292723156028815514494169350e-01, + 8.14148234145536879423689710116613481e-02, 1.70569307751760206622293501491464506e-01, + 6.58861384496479586755412997017070988e-01, 5.05472283170309754584235505965989197e-02, + 8.98905543365938049083152898806802161e-01, 8.39477740995760533721383453929445768e-03, + 7.28492392955404281241000379176061966e-01, 2.63112829634638113421785786284643576e-01 }; + +constexpr RealT tri_deg9_weights[] = { 1.94271592565597667638483965014577269e-01, 1.55655082009548558633478712598807923e-01, + 1.59295477854420506065783548528090548e-01, 6.26694004542781410737096625744186273e-02, + 5.11553513173960625233575971179996460e-02, 8.65670787545787545787545787545787526e-02 }; +constexpr RealT tri_deg9_orbits[] = { 3.33333333333333333333333333333333317e-01, 4.37089591492936637269930364435354971e-01, + 1.25820817014126725460139271129290058e-01, 1.88203535619032730240961280467335542e-01, + 6.23592928761934539518077439065328819e-01, 4.89682519198737627783706924836192818e-01, + 2.06349616025247444325861503276144129e-02, 4.47295133944527098651065899662763588e-02, + 9.10540973211094580269786820067447282e-01, 3.68384120547362836348175987833851049e-02, + 7.41198598784498020690079873523423793e-01, 2.21962989160765695675102527693191078e-01 }; + +constexpr RealT tri_deg10_weights[] = { 1.63486658292571932856237369968355216e-01, 2.67059376262991325511459567981373070e-02, + 9.19159272094894560275758192650956353e-02, 1.27809812792848090865797467525306648e-01, + 6.83692963259188572573831680826845816e-02, 5.05955154145767687780855813656664345e-02 }; +constexpr RealT tri_deg10_orbits[] = { 3.33333333333333333333333333333333317e-01, 3.20553732169435129309845893364897379e-02, + 9.35889253566112974138030821327020524e-01, 1.42161101056564385092162103190958311e-01, + 7.15677797886871229815675793618083377e-01, 3.21812995288835421225097560986048687e-01, + 5.30054118927344028277095673945694069e-01, 1.48132885783820550497806765068257172e-01, + 2.96198894887297676338362694260427776e-02, 6.01233328683459245454742893458687815e-01, + 3.69146781827810986911420837115269408e-01, 2.83676653399384392504357555781301898e-02, + 8.07930600922879065079949902881744115e-01, 1.63701733737182495669614341540125695e-01 }; + +TRIBOL_HOST_DEVICE inline bool GetSymmetricTriangleRuleData( int order, SymmetricTriangleRuleData& rule ) +{ + switch ( order ) { + case 2: + rule = { 0, 1, 0, tri_deg2_weights, tri_deg2_orbits }; + return true; + case 3: + case 4: + rule = { 0, 2, 0, tri_deg4_weights, tri_deg4_orbits }; + return true; + case 5: + rule = { 1, 2, 0, tri_deg5_weights, tri_deg5_orbits }; + return true; + case 6: + rule = { 0, 2, 1, tri_deg6_weights, tri_deg6_orbits }; + return true; + case 7: + rule = { 0, 3, 1, tri_deg7_weights, tri_deg7_orbits }; + return true; + case 8: + rule = { 1, 3, 1, tri_deg8_weights, tri_deg8_orbits }; + return true; + case 9: + rule = { 1, 4, 1, tri_deg9_weights, tri_deg9_orbits }; + return true; + case 10: + rule = { 1, 2, 3, tri_deg10_weights, tri_deg10_orbits }; + return true; + default: + return false; + } +} + +TRIBOL_HOST_DEVICE inline int ExpandSymmetricTriangleRule( const SymmetricTriangleRuleData& rule, RealT* wts, RealT* coords ) +{ + int w_idx = 0; + int c_idx = 0; + int q_idx = 0; + + for ( int orbit = 0; orbit < rule.num_centroid_orbits; ++orbit ) { + const RealT a = rule.orbits[c_idx++]; + wts[q_idx] = symmetric_triangle_weight_scale * rule.weights[w_idx++]; + coords[2 * q_idx] = a; + coords[2 * q_idx + 1] = a; + ++q_idx; + } + + for ( int orbit = 0; orbit < rule.num_edge_orbits; ++orbit ) { + const RealT a = rule.orbits[c_idx++]; + const RealT b = rule.orbits[c_idx++]; + const RealT w = symmetric_triangle_weight_scale * rule.weights[w_idx++]; + + wts[q_idx] = w; + coords[2 * q_idx] = a; + coords[2 * q_idx + 1] = b; + ++q_idx; + + wts[q_idx] = w; + coords[2 * q_idx] = b; + coords[2 * q_idx + 1] = a; + ++q_idx; + + wts[q_idx] = w; + coords[2 * q_idx] = a; + coords[2 * q_idx + 1] = a; + ++q_idx; + } + + for ( int orbit = 0; orbit < rule.num_general_orbits; ++orbit ) { + const RealT a = rule.orbits[c_idx++]; + const RealT b = rule.orbits[c_idx++]; + const RealT c = rule.orbits[c_idx++]; + const RealT w = symmetric_triangle_weight_scale * rule.weights[w_idx++]; + + const RealT xi_eta[6][2] = { { b, c }, { c, b }, { a, c }, { c, a }, { a, b }, { b, a } }; + for ( int i = 0; i < 6; ++i ) { + wts[q_idx] = w; + coords[2 * q_idx] = xi_eta[i][0]; + coords[2 * q_idx + 1] = xi_eta[i][1]; + ++q_idx; + } + } + + return q_idx; +} + +} // namespace detail + +/*! + * \brief Returns the built-in higher-order symmetric triangle quadrature rule. + * + * \param [in] order requested rule order + * \param [out] wts quadrature weights normalized so they sum to 1 on a triangle + * \param [out] coords quadrature coordinates stored as stacked (xi, eta) pairs + * + * \note Orders 2 through 10 are available. Order 3 uses the same minimal rule + * as order 4, matching the PETSc/Witherden-Vincent data set. + */ +TRIBOL_HOST_DEVICE inline int GetCommonPlaneTriangleRule( int order, RealT* wts, RealT* coords ) +{ + detail::SymmetricTriangleRuleData rule; + if ( !detail::GetSymmetricTriangleRuleData( order, rule ) ) { +#ifdef TRIBOL_USE_HOST + SLIC_ERROR( "GetCommonPlaneTriangleRule(): only symmetric triangle integration of order 2-10 is implemented." ); +#endif + return 0; + } + return detail::ExpandSymmetricTriangleRule( rule, wts, coords ); +} + +/*! + * \brief Returns the requested triangle quadrature rule family. + * + * \param [in] order requested rule order + * \param [in] family selector for legacy versus symmetric rule data + * \param [out] wts quadrature weights normalized so they sum to 1 on a triangle + * \param [out] coords quadrature coordinates stored as stacked (xi, eta) pairs + */ +TRIBOL_HOST_DEVICE inline int GetTriangleRule( int order, TriangleQuadratureRuleFamily family, RealT* wts, RealT* coords ) +{ + switch ( family ) { + case TRI_RULE_LEGACY: + return GetLegacyTriangleRule( order, wts, coords ); + case TRI_RULE_SYMMETRIC: + return GetCommonPlaneTriangleRule( order, wts, coords ); + default: +#ifdef TRIBOL_USE_HOST + SLIC_ERROR( "GetTriangleRule(): unsupported triangle rule family." ); #endif return 0; } @@ -264,8 +506,8 @@ TRIBOL_HOST_DEVICE inline void EvalWeakFormIntegralCommonPlaneFullTri( SurfaceCo return; } - constexpr int max_qpts = 6; - RealT rule_wts[max_qpts] = { 0., 0., 0., 0., 0., 0. }; + constexpr int max_qpts = max_symmetric_triangle_qpts; + RealT rule_wts[max_qpts] = { 0. }; RealT rule_coords[2 * max_qpts] = { 0. }; const int num_qpts = GetCommonPlaneTriangleRule( tri_order, rule_wts, rule_coords ); diff --git a/src/tribol/interface/tribol.cpp b/src/tribol/interface/tribol.cpp index 760ebafb..eb7b95a2 100644 --- a/src/tribol/interface/tribol.cpp +++ b/src/tribol/interface/tribol.cpp @@ -111,8 +111,8 @@ void setCommonPlaneIntegrationOptions( IndexT cs_id, PolyInteg rule, int triangl return; } - SLIC_ERROR_ROOT_IF( triangle_order < 2 || triangle_order > 4, - "tribol::setCommonPlaneIntegrationOptions(): triangle quadrature order must be in [2,4]." ); + SLIC_ERROR_ROOT_IF( triangle_order < 2 || triangle_order > 10, + "tribol::setCommonPlaneIntegrationOptions(): triangle quadrature order must be in [2,10]." ); penalty_options.common_plane_rule = rule; penalty_options.common_plane_triangle_order = triangle_order; diff --git a/src/tribol/interface/tribol.hpp b/src/tribol/interface/tribol.hpp index f8c9b26d..cd33b716 100644 --- a/src/tribol/interface/tribol.hpp +++ b/src/tribol/interface/tribol.hpp @@ -70,7 +70,7 @@ void setPenaltyOptions( IndexT cs_id, PenaltyConstraintType pen_enfrc_option, * * \param [in] cs_id coupling scheme id * \param [in] rule polygon integration rule for CommonPlane force integration - * \param [in] triangle_order order of the triangle quadrature used by FULL_TRI_DECOMP + * \param [in] triangle_order order of the triangle quadrature used by FULL_TRI_DECOMP in the range [2,10] * \pre user must register coupling scheme prior to setting CommonPlane integration options */ void setCommonPlaneIntegrationOptions( IndexT cs_id, PolyInteg rule, int triangle_order = 3 ); diff --git a/src/tribol/physics/CommonPlane.cpp b/src/tribol/physics/CommonPlane.cpp index a3b74a99..c3780cc0 100644 --- a/src/tribol/physics/CommonPlane.cpp +++ b/src/tribol/physics/CommonPlane.cpp @@ -576,8 +576,8 @@ int ApplyNormal( CouplingScheme* cs ) mesh2.getFaceVelocities( index2, actual_vf2 ); } - constexpr int max_qpts = 6; - RealT rule_wts[max_qpts] = { 0., 0., 0., 0., 0., 0. }; + constexpr int max_qpts = max_symmetric_triangle_qpts; + RealT rule_wts[max_qpts] = { 0. }; RealT rule_coords[2 * max_qpts] = { 0. }; const int num_qpts = GetCommonPlaneTriangleRule( pen_enfrc_options.common_plane_triangle_order, rule_wts, rule_coords ); @@ -850,8 +850,8 @@ int ApplyTangential( CouplingScheme* mesh1.getFaceVelocities( index1, actual_vf1 ); mesh2.getFaceVelocities( index2, actual_vf2 ); - constexpr int max_qpts = 6; - RealT rule_wts[max_qpts] = { 0., 0., 0., 0., 0., 0. }; + constexpr int max_qpts = max_symmetric_triangle_qpts; + RealT rule_wts[max_qpts] = { 0. }; RealT rule_coords[2 * max_qpts] = { 0. }; const int num_qpts = GetCommonPlaneTriangleRule( pen_enfrc_options.common_plane_triangle_order, rule_wts, rule_coords ); diff --git a/src/tribol/physics/Mortar.cpp b/src/tribol/physics/Mortar.cpp index 91e85f1d..a6d8bb51 100644 --- a/src/tribol/physics/Mortar.cpp +++ b/src/tribol/physics/Mortar.cpp @@ -28,8 +28,8 @@ void ComputeMortarWeights( SurfaceContactElem& elem ) // instantiate integration object IntegPts integ; - // Debug: leave code in for now to call Gauss quadrature on triangle rule - GaussPolyIntTri( elem, integ, 3 ); + // Mortar retains the historical triangle rule for now. + GaussPolyIntTri( elem, integ, 3, TRI_RULE_LEGACY ); // call Taylor-Wingate-Bos integation rule. NOTE: this is not // working. The correct gaps are not being computed. From 2349087ab8228d60db0619020ebd07ad6d99701e Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Fri, 22 May 2026 23:24:49 -0700 Subject: [PATCH 4/6] add pointwise gap eval and higher order quadrature for 2d common plane --- src/tests/tribol_common_plane_penalty.cpp | 115 +++++- src/tests/tribol_enforcement_options.cpp | 6 +- src/tests/tribol_mfem_common_plane.cpp | 8 +- src/tribol/common/Parameters.hpp | 6 +- src/tribol/integ/Integration.hpp | 181 +++++++++- src/tribol/interface/tribol.cpp | 8 +- src/tribol/interface/tribol.hpp | 4 +- src/tribol/physics/CommonPlane.cpp | 416 +++++++++++++++------- src/tribol/utils/TestUtils.cpp | 2 +- src/tribol/utils/TestUtils.hpp | 4 +- 10 files changed, 584 insertions(+), 166 deletions(-) diff --git a/src/tests/tribol_common_plane_penalty.cpp b/src/tests/tribol_common_plane_penalty.cpp index b0be687b..add16f68 100644 --- a/src/tests/tribol_common_plane_penalty.cpp +++ b/src/tests/tribol_common_plane_penalty.cpp @@ -255,6 +255,71 @@ WarpedQuadForceResult runWarpedQuadForceCase( tribol::PolyInteg rule ) return result; } +struct EdgeLocalContactResult { + int err{ -1 }; + tribol::IndexT num_active_pairs{ 0 }; + RealT gap{ 0. }; + RealT total_abs_force{ 0. }; + RealT mesh1_node_force_norm[2] = { 0., 0. }; +}; + +EdgeLocalContactResult runEdgeLocalContactCase( tribol::PolyInteg rule ) +{ + constexpr int numVerts = 2; + + RealT x1[numVerts] = { 1.0, 0.0 }; + RealT y1[numVerts] = { 0.0, 0.0 }; + + RealT x2[numVerts] = { 0.0, 1.0 }; + RealT y2[numVerts] = { -0.2, -0.05 }; + + tribol::IndexT conn1[numVerts] = { 0, 1 }; + tribol::IndexT conn2[numVerts] = { 0, 1 }; + + tribol::registerMesh( 0, 1, numVerts, &conn1[0], (int)( tribol::LINEAR_EDGE ), &x1[0], &y1[0], nullptr, + tribol::MemorySpace::Host ); + tribol::registerMesh( 1, 1, numVerts, &conn2[0], (int)( tribol::LINEAR_EDGE ), &x2[0], &y2[0], nullptr, + tribol::MemorySpace::Host ); + + RealT fx1[numVerts] = { 0., 0. }; + RealT fy1[numVerts] = { 0., 0. }; + RealT fx2[numVerts] = { 0., 0. }; + RealT fy2[numVerts] = { 0., 0. }; + + tribol::registerNodalResponse( 0, &fx1[0], &fy1[0], nullptr ); + tribol::registerNodalResponse( 1, &fx2[0], &fy2[0], nullptr ); + + tribol::setKinematicConstantPenalty( 0, 1. ); + tribol::setKinematicConstantPenalty( 1, 1. ); + + tribol::registerCouplingScheme( 0, 0, 1, tribol::SURFACE_TO_SURFACE, tribol::NO_CASE, tribol::COMMON_PLANE, + tribol::FRICTIONLESS, tribol::PENALTY, tribol::BINNING_GRID, + tribol::ExecutionMode::Sequential ); + + tribol::setPenaltyOptions( 0, tribol::KINEMATIC, tribol::KINEMATIC_CONSTANT ); + tribol::setCommonPlaneIntegrationOptions( 0, rule, 4 ); + tribol::setContactAreaFrac( 0, 1.e-12 ); + + EdgeLocalContactResult result; + RealT dt = 1.; + result.err = tribol::update( 1, 1., dt ); + + tribol::CouplingScheme* couplingScheme = &tribol::CouplingSchemeManager::getInstance().at( 0 ); + result.num_active_pairs = couplingScheme->getNumActivePairs(); + if ( result.num_active_pairs > 0 ) { + result.gap = couplingScheme->getCompGeom().getCommonPlane( 0 ).m_gap; + } + + for ( int i = 0; i < numVerts; ++i ) { + result.total_abs_force += std::abs( fx1[i] ) + std::abs( fy1[i] ) + std::abs( fx2[i] ) + std::abs( fy2[i] ); + result.mesh1_node_force_norm[i] = tribol::magnitude( fx1[i], fy1[i] ); + } + + tribol::finalize(); + + return result; +} + TEST_F( CommonPlaneTest, penetration_gap_check ) { this->m_mesh.mortarMeshId = 0; @@ -311,7 +376,7 @@ TEST_F( CommonPlaneTest, penetration_gap_check ) tribol::finalize(); } -TEST_F( CommonPlaneTest, full_triangle_decomp_quad_execution ) +TEST_F( CommonPlaneTest, multipoint_quad_execution ) { this->m_mesh.mortarMeshId = 0; this->m_mesh.nonmortarMeshId = 1; @@ -337,8 +402,8 @@ TEST_F( CommonPlaneTest, full_triangle_decomp_quad_execution ) tribol::TestControlParameters parameters; parameters.dt = 1.e-3; parameters.const_penalty = 1.0; - parameters.common_plane_rule = tribol::FULL_TRI_DECOMP; - parameters.common_plane_triangle_order = 3; + parameters.common_plane_rule = tribol::MULTI_POINT; + parameters.common_plane_quadrature_order = 3; int err = this->m_mesh.tribolSetupAndUpdate( tribol::COMMON_PLANE, tribol::PENALTY, tribol::FRICTIONLESS, tribol::NO_CASE, false, parameters ); @@ -352,7 +417,7 @@ TEST_F( CommonPlaneTest, full_triangle_decomp_quad_execution ) tribol::finalize(); } -TEST_F( CommonPlaneTest, full_triangle_decomp_triangle_execution ) +TEST_F( CommonPlaneTest, multipoint_triangle_execution ) { this->m_mesh.mortarMeshId = 0; this->m_mesh.nonmortarMeshId = 1; @@ -378,8 +443,8 @@ TEST_F( CommonPlaneTest, full_triangle_decomp_triangle_execution ) tribol::TestControlParameters parameters; parameters.dt = 1.e-3; parameters.const_penalty = 1.0; - parameters.common_plane_rule = tribol::FULL_TRI_DECOMP; - parameters.common_plane_triangle_order = 3; + parameters.common_plane_rule = tribol::MULTI_POINT; + parameters.common_plane_quadrature_order = 3; int err = this->m_mesh.tribolSetupAndUpdate( tribol::COMMON_PLANE, tribol::PENALTY, tribol::FRICTIONLESS, tribol::NO_CASE, false, parameters ); @@ -393,20 +458,46 @@ TEST_F( CommonPlaneTest, full_triangle_decomp_triangle_execution ) tribol::finalize(); } -TEST_F( CommonPlaneTest, full_triangle_decomp_warped_quad_local_contact ) +TEST_F( CommonPlaneTest, multipoint_warped_quad_local_contact ) { const auto single_point = runWarpedQuadForceCase( tribol::SINGLE_POINT ); - const auto full_tri = runWarpedQuadForceCase( tribol::FULL_TRI_DECOMP ); + const auto multi_point = runWarpedQuadForceCase( tribol::MULTI_POINT ); EXPECT_EQ( single_point.err, 0 ); - EXPECT_EQ( full_tri.err, 0 ); + EXPECT_EQ( multi_point.err, 0 ); EXPECT_GT( single_point.gap, 0. ); - EXPECT_GT( full_tri.gap, 0. ); + EXPECT_GT( multi_point.gap, 0. ); EXPECT_NEAR( single_point.total_abs_force, 0., 1.e-12 ); - EXPECT_GT( full_tri.total_abs_force, 1.e-6 ); - EXPECT_NEAR( full_tri.total_force_z, 0., 1.e-12 ); + EXPECT_GT( multi_point.total_abs_force, 1.e-6 ); + EXPECT_NEAR( multi_point.total_force_z, 0., 1.e-12 ); +} + +TEST_F( CommonPlaneTest, multipoint_edge_local_contact ) +{ + const auto single_point = runEdgeLocalContactCase( tribol::SINGLE_POINT ); + const auto multi_point = runEdgeLocalContactCase( tribol::MULTI_POINT ); + + EXPECT_EQ( single_point.err, 0 ); + EXPECT_EQ( multi_point.err, 0 ); + + EXPECT_EQ( single_point.num_active_pairs, 1 ); + EXPECT_EQ( multi_point.num_active_pairs, 1 ); + + EXPECT_LT( single_point.gap, 0. ); + EXPECT_LT( multi_point.gap, 0. ); + + EXPECT_GT( single_point.total_abs_force, 1.e-6 ); + EXPECT_GT( multi_point.total_abs_force, 1.e-6 ); + + const RealT single_point_imbalance = + std::abs( single_point.mesh1_node_force_norm[0] - single_point.mesh1_node_force_norm[1] ); + const RealT multi_point_imbalance = + std::abs( multi_point.mesh1_node_force_norm[0] - multi_point.mesh1_node_force_norm[1] ); + + EXPECT_LT( single_point_imbalance, 1.e-3 ); + EXPECT_GT( multi_point_imbalance, single_point_imbalance + 1.e-3 ); } TEST_F( CommonPlaneTest, separation_gap_check ) diff --git a/src/tests/tribol_enforcement_options.cpp b/src/tests/tribol_enforcement_options.cpp index 0676df8e..11090545 100644 --- a/src/tests/tribol_enforcement_options.cpp +++ b/src/tests/tribol_enforcement_options.cpp @@ -171,13 +171,13 @@ TEST_F( EnforcementOptionsTest, common_plane_integration_options_are_stored ) tribol::setKinematicConstantPenalty( 0, penalty ); tribol::setKinematicConstantPenalty( 1, penalty ); tribol::setPenaltyOptions( 0, tribol::KINEMATIC, tribol::KINEMATIC_CONSTANT ); - tribol::setCommonPlaneIntegrationOptions( 0, tribol::FULL_TRI_DECOMP, 4 ); + tribol::setCommonPlaneIntegrationOptions( 0, tribol::MULTI_POINT, 4 ); tribol::CouplingSchemeManager& csManager = tribol::CouplingSchemeManager::getInstance(); tribol::CouplingScheme* scheme = &csManager.at( 0 ); - EXPECT_EQ( scheme->getEnforcementOptions().penalty_options.common_plane_rule, tribol::FULL_TRI_DECOMP ); - EXPECT_EQ( scheme->getEnforcementOptions().penalty_options.common_plane_triangle_order, 4 ); + EXPECT_EQ( scheme->getEnforcementOptions().penalty_options.common_plane_rule, tribol::MULTI_POINT ); + EXPECT_EQ( scheme->getEnforcementOptions().penalty_options.common_plane_quadrature_order, 4 ); tribol::finalize(); delete mesh; diff --git a/src/tests/tribol_mfem_common_plane.cpp b/src/tests/tribol_mfem_common_plane.cpp index cdabf297..370b2393 100644 --- a/src/tests/tribol_mfem_common_plane.cpp +++ b/src/tests/tribol_mfem_common_plane.cpp @@ -217,13 +217,13 @@ TEST_P( MfemCommonPlaneTest, common_plane ) INSTANTIATE_TEST_SUITE_P( tribol, MfemCommonPlaneTest, testing::Values( std::make_tuple( 1, tribol::KINEMATIC_CONSTANT, tribol::SINGLE_POINT ), - std::make_tuple( 1, tribol::KINEMATIC_CONSTANT, tribol::FULL_TRI_DECOMP ), + std::make_tuple( 1, tribol::KINEMATIC_CONSTANT, tribol::MULTI_POINT ), std::make_tuple( 1, tribol::KINEMATIC_ELEMENT, tribol::SINGLE_POINT ), - std::make_tuple( 1, tribol::KINEMATIC_ELEMENT, tribol::FULL_TRI_DECOMP ), + std::make_tuple( 1, tribol::KINEMATIC_ELEMENT, tribol::MULTI_POINT ), std::make_tuple( 2, tribol::KINEMATIC_CONSTANT, tribol::SINGLE_POINT ), - std::make_tuple( 2, tribol::KINEMATIC_CONSTANT, tribol::FULL_TRI_DECOMP ), + std::make_tuple( 2, tribol::KINEMATIC_CONSTANT, tribol::MULTI_POINT ), std::make_tuple( 2, tribol::KINEMATIC_ELEMENT, tribol::SINGLE_POINT ), - std::make_tuple( 2, tribol::KINEMATIC_ELEMENT, tribol::FULL_TRI_DECOMP ) ) ); + std::make_tuple( 2, tribol::KINEMATIC_ELEMENT, tribol::MULTI_POINT ) ) ); //------------------------------------------------------------------------------ int main( int argc, char* argv[] ) diff --git a/src/tribol/common/Parameters.hpp b/src/tribol/common/Parameters.hpp index 9aadc965..c29d78fe 100644 --- a/src/tribol/common/Parameters.hpp +++ b/src/tribol/common/Parameters.hpp @@ -235,8 +235,8 @@ enum IntNodalFields */ enum PolyInteg { - SINGLE_POINT, ///! Single point integration at centroid of polygon - FULL_TRI_DECOMP, ///! Full integration using triangular decomposition + SINGLE_POINT, ///! Single point integration at centroid of polygon + MULTI_POINT, ///! Multi-point integration over the overlap NUM_INTEG_RULES }; @@ -445,7 +445,7 @@ struct PenaltyEnforcementOptions { KinematicPenaltyCalculation kinematic_calculation; RatePenaltyCalculation rate_calculation; PolyInteg common_plane_rule{ SINGLE_POINT }; - int common_plane_triangle_order{ 3 }; + int common_plane_quadrature_order{ 3 }; bool constraint_type_set{ false }; bool kinematic_calc_set{ false }; diff --git a/src/tribol/integ/Integration.hpp b/src/tribol/integ/Integration.hpp index 30b7f890..f61eb0da 100644 --- a/src/tribol/integ/Integration.hpp +++ b/src/tribol/integ/Integration.hpp @@ -99,6 +99,8 @@ enum TriangleQuadratureRuleFamily /// Maximum number of quadrature points in the built-in symmetric triangle rules. constexpr int max_symmetric_triangle_qpts = 25; +/// Maximum number of quadrature points in the built-in Gauss-Legendre segment rules. +constexpr int max_segment_gauss_legendre_qpts = 10; /*! * @@ -495,21 +497,178 @@ TRIBOL_HOST_DEVICE inline int GetTriangleRule( int order, TriangleQuadratureRule } } -TRIBOL_HOST_DEVICE inline void EvalWeakFormIntegralCommonPlaneFullTri( SurfaceContactElem const& elem, - const int tri_order, RealT* const integ1, - RealT* const integ2 ) +/*! + * \brief Returns a Gauss-Legendre quadrature rule on the unit segment. + * + * \param [in] order requested rule order + * \param [out] wts quadrature weights normalized so they sum to 1 on [0,1] + * \param [out] coords quadrature coordinates on [0,1] + */ +TRIBOL_HOST_DEVICE inline int GetCommonPlaneSegmentRule( int order, RealT* wts, RealT* coords ) { - if ( elem.dim != 3 ) { - RealT cx[3] = { 0., 0., 0. }; - GetCommonPlaneOverlapCentroid( elem, cx ); - AccumulateCommonPlaneIntegralAtPoint( elem, cx, 1.0, integ1, integ2 ); + switch ( order ) { + case 2: + wts[0] = 5.00000000000000000000000000000000000e-01; + wts[1] = 5.00000000000000000000000000000000000e-01; + coords[0] = 2.11324865405187117745425609795414482e-01; + coords[1] = 7.88675134594812882254574390204585518e-01; + return 2; + case 3: + wts[0] = 2.77777777777777777777777777777777778e-01; + wts[1] = 4.44444444444444444444444444444444444e-01; + wts[2] = 2.77777777777777777777777777777777778e-01; + coords[0] = 1.12701665379258311482063373571554511e-01; + coords[1] = 5.00000000000000000000000000000000000e-01; + coords[2] = 8.87298334620741688517936626428445489e-01; + return 3; + case 4: + wts[0] = 1.73927422568726928648300228976864863e-01; + wts[1] = 3.26072577431273071351699771023135137e-01; + wts[2] = 3.26072577431273071351699771023135137e-01; + wts[3] = 1.73927422568726928648300228976864863e-01; + coords[0] = 6.94318442029737123880253935661376383e-02; + coords[1] = 3.30009478207571867549864872986354601e-01; + coords[2] = 6.69990521792428132450135127013645399e-01; + coords[3] = 9.30568155797026287611974606433862362e-01; + return 4; + case 5: + wts[0] = 1.18463442528094543757132020373224693e-01; + wts[1] = 2.39314335249683234020645713783081311e-01; + wts[2] = 2.84444444444444444444444444444444444e-01; + wts[3] = 2.39314335249683234020645713783081311e-01; + wts[4] = 1.18463442528094543757132020373224693e-01; + coords[0] = 4.69100770306680036011865699692305193e-02; + coords[1] = 2.30765344947158454446500534703347781e-01; + coords[2] = 5.00000000000000000000000000000000000e-01; + coords[3] = 7.69234655052841545553499465296652219e-01; + coords[4] = 9.53089922969331996398813430030769481e-01; + return 5; + case 6: + wts[0] = 8.56622461895851725230519499689802902e-02; + wts[1] = 1.80380786524069303841036681987673464e-01; + wts[2] = 2.33956967286345520487813812579846245e-01; + wts[3] = 2.33956967286345520487813812579846245e-01; + wts[4] = 1.80380786524069303841036681987673464e-01; + wts[5] = 8.56622461895851725230519499689802902e-02; + coords[0] = 3.37652428984239962556928330124159244e-02; + coords[1] = 1.69395306766867745483983092697870979e-01; + coords[2] = 3.80690406958401560316832671838599160e-01; + coords[3] = 6.19309593041598439683167328161400840e-01; + coords[4] = 8.30604693233132254516016907302129021e-01; + coords[5] = 9.66234757101576003744307166987584076e-01; + return 6; + case 7: + wts[0] = 6.47424830844348466391116955198397926e-02; + wts[1] = 1.39852695744638333950704650054195659e-01; + wts[2] = 1.90915025252559472475161990594647293e-01; + wts[3] = 2.08979591836734693877551020408163265e-01; + wts[4] = 1.90915025252559472475161990594647293e-01; + wts[5] = 1.39852695744638333950704650054195659e-01; + wts[6] = 6.47424830844348466391116955198397926e-02; + coords[0] = 2.54460438286207377369030550090367355e-02; + coords[1] = 1.29234424311301331233889510925404835e-01; + coords[2] = 2.97483866778698624490873198198455490e-01; + coords[3] = 5.00000000000000000000000000000000000e-01; + coords[4] = 7.02516133221301375509126801801544510e-01; + coords[5] = 8.70765575688698668766110489074595165e-01; + coords[6] = 9.74553956171379262263096944990963264e-01; + return 7; + case 8: + wts[0] = 5.06142681451881693180916895511138059e-02; + wts[1] = 1.11190517226687235272177997268125811e-01; + wts[2] = 1.56853322938943643668981100993387281e-01; + wts[3] = 1.81341891689181001927177820586651362e-01; + wts[4] = 1.81341891689181001927177820586651362e-01; + wts[5] = 1.56853322938943643668981100993387281e-01; + wts[6] = 1.11190517226687235272177997268125811e-01; + wts[7] = 5.06142681451881693180916895511138059e-02; + coords[0] = 1.98550717512318841525047110753824461e-02; + coords[1] = 1.01666761293186647733518599687717188e-01; + coords[2] = 2.37233795041835507091130475405343431e-01; + coords[3] = 4.08282678752175097530261928819908057e-01; + coords[4] = 5.91717321247824902469738071180091943e-01; + coords[5] = 7.62766204958164492908869524594656569e-01; + coords[6] = 8.98333238706813352266481400312282812e-01; + coords[7] = 9.80144928248768115847495288924617554e-01; + return 8; + case 9: + wts[0] = 4.06371941807872005172919364720240956e-02; + wts[1] = 9.03240803474287173109739194798907182e-02; + wts[2] = 1.30305348201467649015160147969872934e-01; + wts[3] = 1.56173538520001468666934369973026078e-01; + wts[4] = 1.65119677500629881523396871610248955e-01; + wts[5] = 1.56173538520001468666934369973026078e-01; + wts[6] = 1.30305348201467649015160147969872934e-01; + wts[7] = 9.03240803474287173109739194798907182e-02; + wts[8] = 4.06371941807872005172919364720240956e-02; + coords[0] = 1.59198802461869216995971690127085065e-02; + coords[1] = 8.19844463366820870961097794728460945e-02; + coords[2] = 1.93314283649704707966974877456187964e-01; + coords[3] = 3.37873288298095542617664572568897352e-01; + coords[4] = 5.00000000000000000000000000000000000e-01; + coords[5] = 6.62126711701904457382335427431102648e-01; + coords[6] = 8.06685716350295292033025122543812036e-01; + coords[7] = 9.18015553663317912903890220527153906e-01; + coords[8] = 9.84080119753813078300402830987291494e-01; + return 9; + case 10: + wts[0] = 3.33356721543440461444643439469389906e-02; + wts[1] = 7.47256745752902979112342760845433925e-02; + wts[2] = 1.09543181257990906041775930711667835e-01; + wts[3] = 1.34633359654998153565044736633326588e-01; + wts[4] = 1.47762112357376435165896479362247582e-01; + wts[5] = 1.47762112357376435165896479362247582e-01; + wts[6] = 1.34633359654998153565044736633326588e-01; + wts[7] = 1.09543181257990906041775930711667835e-01; + wts[8] = 7.47256745752902979112342760845433925e-02; + wts[9] = 3.33356721543440461444643439469389906e-02; + coords[0] = 1.30467357414141598997194531116613714e-02; + coords[1] = 6.74683166555077431721327998540848508e-02; + coords[2] = 1.60295215850487803884800815437504493e-01; + coords[3] = 2.83302302935376372885100587561369794e-01; + coords[4] = 4.25562830509184389676645012342252051e-01; + coords[5] = 5.74437169490815610323354987657747949e-01; + coords[6] = 7.16697697064623627114899412438630206e-01; + coords[7] = 8.39704784149512196115199184562495507e-01; + coords[8] = 9.32531683344492256827867200145915149e-01; + coords[9] = 9.86953264258585840100280546888338629e-01; + return 10; + default: +#ifdef TRIBOL_USE_HOST + SLIC_ERROR( "GetCommonPlaneSegmentRule(): only Gauss-Legendre integration of order 2-10 is implemented." ); +#endif + return 0; + } +} + +TRIBOL_HOST_DEVICE inline void EvalWeakFormIntegralCommonPlaneMultiPoint( SurfaceContactElem const& elem, + const int quadrature_order, + RealT* const integ1, RealT* const integ2 ) +{ + if ( elem.dim == 2 ) { + RealT rule_wts[max_segment_gauss_legendre_qpts] = { 0. }; + RealT rule_coords[max_segment_gauss_legendre_qpts] = { 0. }; + const int num_qpts = GetCommonPlaneSegmentRule( quadrature_order, rule_wts, rule_coords ); + + const RealT x0 = elem.overlapCoords[0]; + const RealT y0 = elem.overlapCoords[1]; + const RealT x1 = elem.overlapCoords[2]; + const RealT y1 = elem.overlapCoords[3]; + const RealT length = magnitude( x1 - x0, y1 - y0 ); + + for ( int qp = 0; qp < num_qpts; ++qp ) { + const RealT s = rule_coords[qp]; + const RealT one_minus_s = 1. - s; + RealT x[3] = { one_minus_s * x0 + s * x1, one_minus_s * y0 + s * y1, 0. }; + AccumulateCommonPlaneIntegralAtPoint( elem, x, length * rule_wts[qp], integ1, integ2 ); + } return; } constexpr int max_qpts = max_symmetric_triangle_qpts; RealT rule_wts[max_qpts] = { 0. }; RealT rule_coords[2 * max_qpts] = { 0. }; - const int num_qpts = GetCommonPlaneTriangleRule( tri_order, rule_wts, rule_coords ); + const int num_qpts = GetCommonPlaneTriangleRule( quadrature_order, rule_wts, rule_coords ); RealT centroid[3]; GetCommonPlaneOverlapCentroid( elem, centroid ); @@ -549,7 +708,7 @@ TRIBOL_HOST_DEVICE inline void EvalWeakFormIntegralCommonPlaneFullTri( SurfaceCo } TRIBOL_HOST_DEVICE inline void EvalWeakFormIntegralCommonPlane( SurfaceContactElem const& elem, const PolyInteg rule, - const int tri_order, RealT* const integ1, + const int quadrature_order, RealT* const integ1, RealT* const integ2 ) { switch ( rule ) { @@ -559,8 +718,8 @@ TRIBOL_HOST_DEVICE inline void EvalWeakFormIntegralCommonPlane( SurfaceContactEl AccumulateCommonPlaneIntegralAtPoint( elem, cx, 1.0, integ1, integ2 ); break; } - case FULL_TRI_DECOMP: - EvalWeakFormIntegralCommonPlaneFullTri( elem, tri_order, integ1, integ2 ); + case MULTI_POINT: + EvalWeakFormIntegralCommonPlaneMultiPoint( elem, quadrature_order, integ1, integ2 ); break; default: #ifdef TRIBOL_USE_HOST diff --git a/src/tribol/interface/tribol.cpp b/src/tribol/interface/tribol.cpp index eb7b95a2..f8640026 100644 --- a/src/tribol/interface/tribol.cpp +++ b/src/tribol/interface/tribol.cpp @@ -97,7 +97,7 @@ void setPenaltyOptions( IndexT cs_id, PenaltyConstraintType pen_enfrc_option, } // end setPenaltyOptions() //------------------------------------------------------------------------------ -void setCommonPlaneIntegrationOptions( IndexT cs_id, PolyInteg rule, int triangle_order ) +void setCommonPlaneIntegrationOptions( IndexT cs_id, PolyInteg rule, int quadrature_order ) { auto cs = CouplingSchemeManager::getInstance().findData( cs_id ); @@ -111,11 +111,11 @@ void setCommonPlaneIntegrationOptions( IndexT cs_id, PolyInteg rule, int triangl return; } - SLIC_ERROR_ROOT_IF( triangle_order < 2 || triangle_order > 10, - "tribol::setCommonPlaneIntegrationOptions(): triangle quadrature order must be in [2,10]." ); + SLIC_ERROR_ROOT_IF( quadrature_order < 2 || quadrature_order > 10, + "tribol::setCommonPlaneIntegrationOptions(): CommonPlane quadrature order must be in [2,10]." ); penalty_options.common_plane_rule = rule; - penalty_options.common_plane_triangle_order = triangle_order; + penalty_options.common_plane_quadrature_order = quadrature_order; } //------------------------------------------------------------------------------ diff --git a/src/tribol/interface/tribol.hpp b/src/tribol/interface/tribol.hpp index cd33b716..991c621e 100644 --- a/src/tribol/interface/tribol.hpp +++ b/src/tribol/interface/tribol.hpp @@ -70,10 +70,10 @@ void setPenaltyOptions( IndexT cs_id, PenaltyConstraintType pen_enfrc_option, * * \param [in] cs_id coupling scheme id * \param [in] rule polygon integration rule for CommonPlane force integration - * \param [in] triangle_order order of the triangle quadrature used by FULL_TRI_DECOMP in the range [2,10] + * \param [in] quadrature_order order of the CommonPlane quadrature used by MULTI_POINT in the range [2,10] * \pre user must register coupling scheme prior to setting CommonPlane integration options */ -void setCommonPlaneIntegrationOptions( IndexT cs_id, PolyInteg rule, int triangle_order = 3 ); +void setCommonPlaneIntegrationOptions( IndexT cs_id, PolyInteg rule, int quadrature_order = 3 ); /*! * \brief Sets the constant kinematic penalty stiffness diff --git a/src/tribol/physics/CommonPlane.cpp b/src/tribol/physics/CommonPlane.cpp index c3780cc0..7c7482b1 100644 --- a/src/tribol/physics/CommonPlane.cpp +++ b/src/tribol/physics/CommonPlane.cpp @@ -250,6 +250,53 @@ TRIBOL_HOST_DEVICE inline bool EvalLinearFaceAtProjectedPoint( const RealT* face return true; } +TRIBOL_HOST_DEVICE inline bool EvalLinearEdgeAtProjectedPoint( const RealT* edge_coords, const RealT x_query[2], + const RealT projection_dir[2], RealT x_edge[3], RealT* phi, + const int value_dim = 0, const RealT* nodal_vals = nullptr, + RealT* values = nullptr ) +{ + const RealT ax = edge_coords[0]; + const RealT ay = edge_coords[1]; + const RealT bx = edge_coords[2]; + const RealT by = edge_coords[3]; + const RealT ex = bx - ax; + const RealT ey = by - ay; + + const RealT det = projection_dir[0] * ey - ex * projection_dir[1]; + constexpr RealT det_tol = 1.e-14; + if ( std::abs( det ) <= det_tol ) { + return false; + } + + const RealT rhs_x = x_query[0] - ax; + const RealT rhs_y = x_query[1] - ay; + RealT edge_param = ( projection_dir[0] * rhs_y - projection_dir[1] * rhs_x ) / det; + + constexpr RealT edge_tol = 1.e-8; + if ( edge_param < -edge_tol || edge_param > 1. + edge_tol ) { + return false; + } + edge_param = std::max( 0., std::min( 1., edge_param ) ); + + phi[0] = 1. - edge_param; + phi[1] = edge_param; + + x_edge[0] = ax + edge_param * ex; + x_edge[1] = ay + edge_param * ey; + x_edge[2] = 0.; + + if ( values != nullptr ) { + initRealArray( values, value_dim, 0. ); + for ( int a = 0; a < 2; ++a ) { + for ( int i = 0; i < value_dim; ++i ) { + values[i] += nodal_vals[i + a * value_dim] * phi[a]; + } + } + } + + return true; +} + TRIBOL_HOST_DEVICE inline void AccumulateContactForce( const MeshData::Viewer& mesh1, const MeshData::Viewer& mesh2, const IndexT index1, const IndexT index2, const int dim, const int num_nodes_per_face, const RealT force_x, @@ -417,8 +464,8 @@ int ApplyNormal( CouplingScheme* cs ) RealT penalty_stiff_per_area{ 0. }; auto& enforcement_options = cs_view.getEnforcementOptions(); const PenaltyEnforcementOptions& pen_enfrc_options = enforcement_options.penalty_options; - const bool use_full_tri = pen_enfrc_options.common_plane_rule == FULL_TRI_DECOMP && cs_view.spatialDimension() == 3; - if ( !use_full_tri && gap > gap_tol ) { + const bool use_multi_point = pen_enfrc_options.common_plane_rule == MULTI_POINT; + if ( !use_multi_point && gap > gap_tol ) { // We are here if we have a pair that passes ALL geometric // filter checks, BUT does not actually violate this method's // gap constraint. @@ -546,20 +593,7 @@ int ApplyNormal( CouplingScheme* cs ) cntctElem.overlapNormal = overlapNormal; cntctElem.overlapArea = plane.m_area; - // create arrays to hold nodal residual weak form integral evaluations - RealT phi1[max_nodes_per_face]; - RealT phi2[max_nodes_per_face]; - initRealArray( phi1, num_nodes_per_face, 0. ); - initRealArray( phi2, num_nodes_per_face, 0. ); - - //////////////////////////////////////////////////////////////////////// - // Integration of contact integrals: integral of shape functions over // - // contact overlap patch // - //////////////////////////////////////////////////////////////////////// - EvalWeakFormIntegralCommonPlane( cntctElem, pen_enfrc_options.common_plane_rule, - pen_enfrc_options.common_plane_triangle_order, phi1, phi2 ); - - if ( use_full_tri ) { + if ( use_multi_point ) { StackArrayT actual_xf1; StackArrayT actual_xf2; mesh1.getFaceCoords( index1, actual_xf1 ); @@ -579,41 +613,107 @@ int ApplyNormal( CouplingScheme* cs ) constexpr int max_qpts = max_symmetric_triangle_qpts; RealT rule_wts[max_qpts] = { 0. }; RealT rule_coords[2 * max_qpts] = { 0. }; - const int num_qpts = GetCommonPlaneTriangleRule( pen_enfrc_options.common_plane_triangle_order, rule_wts, rule_coords ); + bool has_contact = false; - RealT centroid[3]; - GetCommonPlaneOverlapCentroid( cntctElem, centroid ); + if ( dim == 3 ) { + const int num_qpts = + GetCommonPlaneTriangleRule( pen_enfrc_options.common_plane_quadrature_order, rule_wts, rule_coords ); + + RealT centroid[3]; + GetCommonPlaneOverlapCentroid( cntctElem, centroid ); + + RealT xTri[3]; + RealT yTri[3]; + RealT zTri[3]; + + for ( int j = 0; j < numPolyVert; ++j ) { + const int next = ( j == numPolyVert - 1 ) ? 0 : j + 1; + xTri[0] = xVert[dim * j]; + yTri[0] = xVert[dim * j + 1]; + zTri[0] = xVert[dim * j + 2]; + xTri[1] = xVert[dim * next]; + yTri[1] = xVert[dim * next + 1]; + zTri[1] = xVert[dim * next + 2]; + xTri[2] = centroid[0]; + yTri[2] = centroid[1]; + zTri[2] = centroid[2]; + + const RealT area = Area3DTri( xTri, yTri, zTri ); + if ( area <= 0. ) { + continue; + } - RealT xTri[3]; - RealT yTri[3]; - RealT zTri[3]; - bool has_contact = false; + for ( int qp = 0; qp < num_qpts; ++qp ) { + const RealT xi = rule_coords[2 * qp]; + const RealT eta = rule_coords[2 * qp + 1]; + const RealT n0 = 1. - xi - eta; + RealT x_q[3]; + x_q[0] = n0 * xTri[0] + xi * xTri[1] + eta * xTri[2]; + x_q[1] = n0 * yTri[0] + xi * yTri[1] + eta * yTri[2]; + x_q[2] = n0 * zTri[0] + xi * zTri[1] + eta * zTri[2]; + + RealT phi_q1[max_nodes_per_face] = { 0., 0., 0., 0. }; + RealT phi_q2[max_nodes_per_face] = { 0., 0., 0., 0. }; + RealT x_qf1[max_dim]; + RealT x_qf2[max_dim]; + RealT vel_q1[max_dim] = { 0., 0., 0. }; + RealT vel_q2[max_dim] = { 0., 0., 0. }; + + const bool mapped_face1 = + EvalLinearFaceAtProjectedPoint( actual_xf1, num_nodes_per_face, x_q, overlapNormal, x_qf1, phi_q1, + use_rate ? dim : 0, use_rate ? &actual_vf1[0] : nullptr, + use_rate ? vel_q1 : nullptr ); + const bool mapped_face2 = + EvalLinearFaceAtProjectedPoint( actual_xf2, num_nodes_per_face, x_q, overlapNormal, x_qf2, phi_q2, + use_rate ? dim : 0, use_rate ? &actual_vf2[0] : nullptr, + use_rate ? vel_q2 : nullptr ); + if ( !mapped_face1 || !mapped_face2 ) { + continue; + } + + RealT local_gap = + ( x_qf1[0] - x_qf2[0] ) * overlapNormal[0] + ( x_qf1[1] - x_qf2[1] ) * overlapNormal[1] + + ( x_qf1[2] - x_qf2[2] ) * overlapNormal[2]; + if ( local_gap > gap_tol ) { + continue; + } + + has_contact = true; - for ( int j = 0; j < numPolyVert; ++j ) { - const int next = ( j == numPolyVert - 1 ) ? 0 : j + 1; - xTri[0] = xVert[dim * j]; - yTri[0] = xVert[dim * j + 1]; - zTri[0] = xVert[dim * j + 2]; - xTri[1] = xVert[dim * next]; - yTri[1] = xVert[dim * next + 1]; - zTri[1] = xVert[dim * next + 2]; - xTri[2] = centroid[0]; - yTri[2] = centroid[1]; - zTri[2] = centroid[2]; - - const RealT area = Area3DTri( xTri, yTri, zTri ); - if ( area <= 0. ) { - continue; + RealT local_pressure = local_gap * penalty_stiff_per_area; + if ( use_rate && rate_penalty > 0. ) { + RealT local_vel_gap = + ( vel_q1[0] - vel_q2[0] ) * overlapNormal[0] + ( vel_q1[1] - vel_q2[1] ) * overlapNormal[1] + + ( vel_q1[2] - vel_q2[2] ) * overlapNormal[2]; + if ( local_vel_gap <= 0. ) { + local_pressure += local_vel_gap * rate_penalty; + } + } + + const RealT weighted_force = area * rule_wts[qp] * local_pressure; + const RealT force_x = overlapNormal[0] * weighted_force; + const RealT force_y = overlapNormal[1] * weighted_force; + const RealT force_z = overlapNormal[2] * weighted_force; + + AccumulateContactForce( mesh1, mesh2, index1, index2, dim, num_nodes_per_face, force_x, force_y, force_z, + phi_q1, phi_q2 ); + } } + } else { + RealT segment_rule_wts[max_segment_gauss_legendre_qpts] = { 0. }; + RealT segment_rule_coords[max_segment_gauss_legendre_qpts] = { 0. }; + const int num_qpts = GetCommonPlaneSegmentRule( pen_enfrc_options.common_plane_quadrature_order, + segment_rule_wts, segment_rule_coords ); + const RealT x0 = xVert[0]; + const RealT y0 = xVert[1]; + const RealT x1 = xVert[2]; + const RealT y1 = xVert[3]; + const RealT length = magnitude( x1 - x0, y1 - y0 ); for ( int qp = 0; qp < num_qpts; ++qp ) { - const RealT xi = rule_coords[2 * qp]; - const RealT eta = rule_coords[2 * qp + 1]; - const RealT n0 = 1. - xi - eta; - RealT x_q[3]; - x_q[0] = n0 * xTri[0] + xi * xTri[1] + eta * xTri[2]; - x_q[1] = n0 * yTri[0] + xi * yTri[1] + eta * yTri[2]; - x_q[2] = n0 * zTri[0] + xi * zTri[1] + eta * zTri[2]; + const RealT s = segment_rule_coords[qp]; + const RealT one_minus_s = 1. - s; + RealT x_q[2] = { one_minus_s * x0 + s * x1, one_minus_s * y0 + s * y1 }; RealT phi_q1[max_nodes_per_face] = { 0., 0., 0., 0. }; RealT phi_q2[max_nodes_per_face] = { 0., 0., 0., 0. }; @@ -623,19 +723,17 @@ int ApplyNormal( CouplingScheme* cs ) RealT vel_q2[max_dim] = { 0., 0., 0. }; const bool mapped_face1 = - EvalLinearFaceAtProjectedPoint( actual_xf1, num_nodes_per_face, x_q, overlapNormal, x_qf1, phi_q1, - use_rate ? dim : 0, use_rate ? &actual_vf1[0] : nullptr, - use_rate ? vel_q1 : nullptr ); + EvalLinearEdgeAtProjectedPoint( actual_xf1, x_q, overlapNormal, x_qf1, phi_q1, use_rate ? dim : 0, + use_rate ? &actual_vf1[0] : nullptr, use_rate ? vel_q1 : nullptr ); const bool mapped_face2 = - EvalLinearFaceAtProjectedPoint( actual_xf2, num_nodes_per_face, x_q, overlapNormal, x_qf2, phi_q2, - use_rate ? dim : 0, use_rate ? &actual_vf2[0] : nullptr, - use_rate ? vel_q2 : nullptr ); + EvalLinearEdgeAtProjectedPoint( actual_xf2, x_q, overlapNormal, x_qf2, phi_q2, use_rate ? dim : 0, + use_rate ? &actual_vf2[0] : nullptr, use_rate ? vel_q2 : nullptr ); if ( !mapped_face1 || !mapped_face2 ) { continue; } - RealT local_gap = ( x_qf1[0] - x_qf2[0] ) * overlapNormal[0] + ( x_qf1[1] - x_qf2[1] ) * overlapNormal[1] + - ( x_qf1[2] - x_qf2[2] ) * overlapNormal[2]; + RealT local_gap = + ( x_qf1[0] - x_qf2[0] ) * overlapNormal[0] + ( x_qf1[1] - x_qf2[1] ) * overlapNormal[1]; if ( local_gap > gap_tol ) { continue; } @@ -645,20 +743,18 @@ int ApplyNormal( CouplingScheme* cs ) RealT local_pressure = local_gap * penalty_stiff_per_area; if ( use_rate && rate_penalty > 0. ) { RealT local_vel_gap = - ( vel_q1[0] - vel_q2[0] ) * overlapNormal[0] + ( vel_q1[1] - vel_q2[1] ) * overlapNormal[1] + - ( vel_q1[2] - vel_q2[2] ) * overlapNormal[2]; + ( vel_q1[0] - vel_q2[0] ) * overlapNormal[0] + ( vel_q1[1] - vel_q2[1] ) * overlapNormal[1]; if ( local_vel_gap <= 0. ) { local_pressure += local_vel_gap * rate_penalty; } } - const RealT weighted_force = area * rule_wts[qp] * local_pressure; + const RealT weighted_force = length * segment_rule_wts[qp] * local_pressure; const RealT force_x = overlapNormal[0] * weighted_force; const RealT force_y = overlapNormal[1] * weighted_force; - const RealT force_z = overlapNormal[2] * weighted_force; - AccumulateContactForce( mesh1, mesh2, index1, index2, dim, num_nodes_per_face, force_x, force_y, force_z, - phi_q1, phi_q2 ); + AccumulateContactForce( mesh1, mesh2, index1, index2, dim, num_nodes_per_face, force_x, force_y, 0., phi_q1, + phi_q2 ); } } @@ -666,6 +762,19 @@ int ApplyNormal( CouplingScheme* cs ) return; } + // create arrays to hold nodal residual weak form integral evaluations + RealT phi1[max_nodes_per_face]; + RealT phi2[max_nodes_per_face]; + initRealArray( phi1, num_nodes_per_face, 0. ); + initRealArray( phi2, num_nodes_per_face, 0. ); + + //////////////////////////////////////////////////////////////////////// + // Integration of contact integrals: integral of shape functions over // + // contact overlap patch // + //////////////////////////////////////////////////////////////////////// + EvalWeakFormIntegralCommonPlane( cntctElem, pen_enfrc_options.common_plane_rule, + pen_enfrc_options.common_plane_quadrature_order, phi1, phi2 ); + /////////////////////////////////////////////////////////////////////// // Computation of full contact nodal force contributions // // (i.e. premultiplication of contact integrals by normal component, // @@ -730,7 +839,7 @@ int ApplyTangential( CouplingScheme* IndexT index1 = plane.getCpElementId1(); IndexT index2 = plane.getCpElementId2(); - const bool use_full_tri = pen_enfrc_options.common_plane_rule == FULL_TRI_DECOMP && dim == 3; + const bool use_multi_point = pen_enfrc_options.common_plane_rule == MULTI_POINT; // compute the velocity gap and pressure contribution StackArrayT x1; @@ -827,20 +936,7 @@ int ApplyTangential( CouplingScheme* cntctElem.overlapNormal = overlapNormal; cntctElem.overlapArea = plane.m_area; - // create arrays to hold nodal residual weak form integral evaluations - RealT phi1[max_nodes_per_face]; - RealT phi2[max_nodes_per_face]; - initRealArray( phi1, numNodesPerFace1, 0. ); - initRealArray( phi2, numNodesPerFace2, 0. ); - - //////////////////////////////////////////////////////////////////////// - // Integration of contact integrals: integral of shape functions over // - // contact overlap patch // - //////////////////////////////////////////////////////////////////////// - EvalWeakFormIntegralCommonPlane( cntctElem, pen_enfrc_options.common_plane_rule, - pen_enfrc_options.common_plane_triangle_order, phi1, phi2 ); - - if ( use_full_tri ) { + if ( use_multi_point ) { StackArrayT actual_xf1; StackArrayT actual_xf2; StackArrayT actual_vf1; @@ -853,62 +949,125 @@ int ApplyTangential( CouplingScheme* constexpr int max_qpts = max_symmetric_triangle_qpts; RealT rule_wts[max_qpts] = { 0. }; RealT rule_coords[2 * max_qpts] = { 0. }; - const int num_qpts = GetCommonPlaneTriangleRule( pen_enfrc_options.common_plane_triangle_order, rule_wts, rule_coords ); - - RealT centroid[3]; - GetCommonPlaneOverlapCentroid( cntctElem, centroid ); - - RealT xTri[3]; - RealT yTri[3]; - RealT zTri[3]; const RealT visc = 0.5 * ( mesh1.getElementData().m_viscous_damping_coeff + mesh2.getElementData().m_viscous_damping_coeff ); - for ( int j = 0; j < numPolyVert; ++j ) { - const int next = ( j == numPolyVert - 1 ) ? 0 : j + 1; - xTri[0] = xVert[dim * j]; - yTri[0] = xVert[dim * j + 1]; - zTri[0] = xVert[dim * j + 2]; - xTri[1] = xVert[dim * next]; - yTri[1] = xVert[dim * next + 1]; - zTri[1] = xVert[dim * next + 2]; - xTri[2] = centroid[0]; - yTri[2] = centroid[1]; - zTri[2] = centroid[2]; - - const RealT area = Area3DTri( xTri, yTri, zTri ); - if ( area <= 0. ) { - continue; + if ( dim == 3 ) { + const int num_qpts = + GetCommonPlaneTriangleRule( pen_enfrc_options.common_plane_quadrature_order, rule_wts, rule_coords ); + + RealT centroid[3]; + GetCommonPlaneOverlapCentroid( cntctElem, centroid ); + + RealT xTri[3]; + RealT yTri[3]; + RealT zTri[3]; + + for ( int j = 0; j < numPolyVert; ++j ) { + const int next = ( j == numPolyVert - 1 ) ? 0 : j + 1; + xTri[0] = xVert[dim * j]; + yTri[0] = xVert[dim * j + 1]; + zTri[0] = xVert[dim * j + 2]; + xTri[1] = xVert[dim * next]; + yTri[1] = xVert[dim * next + 1]; + zTri[1] = xVert[dim * next + 2]; + xTri[2] = centroid[0]; + yTri[2] = centroid[1]; + zTri[2] = centroid[2]; + + const RealT area = Area3DTri( xTri, yTri, zTri ); + if ( area <= 0. ) { + continue; + } + + for ( int qp = 0; qp < num_qpts; ++qp ) { + const RealT xi = rule_coords[2 * qp]; + const RealT eta = rule_coords[2 * qp + 1]; + const RealT n0 = 1. - xi - eta; + RealT x_q[3]; + x_q[0] = n0 * xTri[0] + xi * xTri[1] + eta * xTri[2]; + x_q[1] = n0 * yTri[0] + xi * yTri[1] + eta * yTri[2]; + x_q[2] = n0 * zTri[0] + xi * zTri[1] + eta * zTri[2]; + + RealT phi_q1[max_nodes_per_face] = { 0., 0., 0., 0. }; + RealT phi_q2[max_nodes_per_face] = { 0., 0., 0., 0. }; + RealT x_qf1[max_dim]; + RealT x_qf2[max_dim]; + RealT vel_q1[max_dim]; + RealT vel_q2[max_dim]; + + const bool mapped_face1 = + EvalLinearFaceAtProjectedPoint( actual_xf1, numNodesPerFace1, x_q, overlapNormal, x_qf1, phi_q1, dim, + &actual_vf1[0], vel_q1 ); + const bool mapped_face2 = + EvalLinearFaceAtProjectedPoint( actual_xf2, numNodesPerFace2, x_q, overlapNormal, x_qf2, phi_q2, dim, + &actual_vf2[0], vel_q2 ); + if ( !mapped_face1 || !mapped_face2 ) { + continue; + } + + RealT local_gap = + ( x_qf1[0] - x_qf2[0] ) * overlapNormal[0] + ( x_qf1[1] - x_qf2[1] ) * overlapNormal[1] + + ( x_qf1[2] - x_qf2[2] ) * overlapNormal[2]; + RealT gap_tol = cs_view.getGapTol( index1, index2 ); + if ( local_gap > gap_tol ) { + continue; + } + + RealT velGap_q[max_dim]; + velGap_q[0] = vel_q1[0] - vel_q2[0]; + velGap_q[1] = vel_q1[1] - vel_q2[1]; + velGap_q[2] = vel_q1[2] - vel_q2[2]; + + RealT velGap_dot_n = + velGap_q[0] * overlapNormal[0] + velGap_q[1] * overlapNormal[1] + velGap_q[2] * overlapNormal[2]; + RealT velGapTan[max_dim]; + velGapTan[0] = velGap_q[0] - velGap_dot_n * overlapNormal[0]; + velGapTan[1] = velGap_q[1] - velGap_dot_n * overlapNormal[1]; + velGapTan[2] = velGap_q[2] - velGap_dot_n * overlapNormal[2]; + + const RealT weighted_force = area * rule_wts[qp] * visc; + const RealT force_x = weighted_force * velGapTan[0]; + const RealT force_y = weighted_force * velGapTan[1]; + const RealT force_z = weighted_force * velGapTan[2]; + + AccumulateContactForce( mesh1, mesh2, index1, index2, dim, numNodesPerFace1, force_x, force_y, force_z, + phi_q1, phi_q2 ); + } } + } else { + RealT segment_rule_wts[max_segment_gauss_legendre_qpts] = { 0. }; + RealT segment_rule_coords[max_segment_gauss_legendre_qpts] = { 0. }; + const int num_qpts = GetCommonPlaneSegmentRule( pen_enfrc_options.common_plane_quadrature_order, + segment_rule_wts, segment_rule_coords ); + const RealT x0 = xVert[0]; + const RealT y0 = xVert[1]; + const RealT x1_seg = xVert[2]; + const RealT y1_seg = xVert[3]; + const RealT length = magnitude( x1_seg - x0, y1_seg - y0 ); for ( int qp = 0; qp < num_qpts; ++qp ) { - const RealT xi = rule_coords[2 * qp]; - const RealT eta = rule_coords[2 * qp + 1]; - const RealT n0 = 1. - xi - eta; - RealT x_q[3]; - x_q[0] = n0 * xTri[0] + xi * xTri[1] + eta * xTri[2]; - x_q[1] = n0 * yTri[0] + xi * yTri[1] + eta * yTri[2]; - x_q[2] = n0 * zTri[0] + xi * zTri[1] + eta * zTri[2]; + const RealT s = segment_rule_coords[qp]; + const RealT one_minus_s = 1. - s; + RealT x_q[2] = { one_minus_s * x0 + s * x1_seg, one_minus_s * y0 + s * y1_seg }; RealT phi_q1[max_nodes_per_face] = { 0., 0., 0., 0. }; RealT phi_q2[max_nodes_per_face] = { 0., 0., 0., 0. }; RealT x_qf1[max_dim]; RealT x_qf2[max_dim]; - RealT vel_q1[max_dim]; - RealT vel_q2[max_dim]; + RealT vel_q1[max_dim] = { 0., 0., 0. }; + RealT vel_q2[max_dim] = { 0., 0., 0. }; const bool mapped_face1 = - EvalLinearFaceAtProjectedPoint( actual_xf1, numNodesPerFace1, x_q, overlapNormal, x_qf1, phi_q1, dim, - &actual_vf1[0], vel_q1 ); + EvalLinearEdgeAtProjectedPoint( actual_xf1, x_q, overlapNormal, x_qf1, phi_q1, dim, &actual_vf1[0], vel_q1 ); const bool mapped_face2 = - EvalLinearFaceAtProjectedPoint( actual_xf2, numNodesPerFace2, x_q, overlapNormal, x_qf2, phi_q2, dim, - &actual_vf2[0], vel_q2 ); + EvalLinearEdgeAtProjectedPoint( actual_xf2, x_q, overlapNormal, x_qf2, phi_q2, dim, &actual_vf2[0], vel_q2 ); if ( !mapped_face1 || !mapped_face2 ) { continue; } - RealT local_gap = ( x_qf1[0] - x_qf2[0] ) * overlapNormal[0] + ( x_qf1[1] - x_qf2[1] ) * overlapNormal[1] + - ( x_qf1[2] - x_qf2[2] ) * overlapNormal[2]; + RealT local_gap = + ( x_qf1[0] - x_qf2[0] ) * overlapNormal[0] + ( x_qf1[1] - x_qf2[1] ) * overlapNormal[1]; RealT gap_tol = cs_view.getGapTol( index1, index2 ); if ( local_gap > gap_tol ) { continue; @@ -917,27 +1076,36 @@ int ApplyTangential( CouplingScheme* RealT velGap_q[max_dim]; velGap_q[0] = vel_q1[0] - vel_q2[0]; velGap_q[1] = vel_q1[1] - vel_q2[1]; - velGap_q[2] = vel_q1[2] - vel_q2[2]; - RealT velGap_dot_n = velGap_q[0] * overlapNormal[0] + velGap_q[1] * overlapNormal[1] + velGap_q[2] * overlapNormal[2]; - RealT velGapTan[max_dim]; - velGapTan[0] = velGap_q[0] - velGap_dot_n * overlapNormal[0]; - velGapTan[1] = velGap_q[1] - velGap_dot_n * overlapNormal[1]; - velGapTan[2] = velGap_q[2] - velGap_dot_n * overlapNormal[2]; + RealT velGap_dot_n = velGap_q[0] * overlapNormal[0] + velGap_q[1] * overlapNormal[1]; + RealT velGapTan_q[max_dim] = { velGap_q[0] - velGap_dot_n * overlapNormal[0], + velGap_q[1] - velGap_dot_n * overlapNormal[1], 0. }; - const RealT weighted_force = area * rule_wts[qp] * visc; - const RealT force_x = weighted_force * velGapTan[0]; - const RealT force_y = weighted_force * velGapTan[1]; - const RealT force_z = weighted_force * velGapTan[2]; + const RealT weighted_force = length * segment_rule_wts[qp] * visc; + const RealT force_x = weighted_force * velGapTan_q[0]; + const RealT force_y = weighted_force * velGapTan_q[1]; - AccumulateContactForce( mesh1, mesh2, index1, index2, dim, numNodesPerFace1, force_x, force_y, force_z, - phi_q1, phi_q2 ); + AccumulateContactForce( mesh1, mesh2, index1, index2, dim, numNodesPerFace1, force_x, force_y, 0., phi_q1, + phi_q2 ); } } return; } + // create arrays to hold nodal residual weak form integral evaluations + RealT phi1[max_nodes_per_face]; + RealT phi2[max_nodes_per_face]; + initRealArray( phi1, numNodesPerFace1, 0. ); + initRealArray( phi2, numNodesPerFace2, 0. ); + + //////////////////////////////////////////////////////////////////////// + // Integration of contact integrals: integral of shape functions over // + // contact overlap patch // + //////////////////////////////////////////////////////////////////////// + EvalWeakFormIntegralCommonPlane( cntctElem, pen_enfrc_options.common_plane_rule, + pen_enfrc_options.common_plane_quadrature_order, phi1, phi2 ); + ///////////////////////////////////////////////////// // Computation of tangential viscous damping force // ///////////////////////////////////////////////////// diff --git a/src/tribol/utils/TestUtils.cpp b/src/tribol/utils/TestUtils.cpp index 22058f8d..9acdc882 100644 --- a/src/tribol/utils/TestUtils.cpp +++ b/src/tribol/utils/TestUtils.cpp @@ -1142,7 +1142,7 @@ int TestMesh::tribolSetupAndUpdate( ContactMethod method, EnforcementMethod enfo // set penalty options after registering coupling scheme setPenaltyOptions( csIndex, constraint_type, pen_calc, rate_calc ); - setCommonPlaneIntegrationOptions( csIndex, params.common_plane_rule, params.common_plane_triangle_order ); + setCommonPlaneIntegrationOptions( csIndex, params.common_plane_rule, params.common_plane_quadrature_order ); } else if ( ( method == SINGLE_MORTAR || method == ALIGNED_MORTAR ) && enforcement == LAGRANGE_MULTIPLIER ) { // note, eval modes and sparse modes not exposed in the interface to this class diff --git a/src/tribol/utils/TestUtils.hpp b/src/tribol/utils/TestUtils.hpp index 3dfe07f9..5b1c269b 100644 --- a/src/tribol/utils/TestUtils.hpp +++ b/src/tribol/utils/TestUtils.hpp @@ -45,7 +45,7 @@ struct TestControlParameters { constant_rate_penalty( false ), percent_rate_penalty( false ), common_plane_rule( SINGLE_POINT ), - common_plane_triangle_order( 3 ), + common_plane_quadrature_order( 3 ), rate_penalty( 1.0 ), rate_penalty_ratio( 0.0 ), const_penalty( 1.0 ), @@ -68,7 +68,7 @@ struct TestControlParameters { bool constant_rate_penalty; bool percent_rate_penalty; PolyInteg common_plane_rule; - int common_plane_triangle_order; + int common_plane_quadrature_order; RealT rate_penalty; RealT rate_penalty_ratio; RealT const_penalty; From af60b3eacb88d9221e3cd6d9c835f95c476e4687 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Fri, 29 May 2026 23:01:26 -0700 Subject: [PATCH 5/6] formatting --- src/tests/tribol_common_plane_penalty.cpp | 2 - src/tests/tribol_iso_integ.cpp | 14 +-- src/tribol/integ/Integration.hpp | 126 ++++++++++++---------- src/tribol/physics/CommonPlane.cpp | 113 +++++++++---------- 4 files changed, 128 insertions(+), 127 deletions(-) diff --git a/src/tests/tribol_common_plane_penalty.cpp b/src/tests/tribol_common_plane_penalty.cpp index add16f68..d6038016 100644 --- a/src/tests/tribol_common_plane_penalty.cpp +++ b/src/tests/tribol_common_plane_penalty.cpp @@ -185,8 +185,6 @@ class CommonPlaneTest : public ::testing::Test { void SetUp() override {} void TearDown() override { this->m_mesh.clear(); } - -protected: }; struct WarpedQuadForceResult { diff --git a/src/tests/tribol_iso_integ.cpp b/src/tests/tribol_iso_integ.cpp index 62fa9d4e..c9b9bd70 100644 --- a/src/tests/tribol_iso_integ.cpp +++ b/src/tests/tribol_iso_integ.cpp @@ -29,22 +29,16 @@ RealT factorial( int n ) return result; } -RealT referenceTriangleMoment( int px, int py ) -{ - return factorial( px ) * factorial( py ) / factorial( px + py + 2 ); -} +RealT referenceTriangleMoment( int px, int py ) { return factorial( px ) * factorial( py ) / factorial( px + py + 2 ); } -RealT normalizedReferenceTriangleMoment( int px, int py ) -{ - return 2. * referenceTriangleMoment( px, py ); -} +RealT normalizedReferenceTriangleMoment( int px, int py ) { return 2. * referenceTriangleMoment( px, py ); } RealT evalTriangleRuleMoment( bool use_legacy, int order, int px, int py ) { RealT wts[tribol::max_symmetric_triangle_qpts] = { 0. }; RealT coords[2 * tribol::max_symmetric_triangle_qpts] = { 0. }; - const int num_qpts = - use_legacy ? tribol::GetLegacyTriangleRule( order, wts, coords ) : tribol::GetCommonPlaneTriangleRule( order, wts, coords ); + const int num_qpts = use_legacy ? tribol::GetLegacyTriangleRule( order, wts, coords ) + : tribol::GetCommonPlaneTriangleRule( order, wts, coords ); RealT value = 0.; for ( int qp = 0; qp < num_qpts; ++qp ) { diff --git a/src/tribol/integ/Integration.hpp b/src/tribol/integ/Integration.hpp index f61eb0da..bb53f752 100644 --- a/src/tribol/integ/Integration.hpp +++ b/src/tribol/integ/Integration.hpp @@ -295,8 +295,7 @@ namespace detail { * area is applied separately, so the imported weights are scaled by 1/2 * during expansion. */ -struct SymmetricTriangleRuleData -{ +struct SymmetricTriangleRuleData { int num_centroid_orbits; int num_edge_orbits; int num_general_orbits; @@ -307,62 +306,77 @@ struct SymmetricTriangleRuleData constexpr RealT symmetric_triangle_weight_scale = 0.5; constexpr RealT tri_deg2_weights[] = { 6.66666666666666666666666666666666635e-01 }; -constexpr RealT tri_deg2_orbits[] = { 1.66666666666666666666666666666666659e-01, 6.66666666666666666666666666666666635e-01 }; +constexpr RealT tri_deg2_orbits[] = { 1.66666666666666666666666666666666659e-01, + 6.66666666666666666666666666666666635e-01 }; -constexpr RealT tri_deg4_weights[] = { 4.46763179356022931390014016866245598e-01, 2.19903487310643735276652649800421061e-01 }; -constexpr RealT tri_deg4_orbits[] = { 4.45948490915964886318329253883051984e-01, 1.08103018168070227363341492233896033e-01, - 9.15762135097707434595714634022014804e-02, 8.16847572980458513080857073195597039e-01 }; +constexpr RealT tri_deg4_weights[] = { 4.46763179356022931390014016866245598e-01, + 2.19903487310643735276652649800421061e-01 }; +constexpr RealT tri_deg4_orbits[] = { + 4.45948490915964886318329253883051984e-01, 1.08103018168070227363341492233896033e-01, + 9.15762135097707434595714634022014804e-02, 8.16847572980458513080857073195597039e-01 }; -constexpr RealT tri_deg5_weights[] = { 4.50000000000000000000000000000000010e-01, 2.51878361089654305191367891000362687e-01, +constexpr RealT tri_deg5_weights[] = { 4.50000000000000000000000000000000010e-01, + 2.51878361089654305191367891000362687e-01, 2.64788305577012361475298775666303977e-01 }; -constexpr RealT tri_deg5_orbits[] = { 3.33333333333333333333333333333333317e-01, 1.01286507323456338800987361915123836e-01, - 7.97426985353087322398025276169752328e-01, 4.70142064105115089770441209513447613e-01, - 5.97158717897698204591175809731048219e-02 }; +constexpr RealT tri_deg5_orbits[] = { + 3.33333333333333333333333333333333317e-01, 1.01286507323456338800987361915123836e-01, + 7.97426985353087322398025276169752328e-01, 4.70142064105115089770441209513447613e-01, + 5.97158717897698204591175809731048219e-02 }; -constexpr RealT tri_deg6_weights[] = { 1.01689812740413633841873618213737963e-01, 2.33572551452758732050579222771158894e-01, +constexpr RealT tri_deg6_weights[] = { 1.01689812740413633841873618213737963e-01, + 2.33572551452758732050579222771158894e-01, 1.65702151236747150387106912840884901e-01 }; -constexpr RealT tri_deg6_orbits[] = { 6.30890144915022283403316028708191300e-02, 8.73821971016995543319336794258361644e-01, - 2.49286745170910421291638553107019076e-01, 5.01426509658179157416722893785961848e-01, - 5.31450498448169473532496716313981651e-02, 6.36502499121398647230142594412049640e-01, - 3.10352451033784405416607733956552146e-01 }; - -constexpr RealT tri_deg7_weights[] = { 3.30901002215842620719558969458348911e-02, 2.55888342460311145565802470369292636e-01, - 1.54173292371972135669643041667482776e-01, 1.11757465806399561679632628842028190e-01 }; -constexpr RealT tri_deg7_orbits[] = { 3.37306485545878487149717263008162317e-02, 9.32538702890824302570056547398367537e-01, - 2.41577382595403558950186769837781999e-01, 5.16845234809192882099626460324436002e-01, - 4.74309692504718234209580735949185780e-01, 5.13806149905635315808385281016284391e-02, - 4.70366446525952333414099753568849895e-02, 7.54280040550053177356239324628119970e-01, - 1.98683314797351589302350700014995040e-01 }; - -constexpr RealT tri_deg8_weights[] = { 2.88631215355574336502182220978129237e-01, 1.90183268534569249587792208777168633e-01, - 2.06434741069436500563583100584258068e-01, 6.49169952463961606218518566835611904e-02, - 5.44606283488699885296893801478178481e-02 }; -constexpr RealT tri_deg8_orbits[] = { 3.33333333333333333333333333333333317e-01, 4.59292588292723156028815514494169350e-01, - 8.14148234145536879423689710116613481e-02, 1.70569307751760206622293501491464506e-01, - 6.58861384496479586755412997017070988e-01, 5.05472283170309754584235505965989197e-02, - 8.98905543365938049083152898806802161e-01, 8.39477740995760533721383453929445768e-03, - 7.28492392955404281241000379176061966e-01, 2.63112829634638113421785786284643576e-01 }; - -constexpr RealT tri_deg9_weights[] = { 1.94271592565597667638483965014577269e-01, 1.55655082009548558633478712598807923e-01, - 1.59295477854420506065783548528090548e-01, 6.26694004542781410737096625744186273e-02, - 5.11553513173960625233575971179996460e-02, 8.65670787545787545787545787545787526e-02 }; -constexpr RealT tri_deg9_orbits[] = { 3.33333333333333333333333333333333317e-01, 4.37089591492936637269930364435354971e-01, - 1.25820817014126725460139271129290058e-01, 1.88203535619032730240961280467335542e-01, - 6.23592928761934539518077439065328819e-01, 4.89682519198737627783706924836192818e-01, - 2.06349616025247444325861503276144129e-02, 4.47295133944527098651065899662763588e-02, - 9.10540973211094580269786820067447282e-01, 3.68384120547362836348175987833851049e-02, - 7.41198598784498020690079873523423793e-01, 2.21962989160765695675102527693191078e-01 }; - -constexpr RealT tri_deg10_weights[] = { 1.63486658292571932856237369968355216e-01, 2.67059376262991325511459567981373070e-02, - 9.19159272094894560275758192650956353e-02, 1.27809812792848090865797467525306648e-01, - 6.83692963259188572573831680826845816e-02, 5.05955154145767687780855813656664345e-02 }; -constexpr RealT tri_deg10_orbits[] = { 3.33333333333333333333333333333333317e-01, 3.20553732169435129309845893364897379e-02, - 9.35889253566112974138030821327020524e-01, 1.42161101056564385092162103190958311e-01, - 7.15677797886871229815675793618083377e-01, 3.21812995288835421225097560986048687e-01, - 5.30054118927344028277095673945694069e-01, 1.48132885783820550497806765068257172e-01, - 2.96198894887297676338362694260427776e-02, 6.01233328683459245454742893458687815e-01, - 3.69146781827810986911420837115269408e-01, 2.83676653399384392504357555781301898e-02, - 8.07930600922879065079949902881744115e-01, 1.63701733737182495669614341540125695e-01 }; +constexpr RealT tri_deg6_orbits[] = { + 6.30890144915022283403316028708191300e-02, 8.73821971016995543319336794258361644e-01, + 2.49286745170910421291638553107019076e-01, 5.01426509658179157416722893785961848e-01, + 5.31450498448169473532496716313981651e-02, 6.36502499121398647230142594412049640e-01, + 3.10352451033784405416607733956552146e-01 }; + +constexpr RealT tri_deg7_weights[] = { + 3.30901002215842620719558969458348911e-02, 2.55888342460311145565802470369292636e-01, + 1.54173292371972135669643041667482776e-01, 1.11757465806399561679632628842028190e-01 }; +constexpr RealT tri_deg7_orbits[] = { + 3.37306485545878487149717263008162317e-02, 9.32538702890824302570056547398367537e-01, + 2.41577382595403558950186769837781999e-01, 5.16845234809192882099626460324436002e-01, + 4.74309692504718234209580735949185780e-01, 5.13806149905635315808385281016284391e-02, + 4.70366446525952333414099753568849895e-02, 7.54280040550053177356239324628119970e-01, + 1.98683314797351589302350700014995040e-01 }; + +constexpr RealT tri_deg8_weights[] = { + 2.88631215355574336502182220978129237e-01, 1.90183268534569249587792208777168633e-01, + 2.06434741069436500563583100584258068e-01, 6.49169952463961606218518566835611904e-02, + 5.44606283488699885296893801478178481e-02 }; +constexpr RealT tri_deg8_orbits[] = { + 3.33333333333333333333333333333333317e-01, 4.59292588292723156028815514494169350e-01, + 8.14148234145536879423689710116613481e-02, 1.70569307751760206622293501491464506e-01, + 6.58861384496479586755412997017070988e-01, 5.05472283170309754584235505965989197e-02, + 8.98905543365938049083152898806802161e-01, 8.39477740995760533721383453929445768e-03, + 7.28492392955404281241000379176061966e-01, 2.63112829634638113421785786284643576e-01 }; + +constexpr RealT tri_deg9_weights[] = { + 1.94271592565597667638483965014577269e-01, 1.55655082009548558633478712598807923e-01, + 1.59295477854420506065783548528090548e-01, 6.26694004542781410737096625744186273e-02, + 5.11553513173960625233575971179996460e-02, 8.65670787545787545787545787545787526e-02 }; +constexpr RealT tri_deg9_orbits[] = { + 3.33333333333333333333333333333333317e-01, 4.37089591492936637269930364435354971e-01, + 1.25820817014126725460139271129290058e-01, 1.88203535619032730240961280467335542e-01, + 6.23592928761934539518077439065328819e-01, 4.89682519198737627783706924836192818e-01, + 2.06349616025247444325861503276144129e-02, 4.47295133944527098651065899662763588e-02, + 9.10540973211094580269786820067447282e-01, 3.68384120547362836348175987833851049e-02, + 7.41198598784498020690079873523423793e-01, 2.21962989160765695675102527693191078e-01 }; + +constexpr RealT tri_deg10_weights[] = { + 1.63486658292571932856237369968355216e-01, 2.67059376262991325511459567981373070e-02, + 9.19159272094894560275758192650956353e-02, 1.27809812792848090865797467525306648e-01, + 6.83692963259188572573831680826845816e-02, 5.05955154145767687780855813656664345e-02 }; +constexpr RealT tri_deg10_orbits[] = { + 3.33333333333333333333333333333333317e-01, 3.20553732169435129309845893364897379e-02, + 9.35889253566112974138030821327020524e-01, 1.42161101056564385092162103190958311e-01, + 7.15677797886871229815675793618083377e-01, 3.21812995288835421225097560986048687e-01, + 5.30054118927344028277095673945694069e-01, 1.48132885783820550497806765068257172e-01, + 2.96198894887297676338362694260427776e-02, 6.01233328683459245454742893458687815e-01, + 3.69146781827810986911420837115269408e-01, 2.83676653399384392504357555781301898e-02, + 8.07930600922879065079949902881744115e-01, 1.63701733737182495669614341540125695e-01 }; TRIBOL_HOST_DEVICE inline bool GetSymmetricTriangleRuleData( int order, SymmetricTriangleRuleData& rule ) { @@ -397,7 +411,8 @@ TRIBOL_HOST_DEVICE inline bool GetSymmetricTriangleRuleData( int order, Symmetri } } -TRIBOL_HOST_DEVICE inline int ExpandSymmetricTriangleRule( const SymmetricTriangleRuleData& rule, RealT* wts, RealT* coords ) +TRIBOL_HOST_DEVICE inline int ExpandSymmetricTriangleRule( const SymmetricTriangleRuleData& rule, RealT* wts, + RealT* coords ) { int w_idx = 0; int c_idx = 0; @@ -482,7 +497,8 @@ TRIBOL_HOST_DEVICE inline int GetCommonPlaneTriangleRule( int order, RealT* wts, * \param [out] wts quadrature weights normalized so they sum to 1 on a triangle * \param [out] coords quadrature coordinates stored as stacked (xi, eta) pairs */ -TRIBOL_HOST_DEVICE inline int GetTriangleRule( int order, TriangleQuadratureRuleFamily family, RealT* wts, RealT* coords ) +TRIBOL_HOST_DEVICE inline int GetTriangleRule( int order, TriangleQuadratureRuleFamily family, RealT* wts, + RealT* coords ) { switch ( family ) { case TRI_RULE_LEGACY: diff --git a/src/tribol/physics/CommonPlane.cpp b/src/tribol/physics/CommonPlane.cpp index 7c7482b1..d185e6be 100644 --- a/src/tribol/physics/CommonPlane.cpp +++ b/src/tribol/physics/CommonPlane.cpp @@ -52,20 +52,17 @@ TRIBOL_HOST_DEVICE inline bool Solve3x3( const RealT A[3][3], const RealT b[3], } const RealT inv_detA = 1. / detA; - x[0] = inv_detA * ( b[0] * ( A[1][1] * A[2][2] - A[1][2] * A[2][1] ) - - A[0][1] * ( b[1] * A[2][2] - A[1][2] * b[2] ) + + x[0] = inv_detA * ( b[0] * ( A[1][1] * A[2][2] - A[1][2] * A[2][1] ) - A[0][1] * ( b[1] * A[2][2] - A[1][2] * b[2] ) + A[0][2] * ( b[1] * A[2][1] - A[1][1] * b[2] ) ); - x[1] = inv_detA * ( A[0][0] * ( b[1] * A[2][2] - A[1][2] * b[2] ) - - b[0] * ( A[1][0] * A[2][2] - A[1][2] * A[2][0] ) + + x[1] = inv_detA * ( A[0][0] * ( b[1] * A[2][2] - A[1][2] * b[2] ) - b[0] * ( A[1][0] * A[2][2] - A[1][2] * A[2][0] ) + A[0][2] * ( A[1][0] * b[2] - b[1] * A[2][0] ) ); - x[2] = inv_detA * ( A[0][0] * ( A[1][1] * b[2] - b[1] * A[2][1] ) - - A[0][1] * ( A[1][0] * b[2] - b[1] * A[2][0] ) + + x[2] = inv_detA * ( A[0][0] * ( A[1][1] * b[2] - b[1] * A[2][1] ) - A[0][1] * ( A[1][0] * b[2] - b[1] * A[2][0] ) + b[0] * ( A[1][0] * A[2][1] - A[1][1] * A[2][0] ) ); return true; } -TRIBOL_HOST_DEVICE inline void AccumulateFaceInterpolation( const RealT* face_coords, const int num_nodes, const RealT* phi, - RealT x_face[3], const int value_dim = 0, +TRIBOL_HOST_DEVICE inline void AccumulateFaceInterpolation( const RealT* face_coords, const int num_nodes, + const RealT* phi, RealT x_face[3], const int value_dim = 0, const RealT* nodal_vals = nullptr, RealT* values = nullptr ) { initRealArray( x_face, max_dim, 0. ); @@ -89,7 +86,8 @@ TRIBOL_HOST_DEVICE inline void AccumulateFaceInterpolation( const RealT* face_co TRIBOL_HOST_DEVICE inline bool EvalLinearFaceAtProjectedPoint( const RealT* face_coords, const int num_nodes, const RealT x_query[3], const RealT projection_dir[3], RealT x_face[3], RealT* phi, const int value_dim = 0, - const RealT* nodal_vals = nullptr, RealT* values = nullptr ) + const RealT* nodal_vals = nullptr, + RealT* values = nullptr ) { // CommonPlane quadrature points lie on the overlap polygon. Evaluate the face // fields at the corresponding on-face point found by projecting along the @@ -104,7 +102,8 @@ TRIBOL_HOST_DEVICE inline bool EvalLinearFaceAtProjectedPoint( const RealT* face RealT n_face[3]; crossProd( e1[0], e1[1], e1[2], e2[0], e2[1], e2[2], n_face[0], n_face[1], n_face[2] ); - const RealT denom = dotProd( n_face[0], n_face[1], n_face[2], projection_dir[0], projection_dir[1], projection_dir[2] ); + const RealT denom = + dotProd( n_face[0], n_face[1], n_face[2], projection_dir[0], projection_dir[1], projection_dir[2] ); constexpr RealT parallel_tol = 1.e-14; if ( std::abs( denom ) <= parallel_tol ) { return false; @@ -113,8 +112,7 @@ TRIBOL_HOST_DEVICE inline bool EvalLinearFaceAtProjectedPoint( const RealT* face const RealT dx0 = x0[0] - x_query[0]; const RealT dy0 = x0[1] - x_query[1]; const RealT dz0 = x0[2] - x_query[2]; - const RealT step = - dotProd( n_face[0], n_face[1], n_face[2], dx0, dy0, dz0 ) / denom; + const RealT step = dotProd( n_face[0], n_face[1], n_face[2], dx0, dy0, dz0 ) / denom; x_face[0] = x_query[0] + step * projection_dir[0]; x_face[1] = x_query[1] + step * projection_dir[1]; @@ -183,7 +181,8 @@ TRIBOL_HOST_DEVICE inline bool EvalLinearFaceAtProjectedPoint( const RealT* face dxdeta[2] += za * dphi_deta; } - RealT residual[3] = { x_face[0] - x_query[0] - s * projection_dir[0], x_face[1] - x_query[1] - s * projection_dir[1], + RealT residual[3] = { x_face[0] - x_query[0] - s * projection_dir[0], + x_face[1] - x_query[1] - s * projection_dir[1], x_face[2] - x_query[2] - s * projection_dir[2] }; const RealT residual_norm = magnitude( residual[0], residual[1], residual[2] ); if ( residual_norm <= residual_tol ) { @@ -222,15 +221,13 @@ TRIBOL_HOST_DEVICE inline bool EvalLinearFaceAtProjectedPoint( const RealT* face phi[3] = 0.25 * ( 1. + xi ) * ( 1. - eta ); AccumulateFaceInterpolation( face_coords, num_nodes, phi, x_face ); - const RealT line_residual = - magnitude( x_face[0] - x_query[0], x_face[1] - x_query[1], x_face[2] - x_query[2] ); - const RealT normal_step = - ( x_face[0] - x_query[0] ) * projection_dir[0] + ( x_face[1] - x_query[1] ) * projection_dir[1] + - ( x_face[2] - x_query[2] ) * projection_dir[2]; - const RealT projection_residual = - magnitude( x_face[0] - x_query[0] - normal_step * projection_dir[0], - x_face[1] - x_query[1] - normal_step * projection_dir[1], - x_face[2] - x_query[2] - normal_step * projection_dir[2] ); + const RealT line_residual = magnitude( x_face[0] - x_query[0], x_face[1] - x_query[1], x_face[2] - x_query[2] ); + const RealT normal_step = ( x_face[0] - x_query[0] ) * projection_dir[0] + + ( x_face[1] - x_query[1] ) * projection_dir[1] + + ( x_face[2] - x_query[2] ) * projection_dir[2]; + const RealT projection_residual = magnitude( x_face[0] - x_query[0] - normal_step * projection_dir[0], + x_face[1] - x_query[1] - normal_step * projection_dir[1], + x_face[2] - x_query[2] - normal_step * projection_dir[2] ); if ( line_residual > 0. && projection_residual > residual_tol * line_residual ) { return false; } @@ -251,8 +248,9 @@ TRIBOL_HOST_DEVICE inline bool EvalLinearFaceAtProjectedPoint( const RealT* face } TRIBOL_HOST_DEVICE inline bool EvalLinearEdgeAtProjectedPoint( const RealT* edge_coords, const RealT x_query[2], - const RealT projection_dir[2], RealT x_edge[3], RealT* phi, - const int value_dim = 0, const RealT* nodal_vals = nullptr, + const RealT projection_dir[2], RealT x_edge[3], + RealT* phi, const int value_dim = 0, + const RealT* nodal_vals = nullptr, RealT* values = nullptr ) { const RealT ax = edge_coords[0]; @@ -601,7 +599,8 @@ int ApplyNormal( CouplingScheme* cs ) const bool use_rate = pen_enfrc_options.constraint_type == KINEMATIC_AND_RATE; const RealT rate_penalty = - use_rate ? ComputeRatePenalty( mesh1, mesh2, penalty_stiff_per_area, pen_enfrc_options.rate_calculation ) : 0.; + use_rate ? ComputeRatePenalty( mesh1, mesh2, penalty_stiff_per_area, pen_enfrc_options.rate_calculation ) + : 0.; StackArrayT actual_vf1; StackArrayT actual_vf2; @@ -659,21 +658,18 @@ int ApplyNormal( CouplingScheme* cs ) RealT vel_q1[max_dim] = { 0., 0., 0. }; RealT vel_q2[max_dim] = { 0., 0., 0. }; - const bool mapped_face1 = - EvalLinearFaceAtProjectedPoint( actual_xf1, num_nodes_per_face, x_q, overlapNormal, x_qf1, phi_q1, - use_rate ? dim : 0, use_rate ? &actual_vf1[0] : nullptr, - use_rate ? vel_q1 : nullptr ); - const bool mapped_face2 = - EvalLinearFaceAtProjectedPoint( actual_xf2, num_nodes_per_face, x_q, overlapNormal, x_qf2, phi_q2, - use_rate ? dim : 0, use_rate ? &actual_vf2[0] : nullptr, - use_rate ? vel_q2 : nullptr ); + const bool mapped_face1 = EvalLinearFaceAtProjectedPoint( + actual_xf1, num_nodes_per_face, x_q, overlapNormal, x_qf1, phi_q1, use_rate ? dim : 0, + use_rate ? &actual_vf1[0] : nullptr, use_rate ? vel_q1 : nullptr ); + const bool mapped_face2 = EvalLinearFaceAtProjectedPoint( + actual_xf2, num_nodes_per_face, x_q, overlapNormal, x_qf2, phi_q2, use_rate ? dim : 0, + use_rate ? &actual_vf2[0] : nullptr, use_rate ? vel_q2 : nullptr ); if ( !mapped_face1 || !mapped_face2 ) { continue; } - RealT local_gap = - ( x_qf1[0] - x_qf2[0] ) * overlapNormal[0] + ( x_qf1[1] - x_qf2[1] ) * overlapNormal[1] + - ( x_qf1[2] - x_qf2[2] ) * overlapNormal[2]; + RealT local_gap = ( x_qf1[0] - x_qf2[0] ) * overlapNormal[0] + ( x_qf1[1] - x_qf2[1] ) * overlapNormal[1] + + ( x_qf1[2] - x_qf2[2] ) * overlapNormal[2]; if ( local_gap > gap_tol ) { continue; } @@ -682,9 +678,9 @@ int ApplyNormal( CouplingScheme* cs ) RealT local_pressure = local_gap * penalty_stiff_per_area; if ( use_rate && rate_penalty > 0. ) { - RealT local_vel_gap = - ( vel_q1[0] - vel_q2[0] ) * overlapNormal[0] + ( vel_q1[1] - vel_q2[1] ) * overlapNormal[1] + - ( vel_q1[2] - vel_q2[2] ) * overlapNormal[2]; + RealT local_vel_gap = ( vel_q1[0] - vel_q2[0] ) * overlapNormal[0] + + ( vel_q1[1] - vel_q2[1] ) * overlapNormal[1] + + ( vel_q1[2] - vel_q2[2] ) * overlapNormal[2]; if ( local_vel_gap <= 0. ) { local_pressure += local_vel_gap * rate_penalty; } @@ -732,8 +728,7 @@ int ApplyNormal( CouplingScheme* cs ) continue; } - RealT local_gap = - ( x_qf1[0] - x_qf2[0] ) * overlapNormal[0] + ( x_qf1[1] - x_qf2[1] ) * overlapNormal[1]; + RealT local_gap = ( x_qf1[0] - x_qf2[0] ) * overlapNormal[0] + ( x_qf1[1] - x_qf2[1] ) * overlapNormal[1]; if ( local_gap > gap_tol ) { continue; } @@ -794,7 +789,8 @@ int ApplyNormal( CouplingScheme* cs ) force_z = overlapNormal[2] * contact_force; } - AccumulateContactForce( mesh1, mesh2, index1, index2, dim, num_nodes_per_face, force_x, force_y, force_z, phi1, phi2 ); + AccumulateContactForce( mesh1, mesh2, index1, index2, dim, num_nodes_per_face, force_x, force_y, force_z, phi1, + phi2 ); // comment out debug logs; too much output during tests. Keep for easy // debugging if needed @@ -949,8 +945,8 @@ int ApplyTangential( CouplingScheme* constexpr int max_qpts = max_symmetric_triangle_qpts; RealT rule_wts[max_qpts] = { 0. }; RealT rule_coords[2 * max_qpts] = { 0. }; - const RealT visc = 0.5 * - ( mesh1.getElementData().m_viscous_damping_coeff + mesh2.getElementData().m_viscous_damping_coeff ); + const RealT visc = + 0.5 * ( mesh1.getElementData().m_viscous_damping_coeff + mesh2.getElementData().m_viscous_damping_coeff ); if ( dim == 3 ) { const int num_qpts = @@ -996,19 +992,16 @@ int ApplyTangential( CouplingScheme* RealT vel_q1[max_dim]; RealT vel_q2[max_dim]; - const bool mapped_face1 = - EvalLinearFaceAtProjectedPoint( actual_xf1, numNodesPerFace1, x_q, overlapNormal, x_qf1, phi_q1, dim, - &actual_vf1[0], vel_q1 ); - const bool mapped_face2 = - EvalLinearFaceAtProjectedPoint( actual_xf2, numNodesPerFace2, x_q, overlapNormal, x_qf2, phi_q2, dim, - &actual_vf2[0], vel_q2 ); + const bool mapped_face1 = EvalLinearFaceAtProjectedPoint( actual_xf1, numNodesPerFace1, x_q, overlapNormal, + x_qf1, phi_q1, dim, &actual_vf1[0], vel_q1 ); + const bool mapped_face2 = EvalLinearFaceAtProjectedPoint( actual_xf2, numNodesPerFace2, x_q, overlapNormal, + x_qf2, phi_q2, dim, &actual_vf2[0], vel_q2 ); if ( !mapped_face1 || !mapped_face2 ) { continue; } - RealT local_gap = - ( x_qf1[0] - x_qf2[0] ) * overlapNormal[0] + ( x_qf1[1] - x_qf2[1] ) * overlapNormal[1] + - ( x_qf1[2] - x_qf2[2] ) * overlapNormal[2]; + RealT local_gap = ( x_qf1[0] - x_qf2[0] ) * overlapNormal[0] + ( x_qf1[1] - x_qf2[1] ) * overlapNormal[1] + + ( x_qf1[2] - x_qf2[2] ) * overlapNormal[2]; RealT gap_tol = cs_view.getGapTol( index1, index2 ); if ( local_gap > gap_tol ) { continue; @@ -1058,16 +1051,15 @@ int ApplyTangential( CouplingScheme* RealT vel_q1[max_dim] = { 0., 0., 0. }; RealT vel_q2[max_dim] = { 0., 0., 0. }; - const bool mapped_face1 = - EvalLinearEdgeAtProjectedPoint( actual_xf1, x_q, overlapNormal, x_qf1, phi_q1, dim, &actual_vf1[0], vel_q1 ); - const bool mapped_face2 = - EvalLinearEdgeAtProjectedPoint( actual_xf2, x_q, overlapNormal, x_qf2, phi_q2, dim, &actual_vf2[0], vel_q2 ); + const bool mapped_face1 = EvalLinearEdgeAtProjectedPoint( actual_xf1, x_q, overlapNormal, x_qf1, phi_q1, dim, + &actual_vf1[0], vel_q1 ); + const bool mapped_face2 = EvalLinearEdgeAtProjectedPoint( actual_xf2, x_q, overlapNormal, x_qf2, phi_q2, dim, + &actual_vf2[0], vel_q2 ); if ( !mapped_face1 || !mapped_face2 ) { continue; } - RealT local_gap = - ( x_qf1[0] - x_qf2[0] ) * overlapNormal[0] + ( x_qf1[1] - x_qf2[1] ) * overlapNormal[1]; + RealT local_gap = ( x_qf1[0] - x_qf2[0] ) * overlapNormal[0] + ( x_qf1[1] - x_qf2[1] ) * overlapNormal[1]; RealT gap_tol = cs_view.getGapTol( index1, index2 ); if ( local_gap > gap_tol ) { continue; @@ -1118,7 +1110,8 @@ int ApplyTangential( CouplingScheme* force_z = visc * velGapTan[2]; } - AccumulateContactForce( mesh1, mesh2, index1, index2, dim, numNodesPerFace1, force_x, force_y, force_z, phi1, phi2 ); + AccumulateContactForce( mesh1, mesh2, index1, index2, dim, numNodesPerFace1, force_x, force_y, force_z, phi1, + phi2 ); } ); return 0; From 5cc845607dfa9c9651dfa586ce3a0ecf5fafc0bd Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Wed, 24 Jun 2026 22:48:13 -0700 Subject: [PATCH 6/6] address review comments --- src/tests/tribol_common_plane_penalty.cpp | 57 ++++++++++++++++------- src/tests/tribol_enforcement_options.cpp | 12 +++-- src/tests/tribol_iso_integ.cpp | 8 ++++ src/tests/tribol_mfem_common_plane.cpp | 3 +- src/tribol/common/Parameters.hpp | 3 +- src/tribol/integ/FE.hpp | 13 +++++- 6 files changed, 70 insertions(+), 26 deletions(-) diff --git a/src/tests/tribol_common_plane_penalty.cpp b/src/tests/tribol_common_plane_penalty.cpp index d6038016..bc11a403 100644 --- a/src/tests/tribol_common_plane_penalty.cpp +++ b/src/tests/tribol_common_plane_penalty.cpp @@ -198,6 +198,8 @@ WarpedQuadForceResult runWarpedQuadForceCase( tribol::PolyInteg rule ) { constexpr int numVerts = 4; + // Mesh 1 is a planar unit quad at z=0. Mesh 2 fully overlaps it in x-y, but is warped + // so that only one corner penetrates the plane. RealT x1[numVerts] = { 0.0, 1.0, 1.0, 0.0 }; RealT y1[numVerts] = { 0.0, 0.0, 1.0, 1.0 }; RealT z1[numVerts] = { 0.0, 0.0, 0.0, 0.0 }; @@ -231,9 +233,11 @@ WarpedQuadForceResult runWarpedQuadForceCase( tribol::PolyInteg rule ) tribol::FRICTIONLESS, tribol::PENALTY, tribol::BINNING_GRID, tribol::ExecutionMode::Sequential ); - tribol::setPenaltyOptions( 0, tribol::KINEMATIC, tribol::KINEMATIC_CONSTANT ); - tribol::setCommonPlaneIntegrationOptions( 0, rule, 3 ); - tribol::setContactAreaFrac( 0, 1.e-12 ); + constexpr tribol::IndexT couplingSchemeId = 0; + constexpr int quadratureOrder = 3; + tribol::setPenaltyOptions( couplingSchemeId, tribol::KINEMATIC, tribol::KINEMATIC_CONSTANT ); + tribol::setCommonPlaneIntegrationOptions( couplingSchemeId, rule, quadratureOrder ); + tribol::setContactAreaFrac( couplingSchemeId, 1.e-12 ); WarpedQuadForceResult result; RealT dt = 1.; @@ -253,7 +257,7 @@ WarpedQuadForceResult runWarpedQuadForceCase( tribol::PolyInteg rule ) return result; } -struct EdgeLocalContactResult { +struct EdgeLocalContactForceResult { int err{ -1 }; tribol::IndexT num_active_pairs{ 0 }; RealT gap{ 0. }; @@ -261,10 +265,12 @@ struct EdgeLocalContactResult { RealT mesh1_node_force_norm[2] = { 0., 0. }; }; -EdgeLocalContactResult runEdgeLocalContactCase( tribol::PolyInteg rule ) +EdgeLocalContactForceResult runEdgeLocalContactCase( tribol::PolyInteg rule ) { constexpr int numVerts = 2; + // Mesh 1 is a horizontal unit edge. Mesh 2 fully overlaps it in x, but is tilted so the + // penetration varies from 0.2 at node 0 to 0.05 at node 1. RealT x1[numVerts] = { 1.0, 0.0 }; RealT y1[numVerts] = { 0.0, 0.0 }; @@ -294,11 +300,13 @@ EdgeLocalContactResult runEdgeLocalContactCase( tribol::PolyInteg rule ) tribol::FRICTIONLESS, tribol::PENALTY, tribol::BINNING_GRID, tribol::ExecutionMode::Sequential ); - tribol::setPenaltyOptions( 0, tribol::KINEMATIC, tribol::KINEMATIC_CONSTANT ); - tribol::setCommonPlaneIntegrationOptions( 0, rule, 4 ); - tribol::setContactAreaFrac( 0, 1.e-12 ); + constexpr tribol::IndexT couplingSchemeId = 0; + constexpr int quadratureOrder = 4; + tribol::setPenaltyOptions( couplingSchemeId, tribol::KINEMATIC, tribol::KINEMATIC_CONSTANT ); + tribol::setCommonPlaneIntegrationOptions( couplingSchemeId, rule, quadratureOrder ); + tribol::setContactAreaFrac( couplingSchemeId, 1.e-12 ); - EdgeLocalContactResult result; + EdgeLocalContactForceResult result; RealT dt = 1.; result.err = tribol::update( 1, 1., dt ); @@ -376,6 +384,8 @@ TEST_F( CommonPlaneTest, penetration_gap_check ) TEST_F( CommonPlaneTest, multipoint_quad_execution ) { + // Planar quad meshes with full face overlap and uniform 0.1 interpenetration. Multi-point + // integration should reproduce the same gap and valid force sense as the single-point rule. this->m_mesh.mortarMeshId = 0; this->m_mesh.nonmortarMeshId = 1; @@ -417,6 +427,8 @@ TEST_F( CommonPlaneTest, multipoint_quad_execution ) TEST_F( CommonPlaneTest, multipoint_triangle_execution ) { + // Planar tet-surface triangle meshes with full overlap and uniform 0.1 interpenetration. + // This exercises triangle overlap integration in the CommonPlane penalty path. this->m_mesh.mortarMeshId = 0; this->m_mesh.nonmortarMeshId = 1; @@ -458,22 +470,27 @@ TEST_F( CommonPlaneTest, multipoint_triangle_execution ) TEST_F( CommonPlaneTest, multipoint_warped_quad_local_contact ) { + // Single-point integration at the overlap centroid misses the local warped-quad penetration. + // Full triangle-decomposition quadrature samples the penetrating region and generates force. const auto single_point = runWarpedQuadForceCase( tribol::SINGLE_POINT ); const auto multi_point = runWarpedQuadForceCase( tribol::MULTI_POINT ); EXPECT_EQ( single_point.err, 0 ); EXPECT_EQ( multi_point.err, 0 ); - EXPECT_GT( single_point.gap, 0. ); - EXPECT_GT( multi_point.gap, 0. ); + constexpr RealT expected_gap = 0.5852486869304208; + EXPECT_NEAR( single_point.gap, expected_gap, 1.e-12 ); + EXPECT_NEAR( multi_point.gap, expected_gap, 1.e-12 ); EXPECT_NEAR( single_point.total_abs_force, 0., 1.e-12 ); - EXPECT_GT( multi_point.total_abs_force, 1.e-6 ); + EXPECT_NEAR( multi_point.total_abs_force, 8.20970073925438e-4, 1.e-12 ); EXPECT_NEAR( multi_point.total_force_z, 0., 1.e-12 ); } TEST_F( CommonPlaneTest, multipoint_edge_local_contact ) { + // Both rules detect the tilted full-overlap edge contact, but multi-point integration resolves + // the nonuniform penetration and therefore produces a larger nodal force imbalance. const auto single_point = runEdgeLocalContactCase( tribol::SINGLE_POINT ); const auto multi_point = runEdgeLocalContactCase( tribol::MULTI_POINT ); @@ -483,19 +500,21 @@ TEST_F( CommonPlaneTest, multipoint_edge_local_contact ) EXPECT_EQ( single_point.num_active_pairs, 1 ); EXPECT_EQ( multi_point.num_active_pairs, 1 ); - EXPECT_LT( single_point.gap, 0. ); - EXPECT_LT( multi_point.gap, 0. ); + constexpr RealT expected_gap = -0.12423774246760939; + EXPECT_NEAR( single_point.gap, expected_gap, 1.e-12 ); + EXPECT_NEAR( multi_point.gap, expected_gap, 1.e-12 ); - EXPECT_GT( single_point.total_abs_force, 1.e-6 ); - EXPECT_GT( multi_point.total_abs_force, 1.e-6 ); + EXPECT_NEAR( single_point.total_abs_force, 0.13126963259866961, 1.e-12 ); + EXPECT_NEAR( multi_point.total_abs_force, 0.13227012255210607, 1.e-12 ); + // Imbalance is the difference between the force norms at the two nodes of mesh 1. const RealT single_point_imbalance = std::abs( single_point.mesh1_node_force_norm[0] - single_point.mesh1_node_force_norm[1] ); const RealT multi_point_imbalance = std::abs( multi_point.mesh1_node_force_norm[0] - multi_point.mesh1_node_force_norm[1] ); - EXPECT_LT( single_point_imbalance, 1.e-3 ); - EXPECT_GT( multi_point_imbalance, single_point_imbalance + 1.e-3 ); + EXPECT_NEAR( single_point_imbalance, 3.37547013399012e-4, 1.e-12 ); + EXPECT_NEAR( multi_point_imbalance, 1.2454070813893655e-2, 1.e-12 ); } TEST_F( CommonPlaneTest, separation_gap_check ) @@ -924,6 +943,8 @@ TEST_F( CommonPlaneTest, common_plane_viscous_tangential_2d ) // by the gap and then divided amongst the edge nodes RealT force_y = 0.5 * gap / numVerts; RealT force_x = visc_coeff * ( vx1[0] - vx2[0] ) / numVerts; + // The analytic force above assumes the nominal edge directions, while the implementation uses + // the CommonPlane normal/tangent for the slightly tilted contact pair. constexpr RealT tol = 5.e-5; for ( int i = 0; i < numVerts; ++i ) { EXPECT_NEAR( fx1[i], -force_x, tol ); diff --git a/src/tests/tribol_enforcement_options.cpp b/src/tests/tribol_enforcement_options.cpp index 11090545..5cd39440 100644 --- a/src/tests/tribol_enforcement_options.cpp +++ b/src/tests/tribol_enforcement_options.cpp @@ -164,20 +164,24 @@ TEST_F( EnforcementOptionsTest, penalty_kinematic_constant_error ) TEST_F( EnforcementOptionsTest, common_plane_integration_options_are_stored ) { + // Verify that user-provided CommonPlane polygon integration options are stored on the + // coupling scheme and are available during penalty enforcement. tribol::TestMesh* mesh = new tribol::TestMesh(); SetupTest( mesh ); + constexpr tribol::IndexT couplingSchemeId = 0; + constexpr int quadratureOrder = 4; RealT penalty = 1.0; tribol::setKinematicConstantPenalty( 0, penalty ); tribol::setKinematicConstantPenalty( 1, penalty ); - tribol::setPenaltyOptions( 0, tribol::KINEMATIC, tribol::KINEMATIC_CONSTANT ); - tribol::setCommonPlaneIntegrationOptions( 0, tribol::MULTI_POINT, 4 ); + tribol::setPenaltyOptions( couplingSchemeId, tribol::KINEMATIC, tribol::KINEMATIC_CONSTANT ); + tribol::setCommonPlaneIntegrationOptions( couplingSchemeId, tribol::MULTI_POINT, quadratureOrder ); tribol::CouplingSchemeManager& csManager = tribol::CouplingSchemeManager::getInstance(); - tribol::CouplingScheme* scheme = &csManager.at( 0 ); + tribol::CouplingScheme* scheme = &csManager.at( couplingSchemeId ); EXPECT_EQ( scheme->getEnforcementOptions().penalty_options.common_plane_rule, tribol::MULTI_POINT ); - EXPECT_EQ( scheme->getEnforcementOptions().penalty_options.common_plane_quadrature_order, 4 ); + EXPECT_EQ( scheme->getEnforcementOptions().penalty_options.common_plane_quadrature_order, quadratureOrder ); tribol::finalize(); delete mesh; diff --git a/src/tests/tribol_iso_integ.cpp b/src/tests/tribol_iso_integ.cpp index c9b9bd70..8f4f7a6f 100644 --- a/src/tests/tribol_iso_integ.cpp +++ b/src/tests/tribol_iso_integ.cpp @@ -31,8 +31,11 @@ RealT factorial( int n ) RealT referenceTriangleMoment( int px, int py ) { return factorial( px ) * factorial( py ) / factorial( px + py + 2 ); } +// Moment integral over the unit-area reference triangle used by the symmetric rules. RealT normalizedReferenceTriangleMoment( int px, int py ) { return 2. * referenceTriangleMoment( px, py ); } +// Evaluate the x^px y^py moment for either the legacy CommonPlane triangle rule or the +// newer symmetric triangle rule. RealT evalTriangleRuleMoment( bool use_legacy, int order, int px, int py ) { RealT wts[tribol::max_symmetric_triangle_qpts] = { 0. }; @@ -275,6 +278,8 @@ TEST_F( IsoIntegTest, nonaffine ) TEST( TriangleRuleTest, legacy_and_symmetric_match_on_shared_orders ) { + // Orders 2 and 4 are supported by both rule implementations and should integrate the + // same low-order reference-triangle moments. for ( int order : { 2, 4 } ) { EXPECT_NEAR( evalTriangleRuleMoment( true, order, 0, 0 ), evalTriangleRuleMoment( false, order, 0, 0 ), 2.e-10 ); EXPECT_NEAR( evalTriangleRuleMoment( true, order, 2, 0 ), evalTriangleRuleMoment( false, order, 2, 0 ), 2.e-10 ); @@ -284,6 +289,8 @@ TEST( TriangleRuleTest, legacy_and_symmetric_match_on_shared_orders ) TEST( TriangleRuleTest, gauss_poly_int_tri_supports_order_10 ) { + // CommonPlane triangle-decomposition integration supports the order-10 symmetric rule + // on triangular overlap facets in 3D. constexpr int dim = 3; constexpr int num_nodes = 3; RealT xyz[dim * num_nodes] = { 0., 0., 0., 1., 0., 0., 0., 1., 0. }; @@ -303,6 +310,7 @@ TEST( TriangleRuleTest, gauss_poly_int_tri_supports_order_10 ) EXPECT_EQ( integ.numIPs, 75 ); EXPECT_NEAR( area, 0.5, 1.e-14 ); + // Integral of x^7 y^3 over the physical unit right triangle. EXPECT_NEAR( moment73, referenceTriangleMoment( 7, 3 ), 1.e-14 ); EXPECT_NEAR( evalTriangleRuleMoment( false, 10, 7, 3 ), normalizedReferenceTriangleMoment( 7, 3 ), 1.e-14 ); } diff --git a/src/tests/tribol_mfem_common_plane.cpp b/src/tests/tribol_mfem_common_plane.cpp index 370b2393..6390ad69 100644 --- a/src/tests/tribol_mfem_common_plane.cpp +++ b/src/tests/tribol_mfem_common_plane.cpp @@ -38,7 +38,8 @@ * * Both the element penalty and a constant penalty are tested, with the constant penalty tuned to match the element * penalty for this case. The contact solve is also exercised with both the legacy single-point CommonPlane integration - * rule and the new full triangle-decomposition rule. + * rule and the new order-3 full triangle-decomposition rule. The contact surfaces are planar, so both rules should + * produce the same force and gap response. * */ class MfemCommonPlaneTest diff --git a/src/tribol/common/Parameters.hpp b/src/tribol/common/Parameters.hpp index c29d78fe..4bfcb09b 100644 --- a/src/tribol/common/Parameters.hpp +++ b/src/tribol/common/Parameters.hpp @@ -236,7 +236,7 @@ enum IntNodalFields enum PolyInteg { SINGLE_POINT, ///! Single point integration at centroid of polygon - MULTI_POINT, ///! Multi-point integration over the overlap + MULTI_POINT, ///! Multi-point integration over a triangle decomposition of the overlap NUM_INTEG_RULES }; @@ -445,6 +445,7 @@ struct PenaltyEnforcementOptions { KinematicPenaltyCalculation kinematic_calculation; RatePenaltyCalculation rate_calculation; PolyInteg common_plane_rule{ SINGLE_POINT }; + ///! Triangle/segment quadrature order used when common_plane_rule is MULTI_POINT; ignored for SINGLE_POINT int common_plane_quadrature_order{ 3 }; bool constraint_type_set{ false }; diff --git a/src/tribol/integ/FE.hpp b/src/tribol/integ/FE.hpp index 494eed29..1aeb594f 100644 --- a/src/tribol/integ/FE.hpp +++ b/src/tribol/integ/FE.hpp @@ -106,9 +106,14 @@ TRIBOL_HOST_DEVICE inline void SegmentBasis( const RealT* const x, const RealT p * \param [in,out] xi (xi,eta) coordinates in parent space * * \pre xA, yA, and zA are pointer to arrays of length, numNodes + * \pre x is a point on the physical edge or face being inverse-mapped * * \note This routine works in 2D or 3D. In 2D, zA is a nullptr and * x[2] is equal to 0. + * \note For a quadrilateral, an off-face point is handled by the same least-squares + * Newton solve, effectively mapping the closest point in the face-normal direction + * for near-planar contact faces. Callers should not rely on off-face behavior for + * points that are far from the surface. * */ TRIBOL_HOST_DEVICE inline void InvIso( const RealT x[3], const RealT* xA, const RealT* yA, const RealT* zA, @@ -457,7 +462,7 @@ TRIBOL_HOST_DEVICE inline int GetNumFaceNodes( int dim, FaceOrderType order_type * \param [in] x pointer to stacked physical coordinates of the face vertices * \param [in] pX x-coordinate of the evaluation point * \param [in] pY y-coordinate of the evaluation point - * \param [in] pZ z-coordinate of the evaluation point + * \param [in] pZ z-coordinate of the evaluation point; ignored for physical edges * \param [in] numNodes number of nodes on the face * \param [in] vertexId local node id whose basis function is to be evaluated * \param [in,out] phi evaluated basis function value @@ -465,6 +470,8 @@ TRIBOL_HOST_DEVICE inline int GetNumFaceNodes( int dim, FaceOrderType order_type * \note For triangles and quadrilaterals this routine inverse-maps the physical * point to parent space and then evaluates the corresponding linear * isoparametric basis function. + * \note The evaluation point is assumed to lie on the physical edge or face. Off-face + * points inherit the implicit projection behavior of InvIso. * */ TRIBOL_HOST_DEVICE inline void EvalBasisOnPhysicalFace( const RealT* const x, const RealT pX, const RealT pY, @@ -511,7 +518,7 @@ TRIBOL_HOST_DEVICE inline void EvalBasisOnPhysicalFace( const RealT* const x, co * \param [in] x pointer to stacked physical coordinates of the face vertices * \param [in] pX x-coordinate of the evaluation point * \param [in] pY y-coordinate of the evaluation point - * \param [in] pZ z-coordinate of the evaluation point + * \param [in] pZ z-coordinate of the evaluation point; ignored for physical edges * \param [in] numNodes number of nodes on the face * \param [in] galerkinDim vector dimension of the nodal coefficients * \param [in] nodeVals stacked nodal coefficients for the Galerkin approximation @@ -519,6 +526,8 @@ TRIBOL_HOST_DEVICE inline void EvalBasisOnPhysicalFace( const RealT* const x, co * * \note This helper is topology-aware and supports the linear segment, * triangle, and quadrilateral basis evaluations used by CommonPlane. + * \note The evaluation point is assumed to lie on the physical edge or face. Off-face + * points inherit the implicit projection behavior of InvIso. * */ TRIBOL_HOST_DEVICE inline void GalerkinEvalOnPhysicalFace( const RealT* const x, const RealT pX, const RealT pY,