diff --git a/src/tests/tribol_common_plane_penalty.cpp b/src/tests/tribol_common_plane_penalty.cpp index 5f264725..2e7e1c25 100644 --- a/src/tests/tribol_common_plane_penalty.cpp +++ b/src/tests/tribol_common_plane_penalty.cpp @@ -185,10 +185,167 @@ class CommonPlaneTest : public ::testing::Test { void SetUp() override {} void TearDown() override { this->m_mesh.clear(); } +}; - protected: +/** Results returned by the warped-quadrilateral force regression helper. */ +struct WarpedQuadForceResult { + int update_error{ -1 }; ///< Return code from tribol::update(). + RealT gap{ 0. }; ///< Geometric gap stored on the CommonPlane pair. + RealT total_absolute_force{ 0. }; ///< Sum of absolute nodal-force components on both surfaces. + RealT total_z_force{ 0. }; ///< Net z-directed force on both surfaces. +}; + +/** + * @brief Run one warped-quadrilateral force case with the requested overlap rule. + * + * @param integration_rule CommonPlane overlap integration rule + * @param quadrature_order Polynomial order requested for multipoint integration + * @return Update status and force diagnostics for the case + */ +WarpedQuadForceResult runWarpedQuadForceCase( tribol::PolyInteg integration_rule, int quadrature_order ) +{ + constexpr int number_of_vertices = 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 first_x_coordinates[number_of_vertices] = { 0.0, 1.0, 1.0, 0.0 }; + RealT first_y_coordinates[number_of_vertices] = { 0.0, 0.0, 1.0, 1.0 }; + RealT first_z_coordinates[number_of_vertices] = { 0.0, 0.0, 0.0, 0.0 }; + + RealT second_x_coordinates[number_of_vertices] = { 0.0, 0.0, 1.0, 1.0 }; + RealT second_y_coordinates[number_of_vertices] = { 0.0, 1.0, 1.0, 0.0 }; + RealT second_z_coordinates[number_of_vertices] = { -0.20, 1.00, 1.00, 1.00 }; + + tribol::IndexT first_connectivity[number_of_vertices] = { 0, 1, 2, 3 }; + tribol::IndexT second_connectivity[number_of_vertices] = { 0, 1, 2, 3 }; + + tribol::registerMesh( 0, 1, number_of_vertices, first_connectivity, static_cast( tribol::LINEAR_QUAD ), + first_x_coordinates, first_y_coordinates, first_z_coordinates, tribol::MemorySpace::Host ); + tribol::registerMesh( 1, 1, number_of_vertices, second_connectivity, static_cast( tribol::LINEAR_QUAD ), + second_x_coordinates, second_y_coordinates, second_z_coordinates, tribol::MemorySpace::Host ); + + RealT first_response_x[number_of_vertices] = { 0., 0., 0., 0. }; + RealT first_response_y[number_of_vertices] = { 0., 0., 0., 0. }; + RealT first_response_z[number_of_vertices] = { 0., 0., 0., 0. }; + RealT second_response_x[number_of_vertices] = { 0., 0., 0., 0. }; + RealT second_response_y[number_of_vertices] = { 0., 0., 0., 0. }; + RealT second_response_z[number_of_vertices] = { 0., 0., 0., 0. }; + + tribol::registerNodalResponse( 0, first_response_x, first_response_y, first_response_z ); + tribol::registerNodalResponse( 1, second_response_x, second_response_y, second_response_z ); + + 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 coupling_scheme_id = 0; + tribol::setPenaltyOptions( coupling_scheme_id, tribol::KINEMATIC, tribol::KINEMATIC_CONSTANT ); + tribol::setCommonPlaneIntegrationOptions( coupling_scheme_id, integration_rule, quadrature_order ); + tribol::setContactAreaFrac( coupling_scheme_id, 1.e-12 ); + + WarpedQuadForceResult result; + RealT timestep = 1.; + result.update_error = tribol::update( 1, 1., timestep ); + + auto& coupling_scheme = tribol::CouplingSchemeManager::getInstance().at( coupling_scheme_id ); + EXPECT_EQ( 1, coupling_scheme.getNumActivePairs() ); + result.gap = coupling_scheme.getCompGeom().getCommonPlane( 0 ).m_gap; + + for ( int node_index = 0; node_index < number_of_vertices; ++node_index ) { + result.total_absolute_force += std::abs( first_response_x[node_index] ) + std::abs( first_response_y[node_index] ) + + std::abs( first_response_z[node_index] ); + result.total_absolute_force += std::abs( second_response_x[node_index] ) + + std::abs( second_response_y[node_index] ) + + std::abs( second_response_z[node_index] ); + result.total_z_force += first_response_z[node_index] + second_response_z[node_index]; + } + + tribol::finalize(); + return result; +} + +/** Results returned by the tilted-edge force regression helper. */ +struct EdgeLocalContactForceResult { + int update_error{ -1 }; ///< Return code from tribol::update(). + tribol::IndexT number_of_active_pairs{ 0 }; ///< Active CommonPlane pair count. + RealT gap{ 0. }; ///< Geometric gap stored on the CommonPlane pair. + RealT total_absolute_force{ 0. }; ///< Sum of absolute nodal-force components. + RealT first_mesh_node_force_magnitudes[2]{}; ///< Force magnitude at each first-surface node. }; +/** + * @brief Run one tilted-edge force case with the requested overlap rule. + * + * @param integration_rule CommonPlane overlap integration rule + * @param quadrature_order Polynomial order requested for multipoint integration + * @return Update status and force diagnostics for the case + */ +EdgeLocalContactForceResult runEdgeLocalContactCase( tribol::PolyInteg integration_rule, int quadrature_order ) +{ + constexpr int number_of_vertices = 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 first_x_coordinates[number_of_vertices] = { 1.0, 0.0 }; + RealT first_y_coordinates[number_of_vertices] = { 0.0, 0.0 }; + + RealT second_x_coordinates[number_of_vertices] = { 0.0, 1.0 }; + RealT second_y_coordinates[number_of_vertices] = { -0.2, -0.05 }; + + tribol::IndexT first_connectivity[number_of_vertices] = { 0, 1 }; + tribol::IndexT second_connectivity[number_of_vertices] = { 0, 1 }; + + tribol::registerMesh( 0, 1, number_of_vertices, first_connectivity, static_cast( tribol::LINEAR_EDGE ), + first_x_coordinates, first_y_coordinates, nullptr, tribol::MemorySpace::Host ); + tribol::registerMesh( 1, 1, number_of_vertices, second_connectivity, static_cast( tribol::LINEAR_EDGE ), + second_x_coordinates, second_y_coordinates, nullptr, tribol::MemorySpace::Host ); + + RealT first_response_x[number_of_vertices] = { 0., 0. }; + RealT first_response_y[number_of_vertices] = { 0., 0. }; + RealT second_response_x[number_of_vertices] = { 0., 0. }; + RealT second_response_y[number_of_vertices] = { 0., 0. }; + + tribol::registerNodalResponse( 0, first_response_x, first_response_y, nullptr ); + tribol::registerNodalResponse( 1, second_response_x, second_response_y, 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 coupling_scheme_id = 0; + tribol::setPenaltyOptions( coupling_scheme_id, tribol::KINEMATIC, tribol::KINEMATIC_CONSTANT ); + tribol::setCommonPlaneIntegrationOptions( coupling_scheme_id, integration_rule, quadrature_order ); + tribol::setContactAreaFrac( coupling_scheme_id, 1.e-12 ); + + EdgeLocalContactForceResult result; + RealT timestep = 1.; + result.update_error = tribol::update( 1, 1., timestep ); + + tribol::CouplingScheme* coupling_scheme = &tribol::CouplingSchemeManager::getInstance().at( coupling_scheme_id ); + result.number_of_active_pairs = coupling_scheme->getNumActivePairs(); + if ( result.number_of_active_pairs > 0 ) { + result.gap = coupling_scheme->getCompGeom().getCommonPlane( 0 ).m_gap; + } + + for ( int node_index = 0; node_index < number_of_vertices; ++node_index ) { + result.total_absolute_force += std::abs( first_response_x[node_index] ) + std::abs( first_response_y[node_index] ) + + std::abs( second_response_x[node_index] ) + + std::abs( second_response_y[node_index] ); + result.first_mesh_node_force_magnitudes[node_index] = + tribol::magnitude( first_response_x[node_index], first_response_y[node_index] ); + } + + tribol::finalize(); + + return result; +} + TEST_F( CommonPlaneTest, penetration_gap_check ) { this->m_mesh.mortarMeshId = 0; @@ -245,6 +402,179 @@ TEST_F( CommonPlaneTest, penetration_gap_check ) tribol::finalize(); } +/** Verify multipoint CommonPlane force integration on planar quadrilateral faces. */ +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(); +} + +/** Verify multipoint CommonPlane force integration on planar triangular faces. */ +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(); +} + +/** Verify that multipoint quadrature resolves localized contact on warped quadrilateral faces. */ +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. + constexpr int lower_quadrature_order = 3; + constexpr int higher_quadrature_order = 6; + constexpr int reference_quadrature_order = 10; + const auto single_point = runWarpedQuadForceCase( tribol::SINGLE_POINT, lower_quadrature_order ); + const auto lower_order_result = runWarpedQuadForceCase( tribol::MULTI_POINT, lower_quadrature_order ); + const auto higher_order_result = runWarpedQuadForceCase( tribol::MULTI_POINT, higher_quadrature_order ); + const auto overintegrated_reference = runWarpedQuadForceCase( tribol::MULTI_POINT, reference_quadrature_order ); + + EXPECT_EQ( single_point.update_error, 0 ); + EXPECT_EQ( lower_order_result.update_error, 0 ); + EXPECT_EQ( higher_order_result.update_error, 0 ); + EXPECT_EQ( overintegrated_reference.update_error, 0 ); + + constexpr RealT expected_gap = 0.5852486869304208; + EXPECT_NEAR( single_point.gap, expected_gap, 1.e-12 ); + EXPECT_NEAR( lower_order_result.gap, expected_gap, 1.e-12 ); + EXPECT_NEAR( higher_order_result.gap, expected_gap, 1.e-12 ); + EXPECT_NEAR( overintegrated_reference.gap, expected_gap, 1.e-12 ); + + // The centroid lies outside the penetrating corner, so the one-point force + // is zero. The pointwise active set makes this a nonsmooth integrand, so + // compare lower- and higher-order errors with an order-ten reference over + // the same fixed LOR overlap polygon. The final checks verify + // equal-and-opposite resultant force. + EXPECT_NEAR( single_point.total_absolute_force, 0., 1.e-12 ); + EXPECT_GT( lower_order_result.total_absolute_force, 0. ); + const RealT lower_order_error = + std::abs( lower_order_result.total_absolute_force - overintegrated_reference.total_absolute_force ); + const RealT higher_order_error = + std::abs( higher_order_result.total_absolute_force - overintegrated_reference.total_absolute_force ); + EXPECT_LT( higher_order_error, lower_order_error ); + EXPECT_LT( higher_order_error, 0.15 * overintegrated_reference.total_absolute_force ); + EXPECT_NEAR( lower_order_result.total_z_force, 0., 1.e-12 ); + EXPECT_NEAR( higher_order_result.total_z_force, 0., 1.e-12 ); + EXPECT_NEAR( overintegrated_reference.total_z_force, 0., 1.e-12 ); +} + +/** Verify that multipoint quadrature resolves the nonuniform force on tilted edges. */ +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. + constexpr int production_quadrature_order = 4; + constexpr int reference_quadrature_order = 10; + const auto single_point = runEdgeLocalContactCase( tribol::SINGLE_POINT, production_quadrature_order ); + const auto multi_point = runEdgeLocalContactCase( tribol::MULTI_POINT, production_quadrature_order ); + const auto overintegrated_reference = runEdgeLocalContactCase( tribol::MULTI_POINT, reference_quadrature_order ); + + EXPECT_EQ( single_point.update_error, 0 ); + EXPECT_EQ( multi_point.update_error, 0 ); + EXPECT_EQ( overintegrated_reference.update_error, 0 ); + + EXPECT_EQ( single_point.number_of_active_pairs, 1 ); + EXPECT_EQ( multi_point.number_of_active_pairs, 1 ); + EXPECT_EQ( overintegrated_reference.number_of_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( overintegrated_reference.gap, expected_gap, 1.e-12 ); + + // A fourth-order segment rule integrates the linear-basis penalty force + // exactly. Compare it with an independent order-ten run over the same fixed + // LOR overlap segment rather than storing implementation-specific totals. + EXPECT_NEAR( single_point.total_absolute_force, overintegrated_reference.total_absolute_force, 1.e-12 ); + EXPECT_NEAR( multi_point.total_absolute_force, overintegrated_reference.total_absolute_force, 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.first_mesh_node_force_magnitudes[0] - single_point.first_mesh_node_force_magnitudes[1] ); + const RealT multi_point_imbalance = + std::abs( multi_point.first_mesh_node_force_magnitudes[0] - multi_point.first_mesh_node_force_magnitudes[1] ); + + const RealT reference_imbalance = std::abs( overintegrated_reference.first_mesh_node_force_magnitudes[0] - + overintegrated_reference.first_mesh_node_force_magnitudes[1] ); + + EXPECT_GT( multi_point_imbalance, single_point_imbalance ); + EXPECT_NEAR( multi_point_imbalance, reference_imbalance, 1.e-12 ); +} + TEST_F( CommonPlaneTest, separation_gap_check ) { this->m_mesh.mortarMeshId = 0; @@ -671,11 +1001,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..0f829400 100644 --- a/src/tests/tribol_enforcement_options.cpp +++ b/src/tests/tribol_enforcement_options.cpp @@ -162,6 +162,32 @@ TEST_F( EnforcementOptionsTest, penalty_kinematic_constant_error ) delete mesh; } +/** Verify that the requested CommonPlane integration rule and order are retained. */ +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..b0640709 100644 --- a/src/tests/tribol_iso_integ.cpp +++ b/src/tests/tribol_iso_integ.cpp @@ -18,6 +18,73 @@ using RealT = tribol::RealT; +namespace { + +/** Return the factorial of a nonnegative integer as a Tribol scalar. */ +RealT Factorial( int value ) +{ + RealT result = 1.; + for ( int factor = 2; factor <= value; ++factor ) { + result *= factor; + } + return result; +} + +/** Return the exact x^p y^q moment on the area-one-half unit right triangle. */ +RealT ReferenceTriangleMoment( int first_exponent, int second_exponent ) +{ + return Factorial( first_exponent ) * Factorial( second_exponent ) / Factorial( first_exponent + second_exponent + 2 ); +} + +/** Return the exact moment normalized to the unit-sum triangle-weight convention. */ +RealT NormalizedReferenceTriangleMoment( int first_exponent, int second_exponent ) +{ + return 2. * ReferenceTriangleMoment( first_exponent, second_exponent ); +} + +/** + * @brief Evaluate one polynomial moment with either CommonPlane triangle-rule family. + * + * @param use_legacy_rule Whether to use the historical rule instead of the symmetric rule + * @param order Requested quadrature order + * @param first_exponent Exponent of the first reference coordinate + * @param second_exponent Exponent of the second reference coordinate + * @return Numerically integrated, unit-sum-normalized moment + */ +RealT EvaluateTriangleRuleMoment( bool use_legacy_rule, int order, int first_exponent, int second_exponent ) +{ + RealT quadrature_weights[tribol::max_symmetric_triangle_qpts] = { 0. }; + RealT reference_coordinates[2 * tribol::max_symmetric_triangle_qpts] = { 0. }; + const int number_of_quadrature_points = + use_legacy_rule ? tribol::GetLegacyTriangleRule( order, quadrature_weights, reference_coordinates ) + : tribol::GetCommonPlaneTriangleRule( order, quadrature_weights, reference_coordinates ); + + RealT value = 0.; + for ( int quadrature_point = 0; quadrature_point < number_of_quadrature_points; ++quadrature_point ) { + value += quadrature_weights[quadrature_point] * + std::pow( reference_coordinates[2 * quadrature_point], first_exponent ) * + std::pow( reference_coordinates[2 * quadrature_point + 1], second_exponent ); + } + return value; +} + +/** Evaluate one monomial moment with a CommonPlane segment quadrature rule. */ +RealT EvaluateSegmentRuleMoment( int order, int exponent ) +{ + RealT quadrature_weights[tribol::max_segment_gauss_legendre_qpts] = { 0. }; + RealT reference_coordinates[tribol::max_segment_gauss_legendre_qpts] = { 0. }; + const int number_of_quadrature_points = + tribol::GetCommonPlaneSegmentRule( order, quadrature_weights, reference_coordinates ); + + RealT value = 0.; + for ( int quadrature_point = 0; quadrature_point < number_of_quadrature_points; ++quadrature_point ) { + value += quadrature_weights[quadrature_point] * std::pow( reference_coordinates[quadrature_point], exponent ); + } + return value; +} + +} // namespace + /*! * Test fixture class with some setup necessary to use the * triangular decomposition of a quadrilateral with integration @@ -242,6 +309,86 @@ TEST_F( IsoIntegTest, nonaffine ) EXPECT_EQ( convrg, true ); } +/** Verify low-order compatibility between the legacy and symmetric triangle rules. */ +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( EvaluateTriangleRuleMoment( true, order, 0, 0 ), EvaluateTriangleRuleMoment( false, order, 0, 0 ), + 2.e-10 ); + EXPECT_NEAR( EvaluateTriangleRuleMoment( true, order, 2, 0 ), EvaluateTriangleRuleMoment( false, order, 2, 0 ), + 2.e-10 ); + EXPECT_NEAR( EvaluateTriangleRuleMoment( true, order, 1, 1 ), EvaluateTriangleRuleMoment( false, order, 1, 1 ), + 2.e-10 ); + } +} + +/** Verify polynomial exactness for every supported CommonPlane segment rule. */ +TEST( SegmentRuleTest, integrates_polynomials_through_each_supported_order ) +{ + // An n-point Gauss-Legendre rule integrates every polynomial through degree + // 2n-1 exactly. Exercising each exposed order also validates its point count. + constexpr RealT integration_tolerance = 2.e-14; + for ( int order = 2; order <= 10; ++order ) { + for ( int exponent = 0; exponent <= 2 * order - 1; ++exponent ) { + const RealT exact_moment = 1. / static_cast( exponent + 1 ); + EXPECT_NEAR( EvaluateSegmentRuleMoment( order, exponent ), exact_moment, integration_tolerance ) + << "quadrature order " << order << ", polynomial exponent " << exponent; + } + } +} + +/** Verify polynomial exactness for every supported symmetric triangle rule. */ +TEST( TriangleRuleTest, integrates_polynomials_through_each_supported_order ) +{ + // The symmetric order-p rule must reproduce every reference-triangle monomial + // whose total degree does not exceed p. This validates every table from 2–10. + constexpr RealT integration_tolerance = 5.e-14; + for ( int order = 2; order <= 10; ++order ) { + for ( int first_exponent = 0; first_exponent <= order; ++first_exponent ) { + for ( int second_exponent = 0; second_exponent <= order - first_exponent; ++second_exponent ) { + const RealT exact_moment = NormalizedReferenceTriangleMoment( first_exponent, second_exponent ); + EXPECT_NEAR( EvaluateTriangleRuleMoment( false, order, first_exponent, second_exponent ), exact_moment, + integration_tolerance ) + << "quadrature order " << order << ", polynomial exponents " << first_exponent << " and " + << second_exponent; + } + } + } +} + +/** Verify polygon-fan integration with the highest supported triangle rule. */ +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 spatial_dimension = 3; + constexpr int number_of_nodes = 3; + RealT coordinates[spatial_dimension * number_of_nodes] = { 0., 0., 0., 1., 0., 0., 0., 1., 0. }; + + tribol::SurfaceContactElem contact_element( spatial_dimension, coordinates, coordinates, coordinates, number_of_nodes, + number_of_nodes, nullptr, nullptr, 0, 0 ); + tribol::IntegPts integration_points; + tribol::GaussPolyIntTri( contact_element, integration_points, 10 ); + + RealT area = 0.; + RealT seventh_third_moment = 0.; + for ( int integration_point = 0; integration_point < integration_points.numIPs; ++integration_point ) { + const RealT x_coordinate = integration_points.xy[spatial_dimension * integration_point]; + const RealT y_coordinate = integration_points.xy[spatial_dimension * integration_point + 1]; + area += integration_points.wts[integration_point]; + seventh_third_moment += + integration_points.wts[integration_point] * std::pow( x_coordinate, 7 ) * std::pow( y_coordinate, 3 ); + } + + EXPECT_EQ( integration_points.numIPs, 75 ); + EXPECT_NEAR( area, 0.5, 1.e-14 ); + // Integral of x^7 y^3 over the physical unit right triangle. + EXPECT_NEAR( seventh_third_moment, ReferenceTriangleMoment( 7, 3 ), 1.e-14 ); + EXPECT_NEAR( EvaluateTriangleRuleMoment( 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 50ccb804..395b5b50 100644 --- a/src/tests/tribol_mfem_common_plane.cpp +++ b/src/tests/tribol_mfem_common_plane.cpp @@ -347,7 +347,9 @@ 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 ) ) ); + std::make_tuple( 2, tribol::KINEMATIC_ELEMENT ), + std::make_tuple( 3, tribol::KINEMATIC_CONSTANT ), + std::make_tuple( 4, tribol::KINEMATIC_CONSTANT ) ) ); /** Verify identity mappings for every supported contact face element type. */ TEST( MfemCommonPlaneParentFaceData, MapsSupportedFaceElementTypesAndRejectsInvalidInput ) @@ -591,9 +593,40 @@ TEST_P( MfemCommonPlaneParentFaceDataTest, MapsRedecomposedQuadrilateralFacesToP second_parent_region_coordinate * lor_factor + first_parent_region_coordinate; } + // Evaluate the native parent basis at the mapped point. Partition of unity + // and reproduction of the LOR-face center prove that the device-side + // basis representation follows MFEM's native parent-node ordering. + tribol::RealT parent_basis_values[tribol::ParentFaceData::max_parent_face_nodes] = { 0.0 }; + validation_result |= + !mesh_view.evaluateParentFaceBasis( face_id, mapped_parent_center, parent_basis_values ) ? 1024 : 0; + tribol::RealT basis_value_sum = 0.0; + const int number_of_parent_nodes = parent_face_data.m_parent_node_counts[face_id]; + for ( int parent_node = 0; parent_node < number_of_parent_nodes; ++parent_node ) { + basis_value_sum += parent_basis_values[parent_node]; + } + validation_result |= std::abs( basis_value_sum - 1.0 ) > comparison_tolerance ? 2048 : 0; + + tribol::RealT parent_position[3] = { 0.0, 0.0, 0.0 }; + mesh_view.evaluateParentFaceFields( face_id, parent_basis_values, parent_position, nullptr ); + + // The two cube boundaries have opposite MFEM face orientations. Their + // known affine coordinate fields provide an independent reproduction + // check without relying on Tribol's derived face-centroid storage. + const tribol::RealT expected_parent_position[3] = { + mesh_id == first_mesh_id ? mapped_parent_center[0] : mapped_parent_center[1], + mesh_id == first_mesh_id ? mapped_parent_center[1] : mapped_parent_center[0], + mesh_id == first_mesh_id ? 1.0 : 1.0 + initial_separation }; + for ( int coordinate_component = 0; coordinate_component < mesh_view.spatialDimension(); + ++coordinate_component ) { + validation_result |= std::abs( parent_position[coordinate_component] - + expected_parent_position[coordinate_component] ) > comparison_tolerance + ? 4096 + : 0; + } + const tribol::RealT exterior_lor_point[2] = { 1.25, 0.5 }; validation_result |= - mesh_view.mapToParentReference( face_id, exterior_lor_point, mapped_parent_center ) ? 256 : 0; + mesh_view.mapToParentReference( face_id, exterior_lor_point, mapped_parent_center ) ? 8192 : 0; face_validation_results_view[face_id] = validation_result; } ); diff --git a/src/tribol/common/Parameters.hpp b/src/tribol/common/Parameters.hpp index cde037af..4bee0f06 100644 --- a/src/tribol/common/Parameters.hpp +++ b/src/tribol/common/Parameters.hpp @@ -244,8 +244,10 @@ 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 the overlap centroid. + FULL_TRI_DECOMP, ///< Legacy name for triangle-decomposition overlap integration. + MULTI_POINT = FULL_TRI_DECOMP, ///< Multipoint integration over the overlap interval or polygon. + AUTO_INTEGRATION, ///< Select the rule and order from the registered parent-face order. NUM_INTEG_RULES }; @@ -453,6 +455,13 @@ struct PenaltyEnforcementOptions { PenaltyConstraintType constraint_type; KinematicPenaltyCalculation kinematic_calculation; RatePenaltyCalculation rate_calculation; + PolyInteg common_plane_rule{ AUTO_INTEGRATION }; ///< CommonPlane overlap integration rule. + /** Triangle or segment quadrature order used by MULTI_POINT and ignored by other rules. */ + int common_plane_quadrature_order{ 3 }; + /** Dimensionless stability limit supplied by the application for its explicit integrator. */ + RealT explicit_integrator_stability_factor{ 0.0 }; + /** Whether the application registered its explicit-integrator stability limit. */ + bool explicit_integrator_stability_factor_set{ false }; 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..e409e367 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,63 @@ 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] face_coordinates pointer to stacked physical coordinates of the face vertices + * \param [in] evaluation_point_x x-coordinate of the evaluation point + * \param [in] evaluation_point_y y-coordinate of the evaluation point + * \param [in] evaluation_point_z z-coordinate of the evaluation point; ignored for physical edges + * \param [in] number_of_nodes number of nodes on the face + * \param [in] basis_index local node id whose basis function is to be evaluated + * \param [out] basis_value 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 face_coordinates, + const RealT evaluation_point_x, const RealT evaluation_point_y, + const RealT evaluation_point_z, const int number_of_nodes, + const int basis_index, RealT& basis_value ) +{ + if ( number_of_nodes == 2 ) { + SegmentBasis( face_coordinates, evaluation_point_x, evaluation_point_y, basis_index, basis_value ); + return; + } + +#ifdef TRIBOL_USE_HOST + SLIC_ERROR_IF( number_of_nodes != 3 && number_of_nodes != 4, + "EvalBasisOnPhysicalFace(): only linear triangle and quadrilateral faces are supported." ); +#endif + + constexpr int max_nodes_per_face = 4; + RealT x_coordinates[max_nodes_per_face] = { 0., 0., 0., 0. }; + RealT y_coordinates[max_nodes_per_face] = { 0., 0., 0., 0. }; + RealT z_coordinates[max_nodes_per_face] = { 0., 0., 0., 0. }; + for ( int node_index = 0; node_index < number_of_nodes; ++node_index ) { + x_coordinates[node_index] = face_coordinates[3 * node_index]; + y_coordinates[node_index] = face_coordinates[3 * node_index + 1]; + z_coordinates[node_index] = face_coordinates[3 * node_index + 2]; + } + + RealT evaluation_point[3] = { evaluation_point_x, evaluation_point_y, evaluation_point_z }; + RealT reference_coordinates[2] = { 0., 0. }; + InvIso( evaluation_point, x_coordinates, y_coordinates, z_coordinates, number_of_nodes, reference_coordinates ); + + if ( number_of_nodes == 4 ) { + LinIsoQuadShapeFunc( reference_coordinates[0], reference_coordinates[1], basis_index, basis_value ); + } else { + LinIsoTriShapeFunc( reference_coordinates[0], reference_coordinates[1], basis_index, basis_value ); + } +} + //----------------------------------------------------------------------------- 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..a13b3627 100644 --- a/src/tribol/integ/Integration.cpp +++ b/src/tribol/integ/Integration.cpp @@ -217,125 +217,58 @@ int NumTWBPointsPerTri( int order ) } //------------------------------------------------------------------------------ -void GaussPolyIntTri( SurfaceContactElem const& elem, IntegPts& integ, int k ) +void GaussPolyIntTri( SurfaceContactElem const& elem, IntegPts& integ, int order, 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; - } - - int parentDim = 2; - - 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; + constexpr int reference_dimension = 2; + RealT quadrature_weights[max_symmetric_triangle_qpts] = { 0. }; + RealT reference_coordinates[reference_dimension * max_symmetric_triangle_qpts] = { 0. }; + const int number_of_triangle_points = GetTriangleRule( order, family, quadrature_weights, reference_coordinates ); + if ( number_of_triangle_points == 0 ) { + SLIC_ERROR( "GaussPolyIntTri: requested triangle integration rule is not available." ); + return; } - - // compute area centroid of polygon - RealT xTri[3] = { 0., 0., 0. }; - RealT yTri[3] = { 0., 0., 0. }; - RealT zTri[3] = { 0., 0., 0. }; - PolyAreaCentroid( elem.overlapCoords, elem.dim, elem.numPolyVert, xTri[2], yTri[2], zTri[2] ); - - // populate xy array - for ( int j = 0; j < elem.numPolyVert; ++j ) { - // group triangle coordinates - int triId = j; - int triIdPlusOne = ( j == ( elem.numPolyVert - 1 ) ) ? 0 : triId + 1; - xTri[0] = elem.overlapCoords[elem.dim * triId]; - yTri[0] = elem.overlapCoords[elem.dim * triId + 1]; - zTri[0] = elem.overlapCoords[elem.dim * triId + 2]; - xTri[1] = elem.overlapCoords[elem.dim * triIdPlusOne]; - yTri[1] = elem.overlapCoords[elem.dim * triIdPlusOne + 1]; - zTri[1] = elem.overlapCoords[elem.dim * triIdPlusOne + 2]; - - // compute area of triangle - RealT area = Area3DTri( xTri, yTri, zTri ); - - for ( int k = 0; k < numTriPoints; ++k ) { - // 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; - - // group parent space ip coordinates - RealT xi[2]; - xi[0] = coords[parentDim * k]; - xi[1] = coords[parentDim * k + 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]; + const int number_of_total_points = number_of_triangle_points * elem.numPolyVert; + integ.initialize( 3, number_of_total_points ); + + // The overlap polygon is convex. Its centroid and each polygon edge form a + // nonoverlapping triangle, so integrating the fan covers the full overlap. + RealT triangle_x_coordinates[3] = { 0., 0., 0. }; + RealT triangle_y_coordinates[3] = { 0., 0., 0. }; + RealT triangle_z_coordinates[3] = { 0., 0., 0. }; + PolyAreaCentroid( elem.overlapCoords, elem.dim, elem.numPolyVert, triangle_x_coordinates[2], + triangle_y_coordinates[2], triangle_z_coordinates[2] ); + + for ( int overlap_vertex = 0; overlap_vertex < elem.numPolyVert; ++overlap_vertex ) { + const int next_overlap_vertex = overlap_vertex == elem.numPolyVert - 1 ? 0 : overlap_vertex + 1; + triangle_x_coordinates[0] = elem.overlapCoords[elem.dim * overlap_vertex]; + triangle_y_coordinates[0] = elem.overlapCoords[elem.dim * overlap_vertex + 1]; + triangle_z_coordinates[0] = elem.overlapCoords[elem.dim * overlap_vertex + 2]; + triangle_x_coordinates[1] = elem.overlapCoords[elem.dim * next_overlap_vertex]; + triangle_y_coordinates[1] = elem.overlapCoords[elem.dim * next_overlap_vertex + 1]; + triangle_z_coordinates[1] = elem.overlapCoords[elem.dim * next_overlap_vertex + 2]; + + const RealT triangle_area = Area3DTri( triangle_x_coordinates, triangle_y_coordinates, triangle_z_coordinates ); + + for ( int integration_point = 0; integration_point < number_of_triangle_points; ++integration_point ) { + // Tribol stores the physical integration weight, so apply the fan-triangle + // area to each unit-sum reference weight. + integ.wts[number_of_triangle_points * overlap_vertex + integration_point] = + triangle_area * quadrature_weights[integration_point]; + + RealT triangle_reference_coordinates[2]; + triangle_reference_coordinates[0] = reference_coordinates[reference_dimension * integration_point]; + triangle_reference_coordinates[1] = reference_coordinates[reference_dimension * integration_point + 1]; + + RealT physical_coordinates[3]; + FwdMapLinTri( triangle_reference_coordinates, triangle_x_coordinates, triangle_y_coordinates, + triangle_z_coordinates, physical_coordinates ); + + const int output_offset = integ.ipDim * ( number_of_triangle_points * overlap_vertex + integration_point ); + integ.xy[output_offset] = physical_coordinates[0]; + integ.xy[output_offset + 1] = physical_coordinates[1]; + integ.xy[output_offset + 2] = physical_coordinates[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..6a84c500 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 ); +/** @brief Selector for the triangle quadrature family used by GaussPolyIntTri(). */ +enum TriangleQuadratureRuleFamily +{ + TRI_RULE_LEGACY, ///< Historical Tribol rules used by mortar integration. + TRI_RULE_SYMMETRIC ///< Symmetric rules used by multipoint CommonPlane integration. +}; + +/** 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 @@ -119,14 +131,16 @@ 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] order 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 order, + TriangleQuadratureRuleFamily family = TRI_RULE_SYMMETRIC ); /*! * @@ -174,92 +188,670 @@ int NumTWBPointsPerTri( int order ); // Implementations //----------------------------------------------------------------------------- -template <> -TRIBOL_HOST_DEVICE inline void EvalWeakFormIntegral( SurfaceContactElem const& elem, - RealT* const integ1, - RealT* const integ2 ) +/** + * @brief Compute the centroid used to decompose a CommonPlane overlap cell. + * + * @param elem Contact element containing the overlap interval or polygon + * @param overlap_centroid Output physical overlap centroid + */ +TRIBOL_HOST_DEVICE inline void GetCommonPlaneOverlapCentroid( SurfaceContactElem const& elem, + RealT overlap_centroid[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. }; + overlap_centroid[0] = 0.; + overlap_centroid[1] = 0.; + overlap_centroid[2] = 0.; + if ( elem.dim == 2 ) { - VertexAvgCentroid( elem.overlapCoords, elem.dim, elem.numPolyVert, cx[0], cx[1], cx[2] ); + VertexAvgCentroid( elem.overlapCoords, elem.dim, elem.numPolyVert, overlap_centroid[0], overlap_centroid[1], + overlap_centroid[2] ); } else { - PolyAreaCentroid( elem.overlapCoords, elem.dim, elem.numPolyVert, cx[0], cx[1], cx[2] ); + PolyAreaCentroid( elem.overlapCoords, elem.dim, elem.numPolyVert, overlap_centroid[0], overlap_centroid[1], + overlap_centroid[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] ); + for ( int orbit = 0; orbit < rule.num_general_orbits; ++orbit ) { + const RealT first_coordinate = rule.orbits[orbit_coordinate_index++]; + const RealT second_coordinate = rule.orbits[orbit_coordinate_index++]; + const RealT third_coordinate = rule.orbits[orbit_coordinate_index++]; + const RealT quadrature_weight = symmetric_triangle_weight_scale * rule.weights[weight_index++]; + + const RealT coordinate_permutations[6][2] = { + { second_coordinate, third_coordinate }, { third_coordinate, second_coordinate }, + { first_coordinate, third_coordinate }, { third_coordinate, first_coordinate }, + { first_coordinate, second_coordinate }, { second_coordinate, first_coordinate } }; + for ( int permutation_index = 0; permutation_index < 6; ++permutation_index ) { + quadrature_weights[quadrature_point_index] = quadrature_weight; + reference_coordinates[2 * quadrature_point_index] = coordinate_permutations[permutation_index][0]; + reference_coordinates[2 * quadrature_point_index + 1] = coordinate_permutations[permutation_index][1]; + ++quadrature_point_index; } - } 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 quadrature_point_index; +} + +} // namespace detail + +/*! + * \brief Returns the built-in higher-order symmetric triangle quadrature rule. + * + * \param [in] order requested rule order + * \param [out] quadrature_weights quadrature weights normalized so they sum to 1 on a triangle + * \param [out] reference_coordinates quadrature coordinates stored as stacked (xi, eta) pairs + * \return number of quadrature points in the selected rule + * + * \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* quadrature_weights, + RealT* reference_coordinates ) +{ + 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, quadrature_weights, reference_coordinates ); +} + +/*! + * \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] quadrature_weights quadrature weights normalized so they sum to 1 on a triangle + * \param [out] reference_coordinates quadrature coordinates stored as stacked (xi, eta) pairs + * \return number of quadrature points in the selected rule + */ +TRIBOL_HOST_DEVICE inline int GetTriangleRule( int order, TriangleQuadratureRuleFamily family, + RealT* quadrature_weights, RealT* reference_coordinates ) +{ + switch ( family ) { + case TRI_RULE_LEGACY: + return GetLegacyTriangleRule( order, quadrature_weights, reference_coordinates ); + case TRI_RULE_SYMMETRIC: + return GetCommonPlaneTriangleRule( order, quadrature_weights, reference_coordinates ); + default: +#ifdef TRIBOL_USE_HOST + SLIC_ERROR( "GetTriangleRule(): unsupported triangle rule family." ); +#endif + return 0; + } +} + +/*! + * \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] + * \return number of quadrature points in the selected rule + */ +TRIBOL_HOST_DEVICE inline int GetCommonPlaneSegmentRule( int order, RealT* wts, RealT* coords ) +{ + 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.29234407200302780068067613359605864e-01; + coords[2] = 2.97077424311301416546967632731488033e-01; + coords[3] = 5.00000000000000000000000000000000000e-01; + coords[4] = 7.02922575688698583453032367268511967e-01; + coords[5] = 8.70765592799697219931932386640394136e-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; + } +} + +/** + * @brief Integrate both CommonPlane face bases with a multipoint rule. + * + * Two-dimensional overlap intervals use Gauss-Legendre quadrature. Three-dimensional + * overlap polygons are decomposed into a nonoverlapping fan about their centroid, + * and every fan triangle uses the requested symmetric triangle rule. + * + * @param elem Contact element containing the overlap interval or polygon + * @param quadrature_order Requested integration order in the supported range [2,10] + * @param first_face_integrals Accumulated first-face basis integrals + * @param second_face_integrals Accumulated second-face basis integrals + */ +TRIBOL_HOST_DEVICE inline void EvalWeakFormIntegralCommonPlaneMultiPoint( SurfaceContactElem const& elem, + const int quadrature_order, + RealT* const first_face_integrals, + RealT* const second_face_integrals ) +{ + if ( elem.dim == 2 ) { + RealT quadrature_weights[max_segment_gauss_legendre_qpts] = { 0. }; + RealT reference_coordinates[max_segment_gauss_legendre_qpts] = { 0. }; + const int number_of_quadrature_points = + GetCommonPlaneSegmentRule( quadrature_order, quadrature_weights, reference_coordinates ); + + const RealT first_endpoint_x = elem.overlapCoords[0]; + const RealT first_endpoint_y = elem.overlapCoords[1]; + const RealT second_endpoint_x = elem.overlapCoords[2]; + const RealT second_endpoint_y = elem.overlapCoords[3]; + const RealT overlap_length = + magnitude( second_endpoint_x - first_endpoint_x, second_endpoint_y - first_endpoint_y ); + + for ( int quadrature_point = 0; quadrature_point < number_of_quadrature_points; ++quadrature_point ) { + const RealT segment_coordinate = reference_coordinates[quadrature_point]; + const RealT first_endpoint_weight = 1. - segment_coordinate; + RealT integration_point[3] = { first_endpoint_weight * first_endpoint_x + segment_coordinate * second_endpoint_x, + first_endpoint_weight * first_endpoint_y + segment_coordinate * second_endpoint_y, + 0. }; + AccumulateCommonPlaneIntegralAtPoint( elem, integration_point, + overlap_length * quadrature_weights[quadrature_point], first_face_integrals, + second_face_integrals ); } + return; } - // 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] ); + RealT quadrature_weights[max_symmetric_triangle_qpts] = { 0. }; + RealT reference_coordinates[2 * max_symmetric_triangle_qpts] = { 0. }; + const int number_of_quadrature_points = + GetCommonPlaneTriangleRule( quadrature_order, quadrature_weights, reference_coordinates ); + + RealT overlap_centroid[3]; + GetCommonPlaneOverlapCentroid( elem, overlap_centroid ); + + RealT triangle_x_coordinates[3]; + RealT triangle_y_coordinates[3]; + RealT triangle_z_coordinates[3]; + + for ( int overlap_vertex = 0; overlap_vertex < elem.numPolyVert; ++overlap_vertex ) { + const int next_overlap_vertex = overlap_vertex == elem.numPolyVert - 1 ? 0 : overlap_vertex + 1; + triangle_x_coordinates[0] = elem.overlapCoords[elem.dim * overlap_vertex]; + triangle_y_coordinates[0] = elem.overlapCoords[elem.dim * overlap_vertex + 1]; + triangle_z_coordinates[0] = elem.overlapCoords[elem.dim * overlap_vertex + 2]; + triangle_x_coordinates[1] = elem.overlapCoords[elem.dim * next_overlap_vertex]; + triangle_y_coordinates[1] = elem.overlapCoords[elem.dim * next_overlap_vertex + 1]; + triangle_z_coordinates[1] = elem.overlapCoords[elem.dim * next_overlap_vertex + 2]; + triangle_x_coordinates[2] = overlap_centroid[0]; + triangle_y_coordinates[2] = overlap_centroid[1]; + triangle_z_coordinates[2] = overlap_centroid[2]; + + const RealT triangle_area = Area3DTri( triangle_x_coordinates, triangle_y_coordinates, triangle_z_coordinates ); + if ( triangle_area <= 0. ) { + continue; + } + + for ( int quadrature_point = 0; quadrature_point < number_of_quadrature_points; ++quadrature_point ) { + const RealT first_triangle_coordinate = reference_coordinates[2 * quadrature_point]; + const RealT second_triangle_coordinate = reference_coordinates[2 * quadrature_point + 1]; + const RealT third_triangle_coordinate = 1. - first_triangle_coordinate - second_triangle_coordinate; + RealT integration_point[3]; + integration_point[0] = third_triangle_coordinate * triangle_x_coordinates[0] + + first_triangle_coordinate * triangle_x_coordinates[1] + + second_triangle_coordinate * triangle_x_coordinates[2]; + integration_point[1] = third_triangle_coordinate * triangle_y_coordinates[0] + + first_triangle_coordinate * triangle_y_coordinates[1] + + second_triangle_coordinate * triangle_y_coordinates[2]; + integration_point[2] = third_triangle_coordinate * triangle_z_coordinates[0] + + first_triangle_coordinate * triangle_z_coordinates[1] + + second_triangle_coordinate * triangle_z_coordinates[2]; + AccumulateCommonPlaneIntegralAtPoint( elem, integration_point, + triangle_area * quadrature_weights[quadrature_point], first_face_integrals, + second_face_integrals ); + } } +} - return; +/** + * @brief Integrate both CommonPlane face bases with the selected overlap rule. + * + * @param elem Contact element containing the overlap interval or polygon + * @param rule CommonPlane overlap integration rule + * @param quadrature_order Requested multipoint integration order + * @param first_face_integrals Accumulated first-face basis integrals + * @param second_face_integrals Accumulated second-face basis integrals + */ +TRIBOL_HOST_DEVICE inline void EvalWeakFormIntegralCommonPlane( SurfaceContactElem const& elem, const PolyInteg rule, + const int quadrature_order, + RealT* const first_face_integrals, + RealT* const second_face_integrals ) +{ + switch ( rule ) { + case SINGLE_POINT: { + RealT overlap_centroid[3] = { 0., 0., 0. }; + GetCommonPlaneOverlapCentroid( elem, overlap_centroid ); + AccumulateCommonPlaneIntegralAtPoint( elem, overlap_centroid, 1.0, first_face_integrals, second_face_integrals ); + break; + } + case MULTI_POINT: + EvalWeakFormIntegralCommonPlaneMultiPoint( elem, quadrature_order, first_face_integrals, second_face_integrals ); + break; + default: +#ifdef TRIBOL_USE_HOST + SLIC_ERROR( "EvalWeakFormIntegralCommonPlane(): unsupported polygon integration rule." ); +#endif + break; + } +} + +/** + * @brief Evaluate legacy one-point CommonPlane weak-form basis integrals. + * + * @param elem Contact element containing the overlap geometry and face coordinates + * @param integ1 First-face basis integral values + * @param integ2 Second-face basis integral values + */ +template <> +TRIBOL_HOST_DEVICE inline void EvalWeakFormIntegral( SurfaceContactElem const& elem, + RealT* const integ1, + RealT* const integ2 ) +{ + RealT overlap_centroid[3] = { 0., 0., 0. }; + GetCommonPlaneOverlapCentroid( elem, overlap_centroid ); + AccumulateCommonPlaneIntegralAtPoint( elem, overlap_centroid, 1.0, integ1, integ2 ); } } // end namespace tribol diff --git a/src/tribol/interface/tribol.cpp b/src/tribol/interface/tribol.cpp index 97f9ad6f..d20a2fcc 100644 --- a/src/tribol/interface/tribol.cpp +++ b/src/tribol/interface/tribol.cpp @@ -97,6 +97,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 629e9bdd..abb11a60 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/mesh/CouplingScheme.cpp b/src/tribol/mesh/CouplingScheme.cpp index 93db4c81..43135312 100644 --- a/src/tribol/mesh/CouplingScheme.cpp +++ b/src/tribol/mesh/CouplingScheme.cpp @@ -1257,6 +1257,10 @@ void CouplingScheme::allocateMethodData() // dynamically allocate method data object for mortar method switch ( this->m_contactMethod ) { + case COMMON_PLANE: { + this->m_methodData = std::make_unique(); + break; + } case ALIGNED_MORTAR: case MORTAR_WEIGHTS: case SINGLE_MORTAR: { diff --git a/src/tribol/mesh/MeshData.cpp b/src/tribol/mesh/MeshData.cpp index 100cd46d..3bbbc30a 100644 --- a/src/tribol/mesh/MeshData.cpp +++ b/src/tribol/mesh/MeshData.cpp @@ -234,6 +234,12 @@ void MeshData::setVelocity( const RealT* vx, const RealT* vy, const RealT* vz ) m_vel = createNodalVector( vx, vy, vz ); } +//------------------------------------------------------------------------------ +void MeshData::setInverseMass( const RealT* inverse_mass_x, const RealT* inverse_mass_y, const RealT* inverse_mass_z ) +{ + m_inverse_mass = createNodalVector( inverse_mass_x, inverse_mass_y, inverse_mass_z ); +} + //------------------------------------------------------------------------------ void MeshData::setResponse( RealT* rx, RealT* ry, RealT* rz ) { m_response = createNodalVector( rx, ry, rz ); } @@ -604,6 +610,7 @@ MeshData::Viewer::Viewer( MeshData& mesh ) m_ref_position( mesh.m_ref_position ), m_disp( mesh.m_disp ), m_vel( mesh.m_vel ), + m_inverse_mass( mesh.m_inverse_mass ), m_response( mesh.m_response ), m_node_n( mesh.m_node_n ), m_connectivity( mesh.m_connectivity ), diff --git a/src/tribol/mesh/MeshData.hpp b/src/tribol/mesh/MeshData.hpp index 416896f1..0b7c67e1 100644 --- a/src/tribol/mesh/MeshData.hpp +++ b/src/tribol/mesh/MeshData.hpp @@ -14,6 +14,7 @@ // Tribol includes #include "tribol/common/ArrayTypes.hpp" +#include "tribol/common/Atomics.hpp" #include "tribol/common/Parameters.hpp" #include "tribol/utils/DataManager.hpp" @@ -106,9 +107,18 @@ struct ParentFaceData { /** Maximum number of vertices on a supported LOR surface element. */ static constexpr int max_lor_face_vertices{ 4 }; + /** Highest native parent-face order supported by device basis evaluation. */ + static constexpr int max_parent_face_order{ 4 }; + + /** Maximum number of native parent-face nodes for supported quadrilateral faces. */ + static constexpr int max_parent_face_nodes{ ( max_parent_face_order + 1 ) * ( max_parent_face_order + 1 ) }; + /** Polynomial order of the native parent coordinate finite element. */ Array1DView m_parent_face_orders; + /** Number of native parent finite-element nodes on each face. */ + Array1DView m_parent_node_counts; + /** Number of vertices used by each child-to-parent reference map. */ Array1DView m_reference_vertex_counts; @@ -120,6 +130,27 @@ struct ParentFaceData { */ Array2DView m_parent_reference_vertex_coordinates; + /** Coefficients that evaluate each native parent-face nodal basis function. */ + Array2DView m_parent_basis_coefficients; + + /** Native parent-face nodal coordinates in node-major ordering. */ + Array2DView m_parent_positions; + + /** Native parent-face nodal velocities in node-major ordering, when registered. */ + Array2DView m_parent_velocities; + + /** Native parent-face inverse diagonal masses in node-major ordering, when registered. */ + Array2DView m_parent_inverse_masses; + + /** Source-rank vector degree-of-freedom identifier for each parent-face node and component. */ + Array2DView m_parent_vector_dof_ids; + + /** Largest source-rank vector degree-of-freedom count represented by these faces. */ + IndexT m_parent_vector_dof_count{ 0 }; + + /** Element-local native parent-face response accumulated by contact. */ + Array2DView m_parent_responses; + /** * @brief Return whether native parent-face mapping data are available. * @@ -132,6 +163,34 @@ struct ParentFaceData { m_parent_reference_vertex_coordinates.shape()[0] == number_of_faces && m_parent_reference_vertex_coordinates.shape()[1] == max_lor_face_vertices * max_reference_dimension; } + + /** + * @brief Return whether native parent-face basis and field data are available. + * + * @return true when all data required for parent-space position and force evaluation are populated + */ + TRIBOL_HOST_DEVICE bool hasParentFields() const + { + return isValid() && !m_parent_node_counts.empty() && !m_parent_basis_coefficients.empty() && + !m_parent_positions.empty() && !m_parent_responses.empty(); + } + + /** + * @brief Return whether native parent-face velocity data are available. + * + * @return true when parent velocities are populated + */ + TRIBOL_HOST_DEVICE bool hasParentVelocity() const { return !m_parent_velocities.empty(); } + + /** + * @brief Return whether native parent-face inverse diagonal masses are available. + * + * @return true when inverse mass values and source degree-of-freedom identifiers are populated + */ + TRIBOL_HOST_DEVICE bool hasParentInverseMass() const + { + return !m_parent_inverse_masses.empty() && !m_parent_vector_dof_ids.empty() && m_parent_vector_dof_count > 0; + } }; class MeshData { @@ -218,6 +277,20 @@ class MeshData { */ TRIBOL_HOST_DEVICE const ParentFaceData& getParentFaceData() const { return m_parent_face_data; } + /** + * @brief Return whether native parent-face basis and field data are registered. + * + * @return true when parent-space position and response evaluation are available + */ + TRIBOL_HOST_DEVICE bool hasParentFaceFields() const { return m_parent_face_data.hasParentFields(); } + + /** + * @brief Return whether native parent-face velocity data are registered. + * + * @return true when parent-space velocity evaluation is available + */ + TRIBOL_HOST_DEVICE bool hasParentFaceVelocity() const { return m_parent_face_data.hasParentVelocity(); } + /** * @brief Map a LOR face reference point to its native parent-face reference point. * @@ -233,6 +306,38 @@ class MeshData { TRIBOL_HOST_DEVICE bool mapToParentReference( IndexT face_id, const RealT* lor_reference_coordinates, RealT* parent_reference_coordinates ) const; + /** + * @brief Evaluate native parent-face nodal basis functions. + * + * @param face_id Tribol surface element identifier + * @param parent_reference_coordinates Native parent-face reference coordinates + * @param basis_values Output array with space for ParentFaceData::max_parent_face_nodes values + * @return true when the face data and reference point are valid + */ + TRIBOL_HOST_DEVICE bool evaluateParentFaceBasis( IndexT face_id, const RealT* parent_reference_coordinates, + RealT* basis_values ) const; + + /** + * @brief Evaluate native parent-face position and optional velocity. + * + * @param face_id Tribol surface element identifier + * @param basis_values Native parent-face basis values + * @param position Output physical position with spatialDimension() components + * @param velocity Optional output velocity with spatialDimension() components + */ + TRIBOL_HOST_DEVICE void evaluateParentFaceFields( IndexT face_id, const RealT* basis_values, RealT* position, + RealT* velocity ) const; + + /** + * @brief Add an element-local force to the native parent-face response. + * + * @param face_id Tribol surface element identifier + * @param basis_values Native parent-face basis values + * @param force Physical force components to scatter + */ + TRIBOL_HOST_DEVICE void addParentFaceResponse( IndexT face_id, const RealT* basis_values, + const RealT* force ) const; + /** * @brief Spatial dimension of the mesh * @@ -322,6 +427,25 @@ class MeshData { */ TRIBOL_HOST_DEVICE const MultiViewArrayView& getVelocity() const { return m_vel; } + /** + * @brief Return whether component-wise inverse diagonal nodal masses are registered. + * + * @return true when inverse mass data are available + */ + TRIBOL_HOST_DEVICE bool hasInverseMass() const { return !m_inverse_mass.empty(); } + + /** + * @brief Return one component of a registered inverse diagonal nodal mass. + * + * @param node_id Contact mesh node identifier + * @param component Vector component identifier + * @return Registered inverse mass value + */ + TRIBOL_HOST_DEVICE RealT getInverseMass( IndexT node_id, int component ) const + { + return m_inverse_mass[component][node_id]; + } + /** * @brief Is the nodal response vector populated? * @@ -484,6 +608,9 @@ class MeshData { /// Array of views of nodal velocity data const MultiViewArrayView m_vel; + /// Array of views of component-wise inverse diagonal nodal mass data + const MultiViewArrayView m_inverse_mass; + /// Array of views of nodal response data const MultiViewArrayView m_response; @@ -584,6 +711,20 @@ class MeshData { */ void setParentFaceData( const ParentFaceData& parent_face_data ) { m_parent_face_data = parent_face_data; } + /** + * @brief Return whether native parent-face basis and field data are registered. + * + * @return true when parent-space position and response evaluation are available + */ + bool hasParentFaceFields() const { return m_parent_face_data.hasParentFields(); } + + /** + * @brief Get native parent-face mapping and field data. + * + * @return Non-owning native parent-face data views + */ + const ParentFaceData& getParentFaceData() const { return m_parent_face_data; } + /** * @brief Marker which can indicate mesh validity * @@ -692,6 +833,18 @@ class MeshData { */ void setVelocity( const RealT* vx, const RealT* vy, const RealT* vz ); + /** + * @brief Set component-wise inverse diagonal nodal masses. + * + * A zero value denotes a constrained velocity degree of freedom with no + * contribution to the explicit contact stability operator. + * + * @param inverse_mass_x Inverse mass for x velocity degrees of freedom + * @param inverse_mass_y Inverse mass for y velocity degrees of freedom + * @param inverse_mass_z Inverse mass for z velocity degrees of freedom in three dimensions + */ + void setInverseMass( const RealT* inverse_mass_x, const RealT* inverse_mass_y, const RealT* inverse_mass_z ); + /** * @brief Is the velocity vector populated? * @@ -700,6 +853,13 @@ class MeshData { */ bool hasVelocity() const { return !m_vel.empty(); } + /** + * @brief Return whether component-wise inverse diagonal nodal masses are registered. + * + * @return true when inverse mass data are available + */ + bool hasInverseMass() const { return !m_inverse_mass.empty(); } + /** * @brief Set the pointers to the nodal response data * @@ -770,6 +930,7 @@ class MeshData { MultiArrayView m_ref_position; ///< Reference coordinates of nodes in mesh MultiArrayView m_disp; ///< Nodal displacements MultiArrayView m_vel; ///< Nodal velocity + MultiArrayView m_inverse_mass; ///< Component-wise inverse diagonal nodal masses MultiArrayView m_response; ///< Nodal responses (forces) ArrayT m_node_n; ///< Outward unit node normals @@ -969,6 +1130,118 @@ TRIBOL_HOST_DEVICE inline bool MeshData::Viewer::mapToParentReference( IndexT fa return true; } +//------------------------------------------------------------------------------ +TRIBOL_HOST_DEVICE inline bool MeshData::Viewer::evaluateParentFaceBasis( IndexT face_id, + const RealT* parent_reference_coordinates, + RealT* basis_values ) const +{ + if ( !hasParentFaceFields() || face_id < 0 || face_id >= numberOfElements() || + parent_reference_coordinates == nullptr || basis_values == nullptr || + face_id >= m_parent_face_data.m_parent_node_counts.size() || + face_id >= m_parent_face_data.m_parent_basis_coefficients.shape()[0] ) { + return false; + } + + const int parent_order = m_parent_face_data.m_parent_face_orders[face_id]; + const int number_of_parent_nodes = m_parent_face_data.m_parent_node_counts[face_id]; + const InterfaceElementType parent_face_geometry = getElementType(); + if ( parent_order < 1 || parent_order > ParentFaceData::max_parent_face_order || number_of_parent_nodes < 1 || + number_of_parent_nodes > ParentFaceData::max_parent_face_nodes ) { + return false; + } + + RealT monomial_values[ParentFaceData::max_parent_face_nodes] = { 0.0 }; + int number_of_monomials = 0; + if ( parent_face_geometry == LINEAR_EDGE ) { + RealT first_coordinate_power = 1.0; + for ( int first_degree = 0; first_degree <= parent_order; ++first_degree ) { + monomial_values[number_of_monomials++] = first_coordinate_power; + first_coordinate_power *= parent_reference_coordinates[0]; + } + } else if ( parent_face_geometry == LINEAR_TRIANGLE ) { + for ( int total_degree = 0; total_degree <= parent_order; ++total_degree ) { + for ( int second_degree = 0; second_degree <= total_degree; ++second_degree ) { + const int first_degree = total_degree - second_degree; + RealT monomial_value = 1.0; + for ( int exponent = 0; exponent < first_degree; ++exponent ) { + monomial_value *= parent_reference_coordinates[0]; + } + for ( int exponent = 0; exponent < second_degree; ++exponent ) { + monomial_value *= parent_reference_coordinates[1]; + } + monomial_values[number_of_monomials++] = monomial_value; + } + } + } else if ( parent_face_geometry == LINEAR_QUAD ) { + for ( int second_degree = 0; second_degree <= parent_order; ++second_degree ) { + for ( int first_degree = 0; first_degree <= parent_order; ++first_degree ) { + RealT monomial_value = 1.0; + for ( int exponent = 0; exponent < first_degree; ++exponent ) { + monomial_value *= parent_reference_coordinates[0]; + } + for ( int exponent = 0; exponent < second_degree; ++exponent ) { + monomial_value *= parent_reference_coordinates[1]; + } + monomial_values[number_of_monomials++] = monomial_value; + } + } + } else { + return false; + } + + if ( number_of_monomials != number_of_parent_nodes ) { + return false; + } + + for ( int parent_node = 0; parent_node < number_of_parent_nodes; ++parent_node ) { + basis_values[parent_node] = 0.0; + for ( int monomial = 0; monomial < number_of_monomials; ++monomial ) { + const int coefficient_index = parent_node * ParentFaceData::max_parent_face_nodes + monomial; + basis_values[parent_node] += + m_parent_face_data.m_parent_basis_coefficients( face_id, coefficient_index ) * monomial_values[monomial]; + } + } + return true; +} + +//------------------------------------------------------------------------------ +TRIBOL_HOST_DEVICE inline void MeshData::Viewer::evaluateParentFaceFields( IndexT face_id, const RealT* basis_values, + RealT* position, RealT* velocity ) const +{ + const int number_of_parent_nodes = m_parent_face_data.m_parent_node_counts[face_id]; + for ( int component = 0; component < spatialDimension(); ++component ) { + position[component] = 0.0; + if ( velocity != nullptr ) { + velocity[component] = 0.0; + } + } + + for ( int parent_node = 0; parent_node < number_of_parent_nodes; ++parent_node ) { + for ( int component = 0; component < spatialDimension(); ++component ) { + const int field_index = parent_node * spatialDimension() + component; + position[component] += basis_values[parent_node] * m_parent_face_data.m_parent_positions( face_id, field_index ); + if ( velocity != nullptr ) { + velocity[component] += + basis_values[parent_node] * m_parent_face_data.m_parent_velocities( face_id, field_index ); + } + } + } +} + +//------------------------------------------------------------------------------ +TRIBOL_HOST_DEVICE inline void MeshData::Viewer::addParentFaceResponse( IndexT face_id, const RealT* basis_values, + const RealT* force ) const +{ + const int number_of_parent_nodes = m_parent_face_data.m_parent_node_counts[face_id]; + for ( int parent_node = 0; parent_node < number_of_parent_nodes; ++parent_node ) { + for ( int component = 0; component < spatialDimension(); ++component ) { + const int field_index = parent_node * spatialDimension() + component; + atomicAdd( &m_parent_face_data.m_parent_responses( face_id, field_index ), + basis_values[parent_node] * force[component] ); + } + } +} + } // end namespace tribol /// \a ostream operator to print a \a MeshData instance to \a os diff --git a/src/tribol/mesh/MethodCouplingData.cpp b/src/tribol/mesh/MethodCouplingData.cpp index c172e0d2..19bc62f5 100644 --- a/src/tribol/mesh/MethodCouplingData.cpp +++ b/src/tribol/mesh/MethodCouplingData.cpp @@ -181,6 +181,102 @@ void MethodData::storeElemBlockJ( ArrayT&& blockJElemIds, const StackArray< } } +//------------------------------------------------------------------------------ +void CommonPlaneContactData::resize( IndexT number_of_pairs, int spatial_dimension, int allocator_id ) +{ + number_of_pairs_ = number_of_pairs; + spatial_dimension_ = spatial_dimension; + row_capacity_ = number_of_pairs * maximum_rows_per_pair; + + pair_row_counts_ = Array1D( number_of_pairs_, number_of_pairs_, allocator_id ); + pair_evaluation_statuses_ = Array1D( number_of_pairs_, number_of_pairs_, allocator_id ); + row_is_valid_ = Array1D( row_capacity_, row_capacity_, allocator_id ); + row_is_active_ = Array1D( row_capacity_, row_capacity_, allocator_id ); + contact_pair_ids_ = Array1D( row_capacity_, row_capacity_, allocator_id ); + first_face_ids_ = Array1D( row_capacity_, row_capacity_, allocator_id ); + second_face_ids_ = Array1D( row_capacity_, row_capacity_, allocator_id ); + first_basis_counts_ = Array1D( row_capacity_, row_capacity_, allocator_id ); + second_basis_counts_ = Array1D( row_capacity_, row_capacity_, allocator_id ); + row_uses_parent_fields_ = Array1D( row_capacity_, row_capacity_, allocator_id ); + integration_points_ = Array2D( { row_capacity_, 3 }, allocator_id ); + first_parent_reference_coordinates_ = Array2D( { row_capacity_, 2 }, allocator_id ); + second_parent_reference_coordinates_ = Array2D( { row_capacity_, 2 }, allocator_id ); + first_positions_ = Array2D( { row_capacity_, 3 }, allocator_id ); + second_positions_ = Array2D( { row_capacity_, 3 }, allocator_id ); + first_velocities_ = Array2D( { row_capacity_, 3 }, allocator_id ); + second_velocities_ = Array2D( { row_capacity_, 3 }, allocator_id ); + normals_ = Array2D( { row_capacity_, 3 }, allocator_id ); + first_basis_values_ = Array2D( { row_capacity_, ParentFaceData::max_parent_face_nodes }, allocator_id ); + second_basis_values_ = Array2D( { row_capacity_, ParentFaceData::max_parent_face_nodes }, allocator_id ); + integration_weights_ = Array1D( row_capacity_, row_capacity_, allocator_id ); + gaps_ = Array1D( row_capacity_, row_capacity_, allocator_id ); + normal_velocity_gaps_ = Array1D( row_capacity_, row_capacity_, allocator_id ); + penalty_stiffnesses_ = Array1D( row_capacity_, row_capacity_, allocator_id ); + rate_penalty_coefficients_ = Array1D( row_capacity_, row_capacity_, allocator_id ); + tangential_viscous_coefficients_ = Array1D( row_capacity_, row_capacity_, allocator_id ); + + pair_row_counts_.fill( 0 ); + pair_evaluation_statuses_.fill( static_cast( CommonPlanePairEvaluationStatus::UNINITIALIZED ) ); + row_is_valid_.fill( 0 ); + row_is_active_.fill( 0 ); + contact_pair_ids_.fill( -1 ); + first_face_ids_.fill( -1 ); + second_face_ids_.fill( -1 ); + first_basis_counts_.fill( 0 ); + second_basis_counts_.fill( 0 ); + row_uses_parent_fields_.fill( 0 ); + integration_points_.fill( 0.0 ); + first_parent_reference_coordinates_.fill( 0.0 ); + second_parent_reference_coordinates_.fill( 0.0 ); + first_positions_.fill( 0.0 ); + second_positions_.fill( 0.0 ); + first_velocities_.fill( 0.0 ); + second_velocities_.fill( 0.0 ); + normals_.fill( 0.0 ); + first_basis_values_.fill( 0.0 ); + second_basis_values_.fill( 0.0 ); + integration_weights_.fill( 0.0 ); + gaps_.fill( 0.0 ); + normal_velocity_gaps_.fill( 0.0 ); + penalty_stiffnesses_.fill( 0.0 ); + rate_penalty_coefficients_.fill( 0.0 ); + tangential_viscous_coefficients_.fill( 0.0 ); +} + +//------------------------------------------------------------------------------ +CommonPlaneContactData::Viewer CommonPlaneContactData::getView() +{ + return { number_of_pairs_, + row_capacity_, + spatial_dimension_, + pair_row_counts_.view(), + pair_evaluation_statuses_.view(), + row_is_valid_.view(), + row_is_active_.view(), + contact_pair_ids_.view(), + first_face_ids_.view(), + second_face_ids_.view(), + first_basis_counts_.view(), + second_basis_counts_.view(), + row_uses_parent_fields_.view(), + integration_points_.view(), + first_parent_reference_coordinates_.view(), + second_parent_reference_coordinates_.view(), + first_positions_.view(), + second_positions_.view(), + first_velocities_.view(), + second_velocities_.view(), + normals_.view(), + first_basis_values_.view(), + second_basis_values_.view(), + integration_weights_.view(), + gaps_.view(), + normal_velocity_gaps_.view(), + penalty_stiffnesses_.view(), + rate_penalty_coefficients_.view(), + tangential_viscous_coefficients_.view() }; +} + /////////////////////////////////////// // // // Routines for the MortarData class // diff --git a/src/tribol/mesh/MethodCouplingData.hpp b/src/tribol/mesh/MethodCouplingData.hpp index 5aa77c07..8d348650 100644 --- a/src/tribol/mesh/MethodCouplingData.hpp +++ b/src/tribol/mesh/MethodCouplingData.hpp @@ -270,6 +270,266 @@ class MethodData { ArrayT, 2> m_blockJ; ///< element Jacobian contributions by block }; +//------------------------------------------------------------------------------ +/** + * @brief Evaluation status for one CommonPlane overlap-cell row batch. + */ +enum class CommonPlanePairEvaluationStatus : int +{ + UNINITIALIZED, ///< The overlap cell has not been evaluated. + VALID, ///< Every requested row was generated successfully. + INVALID_PARENT_DATA, ///< Native parent data are missing or inconsistent between the two faces. + INVALID_PARENT_MAPPING, ///< At least one quadrature point could not be mapped to a parent face. + DEGENERATE_OVERLAP, ///< The accepted overlap cell has no positive integration measure. + INCONSISTENT_NORMAL ///< The CommonPlane normal is invalid or inconsistently oriented. +}; + +/** + * @brief Device-resident CommonPlane quadrature rows shared by explicit operators. + * + * Each active CommonPlane face pair owns a fixed-capacity range of rows. The + * pair-local row index gives a deterministic attempt-local identity without a + * host prefix sum. Only the first entry in each range is used for a one-point + * rule. Multipoint rules use one range entry per segment point or per + * triangle-fan point in a three-dimensional overlap polygon. + */ +class CommonPlaneContactData : public MethodData { + public: + /** Maximum number of overlap polygon vertices supported by CommonPlane geometry. */ + static constexpr int maximum_overlap_vertices{ 10 }; + + /** Maximum number of quadrature points in a supported triangle rule. */ + static constexpr int maximum_triangle_quadrature_points{ 25 }; + + /** Maximum number of rows reserved for one accepted overlap cell. */ + static constexpr int maximum_rows_per_pair{ maximum_overlap_vertices * maximum_triangle_quadrature_points }; + + /** + * @brief Non-owning device views of a CommonPlane quadrature row batch. + */ + struct Viewer { + /** Number of active CommonPlane face pairs represented by this batch. */ + IndexT number_of_pairs{ 0 }; + + /** Number of allocated row slots in this batch. */ + IndexT row_capacity{ 0 }; + + /** Spatial dimension of the coupling scheme. */ + int spatial_dimension{ 0 }; + + /** Number of generated rows for each face pair. */ + Array1DView pair_row_counts; + + /** CommonPlanePairEvaluationStatus value for each face pair. */ + Array1DView pair_evaluation_statuses; + + /** One for a generated row and zero for an unused row slot. */ + Array1DView row_is_valid; + + /** One when the row satisfies the normal contact activation criterion. */ + Array1DView row_is_active; + + /** Stable active-pair identifier for each generated row. */ + Array1DView contact_pair_ids; + + /** First Tribol LOR face identifier for each generated row. */ + Array1DView first_face_ids; + + /** Second Tribol LOR face identifier for each generated row. */ + Array1DView second_face_ids; + + /** Number of basis values on the first field face for each generated row. */ + Array1DView first_basis_counts; + + /** Number of basis values on the second field face for each generated row. */ + Array1DView second_basis_counts; + + /** One when rows scatter directly to native parent-face response storage. */ + Array1DView row_uses_parent_fields; + + /** Physical CommonPlane integration-point coordinates. */ + Array2DView integration_points; + + /** First native parent-face reference coordinates. */ + Array2DView first_parent_reference_coordinates; + + /** Second native parent-face reference coordinates. */ + Array2DView second_parent_reference_coordinates; + + /** First face position evaluated at each integration point. */ + Array2DView first_positions; + + /** Second face position evaluated at each integration point. */ + Array2DView second_positions; + + /** First face velocity evaluated at each integration point. */ + Array2DView first_velocities; + + /** Second face velocity evaluated at each integration point. */ + Array2DView second_velocities; + + /** Consistently oriented CommonPlane unit normal for each row. */ + Array2DView normals; + + /** First face basis values evaluated at each integration point. */ + Array2DView first_basis_values; + + /** Second face basis values evaluated at each integration point. */ + Array2DView second_basis_values; + + /** Physical overlap measure multiplied by the reference quadrature weight. */ + Array1DView integration_weights; + + /** Signed normal gap evaluated from native face positions. */ + Array1DView gaps; + + /** Signed normal relative velocity evaluated from native face velocities. */ + Array1DView normal_velocity_gaps; + + /** Kinematic penalty stiffness per unit overlap measure. */ + Array1DView penalty_stiffnesses; + + /** Normal rate-penalty coefficient per unit overlap measure. */ + Array1DView rate_penalty_coefficients; + + /** Tangential viscous coefficient per unit overlap measure. */ + Array1DView tangential_viscous_coefficients; + + /** + * @brief Return the first row slot assigned to an active face pair. + * + * @param pair_id Active face-pair identifier + * @return First row slot reserved for the face pair + */ + TRIBOL_HOST_DEVICE IndexT pairRowOffset( IndexT pair_id ) const { return pair_id * maximum_rows_per_pair; } + }; + + /** + * @brief Allocate and clear storage for one CommonPlane update. + * + * @param number_of_pairs Number of active CommonPlane face pairs + * @param spatial_dimension Coupling-scheme spatial dimension + * @param allocator_id Umpire allocator identifier used by execution kernels + */ + void resize( IndexT number_of_pairs, int spatial_dimension, int allocator_id ); + + /** + * @brief Return writable non-owning views of every row field. + * + * @return Device-copyable row-batch view + */ + Viewer getView(); + + /** + * @brief Return read-only access to per-pair evaluation statuses. + * + * @return Evaluation-status array + */ + const Array1D& getPairEvaluationStatuses() const { return pair_evaluation_statuses_; } + + /** + * @brief Return read-only access to generated row counts. + * + * @return Per-pair generated-row counts + */ + const Array1D& getPairRowCounts() const { return pair_row_counts_; } + + /** + * @brief Return the allocated number of row slots. + * + * @return Number of allocated row slots + */ + IndexT getRowCapacity() const { return row_capacity_; } + + private: + /** Number of active CommonPlane face pairs represented by this batch. */ + IndexT number_of_pairs_{ 0 }; + + /** Number of allocated row slots in this batch. */ + IndexT row_capacity_{ 0 }; + + /** Spatial dimension of the coupling scheme. */ + int spatial_dimension_{ 0 }; + + /** Number of generated rows for each face pair. */ + Array1D pair_row_counts_; + + /** CommonPlanePairEvaluationStatus value for each face pair. */ + Array1D pair_evaluation_statuses_; + + /** One for generated rows and zero for unused row slots. */ + Array1D row_is_valid_; + + /** One for rows that satisfy the normal contact activation criterion. */ + Array1D row_is_active_; + + /** Stable active-pair identifier for each generated row. */ + Array1D contact_pair_ids_; + + /** First Tribol LOR face identifier for each generated row. */ + Array1D first_face_ids_; + + /** Second Tribol LOR face identifier for each generated row. */ + Array1D second_face_ids_; + + /** Number of basis values on the first field face for each generated row. */ + Array1D first_basis_counts_; + + /** Number of basis values on the second field face for each generated row. */ + Array1D second_basis_counts_; + + /** One when rows scatter directly to native parent-face response storage. */ + Array1D row_uses_parent_fields_; + + /** Physical CommonPlane integration-point coordinates. */ + Array2D integration_points_; + + /** First native parent-face reference coordinates. */ + Array2D first_parent_reference_coordinates_; + + /** Second native parent-face reference coordinates. */ + Array2D second_parent_reference_coordinates_; + + /** First face position evaluated at each integration point. */ + Array2D first_positions_; + + /** Second face position evaluated at each integration point. */ + Array2D second_positions_; + + /** First face velocity evaluated at each integration point. */ + Array2D first_velocities_; + + /** Second face velocity evaluated at each integration point. */ + Array2D second_velocities_; + + /** Consistently oriented CommonPlane unit normal for each row. */ + Array2D normals_; + + /** First face basis values evaluated at each integration point. */ + Array2D first_basis_values_; + + /** Second face basis values evaluated at each integration point. */ + Array2D second_basis_values_; + + /** Physical overlap measure multiplied by the reference quadrature weight. */ + Array1D integration_weights_; + + /** Signed normal gap evaluated from native face positions. */ + Array1D gaps_; + + /** Signed normal relative velocity evaluated from native face velocities. */ + Array1D normal_velocity_gaps_; + + /** Kinematic penalty stiffness per unit overlap measure. */ + Array1D penalty_stiffnesses_; + + /** Normal rate-penalty coefficient per unit overlap measure. */ + Array1D rate_penalty_coefficients_; + + /** Tangential viscous coefficient per unit overlap measure. */ + Array1D tangential_viscous_coefficients_; +}; + //------------------------------------------------------------------------------ class MortarData : public MethodData { public: diff --git a/src/tribol/mesh/MfemData.cpp b/src/tribol/mesh/MfemData.cpp index 80f09e8d..e8b4f2b7 100644 --- a/src/tribol/mesh/MfemData.cpp +++ b/src/tribol/mesh/MfemData.cpp @@ -9,13 +9,12 @@ #ifdef BUILD_REDECOMP +#include #include #include "axom/slic/interface/slic_macros.hpp" #include "shared/infrastructure/Profiling.hpp" -#include "tribol/common/LoopExec.hpp" - #include "redecomp/utils/ArrayUtility.hpp" #include "redecomp/transfer/TransferByElements.hpp" @@ -27,13 +26,98 @@ namespace tribol { namespace { /** First component containing a parent reference-coordinate value. */ -constexpr int parent_reference_coordinate_offset{ 2 }; +constexpr int parent_reference_coordinate_offset{ 3 }; -/** Number of scalar values transferred for one LOR-to-parent mapping record. */ -constexpr int parent_face_mapping_record_size{ parent_reference_coordinate_offset + +/** First component containing a native parent-basis coefficient. */ +constexpr int parent_basis_coefficient_offset{ parent_reference_coordinate_offset + ParentFaceData::max_lor_face_vertices * ParentFaceData::max_reference_dimension }; +/** First component containing a native parent-face nodal coordinate. */ +constexpr int parent_position_offset{ parent_basis_coefficient_offset + + ParentFaceData::max_parent_face_nodes * ParentFaceData::max_parent_face_nodes }; + +/** First component containing a native parent-face nodal velocity. */ +constexpr int parent_velocity_offset{ parent_position_offset + ParentFaceData::max_parent_face_nodes * 3 }; + +/** First component containing a native parent-face inverse diagonal mass. */ +constexpr int parent_inverse_mass_offset{ parent_velocity_offset + ParentFaceData::max_parent_face_nodes * 3 }; + +/** First component containing a parent-mesh vector degree-of-freedom identifier. */ +constexpr int parent_vector_dof_offset{ parent_inverse_mass_offset + ParentFaceData::max_parent_face_nodes * 3 }; + +/** Component containing the source-rank parent vector degree-of-freedom count. */ +constexpr int parent_vector_dof_count_offset{ parent_vector_dof_offset + ParentFaceData::max_parent_face_nodes * 3 }; + +/** Number of scalar values transferred for one LOR-to-parent mapping record. */ +constexpr int parent_face_mapping_record_size{ parent_vector_dof_count_offset + 1 }; + +/** + * @brief Decode an MFEM signed degree-of-freedom index. + * + * @param encoded_degree_of_freedom Signed MFEM degree-of-freedom index + * @param sign Orientation sign encoded in the input index + * @return Nonnegative local degree-of-freedom index + */ +int DecodeDegreeOfFreedom( int encoded_degree_of_freedom, RealT& sign ) +{ + sign = encoded_degree_of_freedom >= 0 ? 1.0 : -1.0; + return encoded_degree_of_freedom >= 0 ? encoded_degree_of_freedom : -1 - encoded_degree_of_freedom; +} + +/** + * @brief Evaluate the complete polynomial monomial basis for a parent face. + * + * Segment terms are ordered by increasing degree. Triangle terms are ordered + * by increasing total degree, then increasing second-coordinate degree. + * Quadrilateral terms are ordered with the first-coordinate degree varying + * fastest. The device evaluator in MeshData uses the same ordering. + * + * @param geometry Native parent-face geometry + * @param order Native parent-face polynomial order + * @param reference_point Native parent-face reference point + * @param monomial_values Output monomial values + */ +void EvaluateParentFaceMonomials( mfem::Geometry::Type geometry, int order, + const mfem::IntegrationPoint& reference_point, mfem::Vector& monomial_values ) +{ + int number_of_monomials = 0; + if ( geometry == mfem::Geometry::SEGMENT ) { + number_of_monomials = order + 1; + } else if ( geometry == mfem::Geometry::TRIANGLE ) { + number_of_monomials = ( order + 1 ) * ( order + 2 ) / 2; + } else if ( geometry == mfem::Geometry::SQUARE ) { + number_of_monomials = ( order + 1 ) * ( order + 1 ); + } else { + SLIC_ERROR_ROOT( "Native parent basis evaluation requires segment, triangle, or quadrilateral faces." ); + } + + monomial_values.SetSize( number_of_monomials ); + int monomial_index = 0; + if ( geometry == mfem::Geometry::SEGMENT ) { + RealT first_coordinate_power = 1.0; + for ( int first_degree = 0; first_degree <= order; ++first_degree ) { + monomial_values[monomial_index++] = first_coordinate_power; + first_coordinate_power *= reference_point.x; + } + } else if ( geometry == mfem::Geometry::TRIANGLE ) { + for ( int total_degree = 0; total_degree <= order; ++total_degree ) { + for ( int second_degree = 0; second_degree <= total_degree; ++second_degree ) { + const int first_degree = total_degree - second_degree; + monomial_values[monomial_index++] = + std::pow( reference_point.x, first_degree ) * std::pow( reference_point.y, second_degree ); + } + } + } else { + for ( int second_degree = 0; second_degree <= order; ++second_degree ) { + for ( int first_degree = 0; first_degree <= order; ++first_degree ) { + monomial_values[monomial_index++] = + std::pow( reference_point.x, first_degree ) * std::pow( reference_point.y, second_degree ); + } + } + } +} + /** Evaluate linear shape functions for a supported contact face geometry. */ void EvaluateLinearFaceShape( mfem::Geometry::Type face_geometry, const RealT* reference_coordinates, RealT* shape_values ) @@ -583,9 +667,16 @@ void ParentRedecompTransfer::RedecompToParent( const mfem::GridFunction& redecom submesh_gridfn_ = 0.0; submesh_redecomp_xfer_.RedecompToSubmesh( redecomp_src, submesh_gridfn_ ); + AddSubmeshToParent( submesh_gridfn_, parent_dst ); +} + +void ParentRedecompTransfer::AddSubmeshToParent( const mfem::Vector& submesh_src, mfem::Vector& parent_dst ) const +{ + const bool use_device = parent_dst.UseDevice(); + // Response is a dual field; scatter local submesh contributions directly instead of using ParSubMesh::Transfer, // which communicates/averages shared DOFs as a primal grid function. - const auto submesh_data = submesh_gridfn_.Read( use_device ); + const auto submesh_data = submesh_src.Read( use_device ); const auto submesh_to_parent_vdof_map = submesh_to_parent_vdof_map_.Read( use_device ); auto parent_data = parent_dst.ReadWrite( use_device ); mfem::forall_switch( use_device, submesh_to_parent_vdof_map_.Size(), [=] MFEM_HOST_DEVICE( int i ) { @@ -819,8 +910,8 @@ bool MfemMeshData::UpdateMfemMeshData( RealT binning_proximity_scale, int n_rank } update_data_ = std::make_unique( submesh_, lor_mesh_.get(), *coords_.GetParentGridFn().ParFESpace(), submesh_xfer_gridfn_, submesh_lor_xfer_.get(), attributes_1_, - attributes_2_, build_parent_face_data_, binning_proximity_scale, - n_ranks, allocator_id_, redecomp_trigger_displacement_, residual_gap ); + attributes_2_, binning_proximity_scale, n_ranks, allocator_id_, + redecomp_trigger_displacement_, residual_gap ); rebuilt = true; } @@ -848,6 +939,11 @@ bool MfemMeshData::UpdateMfemMeshData( RealT binning_proximity_scale, int n_rank if ( velocity_ ) { velocity_->UpdateField( update_data_->vector_xfer_, use_device_ ); } + if ( build_parent_face_data_ ) { + update_data_->BuildParentFaceData( submesh_, lor_mesh_.get(), coords_.GetParentGridFn(), + velocity_ ? &velocity_->GetParentGridFn() : nullptr, + inverse_mass_ ? &inverse_mass_->GetParentGridFn() : nullptr ); + } TRIBOL_MARK_END( "Copy fields to Redecomp mesh" ); if ( rebuilt && elem_thickness_ ) { @@ -923,6 +1019,12 @@ ParentFaceData MfemMeshData::GetMesh2ParentFaceData() const { return GetUpdateDa void MfemMeshData::GetParentResponse( mfem::Vector& r ) const { GetParentRedecompTransfer().RedecompToParent( *redecomp_response_, r ); + if ( build_parent_face_data_ ) { + mfem::Vector submesh_parent_face_response; + GetUpdateData().GetParentFaceResponse( const_cast( submesh_ ), lor_mesh_.get(), + submesh_parent_face_response ); + GetParentRedecompTransfer().AddSubmeshToParent( submesh_parent_face_response, r ); + } } void MfemMeshData::SetParentVelocity( const mfem::ParGridFunction& velocity ) @@ -934,6 +1036,23 @@ void MfemMeshData::SetParentVelocity( const mfem::ParGridFunction& velocity ) } } +void MfemMeshData::SetParentInverseMass( const mfem::ParGridFunction& inverse_mass ) +{ + const mfem::ParFiniteElementSpace& coordinate_space = *coords_.GetParentGridFn().ParFESpace(); + const mfem::ParFiniteElementSpace& inverse_mass_space = *inverse_mass.ParFESpace(); + SLIC_ERROR_ROOT_IF( inverse_mass_space.GetParMesh() != coordinate_space.GetParMesh() || + inverse_mass_space.GetVSize() != coordinate_space.GetVSize() || + inverse_mass_space.GetVDim() != coordinate_space.GetVDim() || + inverse_mass_space.GetOrdering() != coordinate_space.GetOrdering(), + "Parent inverse diagonal mass must use the parent coordinate vector finite-element space." ); + + if ( inverse_mass_ ) { + inverse_mass_->SetParentGridFn( inverse_mass ); + } else { + inverse_mass_ = std::make_unique( inverse_mass ); + } +} + void MfemMeshData::ClearAllPenaltyData() { ClearRatePenaltyData(); @@ -1095,8 +1214,8 @@ MfemMeshData::UpdateData::UpdateData( mfem::ParSubMesh& submesh, mfem::ParMesh* const mfem::ParFiniteElementSpace& parent_fes, mfem::ParGridFunction& submesh_gridfn, SubmeshLORTransfer* submesh_lor_xfer, const std::set& attributes_1, const std::set& attributes_2, - bool build_parent_face_data, RealT binning_proximity_scale, int n_ranks, - int allocator_id, RealT redecomp_trigger_displacement, RealT residual_gap ) + RealT binning_proximity_scale, int n_ranks, int allocator_id, + RealT redecomp_trigger_displacement, RealT residual_gap ) : redecomp_mesh_{ lor_mesh ? redecomp::RedecompMesh( *lor_mesh, @@ -1116,15 +1235,9 @@ MfemMeshData::UpdateData::UpdateData( mfem::ParSubMesh& submesh, mfem::ParMesh* TRIBOL_MARK_FUNCTION; // set element type based on redecomp mesh SetElementData(); - // Keep the element maps in host memory until all host-side MFEM data have - // been gathered for the Tribol surface meshes. + // Build connectivity from host-only MFEM data before copying it to the + // allocator used by Tribol kernels. UpdateConnectivity( attributes_1, attributes_2 ); - if ( build_parent_face_data ) { - // Preserve the native parent-face mapping alongside the redecomposed LOR - // geometry. Later integration-point evaluation uses this mapping without - // transferring physical fields through the LOR finite-element space. - BuildParentFaceData( submesh, lor_mesh, parent_fes ); - } CopyConnectivityToAllocator(); } @@ -1132,17 +1245,31 @@ ParentFaceData MfemMeshData::UpdateData::ParentFaceArrays::GetView() const { ParentFaceData parent_face_data; parent_face_data.m_parent_face_orders = Array1DView( parent_face_orders ); + parent_face_data.m_parent_node_counts = Array1DView( parent_node_counts ); parent_face_data.m_reference_vertex_counts = Array1DView( reference_vertex_counts ); parent_face_data.m_parent_reference_vertex_coordinates = Array2DView( parent_reference_vertex_coordinates ); + parent_face_data.m_parent_basis_coefficients = Array2DView( parent_basis_coefficients ); + parent_face_data.m_parent_positions = Array2DView( parent_positions ); + if ( !parent_velocities.empty() ) { + parent_face_data.m_parent_velocities = Array2DView( parent_velocities ); + } + if ( !parent_inverse_masses.empty() ) { + parent_face_data.m_parent_inverse_masses = Array2DView( parent_inverse_masses ); + parent_face_data.m_parent_vector_dof_ids = Array2DView( parent_vector_dof_ids ); + parent_face_data.m_parent_vector_dof_count = parent_vector_dof_count; + } + parent_face_data.m_parent_responses = Array2DView( const_cast&>( parent_responses ) ); return parent_face_data; } void MfemMeshData::UpdateData::BuildParentFaceData( mfem::ParSubMesh& submesh, mfem::ParMesh* lor_mesh, - const mfem::ParFiniteElementSpace& parent_fes ) + const mfem::ParGridFunction& parent_coordinates, + const mfem::ParGridFunction* parent_velocity, + const mfem::ParGridFunction* parent_inverse_mass ) { mfem::ParMesh& source_mesh = lor_mesh ? *lor_mesh : submesh; - mfem::ParMesh* parent_mesh = parent_fes.GetParMesh(); + const mfem::ParMesh* parent_mesh = submesh.GetParent(); SLIC_ERROR_ROOT_IF( parent_mesh == nullptr, "Parent-face mapping requires an MFEM parallel parent mesh." ); const int reference_dimension = source_mesh.Dimension(); SLIC_ERROR_ROOT_IF( reference_dimension < 1 || reference_dimension > ParentFaceData::max_reference_dimension, @@ -1163,6 +1290,32 @@ void MfemMeshData::UpdateData::BuildParentFaceData( mfem::ParSubMesh& submesh, m source_records = 0.0; redecomp_records = 0.0; + // Parent fields are first restricted to the parent-linked boundary submesh. + // Each LOR element record then carries the complete native parent-face nodal + // data to the rank that owns the corresponding redecomposed contact face. + auto* submesh_coordinates = dynamic_cast( submesh.GetNodes() ); + SLIC_ERROR_ROOT_IF( submesh_coordinates == nullptr, + "Parent-face field evaluation requires a ParGridFunction on the boundary submesh." ); + submesh.Transfer( parent_coordinates, *submesh_coordinates ); + const mfem::ParFiniteElementSpace& submesh_space = *submesh_coordinates->ParFESpace(); + const RealT* submesh_coordinate_values = submesh_coordinates->HostRead(); + + mfem::ParGridFunction submesh_velocity( const_cast( &submesh_space ) ); + const RealT* submesh_velocity_values = nullptr; + if ( parent_velocity != nullptr ) { + submesh_velocity = 0.0; + submesh.Transfer( *parent_velocity, submesh_velocity ); + submesh_velocity_values = submesh_velocity.HostRead(); + } + + mfem::ParGridFunction submesh_inverse_mass( const_cast( &submesh_space ) ); + const RealT* submesh_inverse_mass_values = nullptr; + if ( parent_inverse_mass != nullptr ) { + submesh_inverse_mass = 0.0; + submesh.Transfer( *parent_inverse_mass, submesh_inverse_mass ); + submesh_inverse_mass_values = submesh_inverse_mass.HostRead(); + } + const mfem::CoarseFineTransformations* refinement_transforms = lor_mesh ? &lor_mesh->GetRefinementTransforms() : nullptr; mfem::Array element_dofs; @@ -1181,13 +1334,22 @@ void MfemMeshData::UpdateData::BuildParentFaceData( mfem::ParSubMesh& submesh, m } const IndexT parent_face_id = submesh.GetParentElementIDMap()[parent_submesh_element_id]; - const int parent_face_order = parent_fes.GetBE( parent_face_id )->GetOrder(); + const mfem::FiniteElement& parent_face_element = *submesh_space.GetFE( parent_submesh_element_id ); + const int parent_face_order = parent_face_element.GetOrder(); + const mfem::Geometry::Type parent_face_geometry = parent_face_element.GetGeomType(); + const int number_of_parent_nodes = parent_face_element.GetDof(); const int number_of_reference_vertices = reference_vertices->GetNPoints(); + SLIC_ERROR_ROOT_IF( parent_face_order < 1 || parent_face_order > ParentFaceData::max_parent_face_order, + "Native parent-face evaluation supports polynomial orders one through four." ); + SLIC_ERROR_ROOT_IF( number_of_parent_nodes > ParentFaceData::max_parent_face_nodes, + "Native parent-face evaluation supports at most 25 nodes per face." ); SLIC_ERROR_ROOT_IF( number_of_reference_vertices > ParentFaceData::max_lor_face_vertices, "Parent-face mapping supports at most four vertices per LOR surface face." ); RealT record[parent_face_mapping_record_size] = { 0.0 }; record[0] = static_cast( parent_face_order ); record[1] = static_cast( number_of_reference_vertices ); + record[2] = static_cast( number_of_parent_nodes ); + record[parent_vector_dof_count_offset] = static_cast( parent_coordinates.Size() ); // ParSubMesh elements can use a different reference orientation than their // native parent boundary faces. Match vertices through the parent mesh so @@ -1237,6 +1399,64 @@ void MfemMeshData::UpdateData::BuildParentFaceData( mfem::ParSubMesh& submesh, m } } + // Build a polynomial representation of every nodal basis function from + // the native finite element's own reference nodes. This preserves MFEM's + // local node ordering while keeping evaluation device-callable. + const mfem::IntegrationRule& parent_reference_nodes = parent_face_element.GetNodes(); + SLIC_ERROR_ROOT_IF( parent_reference_nodes.GetNPoints() != number_of_parent_nodes, + "Native parent-face basis construction requires one reference node per degree of freedom." ); + mfem::DenseMatrix parent_vandermonde( number_of_parent_nodes ); + mfem::Vector parent_monomials; + for ( int parent_node = 0; parent_node < number_of_parent_nodes; ++parent_node ) { + EvaluateParentFaceMonomials( parent_face_geometry, parent_face_order, + parent_reference_nodes.IntPoint( parent_node ), parent_monomials ); + SLIC_ERROR_ROOT_IF( parent_monomials.Size() != number_of_parent_nodes, + "Native parent-face polynomial basis size does not match its number of nodes." ); + for ( int monomial = 0; monomial < number_of_parent_nodes; ++monomial ) { + parent_vandermonde( parent_node, monomial ) = parent_monomials[monomial]; + } + } + mfem::DenseMatrixInverse parent_vandermonde_inverse( parent_vandermonde ); + mfem::Vector nodal_value( number_of_parent_nodes ); + mfem::Vector basis_coefficients( number_of_parent_nodes ); + for ( int parent_node = 0; parent_node < number_of_parent_nodes; ++parent_node ) { + nodal_value = 0.0; + nodal_value[parent_node] = 1.0; + parent_vandermonde_inverse.Mult( nodal_value, basis_coefficients ); + for ( int monomial = 0; monomial < number_of_parent_nodes; ++monomial ) { + record[parent_basis_coefficient_offset + parent_node * ParentFaceData::max_parent_face_nodes + monomial] = + basis_coefficients[monomial]; + } + } + + mfem::Array parent_element_dofs; + submesh_space.GetElementDofs( parent_submesh_element_id, parent_element_dofs ); + SLIC_ERROR_ROOT_IF( parent_element_dofs.Size() != number_of_parent_nodes, + "Native parent-face field data do not match the parent finite-element basis." ); + for ( int parent_node = 0; parent_node < number_of_parent_nodes; ++parent_node ) { + RealT degree_of_freedom_sign = 1.0; + const int parent_degree_of_freedom = + DecodeDegreeOfFreedom( parent_element_dofs[parent_node], degree_of_freedom_sign ); + for ( int component = 0; component < parent_coordinates.VectorDim(); ++component ) { + const int submesh_vector_degree_of_freedom = submesh_space.DofToVDof( parent_degree_of_freedom, component ); + const int field_index = parent_node * parent_coordinates.VectorDim() + component; + record[parent_position_offset + field_index] = + degree_of_freedom_sign * submesh_coordinate_values[submesh_vector_degree_of_freedom]; + if ( submesh_velocity_values != nullptr ) { + record[parent_velocity_offset + field_index] = + degree_of_freedom_sign * submesh_velocity_values[submesh_vector_degree_of_freedom]; + } + if ( submesh_inverse_mass_values != nullptr ) { + record[parent_inverse_mass_offset + field_index] = + submesh_inverse_mass_values[submesh_vector_degree_of_freedom]; + RealT parent_degree_of_freedom_sign = 1.0; + const int parent_vector_degree_of_freedom = DecodeDegreeOfFreedom( + vector_xfer_.GetParentVDof( submesh_vector_degree_of_freedom ), parent_degree_of_freedom_sign ); + record[parent_vector_dof_offset + field_index] = static_cast( parent_vector_degree_of_freedom ); + } + } + } + source_space.GetElementDofs( source_element_id, element_dofs ); SLIC_ERROR_ROOT_IF( element_dofs.Size() != 1, "Parent-face mapping requires one scalar degree of freedom per source element." ); @@ -1248,14 +1468,37 @@ void MfemMeshData::UpdateData::BuildParentFaceData( mfem::ParSubMesh& submesh, m redecomp::TransferByElements().TransferToSerial( source_records, redecomp_records ); - auto copy_surface_records = [&]( const Array1D& surface_element_map, ParentFaceArrays& parent_face_arrays ) { + auto copy_surface_records = [&]( const Array1D& surface_element_map, + ParentFaceArrays& parent_face_arrays ) { const IndexT number_of_surface_elements = surface_element_map.size(); ArrayT parent_face_orders_host( number_of_surface_elements ); + ArrayT parent_node_counts_host( number_of_surface_elements ); ArrayT reference_vertex_counts_host( number_of_surface_elements ); ArrayT parent_reference_coordinates_host( { number_of_surface_elements, ParentFaceData::max_lor_face_vertices * ParentFaceData::max_reference_dimension }, getResourceAllocatorID( MemorySpace::Host ) ); parent_reference_coordinates_host.fill( 0.0 ); + ArrayT parent_basis_coefficients_host( + { number_of_surface_elements, ParentFaceData::max_parent_face_nodes * ParentFaceData::max_parent_face_nodes }, + getResourceAllocatorID( MemorySpace::Host ) ); + ArrayT parent_positions_host( + { number_of_surface_elements, ParentFaceData::max_parent_face_nodes * 3 }, + getResourceAllocatorID( MemorySpace::Host ) ); + ArrayT parent_velocities_host( + { number_of_surface_elements, ParentFaceData::max_parent_face_nodes * 3 }, + getResourceAllocatorID( MemorySpace::Host ) ); + ArrayT parent_inverse_masses_host( + { number_of_surface_elements, ParentFaceData::max_parent_face_nodes * 3 }, + getResourceAllocatorID( MemorySpace::Host ) ); + ArrayT parent_vector_dof_ids_host( + { number_of_surface_elements, ParentFaceData::max_parent_face_nodes * 3 }, + getResourceAllocatorID( MemorySpace::Host ) ); + parent_basis_coefficients_host.fill( 0.0 ); + parent_positions_host.fill( 0.0 ); + parent_velocities_host.fill( 0.0 ); + parent_inverse_masses_host.fill( 0.0 ); + parent_vector_dof_ids_host.fill( -1 ); + IndexT maximum_parent_vector_dof_count = 0; for ( IndexT surface_element_id = 0; surface_element_id < number_of_surface_elements; ++surface_element_id ) { const int redecomp_element_id = surface_element_map[surface_element_id]; @@ -1270,22 +1513,164 @@ void MfemMeshData::UpdateData::BuildParentFaceData( mfem::ParSubMesh& submesh, m parent_face_orders_host[surface_element_id] = static_cast( get_record_value( 0 ) ); reference_vertex_counts_host[surface_element_id] = static_cast( get_record_value( 1 ) ); + parent_node_counts_host[surface_element_id] = static_cast( get_record_value( 2 ) ); for ( int coordinate_index = 0; coordinate_index < ParentFaceData::max_lor_face_vertices * ParentFaceData::max_reference_dimension; ++coordinate_index ) { parent_reference_coordinates_host( surface_element_id, coordinate_index ) = get_record_value( parent_reference_coordinate_offset + coordinate_index ); } + for ( int coefficient_index = 0; + coefficient_index < ParentFaceData::max_parent_face_nodes * ParentFaceData::max_parent_face_nodes; + ++coefficient_index ) { + parent_basis_coefficients_host( surface_element_id, coefficient_index ) = + get_record_value( parent_basis_coefficient_offset + coefficient_index ); + } + for ( int field_index = 0; field_index < ParentFaceData::max_parent_face_nodes * 3; ++field_index ) { + parent_positions_host( surface_element_id, field_index ) = + get_record_value( parent_position_offset + field_index ); + parent_velocities_host( surface_element_id, field_index ) = + get_record_value( parent_velocity_offset + field_index ); + parent_inverse_masses_host( surface_element_id, field_index ) = + get_record_value( parent_inverse_mass_offset + field_index ); + parent_vector_dof_ids_host( surface_element_id, field_index ) = + static_cast( get_record_value( parent_vector_dof_offset + field_index ) ); + } + const IndexT parent_vector_dof_count = static_cast( get_record_value( parent_vector_dof_count_offset ) ); + maximum_parent_vector_dof_count = std::max( maximum_parent_vector_dof_count, parent_vector_dof_count ); } parent_face_arrays.parent_face_orders = Array1D( parent_face_orders_host, allocator_id_ ); + parent_face_arrays.parent_node_counts = Array1D( parent_node_counts_host, allocator_id_ ); parent_face_arrays.reference_vertex_counts = Array1D( reference_vertex_counts_host, allocator_id_ ); parent_face_arrays.parent_reference_vertex_coordinates = Array2D( parent_reference_coordinates_host, allocator_id_ ); + parent_face_arrays.parent_basis_coefficients = Array2D( parent_basis_coefficients_host, allocator_id_ ); + parent_face_arrays.parent_positions = Array2D( parent_positions_host, allocator_id_ ); + if ( parent_velocity != nullptr ) { + parent_face_arrays.parent_velocities = Array2D( parent_velocities_host, allocator_id_ ); + } else { + parent_face_arrays.parent_velocities = Array2D(); + } + if ( parent_inverse_mass != nullptr ) { + parent_face_arrays.parent_inverse_masses = Array2D( parent_inverse_masses_host, allocator_id_ ); + parent_face_arrays.parent_vector_dof_ids = Array2D( parent_vector_dof_ids_host, allocator_id_ ); + parent_face_arrays.parent_vector_dof_count = maximum_parent_vector_dof_count; + } else { + parent_face_arrays.parent_inverse_masses = Array2D(); + parent_face_arrays.parent_vector_dof_ids = Array2D(); + parent_face_arrays.parent_vector_dof_count = 0; + } + parent_face_arrays.parent_responses = + Array2D( { number_of_surface_elements, ParentFaceData::max_parent_face_nodes * 3 }, allocator_id_ ); + parent_face_arrays.parent_responses.fill( 0.0 ); + }; + + // Connectivity is stored in the requested execution memory space, while + // MFEM element access remains host-only. Mirror only the two element maps + // needed to unpack the transferred records. + const Array1D first_surface_element_map( elem_map_1_ ); + const Array1D second_surface_element_map( elem_map_2_ ); + copy_surface_records( first_surface_element_map, parent_face_data_1_ ); + copy_surface_records( second_surface_element_map, parent_face_data_2_ ); +} + +void MfemMeshData::UpdateData::GetParentFaceResponse( mfem::ParSubMesh& submesh, mfem::ParMesh* lor_mesh, + mfem::Vector& submesh_response ) const +{ + mfem::ParMesh& source_mesh = lor_mesh ? *lor_mesh : submesh; + const int response_components = ParentFaceData::max_parent_face_nodes * 3; + mfem::L2_FECollection response_collection( 0, source_mesh.Dimension() ); + mfem::FiniteElementSpace redecomp_response_space( const_cast( &redecomp_mesh_ ), + &response_collection, response_components, + mfem::Ordering::byNODES ); + mfem::ParFiniteElementSpace source_response_space( &source_mesh, &response_collection, response_components, + mfem::Ordering::byNODES ); + mfem::GridFunction redecomp_response( &redecomp_response_space ); + mfem::ParGridFunction source_response( &source_response_space ); + redecomp_response.UseDevice( false ); + source_response.UseDevice( false ); + redecomp_response = 0.0; + source_response = 0.0; + + // Contact kernels accumulate one parent-basis force vector per Tribol face. + // Place those vectors back on their redecomp elements before the element-wise + // reverse transfer returns each owned contribution to its source LOR face. + RealT* redecomp_response_values = redecomp_response.HostReadWrite(); + auto copy_surface_response = [&]( const Array1D& surface_element_map, + const ParentFaceArrays& parent_face_arrays ) { + ArrayT parent_response_host( parent_face_arrays.parent_responses ); + mfem::Array response_dofs; + for ( IndexT surface_element_id = 0; surface_element_id < surface_element_map.size(); ++surface_element_id ) { + const int redecomp_element_id = surface_element_map[surface_element_id]; + redecomp_response_space.GetElementDofs( redecomp_element_id, response_dofs ); + SLIC_ERROR_ROOT_IF( response_dofs.Size() != 1, + "Parent-face response transfer requires one scalar degree of freedom per element." ); + for ( int component = 0; component < response_components; ++component ) { + const int response_vector_degree_of_freedom = redecomp_response_space.DofToVDof( response_dofs[0], component ); + redecomp_response_values[response_vector_degree_of_freedom] += + parent_response_host( surface_element_id, component ); + } + } }; + // MFEM element degree-of-freedom lookup is host-only, so mirror the element + // maps before placing device-accumulated parent responses on redecomp faces. + const Array1D first_surface_element_map( elem_map_1_ ); + const Array1D second_surface_element_map( elem_map_2_ ); + copy_surface_response( first_surface_element_map, parent_face_data_1_ ); + copy_surface_response( second_surface_element_map, parent_face_data_2_ ); + + redecomp::TransferByElements().TransferToParallel( redecomp_response, source_response ); + + const mfem::ParFiniteElementSpace& submesh_space = + *dynamic_cast( submesh.GetNodes() )->ParFESpace(); + submesh_response.SetSize( submesh_space.GetVSize() ); + submesh_response.UseDevice( false ); + submesh_response = 0.0; + RealT* submesh_response_values = submesh_response.HostReadWrite(); + const RealT* source_response_values = source_response.HostRead(); + const mfem::CoarseFineTransformations* refinement_transforms = + lor_mesh ? &lor_mesh->GetRefinementTransforms() : nullptr; + mfem::Array source_response_dofs; + mfem::Array parent_element_dofs; + for ( int source_element_id = 0; source_element_id < source_mesh.GetNE(); ++source_element_id ) { + const int parent_submesh_element_id = + refinement_transforms ? refinement_transforms->embeddings[source_element_id].parent : source_element_id; + source_response_space.GetElementDofs( source_element_id, source_response_dofs ); + submesh_space.GetElementDofs( parent_submesh_element_id, parent_element_dofs ); + SLIC_ERROR_ROOT_IF( + source_response_dofs.Size() != 1 || parent_element_dofs.Size() > ParentFaceData::max_parent_face_nodes, + "Parent-face response transfer encountered an unsupported finite-element layout." ); + + for ( int parent_node = 0; parent_node < parent_element_dofs.Size(); ++parent_node ) { + RealT degree_of_freedom_sign = 1.0; + const int parent_degree_of_freedom = + DecodeDegreeOfFreedom( parent_element_dofs[parent_node], degree_of_freedom_sign ); + for ( int component = 0; component < submesh_space.GetVDim(); ++component ) { + const int source_component = parent_node * submesh_space.GetVDim() + component; + const int source_vector_degree_of_freedom = + source_response_space.DofToVDof( source_response_dofs[0], source_component ); + const int submesh_vector_degree_of_freedom = submesh_space.DofToVDof( parent_degree_of_freedom, component ); + submesh_response_values[submesh_vector_degree_of_freedom] += + degree_of_freedom_sign * source_response_values[source_vector_degree_of_freedom]; + } + } + } - copy_surface_records( elem_map_1_, parent_face_data_1_ ); - copy_surface_records( elem_map_2_, parent_face_data_2_ ); + // Sum shared native degrees of freedom and retain values only on their owner + // ranks, matching the dual-vector convention used by existing MFEM response transfer. + mfem::Vector local_response_copy( submesh_response ); + mfem::ParGridFunction local_response_grid_function( const_cast( &submesh_space ), + local_response_copy ); + std::unique_ptr true_response( local_response_grid_function.ParallelAssemble() ); + submesh_space.Dof_TrueDof_Matrix()->Mult( *true_response, submesh_response ); + RealT* owned_response_values = submesh_response.HostReadWrite(); + for ( int local_degree_of_freedom = 0; local_degree_of_freedom < submesh_space.GetVSize(); + ++local_degree_of_freedom ) { + if ( submesh_space.GetLocalTDofNumber( local_degree_of_freedom ) < 0 ) { + owned_response_values[local_degree_of_freedom] = 0.0; + } + } } void MfemMeshData::UpdateData::UpdateConnectivity( const std::set& attributes_1, diff --git a/src/tribol/mesh/MfemData.hpp b/src/tribol/mesh/MfemData.hpp index 7a5300b6..8332f4a8 100644 --- a/src/tribol/mesh/MfemData.hpp +++ b/src/tribol/mesh/MfemData.hpp @@ -295,6 +295,17 @@ class ParentRedecompTransfer { */ void RedecompToParent( const mfem::GridFunction& redecomp_src, mfem::Vector& parent_dst ) const; + /** + * @brief Add a boundary-submesh dual vector to its parent-mesh vector. + * + * Shared submesh degrees of freedom must already follow MFEM's owned-value + * convention before this local signed scatter. + * + * @param [in] submesh_src Boundary-submesh dual vector + * @param [in,out] parent_dst Parent-mesh vector receiving the contribution + */ + void AddSubmeshToParent( const mfem::Vector& submesh_src, mfem::Vector& parent_dst ) const; + /** * @brief Get the parent-linked boundary submesh finite element space * associated with this transfer object @@ -303,6 +314,17 @@ class ParentRedecompTransfer { */ const mfem::ParFiniteElementSpace& GetSubmeshFESpace() const { return *submesh_gridfn_.ParFESpace(); } + /** + * @brief Return the parent-mesh vector degree of freedom for a submesh vector degree of freedom. + * + * The returned MFEM index may encode an orientation sign. Callers that only + * need a degree-of-freedom identity must decode the sign before using it. + * + * @param submesh_vector_dof Boundary-submesh vector degree of freedom + * @return Corresponding signed parent-mesh vector degree of freedom + */ + int GetParentVDof( int submesh_vector_dof ) const { return submesh_to_parent_vdof_map_[submesh_vector_dof]; } + /** * @brief Returns finite element space on the redecomp mesh associated with * this transfer object @@ -781,6 +803,20 @@ class MfemMeshData { */ bool HasVelocity() const { return velocity_ != nullptr; } + /** + * @brief Add or replace the parent inverse diagonal mass field. + * + * @param inverse_mass Component-wise inverse diagonal mass in the parent velocity space + */ + void SetParentInverseMass( const mfem::ParGridFunction& inverse_mass ); + + /** + * @brief Determine whether an inverse diagonal mass field is registered. + * + * @return true when the parent inverse diagonal mass field is available + */ + bool HasInverseMass() const { return inverse_mass_ != nullptr; } + /** * @brief Get pointers to component arrays of the velocity on the RedecompMesh * @@ -1154,12 +1190,36 @@ class MfemMeshData { /** Polynomial order of each native parent coordinate face. */ Array1D parent_face_orders; + /** Number of native parent finite-element nodes on each face. */ + Array1D parent_node_counts; + /** Number of vertices defining each child-to-parent reference map. */ Array1D reference_vertex_counts; /** Parent reference coordinates at the LOR face vertices. */ Array2D parent_reference_vertex_coordinates; + /** Coefficients that evaluate the native parent nodal basis. */ + Array2D parent_basis_coefficients; + + /** Native parent-face nodal coordinates in node-major ordering. */ + Array2D parent_positions; + + /** Native parent-face nodal velocities in node-major ordering. */ + Array2D parent_velocities; + + /** Native parent-face inverse diagonal masses in node-major ordering. */ + Array2D parent_inverse_masses; + + /** Parent-mesh vector degree-of-freedom identifiers in node-major ordering. */ + Array2D parent_vector_dof_ids; + + /** Largest source-rank vector degree-of-freedom count represented by these faces. */ + IndexT parent_vector_dof_count{ 0 }; + + /** Element-local native parent-face response accumulated by contact. */ + Array2D parent_responses; + /** * @brief Create non-owning views of the mapping arrays. * @@ -1179,7 +1239,6 @@ class MfemMeshData { * @param submesh_lor_xfer Submesh to LOR grid function transfer object (if using LOR; nullptr otherwise) * @param attributes_1 Set of boundary attributes identifying elements in the first Tribol registered mesh * @param attributes_2 Set of boundary attributes identifying elements in the second Tribol registered mesh - * @param build_parent_face_data Whether native parent-face data are needed by the contact method * @param binning_proximity_scale Element length multiplier for coarse binning and proximity detection inclusion. * This is needed to size the ghost element layer in the redecomp mesh. * @param n_ranks Number of ranks in the parallel decomposition @@ -1190,9 +1249,8 @@ class MfemMeshData { */ UpdateData( mfem::ParSubMesh& submesh, mfem::ParMesh* lor_mesh, const mfem::ParFiniteElementSpace& parent_fes, mfem::ParGridFunction& submesh_gridfn, SubmeshLORTransfer* submesh_lor_xfer, - const std::set& attributes_1, const std::set& attributes_2, bool build_parent_face_data, - RealT binning_proximity_scale, int n_ranks, int allocator_id, RealT redecomp_trigger_displacement, - RealT residual_gap ); + const std::set& attributes_1, const std::set& attributes_2, RealT binning_proximity_scale, + int n_ranks, int allocator_id, RealT redecomp_trigger_displacement, RealT residual_gap ); /** * @brief Redecomposed boundary element mesh @@ -1250,6 +1308,8 @@ class MfemMeshData { int allocator_id_; private: + friend class MfemMeshData; + /** * @brief Builds connectivity arrays and redecomp mesh to Tribol registered * mesh element maps @@ -1269,10 +1329,24 @@ class MfemMeshData { * * @param submesh Parent-linked contact boundary submesh * @param lor_mesh Optional low-order-refined contact mesh - * @param parent_fes Native parent coordinate finite-element space + * @param parent_coordinates Native parent coordinate field + * @param parent_velocity Optional native parent velocity field + * @param parent_inverse_mass Optional native parent inverse diagonal mass field */ void BuildParentFaceData( mfem::ParSubMesh& submesh, mfem::ParMesh* lor_mesh, - const mfem::ParFiniteElementSpace& parent_fes ); + const mfem::ParGridFunction& parent_coordinates, + const mfem::ParGridFunction* parent_velocity, + const mfem::ParGridFunction* parent_inverse_mass ); + + /** + * @brief Accumulate element-local parent-face response onto the boundary submesh. + * + * @param submesh Parent-linked contact boundary submesh + * @param lor_mesh Optional low-order-refined contact mesh + * @param submesh_response Boundary-submesh response receiving contact contributions + */ + void GetParentFaceResponse( mfem::ParSubMesh& submesh, mfem::ParMesh* lor_mesh, + mfem::Vector& submesh_response ) const; /** * @brief Sets the number of vertices per element and the element type for the redecomp mesh @@ -1378,6 +1452,11 @@ class MfemMeshData { */ std::unique_ptr velocity_; + /** + * @brief Contains inverse diagonal mass data if registered; nullptr otherwise + */ + std::unique_ptr inverse_mass_; + /** * @brief Kinematic constant contact penalty for the first Tribol registered mesh */ diff --git a/src/tribol/physics/CommonPlane.cpp b/src/tribol/physics/CommonPlane.cpp index 66ec2190..3e2e8fe7 100644 --- a/src/tribol/physics/CommonPlane.cpp +++ b/src/tribol/physics/CommonPlane.cpp @@ -15,550 +15,1184 @@ #include "tribol/integ/FE.hpp" #include "tribol/utils/Math.hpp" +#include + namespace tribol { -TRIBOL_HOST_DEVICE RealT ComputeGapRatePressure( CommonPlanePair& plane, const MeshData::Viewer& m1, - const MeshData::Viewer& m2, RealT element_penalty, - RatePenaltyCalculation rate_calc ) +namespace { + +constexpr int max_dim = 3; +constexpr int max_nodes_per_face = ParentFaceData::max_parent_face_nodes; + +/** + * @brief Compute the equivalent normal rate-penalty coefficient for a face pair. + * + * @param first_mesh First contact surface mesh + * @param second_mesh Second contact surface mesh + * @param penalty_stiffness Equivalent kinematic penalty stiffness + * @param rate_calculation Registered rate-penalty calculation + * @return Equivalent normal rate coefficient per unit overlap measure + */ +TRIBOL_HOST_DEVICE inline RealT ComputeRatePenalty( const MeshData::Viewer& first_mesh, + const MeshData::Viewer& second_mesh, RealT penalty_stiffness, + RatePenaltyCalculation rate_calculation ) { - auto fId1 = plane.getCpElementId1(); - auto fId2 = plane.getCpElementId2(); - - const auto dim = plane.m_dim; - - // compute the correct rate_penalty - RealT rate_penalty = 0.; - switch ( rate_calc ) { + switch ( rate_calculation ) { 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 * ( first_mesh.getElementData().m_rate_penalty_stiffness + + second_mesh.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 penalty_stiffness * 0.5 * + ( first_mesh.getElementData().m_rate_percent_stiffness + + second_mesh.getElementData().m_rate_percent_stiffness ); } default: - // no-op, quiet compiler - break; - } // end switch on rate_calc - - // compute the velocity gap and pressure contribution - constexpr int max_dim = 3; - constexpr int max_nodes_per_elem = 4; - - 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; - auto numNodesPerFace2 = m2.numberOfNodesPerElement(); - plane.getFace2Coords( x2, numNodesPerFace2 ); // get avg face coords off the contact plane - m2.getFaceVelocities( fId2, v2 ); - - ////////////////////////////////////////////////////////// - // compute velocity Galerkin approximation at projected // - // overlap centroid // - ////////////////////////////////////////////////////////// - RealT vel_f1[max_dim]; - RealT vel_f2[max_dim]; - initRealArray( vel_f1, dim, 0. ); - initRealArray( vel_f2, dim, 0. ); - - // interpolate nodal velocity at overlap centroid as projected - // onto face 1 - RealT cXf1 = plane.m_cXf1; - RealT cYf1 = plane.m_cYf1; - RealT cZf1 = ( dim == 3 ) ? plane.m_cZf1 : 0.; - GalerkinEval( x1, cXf1, cYf1, cZf1, LINEAR, PHYSICAL, dim, dim, v1, vel_f1 ); - - // interpolate nodal velocity at overlap centroid as projected - // onto face 2 - RealT cXf2 = plane.m_cXf2; - RealT cYf2 = plane.m_cYf2; - RealT cZf2 = ( dim == 3 ) ? plane.m_cZf2 : 0.; - GalerkinEval( x2, cXf2, cYf2, cZf2, LINEAR, PHYSICAL, dim, dim, v2, vel_f2 ); - - // compute velocity gap vector - RealT velGap[max_dim]; - velGap[0] = vel_f1[0] - vel_f2[0]; - velGap[1] = vel_f1[1] - vel_f2[1]; - if ( dim == 3 ) { - velGap[2] = vel_f1[2] - vel_f2[2]; + return 0.; } +} - // compute velocity gap scalar - plane.m_velGap = 0.; - plane.m_velGap += velGap[0] * plane.m_nX; - plane.m_velGap += velGap[1] * plane.m_nY; - if ( dim == 3 ) { - plane.m_velGap += velGap[2] * plane.m_nZ; +/** + * @brief Select multipoint integration for one CommonPlane face pair. + * + * Explicit user selections are preserved. Automatic selection retains the + * legacy single-point rule for ordinary linear Tribol meshes and selects + * multipoint integration whenever native MFEM parent-face data are available. + * + * @param penalty_options CommonPlane penalty options + * @param first_mesh First contact surface mesh + * @param second_mesh Second contact surface mesh + * @return true when the overlap must use multipoint integration + */ +TRIBOL_HOST_DEVICE inline bool UseMultipleIntegrationPoints( const PenaltyEnforcementOptions& penalty_options, + const MeshData::Viewer& first_mesh, + const MeshData::Viewer& second_mesh ) +{ + if ( penalty_options.common_plane_rule == MULTI_POINT ) { + return true; + } + if ( penalty_options.common_plane_rule == SINGLE_POINT ) { + return false; } + return first_mesh.hasParentFaceFields() || second_mesh.hasParentFaceFields(); +} - // check the gap rate sense. - // (v1-v2) * \nu < 0 : velocities lead to more interpenetration; - // note, \nu is in direction of face_2 outward unit normal - // TODO consider a velocity gap tolerance. Checking this against - // 0. actually smoothed out contact behavior in contact problem 1 - // for certain percent rate penalties. - if ( plane.m_velGap <= 0. ) // TODO do we want = or just parent_order ) { + parent_order = second_mesh.getParentFaceData().m_parent_face_orders[second_face_id]; + } + const int requested_order = 2 * parent_order; + return requested_order < 2 ? 2 : requested_order > 10 ? 10 : requested_order; +} -} // end ComputeGapRatePressure() +/** + * @brief Solve a three-by-three linear system using Cramer's rule. + * + * @param matrix System matrix + * @param right_hand_side System right-hand-side vector + * @param solution Computed solution vector + * @return true when the matrix determinant is sufficiently far from zero + */ +TRIBOL_HOST_DEVICE inline bool SolveThreeByThreeLinearSystem( const RealT matrix[3][3], const RealT right_hand_side[3], + RealT solution[3] ) +{ + const RealT determinant = matrix[0][0] * ( matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1] ) - + matrix[0][1] * ( matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0] ) + + matrix[0][2] * ( matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0] ); + constexpr RealT determinant_tolerance = 1.e-15; + if ( std::abs( determinant ) <= determinant_tolerance ) { + return false; + } -//------------------------------------------------------------------------------ -template <> -int ApplyNormal( CouplingScheme* cs ) + const RealT inverse_determinant = 1. / determinant; + solution[0] = inverse_determinant * + ( right_hand_side[0] * ( matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1] ) - + matrix[0][1] * ( right_hand_side[1] * matrix[2][2] - matrix[1][2] * right_hand_side[2] ) + + matrix[0][2] * ( right_hand_side[1] * matrix[2][1] - matrix[1][1] * right_hand_side[2] ) ); + solution[1] = inverse_determinant * + ( matrix[0][0] * ( right_hand_side[1] * matrix[2][2] - matrix[1][2] * right_hand_side[2] ) - + right_hand_side[0] * ( matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0] ) + + matrix[0][2] * ( matrix[1][0] * right_hand_side[2] - right_hand_side[1] * matrix[2][0] ) ); + solution[2] = + inverse_determinant * ( matrix[0][0] * ( matrix[1][1] * right_hand_side[2] - right_hand_side[1] * matrix[2][1] ) - + matrix[0][1] * ( matrix[1][0] * right_hand_side[2] - right_hand_side[1] * matrix[2][0] ) + + right_hand_side[0] * ( matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0] ) ); + return true; +} + +/** + * @brief Evaluate physical position and an optional nodal field from basis values. + * + * @param face_coordinates Node-major physical face coordinates + * @param number_of_nodes Number of face nodes + * @param spatial_dimension Number of coordinate components per face node + * @param basis_values Basis values at the evaluation point + * @param face_position Evaluated physical position + * @param value_dimension Number of components in the optional nodal field + * @param nodal_values Optional node-major field values + * @param interpolated_values Optional evaluated field values + */ +TRIBOL_HOST_DEVICE inline void EvaluateFaceFieldsFromBasis( const RealT* face_coordinates, int number_of_nodes, + int spatial_dimension, const RealT* basis_values, + RealT face_position[3], int value_dimension = 0, + const RealT* nodal_values = nullptr, + RealT* interpolated_values = nullptr ) { - /////////////////////////////// - // loop over interface pairs // - /////////////////////////////// - ArrayT err_data( { 0 }, cs->getAllocatorId() ); - ArrayViewT err = err_data; - ArrayT neg_thickness_data( { false }, cs->getAllocatorId() ); - ArrayViewT neg_thickness = neg_thickness_data; - auto cs_view = cs->getView(); - const auto num_pairs = cs->getNumActivePairs(); - forAllExec( cs->getExecutionMode(), num_pairs, [cs_view, err, neg_thickness] TRIBOL_HOST_DEVICE( IndexT i ) { - auto& cg_view = cs_view.getCompGeomView(); - auto& plane = cg_view.getCommonPlane( i ); - - auto& mesh1 = cs_view.getMesh1View(); - auto& mesh2 = cs_view.getMesh2View(); - - // get pair indices - IndexT index1 = plane.getCpElementId1(); - IndexT index2 = plane.getCpElementId2(); - - RealT gap = plane.m_gap; - RealT A = plane.m_area; // face-pair overlap area - - // don't proceed for gaps that don't violate the constraints. This check - // 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; + initRealArray( face_position, max_dim, 0. ); + if ( interpolated_values != nullptr ) { + initRealArray( interpolated_values, value_dimension, 0. ); + } + + for ( int node = 0; node < number_of_nodes; ++node ) { + for ( int component = 0; component < spatial_dimension; ++component ) { + face_position[component] += face_coordinates[spatial_dimension * node + component] * basis_values[node]; } - // debug force sums - // RealT dbg_sum_force1 {0.}; - // RealT dbg_sum_force2 {0.}; - ///////////////////////////////////////////// - // kinematic penalty stiffness calculation // - ///////////////////////////////////////////// - RealT penalty_stiff_per_area{ 0. }; - auto& enforcement_options = cs_view.getEnforcementOptions(); - const PenaltyEnforcementOptions& pen_enfrc_options = enforcement_options.penalty_options; - RealT pen_scale1 = mesh1.getElementData().m_penalty_scale; - RealT pen_scale2 = mesh2.getElementData().m_penalty_scale; - switch ( pen_enfrc_options.kinematic_calculation ) { - case KINEMATIC_CONSTANT: { - // pre-multiply each spring stiffness by each mesh's penalty scale - auto stiffness1 = pen_scale1 * mesh1.getElementData().m_penalty_stiffness; - auto stiffness2 = pen_scale2 * mesh2.getElementData().m_penalty_stiffness; - // compute the equivalent contact penalty spring stiffness per area - penalty_stiff_per_area = ComputePenaltyStiffnessPerArea( stiffness1, stiffness2 ); - break; + if ( interpolated_values != nullptr ) { + for ( int component = 0; component < value_dimension; ++component ) { + interpolated_values[component] += nodal_values[component + node * value_dimension] * basis_values[node]; } - case KINEMATIC_ELEMENT: { - // add tiny_length to element thickness to avoid division by zero - auto t1 = mesh1.getElementData().m_thickness[index1] + pen_enfrc_options.tiny_length; - auto t2 = mesh2.getElementData().m_thickness[index2] + pen_enfrc_options.tiny_length; - - if ( t1 < 0. || t2 < 0. ) { - neg_thickness[0] = true; - err[0] = 1; - } + } + } +} - // compute each element spring stiffness. Pre-multiply the material modulus - // (i.e. material stiffness) by each mesh's penalty scale - auto stiffness1 = pen_scale1 * mesh1.getElementData().m_mat_mod[index1] / t1; - auto stiffness2 = pen_scale2 * mesh2.getElementData().m_mat_mod[index2] / t2; - // compute the equivalent contact penalty spring stiffness per area - penalty_stiff_per_area = ComputePenaltyStiffnessPerArea( stiffness1, stiffness2 ); - break; +/** + * @brief Evaluate a linear face at a point projected along the CommonPlane normal. + * + * @param face_coordinates Node-major physical face coordinates + * @param number_of_nodes Number of linear triangle or quadrilateral nodes + * @param query_point Physical CommonPlane point to project + * @param projection_direction Direction used by the projection equation + * @param face_position Evaluated physical position on the face + * @param basis_values Evaluated linear face basis values + * @param field_dimension Number of components in the optional nodal field + * @param nodal_values Optional node-major field values + * @param interpolated_values Optional evaluated field values + * @param lor_reference_coordinates Optional evaluated LOR reference coordinates + * @return true when the projection converges inside the face tolerance + */ +TRIBOL_HOST_DEVICE inline bool EvaluateLinearFaceAtProjectedPoint( + const RealT* face_coordinates, const int number_of_nodes, const RealT query_point[3], + const RealT projection_direction[3], RealT face_position[3], RealT* basis_values, const int field_dimension = 0, + const RealT* nodal_values = nullptr, RealT* interpolated_values = nullptr, + RealT* lor_reference_coordinates = 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( basis_values, number_of_nodes, 0. ); + + if ( number_of_nodes == 3 ) { + const RealT first_vertex[3] = { face_coordinates[0], face_coordinates[1], face_coordinates[2] }; + const RealT first_edge[3] = { face_coordinates[3] - first_vertex[0], face_coordinates[4] - first_vertex[1], + face_coordinates[5] - first_vertex[2] }; + const RealT second_edge[3] = { face_coordinates[6] - first_vertex[0], face_coordinates[7] - first_vertex[1], + face_coordinates[8] - first_vertex[2] }; + + RealT face_normal[3]; + crossProd( first_edge[0], first_edge[1], first_edge[2], second_edge[0], second_edge[1], second_edge[2], + face_normal[0], face_normal[1], face_normal[2] ); + + const RealT projection_denominator = + dotProd( face_normal[0], face_normal[1], face_normal[2], projection_direction[0], projection_direction[1], + projection_direction[2] ); + constexpr RealT parallel_tolerance = 1.e-14; + if ( std::abs( projection_denominator ) <= parallel_tolerance ) { + return false; + } + + const RealT query_to_first_vertex[3] = { first_vertex[0] - query_point[0], first_vertex[1] - query_point[1], + first_vertex[2] - query_point[2] }; + const RealT projection_distance = dotProd( face_normal[0], face_normal[1], face_normal[2], query_to_first_vertex[0], + query_to_first_vertex[1], query_to_first_vertex[2] ) / + projection_denominator; + + face_position[0] = query_point[0] + projection_distance * projection_direction[0]; + face_position[1] = query_point[1] + projection_distance * projection_direction[1]; + face_position[2] = query_point[2] + projection_distance * projection_direction[2]; + + const RealT first_vertex_to_position[3] = { face_position[0] - first_vertex[0], face_position[1] - first_vertex[1], + face_position[2] - first_vertex[2] }; + const RealT first_edge_squared_norm = + dotProd( first_edge[0], first_edge[1], first_edge[2], first_edge[0], first_edge[1], first_edge[2] ); + const RealT edge_inner_product = + dotProd( first_edge[0], first_edge[1], first_edge[2], second_edge[0], second_edge[1], second_edge[2] ); + const RealT second_edge_squared_norm = + dotProd( second_edge[0], second_edge[1], second_edge[2], second_edge[0], second_edge[1], second_edge[2] ); + const RealT first_edge_position_projection = + dotProd( first_edge[0], first_edge[1], first_edge[2], first_vertex_to_position[0], first_vertex_to_position[1], + first_vertex_to_position[2] ); + const RealT second_edge_position_projection = + dotProd( second_edge[0], second_edge[1], second_edge[2], first_vertex_to_position[0], + first_vertex_to_position[1], first_vertex_to_position[2] ); + + const RealT gram_determinant = + first_edge_squared_norm * second_edge_squared_norm - edge_inner_product * edge_inner_product; + constexpr RealT gram_tolerance = 1.e-15; + if ( std::abs( gram_determinant ) <= gram_tolerance ) { + return false; + } + + const RealT xi = ( second_edge_squared_norm * first_edge_position_projection - + edge_inner_product * second_edge_position_projection ) / + gram_determinant; + const RealT eta = ( first_edge_squared_norm * second_edge_position_projection - + edge_inner_product * first_edge_position_projection ) / + gram_determinant; + if ( lor_reference_coordinates != nullptr ) { + lor_reference_coordinates[0] = xi; + lor_reference_coordinates[1] = eta; + } + basis_values[0] = 1. - xi - eta; + basis_values[1] = xi; + basis_values[2] = eta; + } else if ( number_of_nodes == 4 ) { + constexpr int maximum_iterations = 25; + constexpr RealT step_tolerance = 1.e-12; + constexpr RealT residual_tolerance = 1.e-12; + constexpr RealT reference_coordinate_tolerance = 1.e-8; + + RealT xi = 0.; + RealT eta = 0.; + RealT center_basis_values[max_nodes_per_face] = { 0.25, 0.25, 0.25, 0.25 }; + RealT center_position[3]; + EvaluateFaceFieldsFromBasis( face_coordinates, number_of_nodes, max_dim, center_basis_values, center_position ); + RealT projection_distance = ( center_position[0] - query_point[0] ) * projection_direction[0] + + ( center_position[1] - query_point[1] ) * projection_direction[1] + + ( center_position[2] - query_point[2] ) * projection_direction[2]; + bool converged = false; + + for ( int iteration = 0; iteration < maximum_iterations; ++iteration ) { + const RealT xi_node_signs[4] = { 1., -1., -1., 1. }; + const RealT eta_node_signs[4] = { 1., 1., -1., -1. }; + RealT position_derivative_xi[3] = { 0., 0., 0. }; + RealT position_derivative_eta[3] = { 0., 0., 0. }; + initRealArray( basis_values, number_of_nodes, 0. ); + initRealArray( face_position, max_dim, 0. ); + + for ( int node_index = 0; node_index < number_of_nodes; ++node_index ) { + basis_values[node_index] = + 0.25 * ( 1. + xi_node_signs[node_index] * xi ) * ( 1. + eta_node_signs[node_index] * eta ); + const RealT basis_derivative_xi = 0.25 * xi_node_signs[node_index] * ( 1. + eta_node_signs[node_index] * eta ); + const RealT basis_derivative_eta = 0.25 * eta_node_signs[node_index] * ( 1. + xi_node_signs[node_index] * xi ); + + const RealT node_x = face_coordinates[3 * node_index]; + const RealT node_y = face_coordinates[3 * node_index + 1]; + const RealT node_z = face_coordinates[3 * node_index + 2]; + + face_position[0] += node_x * basis_values[node_index]; + face_position[1] += node_y * basis_values[node_index]; + face_position[2] += node_z * basis_values[node_index]; + + position_derivative_xi[0] += node_x * basis_derivative_xi; + position_derivative_xi[1] += node_y * basis_derivative_xi; + position_derivative_xi[2] += node_z * basis_derivative_xi; + + position_derivative_eta[0] += node_x * basis_derivative_eta; + position_derivative_eta[1] += node_y * basis_derivative_eta; + position_derivative_eta[2] += node_z * basis_derivative_eta; } - default: - // no-op, quiet compiler - break; - } // end switch on kinematic penalty calculation option - - //////////////////////////////////////////////////// - // Compute contact pressure(s) on current overlap // - //////////////////////////////////////////////////// - - // compute total pressure based on constraint type - RealT totalPressure = 0.; - RealT residual_gap = cs_view.getParameters().residual_gap; - plane.m_pressure = ( gap - residual_gap ) * penalty_stiff_per_area; // kinematic contribution - switch ( pen_enfrc_options.constraint_type ) { - case KINEMATIC_AND_RATE: { - // kinematic contribution - totalPressure += plane.m_pressure; - // add gap-rate contribution - totalPressure += - ComputeGapRatePressure( plane, mesh1, mesh2, penalty_stiff_per_area, pen_enfrc_options.rate_calculation ); + + RealT residual[3] = { face_position[0] - query_point[0] - projection_distance * projection_direction[0], + face_position[1] - query_point[1] - projection_distance * projection_direction[1], + face_position[2] - query_point[2] - projection_distance * projection_direction[2] }; + const RealT residual_norm = magnitude( residual[0], residual[1], residual[2] ); + if ( residual_norm <= residual_tolerance ) { + converged = true; break; } - case KINEMATIC: - // kinematic gap pressure contribution only - totalPressure += plane.m_pressure; - break; - default: - // no-op + + RealT jacobian[3][3] = { { position_derivative_xi[0], position_derivative_eta[0], -projection_direction[0] }, + { position_derivative_xi[1], position_derivative_eta[1], -projection_direction[1] }, + { position_derivative_xi[2], position_derivative_eta[2], -projection_direction[2] } }; + RealT newton_right_hand_side[3] = { -residual[0], -residual[1], -residual[2] }; + RealT newton_increment[3]; + if ( !SolveThreeByThreeLinearSystem( jacobian, newton_right_hand_side, newton_increment ) ) { + return false; + } + + xi += newton_increment[0]; + eta += newton_increment[1]; + projection_distance += newton_increment[2]; + + const RealT step_norm = magnitude( newton_increment[0], newton_increment[1], newton_increment[2] ); + if ( step_norm <= step_tolerance ) { + converged = true; break; - } // end switch on registered penalty enforcement option - - // debug prints. Comment out for now, but keep for future common plane - // debugging - // SLIC_DEBUG("gap: " << gap); - // SLIC_DEBUG("area: " << A); - // SLIC_DEBUG("penalty stiffness: " << penalty_stiff_per_area); - // SLIC_DEBUG("pressure: " << cpManager.m_pressure[ cpID ]); - - /////////////////////////////////////////// - // create surface contact element struct // - /////////////////////////////////////////// - - // 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]; - int dim = cs_view.spatialDimension(); - int num_nodes_per_face = mesh1.numberOfNodesPerElement(); - initRealArray( xf1, dim * num_nodes_per_face, 0. ); - initRealArray( xf2, dim * num_nodes_per_face, 0. ); - // initialize assuming 2d - auto xVert_size = 4; - auto numPolyVert = 2; - // update if we are in 3d - if ( dim == 3 ) { - numPolyVert = plane.m_numPolyVert; - xVert_size = 3 * numPolyVert; + } } - initRealArray( xVert, xVert_size, 0. ); - - // get current configuration, physical coordinates of each face - plane.getFace1Coords( &xf1[0], num_nodes_per_face ); - plane.getFace2Coords( &xf2[0], num_nodes_per_face ); - - // construct array of polygon overlap vertex coordinates - plane.getOverlapVertices( &xVert[0] ); - - // instantiate surface contact element struct. Note, this is done with current - // configuration face coordinates (i.e. NOT on the contact plane) and overlap - // coordinates ON the contact plane. The surface contact element does not need - // to be used this way, but the developer should do the book-keeping. - SurfaceContactElem cntctElem( dim, xf1, xf2, xVert, num_nodes_per_face, numPolyVert, &mesh1, &mesh2, index1, - index2 ); - - // set SurfaceContactElem face normals and overlap normal - RealT faceNormal1[max_dim]; - RealT faceNormal2[max_dim]; - RealT overlapNormal[max_dim]; - - mesh1.getFaceNormal( index1, faceNormal1 ); - mesh2.getFaceNormal( index2, faceNormal2 ); - overlapNormal[0] = plane.m_nX; - overlapNormal[1] = plane.m_nY; - if ( dim == 3 ) { - overlapNormal[2] = plane.m_nZ; + + if ( !converged || xi < -1. - reference_coordinate_tolerance || xi > 1. + reference_coordinate_tolerance || + eta < -1. - reference_coordinate_tolerance || eta > 1. + reference_coordinate_tolerance ) { + return false; } - cntctElem.faceNormal1 = faceNormal1; - cntctElem.faceNormal2 = faceNormal2; - 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 // - //////////////////////////////////////////////////////////////////////// - EvalWeakFormIntegral( cntctElem, phi1, phi2 ); - - /////////////////////////////////////////////////////////////////////// - // Computation of full contact nodal force contributions // - // (i.e. premultiplication of contact integrals by normal component, // - // contact pressure, and overlap area) // - /////////////////////////////////////////////////////////////////////// - - // RealT phi_sum_1 = 0.; - // RealT phi_sum_2 = 0.; - - // compute contact force (spring force) - RealT contact_force = totalPressure * A; - - RealT force_x = overlapNormal[0] * contact_force; - RealT force_y = overlapNormal[1] * contact_force; - RealT force_z = 0.; - if ( dim == 3 ) { - force_z = overlapNormal[2] * contact_force; + initRealArray( basis_values, number_of_nodes, 0. ); + basis_values[0] = 0.25 * ( 1. + xi ) * ( 1. + eta ); + basis_values[1] = 0.25 * ( 1. - xi ) * ( 1. + eta ); + basis_values[2] = 0.25 * ( 1. - xi ) * ( 1. - eta ); + basis_values[3] = 0.25 * ( 1. + xi ) * ( 1. - eta ); + if ( lor_reference_coordinates != nullptr ) { + // MFEM's square reference element is [0,1]^2. The Newton solve above + // uses the equivalent [-1,1]^2 coordinates aligned with Tribol's local + // quadrilateral vertex ordering. + lor_reference_coordinates[0] = 0.5 * ( 1.0 - xi ); + lor_reference_coordinates[1] = 0.5 * ( 1.0 - eta ); + } + EvaluateFaceFieldsFromBasis( face_coordinates, number_of_nodes, max_dim, basis_values, face_position ); + + const RealT normal_distance = ( face_position[0] - query_point[0] ) * projection_direction[0] + + ( face_position[1] - query_point[1] ) * projection_direction[1] + + ( face_position[2] - query_point[2] ) * projection_direction[2]; + const RealT projection_residual = + magnitude( face_position[0] - query_point[0] - normal_distance * projection_direction[0], + face_position[1] - query_point[1] - normal_distance * projection_direction[1], + face_position[2] - query_point[2] - normal_distance * projection_direction[2] ); + if ( projection_residual > residual_tolerance ) { + return false; } + } else { + return false; + } - ////////////////////////////////////////////////////// - // 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 ); + if ( interpolated_values != nullptr ) { + initRealArray( interpolated_values, field_dimension, 0. ); + for ( int node_index = 0; node_index < number_of_nodes; ++node_index ) { + for ( int component = 0; component < field_dimension; ++component ) { + interpolated_values[component] += + nodal_values[component + node_index * field_dimension] * basis_values[node_index]; } - } // end for loop over face nodes - - // comment out debug logs; too much output during tests. Keep for easy - // debugging if needed - // SLIC_DEBUG("force sum, side 1, pair " << kp << ": " << -dbg_sum_force1 ); - // SLIC_DEBUG("force sum, side 2, pair " << kp << ": " << dbg_sum_force2 ); - // SLIC_DEBUG("phi 1 sum: " << phi_sum_1 ); - // SLIC_DEBUG("phi 2 sum: " << phi_sum_2 ); - } ); + } + } - ArrayT neg_thickness_host( neg_thickness_data ); - SLIC_DEBUG_IF( neg_thickness_host[0], - "ApplyNormal: negative element thicknesses encountered." ); + return true; +} - ArrayT err_host( err_data ); - return err_host[0]; +/** + * @brief Evaluate a linear edge at a point projected along the CommonPlane normal. + * + * @param edge_coordinates Node-major physical edge coordinates + * @param query_point Physical CommonPlane point to project + * @param projection_direction Direction used by the projection equation + * @param edge_position Evaluated physical position on the edge + * @param basis_values Evaluated linear edge basis values + * @param field_dimension Number of components in the optional nodal field + * @param nodal_values Optional node-major field values + * @param interpolated_values Optional evaluated field values + * @param lor_reference_coordinates Optional evaluated LOR reference coordinate + * @return true when the projected point lies inside the edge tolerance + */ +TRIBOL_HOST_DEVICE inline bool EvaluateLinearEdgeAtProjectedPoint( + const RealT* edge_coordinates, const RealT query_point[2], const RealT projection_direction[2], + RealT edge_position[3], RealT* basis_values, const int field_dimension = 0, const RealT* nodal_values = nullptr, + RealT* interpolated_values = nullptr, RealT* lor_reference_coordinates = nullptr ) +{ + const RealT first_endpoint_x = edge_coordinates[0]; + const RealT first_endpoint_y = edge_coordinates[1]; + const RealT second_endpoint_x = edge_coordinates[2]; + const RealT second_endpoint_y = edge_coordinates[3]; + const RealT edge_direction_x = second_endpoint_x - first_endpoint_x; + const RealT edge_direction_y = second_endpoint_y - first_endpoint_y; + + const RealT projection_determinant = + projection_direction[0] * edge_direction_y - edge_direction_x * projection_direction[1]; + constexpr RealT determinant_tolerance = 1.e-14; + if ( std::abs( projection_determinant ) <= determinant_tolerance ) { + return false; + } -} // end ApplyNormal() + const RealT query_offset_x = query_point[0] - first_endpoint_x; + const RealT query_offset_y = query_point[1] - first_endpoint_y; + RealT edge_parameter = + ( projection_direction[0] * query_offset_y - projection_direction[1] * query_offset_x ) / projection_determinant; -//------------------------------------------------------------------------------ -template <> -int ApplyTangential( CouplingScheme* cs ) -{ - /////////////////////////////// - // loop over interface pairs // - /////////////////////////////// - auto cs_view = cs->getView(); - const auto num_pairs = cs->getNumActivePairs(); - forAllExec( cs->getExecutionMode(), num_pairs, [cs_view] TRIBOL_HOST_DEVICE( IndexT i ) { - auto& cg_view = cs_view.getCompGeomView(); - auto& plane = cg_view.getCommonPlane( i ); - - if ( !plane.m_inContact ) { - return; - } + constexpr RealT edge_tolerance = 1.e-8; + if ( edge_parameter < -edge_tolerance || edge_parameter > 1. + edge_tolerance ) { + return false; + } + edge_parameter = std::max( 0., std::min( 1., edge_parameter ) ); + if ( lor_reference_coordinates != nullptr ) { + lor_reference_coordinates[0] = edge_parameter; + } - const auto dim = plane.m_dim; - auto& mesh1 = cs_view.getMesh1View(); - auto& mesh2 = cs_view.getMesh2View(); - - // 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; - - 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; - auto numNodesPerFace2 = mesh2.numberOfNodesPerElement(); - plane.getFace2Coords( x2, numNodesPerFace2 ); // get avg face coords off the contact plane - mesh2.getFaceVelocities( index2, v2 ); - - ////////////////////////////////////////////////////////// - // compute velocity Galerkin approximation at projected // - // overlap centroid // - ////////////////////////////////////////////////////////// - RealT vel_f1[max_dim]; - RealT vel_f2[max_dim]; - initRealArray( vel_f1, dim, 0. ); - initRealArray( vel_f2, dim, 0. ); - - // interpolate nodal velocity at overlap centroid as projected - // onto face 1 - RealT cXf1 = plane.m_cXf1; - RealT cYf1 = plane.m_cYf1; - RealT cZf1 = ( dim == 3 ) ? plane.m_cZf1 : 0.; - GalerkinEval( x1, cXf1, cYf1, cZf1, LINEAR, PHYSICAL, dim, dim, v1, vel_f1 ); - - // interpolate nodal velocity at overlap centroid as projected - // onto face 2 - RealT cXf2 = plane.m_cXf2; - RealT cYf2 = plane.m_cYf2; - RealT cZf2 = ( dim == 3 ) ? plane.m_cZf2 : 0.; - GalerkinEval( x2, cXf2, cYf2, cZf2, LINEAR, PHYSICAL, dim, dim, v2, vel_f2 ); - - // compute velocity gap vector - RealT velGap[max_dim]; - velGap[0] = vel_f1[0] - vel_f2[0]; - velGap[1] = vel_f1[1] - vel_f2[1]; - if ( dim == 3 ) { - velGap[2] = vel_f1[2] - vel_f2[2]; - } + basis_values[0] = 1. - edge_parameter; + basis_values[1] = edge_parameter; - // subtract off the common-plane normal component of the velocity gap - RealT velGap_dot_n = velGap[0] * plane.m_nX + velGap[1] * plane.m_nY; - if ( dim == 3 ) { - velGap_dot_n += velGap[2] * plane.m_nZ; + edge_position[0] = first_endpoint_x + edge_parameter * edge_direction_x; + edge_position[1] = first_endpoint_y + edge_parameter * edge_direction_y; + edge_position[2] = 0.; + + if ( interpolated_values != nullptr ) { + initRealArray( interpolated_values, field_dimension, 0. ); + for ( int node_index = 0; node_index < 2; ++node_index ) { + for ( int component = 0; component < field_dimension; ++component ) { + interpolated_values[component] += + nodal_values[component + node_index * field_dimension] * basis_values[node_index]; + } } - RealT velGapTan[max_dim]; - velGapTan[0] = velGap[0] - velGap_dot_n * plane.m_nX; - velGapTan[1] = velGap[1] - velGap_dot_n * plane.m_nY; - if ( dim == 3 ) { - velGapTan[2] = velGap[2] - velGap_dot_n * plane.m_nZ; + } + + return true; +} + +/** + * @brief Evaluate native parent-face fields at a projected CommonPlane point. + * + * The physical point is first projected to the linear LOR face. Its LOR + * reference coordinates are mapped through the stored LOR-to-parent mapping, + * after which the native parent basis evaluates position and optional velocity. + * + * @param mesh Contact surface mesh + * @param face_id Tribol LOR face identifier + * @param lor_face_coordinates Physical coordinates of the LOR face vertices + * @param physical_point CommonPlane quadrature point + * @param projection_direction CommonPlane projection direction + * @param evaluate_velocity Whether velocity values are required + * @param parent_basis_values Native parent-face basis values + * @param parent_position Native parent-face position at the mapped point + * @param parent_velocity Native parent-face velocity at the mapped point + * @param mapped_parent_reference_coordinates Native parent reference coordinates at the mapped point + * @return zero when both LOR and native parent mappings succeed; nonzero otherwise + */ +TRIBOL_HOST_DEVICE inline int EvaluateParentFaceAtProjectedPoint( + const MeshData::Viewer& mesh, IndexT face_id, const RealT* lor_face_coordinates, const RealT* physical_point, + const RealT* projection_direction, bool evaluate_velocity, RealT* parent_basis_values, RealT* parent_position, + RealT* parent_velocity, RealT* mapped_parent_reference_coordinates ) +{ + RealT lor_basis_values[4] = { 0.0, 0.0, 0.0, 0.0 }; + RealT projected_lor_position[3] = { 0.0, 0.0, 0.0 }; + RealT lor_reference_coordinates[2] = { 0.0, 0.0 }; + const bool mapped_to_lor = + mesh.spatialDimension() == 2 + ? EvaluateLinearEdgeAtProjectedPoint( lor_face_coordinates, physical_point, projection_direction, + projected_lor_position, lor_basis_values, 0, nullptr, nullptr, + lor_reference_coordinates ) + : EvaluateLinearFaceAtProjectedPoint( lor_face_coordinates, mesh.numberOfNodesPerElement(), physical_point, + projection_direction, projected_lor_position, lor_basis_values, 0, + nullptr, nullptr, lor_reference_coordinates ); + if ( !mapped_to_lor ) { + return 1; + } + + RealT parent_reference_coordinates[2] = { 0.0, 0.0 }; + if ( !mesh.mapToParentReference( face_id, lor_reference_coordinates, parent_reference_coordinates ) ) { + return 2; + } + if ( !mesh.evaluateParentFaceBasis( face_id, parent_reference_coordinates, parent_basis_values ) ) { + return 3; + } + + mapped_parent_reference_coordinates[0] = parent_reference_coordinates[0]; + mapped_parent_reference_coordinates[1] = parent_reference_coordinates[1]; + + mesh.evaluateParentFaceFields( face_id, parent_basis_values, parent_position, + evaluate_velocity ? parent_velocity : nullptr ); + return 0; +} + +/** + * @brief Scatter an equal-and-opposite force through linear contact-face bases. + * + * @param first_mesh First contact surface mesh + * @param second_mesh Second contact surface mesh + * @param first_face_id First contact-surface face identifier + * @param second_face_id Second contact-surface face identifier + * @param spatial_dimension Contact problem dimension + * @param number_of_nodes_per_face Number of nodes on each linear face + * @param force_x Integrated x force applied to the second face + * @param force_y Integrated y force applied to the second face + * @param force_z Integrated z force applied to the second face + * @param first_basis_values Basis values on the first face + * @param second_basis_values Basis values on the second face + */ +TRIBOL_HOST_DEVICE inline void AccumulateContactForce( + const MeshData::Viewer& first_mesh, const MeshData::Viewer& second_mesh, const IndexT first_face_id, + const IndexT second_face_id, const int spatial_dimension, const int number_of_nodes_per_face, const RealT force_x, + const RealT force_y, const RealT force_z, const RealT* first_basis_values, const RealT* second_basis_values ) +{ + for ( IndexT basis_index = 0; basis_index < number_of_nodes_per_face; ++basis_index ) { + const IndexT first_node_id = first_mesh.getGlobalNodeId( first_face_id, basis_index ); + const IndexT second_node_id = second_mesh.getGlobalNodeId( second_face_id, basis_index ); + + const RealT first_nodal_force_x = force_x * first_basis_values[basis_index]; + const RealT first_nodal_force_y = force_y * first_basis_values[basis_index]; + const RealT first_nodal_force_z = force_z * first_basis_values[basis_index]; + + const RealT second_nodal_force_x = force_x * second_basis_values[basis_index]; + const RealT second_nodal_force_y = force_y * second_basis_values[basis_index]; + const RealT second_nodal_force_z = force_z * second_basis_values[basis_index]; + + tribol::atomicAdd( &first_mesh.getResponse()[0][first_node_id], -first_nodal_force_x ); + tribol::atomicAdd( &second_mesh.getResponse()[0][second_node_id], second_nodal_force_x ); + + tribol::atomicAdd( &first_mesh.getResponse()[1][first_node_id], -first_nodal_force_y ); + tribol::atomicAdd( &second_mesh.getResponse()[1][second_node_id], second_nodal_force_y ); + + if ( spatial_dimension == 3 ) { + tribol::atomicAdd( &first_mesh.getResponse()[2][first_node_id], -first_nodal_force_z ); + tribol::atomicAdd( &second_mesh.getResponse()[2][second_node_id], second_nodal_force_z ); } + } +} - // setup the contact element struct for purposes of evaluating basis functions on overlap - // initialize assuming 2d - RealT xVert[max_dim * max_nodes_per_overlap]; - auto xVert_size = 4; - auto numPolyVert = 2; - // update if we are in 3d - if ( dim == 3 ) { - numPolyVert = plane.m_numPolyVert; - xVert_size = 3 * numPolyVert; +/** + * @brief Scatter an equal-and-opposite force through native parent-face bases. + * + * @param mesh1 First contact surface mesh + * @param mesh2 Second contact surface mesh + * @param face_id1 First Tribol face identifier + * @param face_id2 Second Tribol face identifier + * @param force Force applied to the second face + * @param parent_basis_values1 Native basis values on the first face + * @param parent_basis_values2 Native basis values on the second face + */ +TRIBOL_HOST_DEVICE inline void AccumulateParentContactForce( const MeshData::Viewer& mesh1, + const MeshData::Viewer& mesh2, IndexT face_id1, + IndexT face_id2, const RealT* force, + const RealT* parent_basis_values1, + const RealT* parent_basis_values2 ) +{ + RealT opposite_force[max_dim] = { -force[0], -force[1], -force[2] }; + mesh1.addParentFaceResponse( face_id1, parent_basis_values1, opposite_force ); + mesh2.addParentFaceResponse( face_id2, parent_basis_values2, force ); +} + +/** + * @brief Compute the CommonPlane penalty stiffness for one LOR face pair. + * + * @param first_mesh First contact surface mesh + * @param second_mesh Second contact surface mesh + * @param first_face_id First Tribol face identifier + * @param second_face_id Second Tribol face identifier + * @param penalty_options Registered penalty options + * @param penalty_stiffness Computed equivalent stiffness per unit overlap measure + * @return true when the registered penalty data are physically admissible + */ +TRIBOL_HOST_DEVICE inline bool ComputePairPenaltyStiffness( const MeshData::Viewer& first_mesh, + const MeshData::Viewer& second_mesh, IndexT first_face_id, + IndexT second_face_id, + const PenaltyEnforcementOptions& penalty_options, + RealT& penalty_stiffness ) +{ + const RealT first_penalty_scale = first_mesh.getElementData().m_penalty_scale; + const RealT second_penalty_scale = second_mesh.getElementData().m_penalty_scale; + + switch ( penalty_options.kinematic_calculation ) { + case KINEMATIC_CONSTANT: { + const RealT first_stiffness = first_penalty_scale * first_mesh.getElementData().m_penalty_stiffness; + const RealT second_stiffness = second_penalty_scale * second_mesh.getElementData().m_penalty_stiffness; + penalty_stiffness = ComputePenaltyStiffnessPerArea( first_stiffness, second_stiffness ); + return true; } - initRealArray( xVert, xVert_size, 0. ); - - // construct array of polygon overlap vertex coordinates - plane.getOverlapVertices( &xVert[0] ); - - // instantiate surface contact element struct. Note, this is done with current - // configuration face coordinates (i.e. NOT on the contact plane) and overlap - // coordinates ON the contact plane. The surface contact element does not need - // to be used this way, but the developer should do the book-keeping. - SurfaceContactElem cntctElem( dim, x1, x2, xVert, numNodesPerFace1, numPolyVert, &mesh1, &mesh2, index1, index2 ); - - // set SurfaceContactElem face normals and overlap normal - RealT faceNormal1[max_dim]; - RealT faceNormal2[max_dim]; - RealT overlapNormal[max_dim]; - - mesh1.getFaceNormal( index1, faceNormal1 ); - mesh2.getFaceNormal( index2, faceNormal2 ); - overlapNormal[0] = plane.m_nX; - overlapNormal[1] = plane.m_nY; - if ( dim == 3 ) { - overlapNormal[2] = plane.m_nZ; + case KINEMATIC_ELEMENT: { + const RealT first_thickness = + first_mesh.getElementData().m_thickness[first_face_id] + penalty_options.tiny_length; + const RealT second_thickness = + second_mesh.getElementData().m_thickness[second_face_id] + penalty_options.tiny_length; + if ( first_thickness <= 0.0 || second_thickness <= 0.0 ) { + return false; + } + const RealT first_stiffness = + first_penalty_scale * first_mesh.getElementData().m_mat_mod[first_face_id] / first_thickness; + const RealT second_stiffness = + second_penalty_scale * second_mesh.getElementData().m_mat_mod[second_face_id] / second_thickness; + penalty_stiffness = ComputePenaltyStiffnessPerArea( first_stiffness, second_stiffness ); + return true; } + default: + penalty_stiffness = 0.0; + return false; + } +} - cntctElem.faceNormal1 = faceNormal1; - cntctElem.faceNormal2 = faceNormal2; - cntctElem.overlapNormal = overlapNormal; - 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]; - initRealArray( phi1, numNodesPerFace1, 0. ); - initRealArray( phi2, numNodesPerFace2, 0. ); - - //////////////////////////////////////////////////////////////////////// - // Integration of contact integrals: integral of shape functions over // - // contact overlap patch // - //////////////////////////////////////////////////////////////////////// - EvalWeakFormIntegral( cntctElem, phi1, phi2 ); - - ///////////////////////////////////////////////////// - // Computation of tangential viscous damping force // - ///////////////////////////////////////////////////// - RealT visc = - 0.5 * ( mesh1.getElementData().m_viscous_damping_coeff + mesh2.getElementData().m_viscous_damping_coeff ); - RealT force_x = visc * velGapTan[0]; - RealT force_y = visc * velGapTan[1]; - RealT force_z = 0.; - if ( dim == 3 ) { - force_z = visc * velGapTan[2]; +/** + * @brief Store one fully evaluated CommonPlane quadrature row. + * + * The helper maps the physical integration point to each LOR face, maps the + * resulting reference coordinates to native parent space when available, and + * stores all field and basis values needed by downstream explicit operators. + * + * @param rows Writable CommonPlane row-batch view + * @param row_id Stable row slot assigned to this integration point + * @param contact_pair_id Active CommonPlane pair identifier + * @param first_mesh First contact surface mesh + * @param second_mesh Second contact surface mesh + * @param first_face_id First Tribol face identifier + * @param second_face_id Second Tribol face identifier + * @param first_face_coordinates First physical LOR face coordinates + * @param second_face_coordinates Second physical LOR face coordinates + * @param first_projected_face_coordinates First LOR face projected onto the CommonPlane + * @param second_projected_face_coordinates Second LOR face projected onto the CommonPlane + * @param first_face_velocities First LOR face velocities when requested + * @param second_face_velocities Second LOR face velocities when requested + * @param integration_point Physical CommonPlane integration point + * @param normal Consistently oriented CommonPlane unit normal + * @param integration_weight Physical integration measure for this row + * @param penalty_stiffness Equivalent normal penalty stiffness per unit measure + * @param rate_penalty_coefficient Normal rate-penalty coefficient per unit measure + * @param tangential_viscous_coefficient Tangential viscous coefficient per unit measure + * @param gap_tolerance Normal contact activation tolerance + * @param evaluate_velocity Whether velocity fields are required by any consumer + * @param use_parent_fields Whether fields and responses use native parent faces + * @return zero when both face mappings and field evaluations succeed; nonzero otherwise + */ +TRIBOL_HOST_DEVICE inline int StoreCommonPlaneContactRow( + CommonPlaneContactData::Viewer rows, IndexT row_id, IndexT contact_pair_id, const MeshData::Viewer& first_mesh, + const MeshData::Viewer& second_mesh, IndexT first_face_id, IndexT second_face_id, + const RealT* first_face_coordinates, const RealT* second_face_coordinates, const RealT* first_face_velocities, + const RealT* second_face_velocities, const RealT* first_projected_face_coordinates, + const RealT* second_projected_face_coordinates, const RealT* integration_point, const RealT* normal, + RealT integration_weight, RealT penalty_stiffness, RealT rate_penalty_coefficient, + RealT tangential_viscous_coefficient, RealT gap_tolerance, bool evaluate_velocity, bool use_parent_fields ) +{ + RealT first_basis_values[max_nodes_per_face] = { 0.0 }; + RealT second_basis_values[max_nodes_per_face] = { 0.0 }; + RealT first_position[max_dim] = { 0.0, 0.0, 0.0 }; + RealT second_position[max_dim] = { 0.0, 0.0, 0.0 }; + RealT first_velocity[max_dim] = { 0.0, 0.0, 0.0 }; + RealT second_velocity[max_dim] = { 0.0, 0.0, 0.0 }; + RealT first_reference_coordinates[ParentFaceData::max_reference_dimension] = { 0.0, 0.0 }; + RealT second_reference_coordinates[ParentFaceData::max_reference_dimension] = { 0.0, 0.0 }; + + const int spatial_dimension = rows.spatial_dimension; + const int first_basis_count = use_parent_fields ? first_mesh.getParentFaceData().m_parent_node_counts[first_face_id] + : first_mesh.numberOfNodesPerElement(); + const int second_basis_count = use_parent_fields + ? second_mesh.getParentFaceData().m_parent_node_counts[second_face_id] + : second_mesh.numberOfNodesPerElement(); + + int first_mapping_status = 0; + int second_mapping_status = 0; + if ( use_parent_fields ) { + first_mapping_status = EvaluateParentFaceAtProjectedPoint( + first_mesh, first_face_id, first_projected_face_coordinates, integration_point, normal, evaluate_velocity, + first_basis_values, first_position, first_velocity, first_reference_coordinates ); + second_mapping_status = EvaluateParentFaceAtProjectedPoint( + second_mesh, second_face_id, second_projected_face_coordinates, integration_point, normal, evaluate_velocity, + second_basis_values, second_position, second_velocity, second_reference_coordinates ); + } else if ( spatial_dimension == 2 ) { + first_mapping_status = + EvaluateLinearEdgeAtProjectedPoint( first_projected_face_coordinates, integration_point, normal, first_position, + first_basis_values, 0, nullptr, nullptr, first_reference_coordinates ) + ? 0 + : 1; + second_mapping_status = EvaluateLinearEdgeAtProjectedPoint( second_projected_face_coordinates, integration_point, + normal, second_position, second_basis_values, 0, + nullptr, nullptr, second_reference_coordinates ) + ? 0 + : 1; + } else { + first_mapping_status = EvaluateLinearFaceAtProjectedPoint( + first_projected_face_coordinates, first_basis_count, integration_point, normal, + first_position, first_basis_values, 0, nullptr, nullptr, first_reference_coordinates ) + ? 0 + : 1; + second_mapping_status = + EvaluateLinearFaceAtProjectedPoint( second_projected_face_coordinates, second_basis_count, integration_point, + normal, second_position, second_basis_values, 0, nullptr, nullptr, + second_reference_coordinates ) + ? 0 + : 1; + } + + if ( first_mapping_status != 0 || second_mapping_status != 0 ) { + return first_mapping_status != 0 ? first_mapping_status : 10 + second_mapping_status; + } + + if ( !use_parent_fields ) { + EvaluateFaceFieldsFromBasis( first_face_coordinates, first_basis_count, spatial_dimension, first_basis_values, + first_position, evaluate_velocity ? spatial_dimension : 0, + evaluate_velocity ? first_face_velocities : nullptr, + evaluate_velocity ? first_velocity : nullptr ); + EvaluateFaceFieldsFromBasis( second_face_coordinates, second_basis_count, spatial_dimension, second_basis_values, + second_position, evaluate_velocity ? spatial_dimension : 0, + evaluate_velocity ? second_face_velocities : nullptr, + evaluate_velocity ? second_velocity : nullptr ); + } + + rows.row_is_valid[row_id] = 1; + rows.contact_pair_ids[row_id] = contact_pair_id; + rows.first_face_ids[row_id] = first_face_id; + rows.second_face_ids[row_id] = second_face_id; + rows.first_basis_counts[row_id] = first_basis_count; + rows.second_basis_counts[row_id] = second_basis_count; + rows.row_uses_parent_fields[row_id] = use_parent_fields ? 1 : 0; + rows.integration_weights[row_id] = integration_weight; + rows.penalty_stiffnesses[row_id] = penalty_stiffness; + rows.rate_penalty_coefficients[row_id] = rate_penalty_coefficient; + rows.tangential_viscous_coefficients[row_id] = tangential_viscous_coefficient; + + RealT normal_gap = 0.0; + RealT normal_velocity_gap = 0.0; + for ( int component = 0; component < max_dim; ++component ) { + rows.integration_points( row_id, component ) = component < spatial_dimension ? integration_point[component] : 0.0; + rows.first_positions( row_id, component ) = first_position[component]; + rows.second_positions( row_id, component ) = second_position[component]; + rows.first_velocities( row_id, component ) = first_velocity[component]; + rows.second_velocities( row_id, component ) = second_velocity[component]; + rows.normals( row_id, component ) = component < spatial_dimension ? normal[component] : 0.0; + if ( component < spatial_dimension ) { + normal_gap += ( first_position[component] - second_position[component] ) * normal[component]; + normal_velocity_gap += ( first_velocity[component] - second_velocity[component] ) * normal[component]; } + } + for ( int reference_component = 0; reference_component < ParentFaceData::max_reference_dimension; + ++reference_component ) { + rows.first_parent_reference_coordinates( row_id, reference_component ) = + first_reference_coordinates[reference_component]; + rows.second_parent_reference_coordinates( row_id, reference_component ) = + second_reference_coordinates[reference_component]; + } + for ( int basis_index = 0; basis_index < first_basis_count; ++basis_index ) { + rows.first_basis_values( row_id, basis_index ) = first_basis_values[basis_index]; + } + for ( int basis_index = 0; basis_index < second_basis_count; ++basis_index ) { + rows.second_basis_values( row_id, basis_index ) = second_basis_values[basis_index]; + } - ////////////////////////////////////////////////////// - // 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 ); + rows.gaps[row_id] = normal_gap; + rows.normal_velocity_gaps[row_id] = normal_velocity_gap; + rows.row_is_active[row_id] = normal_gap <= gap_tolerance ? 1 : 0; + return 0; +} - 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]; +/** + * @brief Scatter one vector force from a shared CommonPlane row. + * + * @param rows CommonPlane row-batch view + * @param row_id Row carrying basis and face identifiers + * @param first_mesh First contact surface mesh + * @param second_mesh Second contact surface mesh + * @param force_on_second_face Integrated force applied to the second face + */ +TRIBOL_HOST_DEVICE inline void ScatterCommonPlaneRowForce( const CommonPlaneContactData::Viewer& rows, IndexT row_id, + const MeshData::Viewer& first_mesh, + const MeshData::Viewer& second_mesh, + const RealT* force_on_second_face ) +{ + const RealT* first_basis_values = &rows.first_basis_values( row_id, 0 ); + const RealT* second_basis_values = &rows.second_basis_values( row_id, 0 ); + if ( rows.row_uses_parent_fields[row_id] != 0 ) { + AccumulateParentContactForce( first_mesh, second_mesh, rows.first_face_ids[row_id], rows.second_face_ids[row_id], + force_on_second_face, first_basis_values, second_basis_values ); + return; + } - 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]; + AccumulateContactForce( first_mesh, second_mesh, rows.first_face_ids[row_id], rows.second_face_ids[row_id], + rows.spatial_dimension, rows.first_basis_counts[row_id], force_on_second_face[0], + force_on_second_face[1], force_on_second_face[2], first_basis_values, second_basis_values ); +} - // 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 ); +} // namespace + +/** + * @brief Build the shared CommonPlane quadrature row batch for one update. + * + * One execution thread owns each accepted overlap cell. Rows are written into + * a deterministic pair-local range, which keeps generation device-resident and + * gives every row a stable identity without a host prefix sum. + * + * @param cs CommonPlane coupling scheme + * @return zero on success and nonzero when row generation detects invalid data + */ +int BuildCommonPlaneContactRows( CouplingScheme* cs ) +{ + auto* common_plane_data = static_cast( cs->getMethodData() ); + SLIC_ERROR_ROOT_IF( common_plane_data == nullptr, + "BuildCommonPlaneContactRows(): CommonPlane row storage is unavailable." ); + + const IndexT number_of_pairs = cs->getNumActivePairs(); + const int spatial_dimension = cs->spatialDimension(); + common_plane_data->resize( number_of_pairs, spatial_dimension, cs->getAllocatorId() ); + if ( number_of_pairs == 0 ) { + return 0; + } - tribol::atomicAdd( &mesh1.getResponse()[1][node0], -nodal_force_y1 ); - tribol::atomicAdd( &mesh2.getResponse()[1][node1], nodal_force_y2 ); + const CommonPlaneContactData::Viewer rows = common_plane_data->getView(); + const CouplingScheme::Viewer coupling_scheme = cs->getView(); + const bool tangential_velocity_is_required = cs->getContactModel() == VISCOUS_TANGENTIAL; + Array1D evaluation_error_data( { 0 }, cs->getAllocatorId() ); + Array1DView evaluation_error = evaluation_error_data.view(); + forAllExec( + cs->getExecutionMode(), number_of_pairs, + [rows, coupling_scheme, spatial_dimension, tangential_velocity_is_required, + evaluation_error] TRIBOL_HOST_DEVICE( IndexT contact_pair_id ) { + CommonPlanePair& common_plane = coupling_scheme.getCompGeomView().getCommonPlane( contact_pair_id ); + const MeshData::Viewer& first_mesh = coupling_scheme.getMesh1View(); + const MeshData::Viewer& second_mesh = coupling_scheme.getMesh2View(); + const PenaltyEnforcementOptions& penalty_options = coupling_scheme.getEnforcementOptions().penalty_options; + const IndexT first_face_id = common_plane.getCpElementId1(); + const IndexT second_face_id = common_plane.getCpElementId2(); + const IndexT first_row_id = rows.pairRowOffset( contact_pair_id ); + + rows.pair_row_counts[contact_pair_id] = 0; + common_plane.m_inContact = false; + common_plane.m_pressure = 0.0; + common_plane.m_ratePressure = 0.0; + common_plane.m_velGap = 0.0; + + const bool first_face_uses_parent_fields = first_mesh.hasParentFaceFields(); + const bool second_face_uses_parent_fields = second_mesh.hasParentFaceFields(); + if ( first_face_uses_parent_fields != second_face_uses_parent_fields ) { + rows.pair_evaluation_statuses[contact_pair_id] = + static_cast( CommonPlanePairEvaluationStatus::INVALID_PARENT_DATA ); + tribol::atomicMax( &evaluation_error[0], 1 ); + return; + } + const bool use_parent_fields = first_face_uses_parent_fields && second_face_uses_parent_fields; + const bool normal_rate_velocity_is_required = penalty_options.constraint_type == KINEMATIC_AND_RATE; + const bool evaluate_velocity = normal_rate_velocity_is_required || tangential_velocity_is_required; + const bool first_velocity_is_available = + use_parent_fields ? first_mesh.hasParentFaceVelocity() : first_mesh.hasVelocity(); + const bool second_velocity_is_available = + use_parent_fields ? second_mesh.hasParentFaceVelocity() : second_mesh.hasVelocity(); + if ( evaluate_velocity && ( !first_velocity_is_available || !second_velocity_is_available ) ) { + rows.pair_evaluation_statuses[contact_pair_id] = + static_cast( CommonPlanePairEvaluationStatus::INVALID_PARENT_DATA ); + tribol::atomicMax( &evaluation_error[0], 1 ); + return; + } + + RealT common_plane_normal[max_dim] = { common_plane.m_nX, common_plane.m_nY, + spatial_dimension == 3 ? common_plane.m_nZ : 0.0 }; + const RealT normal_magnitude = + spatial_dimension == 3 ? magnitude( common_plane_normal[0], common_plane_normal[1], common_plane_normal[2] ) + : magnitude( common_plane_normal[0], common_plane_normal[1] ); + if ( normal_magnitude <= 1.e-12 || std::abs( normal_magnitude - 1.0 ) > 1.e-6 ) { + rows.pair_evaluation_statuses[contact_pair_id] = + static_cast( CommonPlanePairEvaluationStatus::INCONSISTENT_NORMAL ); + tribol::atomicMax( &evaluation_error[0], 1 ); + return; + } + + RealT penalty_stiffness = 0.0; + if ( !ComputePairPenaltyStiffness( first_mesh, second_mesh, first_face_id, second_face_id, penalty_options, + penalty_stiffness ) ) { + rows.pair_evaluation_statuses[contact_pair_id] = + static_cast( CommonPlanePairEvaluationStatus::INVALID_PARENT_DATA ); + tribol::atomicMax( &evaluation_error[0], 1 ); + return; + } + const RealT rate_penalty_coefficient = + normal_rate_velocity_is_required + ? ComputeRatePenalty( first_mesh, second_mesh, penalty_stiffness, penalty_options.rate_calculation ) + : 0.0; + const RealT tangential_viscous_coefficient = + tangential_velocity_is_required ? 0.5 * ( first_mesh.getElementData().m_viscous_damping_coeff + + second_mesh.getElementData().m_viscous_damping_coeff ) + : 0.0; + const RealT gap_tolerance = coupling_scheme.getGapTol( first_face_id, second_face_id ); + + StackArrayT first_face_coordinates; + StackArrayT second_face_coordinates; + StackArrayT first_projected_face_coordinates; + StackArrayT second_projected_face_coordinates; + StackArrayT first_face_velocities; + StackArrayT second_face_velocities; + first_mesh.getFaceCoords( first_face_id, first_face_coordinates ); + second_mesh.getFaceCoords( second_face_id, second_face_coordinates ); + common_plane.getFace1ProjectedCoords( first_projected_face_coordinates, first_mesh.numberOfNodesPerElement() ); + common_plane.getFace2ProjectedCoords( second_projected_face_coordinates, + second_mesh.numberOfNodesPerElement() ); + if ( evaluate_velocity && !use_parent_fields ) { + first_mesh.getFaceVelocities( first_face_id, first_face_velocities ); + second_mesh.getFaceVelocities( second_face_id, second_face_velocities ); + } + + RealT overlap_vertices[max_dim * CommonPlaneContactData::maximum_overlap_vertices] = { 0.0 }; + common_plane.getOverlapVertices( overlap_vertices ); + const int number_of_overlap_vertices = spatial_dimension == 2 ? 2 : common_plane.m_numPolyVert; + const bool use_multiple_points = UseMultipleIntegrationPoints( penalty_options, first_mesh, second_mesh ); + const int quadrature_order = + GetCommonPlaneQuadratureOrder( penalty_options, first_mesh, second_mesh, first_face_id, second_face_id ); + int generated_row_count = 0; + bool has_active_row = false; + bool mapping_failed = false; + RealT integrated_measure = 0.0; + + if ( !use_multiple_points ) { + const RealT integration_point[max_dim] = { common_plane.m_cX, common_plane.m_cY, + spatial_dimension == 3 ? common_plane.m_cZ : 0.0 }; + const int row_evaluation_status = StoreCommonPlaneContactRow( + rows, first_row_id, contact_pair_id, first_mesh, second_mesh, first_face_id, second_face_id, + first_face_coordinates, second_face_coordinates, first_face_velocities, second_face_velocities, + first_projected_face_coordinates, second_projected_face_coordinates, integration_point, + common_plane_normal, common_plane.m_area, penalty_stiffness, rate_penalty_coefficient, + tangential_viscous_coefficient, gap_tolerance, evaluate_velocity, use_parent_fields ); + mapping_failed = row_evaluation_status != 0; + if ( row_evaluation_status == 0 ) { + generated_row_count = 1; + integrated_measure = common_plane.m_area; + has_active_row = rows.row_is_active[first_row_id] != 0; + } + } else if ( spatial_dimension == 2 ) { + RealT quadrature_weights[max_segment_gauss_legendre_qpts] = { 0.0 }; + RealT quadrature_coordinates[max_segment_gauss_legendre_qpts] = { 0.0 }; + const int number_of_quadrature_points = + GetCommonPlaneSegmentRule( quadrature_order, quadrature_weights, quadrature_coordinates ); + const RealT first_endpoint_x = overlap_vertices[0]; + const RealT first_endpoint_y = overlap_vertices[1]; + const RealT second_endpoint_x = overlap_vertices[2]; + const RealT second_endpoint_y = overlap_vertices[3]; + const RealT overlap_length = + magnitude( second_endpoint_x - first_endpoint_x, second_endpoint_y - first_endpoint_y ); + + for ( int quadrature_point = 0; quadrature_point < number_of_quadrature_points; ++quadrature_point ) { + const RealT segment_coordinate = quadrature_coordinates[quadrature_point]; + const RealT integration_point[max_dim] = { + ( 1.0 - segment_coordinate ) * first_endpoint_x + segment_coordinate * second_endpoint_x, + ( 1.0 - segment_coordinate ) * first_endpoint_y + segment_coordinate * second_endpoint_y, 0.0 }; + const RealT integration_weight = overlap_length * quadrature_weights[quadrature_point]; + const IndexT row_id = first_row_id + generated_row_count; + const int row_evaluation_status = StoreCommonPlaneContactRow( + rows, row_id, contact_pair_id, first_mesh, second_mesh, first_face_id, second_face_id, + first_face_coordinates, second_face_coordinates, first_face_velocities, second_face_velocities, + first_projected_face_coordinates, second_projected_face_coordinates, integration_point, + common_plane_normal, integration_weight, penalty_stiffness, rate_penalty_coefficient, + tangential_viscous_coefficient, gap_tolerance, evaluate_velocity, use_parent_fields ); + if ( row_evaluation_status != 0 ) { + mapping_failed = true; + break; + } + integrated_measure += integration_weight; + has_active_row = has_active_row || rows.row_is_active[row_id] != 0; + ++generated_row_count; + } + } else { + RealT quadrature_weights[max_symmetric_triangle_qpts] = { 0.0 }; + RealT quadrature_coordinates[2 * max_symmetric_triangle_qpts] = { 0.0 }; + const int number_of_quadrature_points = + GetCommonPlaneTriangleRule( quadrature_order, quadrature_weights, quadrature_coordinates ); + const RealT overlap_centroid[max_dim] = { common_plane.m_cX, common_plane.m_cY, common_plane.m_cZ }; + + // CommonPlane overlap polygons are convex and consistently ordered. + // A fan about the polygon centroid therefore creates nonoverlapping + // integration triangles while preserving the LOR overlap measure. + for ( int overlap_vertex = 0; overlap_vertex < number_of_overlap_vertices; ++overlap_vertex ) { + const int next_overlap_vertex = overlap_vertex + 1 == number_of_overlap_vertices ? 0 : overlap_vertex + 1; + RealT triangle_x[3] = { overlap_vertices[spatial_dimension * overlap_vertex], + overlap_vertices[spatial_dimension * next_overlap_vertex], overlap_centroid[0] }; + RealT triangle_y[3] = { overlap_vertices[spatial_dimension * overlap_vertex + 1], + overlap_vertices[spatial_dimension * next_overlap_vertex + 1], + overlap_centroid[1] }; + RealT triangle_z[3] = { overlap_vertices[spatial_dimension * overlap_vertex + 2], + overlap_vertices[spatial_dimension * next_overlap_vertex + 2], + overlap_centroid[2] }; + const RealT triangle_area = Area3DTri( triangle_x, triangle_y, triangle_z ); + if ( triangle_area <= 0.0 ) { + continue; + } + + for ( int quadrature_point = 0; quadrature_point < number_of_quadrature_points; ++quadrature_point ) { + const RealT first_triangle_coordinate = quadrature_coordinates[2 * quadrature_point]; + const RealT second_triangle_coordinate = quadrature_coordinates[2 * quadrature_point + 1]; + const RealT centroid_coordinate = 1.0 - first_triangle_coordinate - second_triangle_coordinate; + const RealT integration_point[max_dim] = { + centroid_coordinate * triangle_x[0] + first_triangle_coordinate * triangle_x[1] + + second_triangle_coordinate * triangle_x[2], + centroid_coordinate * triangle_y[0] + first_triangle_coordinate * triangle_y[1] + + second_triangle_coordinate * triangle_y[2], + centroid_coordinate * triangle_z[0] + first_triangle_coordinate * triangle_z[1] + + second_triangle_coordinate * triangle_z[2] }; + const RealT integration_weight = triangle_area * quadrature_weights[quadrature_point]; + const IndexT row_id = first_row_id + generated_row_count; + const int row_evaluation_status = StoreCommonPlaneContactRow( + rows, row_id, contact_pair_id, first_mesh, second_mesh, first_face_id, second_face_id, + first_face_coordinates, second_face_coordinates, first_face_velocities, second_face_velocities, + first_projected_face_coordinates, second_projected_face_coordinates, integration_point, + common_plane_normal, integration_weight, penalty_stiffness, rate_penalty_coefficient, + tangential_viscous_coefficient, gap_tolerance, evaluate_velocity, use_parent_fields ); + if ( row_evaluation_status != 0 ) { + mapping_failed = true; + break; + } + integrated_measure += integration_weight; + has_active_row = has_active_row || rows.row_is_active[row_id] != 0; + ++generated_row_count; + } + if ( mapping_failed ) { + break; + } + } + } - // 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 ); + if ( mapping_failed ) { + for ( int generated_row = 0; generated_row < generated_row_count; ++generated_row ) { + rows.row_is_valid[first_row_id + generated_row] = 0; + rows.row_is_active[first_row_id + generated_row] = 0; + } + rows.pair_evaluation_statuses[contact_pair_id] = + static_cast( CommonPlanePairEvaluationStatus::INVALID_PARENT_MAPPING ); + tribol::atomicMax( &evaluation_error[0], 1 ); + return; + } + if ( integrated_measure <= 0.0 || generated_row_count == 0 ) { + rows.pair_evaluation_statuses[contact_pair_id] = + static_cast( CommonPlanePairEvaluationStatus::DEGENERATE_OVERLAP ); + tribol::atomicMax( &evaluation_error[0], 1 ); + return; + } + + rows.pair_row_counts[contact_pair_id] = generated_row_count; + rows.pair_evaluation_statuses[contact_pair_id] = static_cast( CommonPlanePairEvaluationStatus::VALID ); + common_plane.m_inContact = has_active_row; + } ); + + Array1D evaluation_error_host( evaluation_error_data ); + if ( evaluation_error_host[0] != 0 ) { + Array1D pair_statuses_host( common_plane_data->getPairEvaluationStatuses() ); + for ( IndexT contact_pair_id = 0; contact_pair_id < number_of_pairs; ++contact_pair_id ) { + if ( pair_statuses_host[contact_pair_id] != static_cast( CommonPlanePairEvaluationStatus::VALID ) ) { + SLIC_DEBUG( "BuildCommonPlaneContactRows(): row generation failed for active pair " + << contact_pair_id << " with status " << pair_statuses_host[contact_pair_id] << "." ); } - } // end for loop over face nodes + } + } + return evaluation_error_host[0]; +} + +//------------------------------------------------------------------------------ +template <> +int ApplyNormal( CouplingScheme* cs ) +{ + auto* common_plane_data = static_cast( cs->getMethodData() ); + SLIC_ERROR_ROOT_IF( common_plane_data == nullptr, + "ApplyNormal(): CommonPlane row storage is unavailable." ); + + const int row_build_error = BuildCommonPlaneContactRows( cs ); + if ( row_build_error != 0 ) { + return row_build_error; + } + + const CommonPlaneContactData::Viewer rows = common_plane_data->getView(); + const CouplingScheme::Viewer coupling_scheme = cs->getView(); + const RealT residual_gap = coupling_scheme.getParameters().residual_gap; + const PenaltyConstraintType constraint_type = coupling_scheme.getEnforcementOptions().penalty_options.constraint_type; + + // These per-pair reductions preserve the existing contact-plane diagnostics, + // but derive them from the same pointwise values used by the force scatter. + Array1D integrated_kinematic_pressure( rows.number_of_pairs, rows.number_of_pairs, cs->getAllocatorId() ); + Array1D integrated_rate_pressure( rows.number_of_pairs, rows.number_of_pairs, cs->getAllocatorId() ); + Array1D integrated_normal_velocity( rows.number_of_pairs, rows.number_of_pairs, cs->getAllocatorId() ); + Array1D active_measure( rows.number_of_pairs, rows.number_of_pairs, cs->getAllocatorId() ); + integrated_kinematic_pressure.fill( 0.0 ); + integrated_rate_pressure.fill( 0.0 ); + integrated_normal_velocity.fill( 0.0 ); + active_measure.fill( 0.0 ); + + const Array1DView integrated_kinematic_pressure_view = integrated_kinematic_pressure.view(); + const Array1DView integrated_rate_pressure_view = integrated_rate_pressure.view(); + const Array1DView integrated_normal_velocity_view = integrated_normal_velocity.view(); + const Array1DView active_measure_view = active_measure.view(); + + forAllExec( cs->getExecutionMode(), rows.row_capacity, + [rows, coupling_scheme, residual_gap, constraint_type, integrated_kinematic_pressure_view, + integrated_rate_pressure_view, integrated_normal_velocity_view, + active_measure_view] TRIBOL_HOST_DEVICE( IndexT row_id ) { + if ( rows.row_is_valid[row_id] == 0 || rows.row_is_active[row_id] == 0 ) { + return; + } + + const RealT kinematic_pressure = + ( rows.gaps[row_id] - residual_gap ) * rows.penalty_stiffnesses[row_id]; + RealT rate_pressure = 0.0; + if ( constraint_type == KINEMATIC_AND_RATE && rows.normal_velocity_gaps[row_id] <= 0.0 ) { + rate_pressure = rows.normal_velocity_gaps[row_id] * rows.rate_penalty_coefficients[row_id]; + } + const RealT applied_pressure = kinematic_pressure + rate_pressure; + const RealT weighted_pressure = rows.integration_weights[row_id] * applied_pressure; + RealT force_on_second_face[max_dim] = { 0.0, 0.0, 0.0 }; + for ( int component = 0; component < rows.spatial_dimension; ++component ) { + force_on_second_face[component] = rows.normals( row_id, component ) * weighted_pressure; + } + + ScatterCommonPlaneRowForce( rows, row_id, coupling_scheme.getMesh1View(), + coupling_scheme.getMesh2View(), force_on_second_face ); + + const IndexT contact_pair_id = rows.contact_pair_ids[row_id]; + tribol::atomicAdd( &integrated_kinematic_pressure_view[contact_pair_id], + rows.integration_weights[row_id] * kinematic_pressure ); + tribol::atomicAdd( &integrated_rate_pressure_view[contact_pair_id], + rows.integration_weights[row_id] * rate_pressure ); + tribol::atomicAdd( &integrated_normal_velocity_view[contact_pair_id], + rows.integration_weights[row_id] * rows.normal_velocity_gaps[row_id] ); + tribol::atomicAdd( &active_measure_view[contact_pair_id], rows.integration_weights[row_id] ); + } ); + + forAllExec( cs->getExecutionMode(), rows.number_of_pairs, + [coupling_scheme, integrated_kinematic_pressure_view, integrated_rate_pressure_view, + integrated_normal_velocity_view, active_measure_view] TRIBOL_HOST_DEVICE( IndexT contact_pair_id ) { + CommonPlanePair& common_plane = coupling_scheme.getCompGeomView().getCommonPlane( contact_pair_id ); + const RealT measure = active_measure_view[contact_pair_id]; + if ( measure > 0.0 ) { + common_plane.m_pressure = integrated_kinematic_pressure_view[contact_pair_id] / measure; + common_plane.m_ratePressure = integrated_rate_pressure_view[contact_pair_id] / measure; + common_plane.m_velGap = integrated_normal_velocity_view[contact_pair_id] / measure; + } else { + common_plane.m_pressure = 0.0; + common_plane.m_ratePressure = 0.0; + common_plane.m_velGap = 0.0; + } + } ); + + return 0; +} + +//------------------------------------------------------------------------------ +template <> +int ApplyTangential( CouplingScheme* cs ) +{ + auto* common_plane_data = static_cast( cs->getMethodData() ); + SLIC_ERROR_ROOT_IF( common_plane_data == nullptr, + "ApplyTangential(): CommonPlane row storage is " + "unavailable." ); + + const CommonPlaneContactData::Viewer rows = common_plane_data->getView(); + const CouplingScheme::Viewer coupling_scheme = cs->getView(); + forAllExec( cs->getExecutionMode(), rows.row_capacity, [rows, coupling_scheme] TRIBOL_HOST_DEVICE( IndexT row_id ) { + if ( rows.row_is_valid[row_id] == 0 || rows.row_is_active[row_id] == 0 ) { + return; + } + + RealT tangential_velocity[max_dim] = { 0.0, 0.0, 0.0 }; + for ( int component = 0; component < rows.spatial_dimension; ++component ) { + const RealT relative_velocity = + rows.first_velocities( row_id, component ) - rows.second_velocities( row_id, component ); + tangential_velocity[component] = + relative_velocity - rows.normal_velocity_gaps[row_id] * rows.normals( row_id, component ); + } + + const RealT force_scale = rows.integration_weights[row_id] * rows.tangential_viscous_coefficients[row_id]; + RealT force_on_second_face[max_dim] = { force_scale * tangential_velocity[0], force_scale * tangential_velocity[1], + force_scale * tangential_velocity[2] }; + ScatterCommonPlaneRowForce( rows, row_id, coupling_scheme.getMesh1View(), coupling_scheme.getMesh2View(), + force_on_second_face ); } ); return 0; diff --git a/src/tribol/physics/CommonPlane.hpp b/src/tribol/physics/CommonPlane.hpp index a7b0ef8c..ecc5efe2 100644 --- a/src/tribol/physics/CommonPlane.hpp +++ b/src/tribol/physics/CommonPlane.hpp @@ -43,6 +43,19 @@ TRIBOL_HOST_DEVICE inline RealT ComputePenaltyStiffnessPerArea( const RealT K1_o } // end ComputePenaltyStiffnessPerArea +/** + * @brief Build the shared CommonPlane quadrature rows for the current overlap cells. + * + * The generated rows contain mapped parent-face coordinates, field values, + * basis values, integration weights, normals, gaps, and penalty coefficients. + * Explicit force, damping, diagnostics, and downstream operators consume this + * same batch. + * + * @param [in,out] cs CommonPlane coupling scheme that owns the row batch + * @return zero on success and nonzero if row generation fails + */ +int BuildCommonPlaneContactRows( CouplingScheme* cs ); + /*! * * \brief routine to apply interface physics in the direction normal to the interface diff --git a/src/tribol/physics/Mortar.cpp b/src/tribol/physics/Mortar.cpp index 1d99580e..426c63c5 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..7fd3978d 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; ///< CommonPlane overlap integration rule used by the test. + int common_plane_quadrature_order; ///< Explicit CommonPlane multipoint integration order used by the test. RealT rate_penalty; RealT rate_penalty_ratio; RealT const_penalty;