Add higher-order CommonPlane quadrature rows

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.

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<RealT>( 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;

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;
     } );
 

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 };