diff --git a/src/tests/tribol_common_plane_penalty.cpp b/src/tests/tribol_common_plane_penalty.cpp index 5f264725..bc11a403 100644 --- a/src/tests/tribol_common_plane_penalty.cpp +++ b/src/tests/tribol_common_plane_penalty.cpp @@ -185,10 +185,147 @@ class CommonPlaneTest : public ::testing::Test { void SetUp() override {} void TearDown() override { this->m_mesh.clear(); } +}; - 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; + + // 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 }; + + 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 ); + + 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.; + 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; +} + +struct EdgeLocalContactForceResult { + 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. }; }; +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 }; + + 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 ); + + 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 ); + + EdgeLocalContactForceResult 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; @@ -245,6 +382,141 @@ TEST_F( CommonPlaneTest, penetration_gap_check ) tribol::finalize(); } +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; + + 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::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 ); + + 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, 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; + + 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::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 ); + + 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, 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 ); + + 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_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 ); + + 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 ); + + 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_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_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 ) { this->m_mesh.mortarMeshId = 0; @@ -671,11 +943,14 @@ 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, 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..5cd39440 100644 --- a/src/tests/tribol_enforcement_options.cpp +++ b/src/tests/tribol_enforcement_options.cpp @@ -162,6 +162,31 @@ TEST_F( EnforcementOptionsTest, penalty_kinematic_constant_error ) delete mesh; } +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( couplingSchemeId, tribol::KINEMATIC, tribol::KINEMATIC_CONSTANT ); + tribol::setCommonPlaneIntegrationOptions( couplingSchemeId, tribol::MULTI_POINT, quadratureOrder ); + + tribol::CouplingSchemeManager& csManager = tribol::CouplingSchemeManager::getInstance(); + 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, quadratureOrder ); + + tribol::finalize(); + delete mesh; +} + TEST_F( EnforcementOptionsTest, penalty_kinematic_element_error ) { // Setup boiler plate test data etc. diff --git a/src/tests/tribol_iso_integ.cpp b/src/tests/tribol_iso_integ.cpp index d34fdb23..8f4f7a6f 100644 --- a/src/tests/tribol_iso_integ.cpp +++ b/src/tests/tribol_iso_integ.cpp @@ -18,6 +18,40 @@ 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 ); } + +// 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. }; + 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 +276,45 @@ TEST_F( IsoIntegTest, nonaffine ) EXPECT_EQ( convrg, true ); } +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 ); + EXPECT_NEAR( evalTriangleRuleMoment( true, order, 1, 1 ), evalTriangleRuleMoment( false, order, 1, 1 ), 2.e-10 ); + } +} + +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. }; + + 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 ); + // 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 ); +} + int main( int argc, char* argv[] ) { int result = 0; diff --git a/src/tests/tribol_mfem_common_plane.cpp b/src/tests/tribol_mfem_common_plane.cpp index d6855d71..6390ad69 100644 --- a/src/tests/tribol_mfem_common_plane.cpp +++ b/src/tests/tribol_mfem_common_plane.cpp @@ -37,10 +37,13 @@ * 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 order-3 full triangle-decomposition rule. The contact surfaces are planar, so both rules should + * produce the same force and gap response. * */ -class MfemCommonPlaneTest : public testing::TestWithParam> { +class MfemCommonPlaneTest + : public testing::TestWithParam> { protected: tribol::RealT max_disp_; void SetUp() override @@ -166,6 +169,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 +217,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::MULTI_POINT ), + std::make_tuple( 1, tribol::KINEMATIC_ELEMENT, tribol::SINGLE_POINT ), + 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::MULTI_POINT ), + std::make_tuple( 2, tribol::KINEMATIC_ELEMENT, tribol::SINGLE_POINT ), + 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 6bb52366..4bfcb09b 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 a triangle decomposition of the overlap NUM_INTEG_RULES }; @@ -444,6 +444,9 @@ struct PenaltyEnforcementOptions { PenaltyConstraintType constraint_type; 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 }; bool kinematic_calc_set{ false }; diff --git a/src/tribol/integ/FE.hpp b/src/tribol/integ/FE.hpp index 56e9dc8d..050e1f46 100644 --- a/src/tribol/integ/FE.hpp +++ b/src/tribol/integ/FE.hpp @@ -106,13 +106,18 @@ 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. * */ -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 +314,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 +343,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 +364,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 +424,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 +453,103 @@ 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; 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 + * + * \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, + 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; 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 + * \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. + * \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, + 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.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 8f9a57eb..bb53f752 100644 --- a/src/tribol/integ/Integration.hpp +++ b/src/tribol/integ/Integration.hpp @@ -90,6 +90,18 @@ 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; +/// Maximum number of quadrature points in the built-in Gauss-Legendre segment rules. +constexpr int max_segment_gauss_legendre_qpts = 10; + /*! * * \brief Populates the integration points and weights on the IntegPts object @@ -120,13 +132,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 ); /*! * @@ -174,92 +188,571 @@ 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] ); + 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 ); } - } 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] ); + 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( quadrature_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 quadrature_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 MULTI_POINT: + EvalWeakFormIntegralCommonPlaneMultiPoint( elem, quadrature_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..f8640026 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 quadrature_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( 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_quadrature_order = quadrature_order; +} + //------------------------------------------------------------------------------ void setKinematicConstantPenalty( IndexT mesh_id, RealT k ) { diff --git a/src/tribol/interface/tribol.hpp b/src/tribol/interface/tribol.hpp index 36e6fa03..991c621e 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] 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 quadrature_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 5278e3b5..b555eaf9 100644 --- a/src/tribol/physics/CommonPlane.cpp +++ b/src/tribol/physics/CommonPlane.cpp @@ -17,48 +17,342 @@ 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 bool Solve3x3( const RealT A[3][3], const RealT b[3], RealT x[3] ) +{ + 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; + } + + 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. ); + if ( values != nullptr ) { + initRealArray( values, value_dim, 0. ); + } + + for ( int a = 0; a < num_nodes; ++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; + } + + 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; + } + + 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 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, + 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]; + + // accumulate contributions in host code's registered nodal force arrays + tribol::atomicAdd( &mesh1.getResponse()[0][node0], -nodal_force_x1 ); + tribol::atomicAdd( &mesh2.getResponse()[0][node1], nodal_force_x2 ); + + tribol::atomicAdd( &mesh1.getResponse()[1][node0], -nodal_force_y1 ); + tribol::atomicAdd( &mesh2.getResponse()[1][node1], nodal_force_y2 ); + + // there is no z component for 2D + if ( dim == 3 ) { + tribol::atomicAdd( &mesh1.getResponse()[2][node0], -nodal_force_z1 ); + tribol::atomicAdd( &mesh2.getResponse()[2][node1], nodal_force_z2 ); + } } // 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(); - StackArrayT x1; - StackArrayT v1; + const auto dim = plane.m_dim; + const RealT rate_penalty = ComputeRatePenalty( m1, m2, element_penalty, rate_calc ); + if ( rate_penalty == 0. ) { + return 0.; + } + + // 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 ); @@ -77,14 +371,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]; @@ -149,14 +443,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.}; @@ -166,6 +452,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_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. + 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 ) { @@ -237,9 +532,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]; @@ -289,6 +581,172 @@ int ApplyNormal( CouplingScheme* cs ) cntctElem.overlapNormal = overlapNormal; cntctElem.overlapArea = plane.m_area; + if ( use_multi_point ) { + 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 = max_symmetric_triangle_qpts; + RealT rule_wts[max_qpts] = { 0. }; + RealT rule_coords[2 * max_qpts] = { 0. }; + bool has_contact = false; + + 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] = { 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; + + 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 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. }; + 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 = + 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 = + 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]; + 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]; + if ( local_vel_gap <= 0. ) { + local_pressure += local_vel_gap * rate_penalty; + } + } + + 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; + + AccumulateContactForce( mesh1, mesh2, index1, index2, dim, num_nodes_per_face, force_x, force_y, 0., phi_q1, + phi_q2 ); + } + } + + plane.m_inContact = has_contact; + return; + } + // create arrays to hold nodal residual weak form integral evaluations RealT phi1[max_nodes_per_face]; RealT phi2[max_nodes_per_face]; @@ -299,7 +757,8 @@ 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_quadrature_order, phi1, phi2 ); /////////////////////////////////////////////////////////////////////// // Computation of full contact nodal force contributions // @@ -320,50 +779,8 @@ 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 - tribol::atomicAdd( &mesh1.getResponse()[0][node0], -nodal_force_x1 ); - tribol::atomicAdd( &mesh2.getResponse()[0][node1], nodal_force_x2 ); - - tribol::atomicAdd( &mesh1.getResponse()[1][node0], -nodal_force_y1 ); - tribol::atomicAdd( &mesh2.getResponse()[1][node1], nodal_force_y2 ); - - // there is no z component for 2D - if ( dim == 3 ) { - tribol::atomicAdd( &mesh1.getResponse()[2][node0], -nodal_force_z1 ); - tribol::atomicAdd( &mesh2.getResponse()[2][node1], nodal_force_z2 ); - } - } // 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 @@ -402,24 +819,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_multi_point = pen_enfrc_options.common_plane_rule == MULTI_POINT; - 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 ); @@ -438,14 +854,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]; @@ -506,9 +922,162 @@ int ApplyTangential( CouplingScheme* cntctElem.overlapNormal = overlapNormal; cntctElem.overlapArea = plane.m_area; + if ( use_multi_point ) { + 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 = 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 ); + + 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 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] = { 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 ); + 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 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]; + + 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 = 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, 0., phi_q1, + phi_q2 ); + } + } + + return; + } + // 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. ); @@ -516,7 +1085,8 @@ 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_quadrature_order, phi1, phi2 ); ///////////////////////////////////////////////////// // Computation of tangential viscous damping force // @@ -530,34 +1100,8 @@ 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 - tribol::atomicAdd( &mesh1.getResponse()[0][node0], -nodal_force_x1 ); - tribol::atomicAdd( &mesh2.getResponse()[0][node1], nodal_force_x2 ); - - tribol::atomicAdd( &mesh1.getResponse()[1][node0], -nodal_force_y1 ); - tribol::atomicAdd( &mesh2.getResponse()[1][node1], nodal_force_y2 ); - - // there is no z component for 2D - if ( dim == 3 ) { - tribol::atomicAdd( &mesh1.getResponse()[2][node0], -nodal_force_z1 ); - tribol::atomicAdd( &mesh2.getResponse()[2][node1], nodal_force_z2 ); - } - } // 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/physics/Mortar.cpp b/src/tribol/physics/Mortar.cpp index 03303af1..3179bb33 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. diff --git a/src/tribol/utils/TestUtils.cpp b/src/tribol/utils/TestUtils.cpp index f09bd38c..4e9e2ada 100644 --- a/src/tribol/utils/TestUtils.cpp +++ b/src/tribol/utils/TestUtils.cpp @@ -1116,6 +1116,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_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 52272275..5b1c269b 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_quadrature_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_quadrature_order; RealT rate_penalty; RealT rate_penalty_ratio; RealT const_penalty;