From 7492e7c4d3a374d4fabe7c8d69e13b9d9f755b20 Mon Sep 17 00:00:00 2001 From: EB Chin Date: Thu, 20 Aug 2026 15:30:07 -0700 Subject: [PATCH 1/7] fix a leak in enzyme AD code with stack memory --- src/tribol/physics/EnergyMortar.cpp | 188 +++++++++++++++------------- 1 file changed, 100 insertions(+), 88 deletions(-) diff --git a/src/tribol/physics/EnergyMortar.cpp b/src/tribol/physics/EnergyMortar.cpp index 52b4a0cd..6bbe82e3 100644 --- a/src/tribol/physics/EnergyMortar.cpp +++ b/src/tribol/physics/EnergyMortar.cpp @@ -44,9 +44,8 @@ TRIBOL_ENZYME_INLINE void find_normal( const double* coord1, const double* coord } // Gets the respective gauss-legendre nodes dependant on quadrature order -TRIBOL_ENZYME_INLINE void determine_legendre_nodes( int N, std::array& x ) +TRIBOL_ENZYME_INLINE void determine_legendre_nodes_raw( int N, double* x ) { - // x.resize( N ); if ( N == 1 ) { x[0] = 0.0; } else if ( N == 2 ) { @@ -79,9 +78,8 @@ TRIBOL_ENZYME_INLINE void determine_legendre_nodes( int N, std::array } // Gets the respective gauss-legendre weights dependant on quadrature order -TRIBOL_ENZYME_INLINE void determine_legendre_weights( int N, std::array& W ) +TRIBOL_ENZYME_INLINE void determine_legendre_weights_raw( int N, double* W ) { - // W.resize( N ); if ( N == 1 ) { W[0] = 2.0; } else if ( N == 2 ) { @@ -107,6 +105,24 @@ TRIBOL_ENZYME_INLINE void determine_legendre_weights( int N, std::arrayqp[i] = 0.5 * ( xi_max - xi_min ) * qpoints[i] + 0.5 * ( xi_max + xi_min ); + out->w[i] = weights[i] * J; + } +} + // Map a point from the 1D parent segment coordinate to physical coordinates. // Parametric space: [-0.5, 0.5] TRIBOL_ENZYME_INLINE void iso_map( const double* coord1, const double* coord2, double xi, double* mapped_coord ) @@ -191,21 +207,16 @@ TRIBOL_ENZYME_INLINE void get_projections( const double* A0, const double* A1, c const double dyA = A1[1] - A0[1]; const double len2A = dxA * dxA + dyA * dyA; - const double* B_endpoints[2] = { B0, B1 }; + double q0[2] = { 0.0, 0.0 }; + find_intersection( A0, A1, B0, nB, q0 ); + // Convert the physical projection point on A to the local coordinate xi in [-0.5, 0.5]. + const double alphaA0 = ( ( q0[0] - A0[0] ) * dxA + ( q0[1] - A0[1] ) * dyA ) / len2A; + const double xi0 = alphaA0 - 0.5; - double xi0 = 0.0, xi1 = 0.0; - for ( int i = 0; i < 2; ++i ) { - double q[2] = { 0.0, 0.0 }; - find_intersection( A0, A1, B_endpoints[i], nB, q ); - // Convert the physical projection point on A to the local coordinate xi in [-0.5, 0.5]. - const double alphaA = ( ( q[0] - A0[0] ) * dxA + ( q[1] - A0[1] ) * dyA ) / len2A; - const double xiA = alphaA - 0.5; - - if ( i == 0 ) - xi0 = xiA; - else - xi1 = xiA; - } + double q1[2] = { 0.0, 0.0 }; + find_intersection( A0, A1, B1, nB, q1 ); + const double alphaA1 = ( ( q1[0] - A0[0] ) * dxA + ( q1[1] - A0[1] ) * dyA ) / len2A; + const double xi1 = alphaA1 - 0.5; double xi_min = std::min( xi0, xi1 ); double xi_max = std::max( xi0, xi1 ); @@ -214,6 +225,57 @@ TRIBOL_ENZYME_INLINE void get_projections( const double* A0, const double* A1, c projections[1] = xi_max; } +TRIBOL_ENZYME_INLINE void bounds_from_projections_raw( const double* projections, double del, double* bounds ) +{ + double xi_min = ( projections[0] < projections[1] ) ? projections[0] : projections[1]; + double xi_max = ( projections[0] > projections[1] ) ? projections[0] : projections[1]; + + if ( xi_max < -0.5 - del ) { + xi_max = -0.5 - del; + } + if ( xi_min > 0.5 + del ) { + xi_min = 0.5 + del; + } + if ( xi_min < -0.5 - del ) { + xi_min = -0.5 - del; + } + if ( xi_max > 0.5 + del ) { + xi_max = 0.5 + del; + } + + bounds[0] = xi_min; + bounds[1] = xi_max; +} + +TRIBOL_ENZYME_INLINE double smooth_bound_value_raw( double bound, double del ) +{ + const double xi = bound + 0.5; + double xi_hat = 0.0; + + if ( del == 0.0 ) { + xi_hat = xi; + } else if ( 0.0 - del <= xi && xi <= del ) { + xi_hat = ( 1.0 / ( 4.0 * del ) ) * ( xi * xi ) + 0.5 * xi + del / 4.0; + } else if ( ( 1.0 - del ) <= xi && xi <= 1.0 + del ) { + const double b = -1.0 / ( 4.0 * del ); + const double c = 0.5 + 1.0 / ( 2.0 * del ); + const double one_minus_del = 1.0 - del; + const double d = 1.0 - del + ( 1.0 / ( 4.0 * del ) ) * one_minus_del * one_minus_del - 0.5 * one_minus_del - + one_minus_del / ( 2.0 * del ); + xi_hat = b * xi * xi + c * xi + d; + } else if ( del <= xi && xi <= ( 1.0 - del ) ) { + xi_hat = xi; + } + + return xi_hat - 0.5; +} + +TRIBOL_ENZYME_INLINE void smooth_bounds_raw( const double* bounds, double del, double* smooth_bounds ) +{ + smooth_bounds[0] = smooth_bound_value_raw( bounds[0], del ); + smooth_bounds[1] = smooth_bound_value_raw( bounds[1], del ); +} + // Integrate the nodal smoothed gap and tributary area contributions over edge A. // The quadrature rule is supplied through gp, allowing this kernel to be reused // for both fixed-quadrature and geometry-dependent quadrature paths. @@ -412,13 +474,14 @@ static void kernel_out_enzyme( const double* x, const void* kp_void, double* out double projs[2] = { 0 }; get_projections( A0, A1, B0, B1, projs ); - std::array projections = { projs[0], projs[1] }; // Recompute the integration bounds and quadrature from the current geometry. - auto bounds = ContactSmoothing::bounds_from_projections( projections, kp->del ); - auto xi_bounds = ContactSmoothing::smooth_bounds( bounds, kp->del ); - - auto qp = EnergyMortarCalculator::compute_quadrature( xi_bounds, kp->N ); + double bounds[2]; + bounds_from_projections_raw( projs, kp->del, bounds ); + double xi_bounds[2]; + smooth_bounds_raw( bounds, kp->del, xi_bounds ); + QuadPoints qp; + compute_quadrature_raw( xi_bounds, kp->N, &qp ); Gparams gp; for ( std::size_t i = 0; i < qp.qp.size(); ++i ) { @@ -486,10 +549,12 @@ TRIBOL_ENZYME_INLINE void qp_penalty_kernel( const double* x, const KernelParams double projs[2] = { 0.0, 0.0 }; get_projections( A0, A1, B0, B1, projs ); - const std::array projections{ projs[0], projs[1] }; - const auto bounds = ContactSmoothing::bounds_from_projections( projections, kp->del ); - const auto xi_bounds = ContactSmoothing::smooth_bounds( bounds, kp->del ); - const auto qp = EnergyMortarCalculator::compute_quadrature( xi_bounds, kp->N ); + double bounds[2]; + bounds_from_projections_raw( projs, kp->del, bounds ); + double xi_bounds[2]; + smooth_bounds_raw( bounds, kp->del, xi_bounds ); + QuadPoints qp; + compute_quadrature_raw( xi_bounds, kp->N, &qp ); double nB[2]; find_normal( B0, B1, nB ); @@ -634,24 +699,10 @@ std::array EnergyMortarCalculator::projections( const InterfacePair& TRIBOL_ENZYME_INLINE std::array ContactSmoothing::bounds_from_projections( const std::array& proj, double del ) { - double xi_min = std::min( proj[0], proj[1] ); - double xi_max = std::max( proj[0], proj[1] ); - - // Limit the integration interval to the extended range [-0.5 - del, 0.5 + del]. - if ( xi_max < -0.5 - del ) { - xi_max = -0.5 - del; - } - if ( xi_min > 0.5 + del ) { - xi_min = 0.5 + del; - } - if ( xi_min < -0.5 - del ) { - xi_min = -0.5 - del; - } - if ( xi_max > 0.5 + del ) { - xi_max = 0.5 + del; - } - - return { xi_min, xi_max }; + const double projections[2] = { proj[0], proj[1] }; + double bounds[2]; + bounds_from_projections_raw( projections, del, bounds ); + return { bounds[0], bounds[1] }; } // Smooth the integration bounds using a C1 ramp near the ends of edge A. @@ -663,32 +714,8 @@ TRIBOL_ENZYME_INLINE std::array ContactSmoothing::smooth_bounds( cons double del ) { std::array smooth_bounds; - for ( int i = 0; i < 2; ++i ) { - double xi = 0.0; - double xi_hat = 0.0; - - // Shift from the local coordinate interval [-0.5, 0.5] to [0, 1]. - xi = bounds[i] + 0.5; - if ( del == 0.0 ) { - xi_hat = xi; - } else { - // Apply quadratic ramps near the endpoints and leave the interior unchanged. - if ( 0.0 - del <= xi && xi <= del ) { - xi_hat = ( 1.0 / ( 4 * del ) ) * ( xi * xi ) + 0.5 * xi + del / 4.0; - } else if ( ( 1.0 - del ) <= xi && xi <= 1.0 + del ) { - double b = -1.0 / ( 4.0 * del ); - double c = 0.5 + 1.0 / ( 2.0 * del ); - double d = 1.0 - del + ( 1.0 / ( 4.0 * del ) ) * pow( 1.0 - del, 2 ) - 0.5 * ( 1.0 - del ) - - ( 1.0 - del ) / ( 2.0 * del ); - - xi_hat = b * xi * xi + c * xi + d; - } else if ( del <= xi && xi <= ( 1.0 - del ) ) { - xi_hat = xi; - } - } - // Shift the smoothed coordinate back to [-0.5, 0.5]. - smooth_bounds[i] = xi_hat - 0.5; - } + const double bounds_raw[2] = { bounds[0], bounds[1] }; + smooth_bounds_raw( bounds_raw, del, smooth_bounds.data() ); return smooth_bounds; } @@ -698,23 +725,8 @@ TRIBOL_ENZYME_INLINE QuadPoints EnergyMortarCalculator::compute_quadrature( cons int N ) { QuadPoints out; - - std::array qpoints; - std::array weights; - - determine_legendre_nodes( N, qpoints ); - determine_legendre_weights( N, weights ); - - const double xi_min = xi_bounds[0]; - const double xi_max = xi_bounds[1]; - // Map the reference quadrature rule to [xi_min, xi_max]. - const double J = 0.5 * ( xi_max - xi_min ); - - for ( int i = 0; i < N; ++i ) { - out.qp[i] = 0.5 * ( xi_max - xi_min ) * qpoints[i] + 0.5 * ( xi_max + xi_min ); - out.w[i] = weights[i] * J; - } - + const double bounds[2] = { xi_bounds[0], xi_bounds[1] }; + compute_quadrature_raw( bounds, N, &out ); return out; } From 2a67f80aba15d2c6a256016adb2cfb1a3801599d Mon Sep 17 00:00:00 2001 From: EB Chin Date: Thu, 20 Aug 2026 16:03:05 -0700 Subject: [PATCH 2/7] fix another leak --- src/tribol/physics/EnergyMortar.cpp | 49 ++++++++++++++++++++--------- 1 file changed, 34 insertions(+), 15 deletions(-) diff --git a/src/tribol/physics/EnergyMortar.cpp b/src/tribol/physics/EnergyMortar.cpp index 6bbe82e3..ee2dca27 100644 --- a/src/tribol/physics/EnergyMortar.cpp +++ b/src/tribol/physics/EnergyMortar.cpp @@ -540,6 +540,37 @@ void d2_kernel( const double* x, const KernelParams* kp, double* H ) } } +TRIBOL_ENZYME_INLINE double qp_penalty_kernel_qp_energy( double xiA, double w, const double* A0, + const double* A1, const double* B0, const double* B1, + const double* nB, double eta, double penalty, double J ) +{ + const double N1 = 0.5 - xiA; + const double N2 = 0.5 + xiA; + const double x1x = N1 * A0[0] + N2 * A1[0]; + const double x1y = N1 * A0[1] + N2 * A1[1]; + + const double tBx = B1[0] - B0[0]; + const double tBy = B1[1] - B0[1]; + const double dxB = x1x - B0[0]; + const double dyB = x1y - B0[1]; + const double det = tBx * nB[1] - tBy * nB[0]; + + double x2x = x1x; + double x2y = x1y; + if ( std::abs( det ) >= 1e-12 ) { + const double alpha = ( dxB * nB[1] - dyB * nB[0] ) / det; + x2x = B0[0] + alpha * tBx; + x2y = B0[1] + alpha * tBy; + } + + const double dx = x1x - x2x; + const double dy = x1y - x2y; + const double gn = -( dx * nB[0] + dy * nB[1] ); + const double gap = gn * eta; + + return gap < 0.0 ? 0.5 * penalty * gap * gap * w * J : 0.0; +} + TRIBOL_ENZYME_INLINE void qp_penalty_kernel( const double* x, const KernelParams* kp, double* energy ) { double A0[2] = { x[0], x[1] }; @@ -568,21 +599,9 @@ TRIBOL_ENZYME_INLINE void qp_penalty_kernel( const double* x, const KernelParams const double J = line_jacobian( A0, A1 ); double value = 0.0; - for ( int i = 0; i < kp->N; ++i ) { - double x1[2]; - iso_map( A0, A1, qp.qp[i], x1 ); - - double x2[2]; - find_intersection( B0, B1, x1, nB, x2 ); - - const double dx = x1[0] - x2[0]; - const double dy = x1[1] - x2[1]; - const double gn = -( dx * nB[0] + dy * nB[1] ); - const double gap = gn * eta; - if ( gap < 0.0 ) { - value += 0.5 * kp->k * gap * gap * qp.w[i] * J; - } - } + value += qp_penalty_kernel_qp_energy( qp.qp[0], qp.w[0], A0, A1, B0, B1, nB, eta, kp->k, J ); + value += qp_penalty_kernel_qp_energy( qp.qp[1], qp.w[1], A0, A1, B0, B1, nB, eta, kp->k, J ); + value += qp_penalty_kernel_qp_energy( qp.qp[2], qp.w[2], A0, A1, B0, B1, nB, eta, kp->k, J ); *energy = value; } From d498f476ca71a2c894e84e9818115d5beb88bcc4 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Thu, 20 Aug 2026 16:13:07 -0700 Subject: [PATCH 3/7] formatting --- src/tribol/physics/EnergyMortar.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/tribol/physics/EnergyMortar.cpp b/src/tribol/physics/EnergyMortar.cpp index ee2dca27..69b0e557 100644 --- a/src/tribol/physics/EnergyMortar.cpp +++ b/src/tribol/physics/EnergyMortar.cpp @@ -540,9 +540,9 @@ void d2_kernel( const double* x, const KernelParams* kp, double* H ) } } -TRIBOL_ENZYME_INLINE double qp_penalty_kernel_qp_energy( double xiA, double w, const double* A0, - const double* A1, const double* B0, const double* B1, - const double* nB, double eta, double penalty, double J ) +TRIBOL_ENZYME_INLINE double qp_penalty_kernel_qp_energy( double xiA, double w, const double* A0, const double* A1, + const double* B0, const double* B1, const double* nB, + double eta, double penalty, double J ) { const double N1 = 0.5 - xiA; const double N2 = 0.5 + xiA; From 71f9a54bd7139eb8edbe5a8f0c021ba95ee81cc3 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Fri, 21 Aug 2026 10:47:04 -0700 Subject: [PATCH 4/7] fix naming, make delta from develop smaller --- src/tribol/physics/EnergyMortar.cpp | 133 ++++++++++++++++++---------- 1 file changed, 87 insertions(+), 46 deletions(-) diff --git a/src/tribol/physics/EnergyMortar.cpp b/src/tribol/physics/EnergyMortar.cpp index 69b0e557..92e17130 100644 --- a/src/tribol/physics/EnergyMortar.cpp +++ b/src/tribol/physics/EnergyMortar.cpp @@ -44,7 +44,7 @@ TRIBOL_ENZYME_INLINE void find_normal( const double* coord1, const double* coord } // Gets the respective gauss-legendre nodes dependant on quadrature order -TRIBOL_ENZYME_INLINE void determine_legendre_nodes_raw( int N, double* x ) +TRIBOL_ENZYME_INLINE void determine_legendre_nodes( int N, double* x ) { if ( N == 1 ) { x[0] = 0.0; @@ -78,7 +78,7 @@ TRIBOL_ENZYME_INLINE void determine_legendre_nodes_raw( int N, double* x ) } // Gets the respective gauss-legendre weights dependant on quadrature order -TRIBOL_ENZYME_INLINE void determine_legendre_weights_raw( int N, double* W ) +TRIBOL_ENZYME_INLINE void determine_legendre_weights( int N, double* W ) { if ( N == 1 ) { W[0] = 2.0; @@ -105,13 +105,13 @@ TRIBOL_ENZYME_INLINE void determine_legendre_weights_raw( int N, double* W ) } } -TRIBOL_ENZYME_INLINE void compute_quadrature_raw( const double* xi_bounds, int N, QuadPoints* out ) +TRIBOL_ENZYME_INLINE void compute_quadrature( const double* xi_bounds, int N, QuadPoints* out ) { double qpoints[3] = { 0.0, 0.0, 0.0 }; double weights[3] = { 0.0, 0.0, 0.0 }; - determine_legendre_nodes_raw( N, qpoints ); - determine_legendre_weights_raw( N, weights ); + determine_legendre_nodes( N, qpoints ); + determine_legendre_weights( N, weights ); const double xi_min = xi_bounds[0]; const double xi_max = xi_bounds[1]; @@ -225,7 +225,7 @@ TRIBOL_ENZYME_INLINE void get_projections( const double* A0, const double* A1, c projections[1] = xi_max; } -TRIBOL_ENZYME_INLINE void bounds_from_projections_raw( const double* projections, double del, double* bounds ) +TRIBOL_ENZYME_INLINE void bounds_from_projections( const double* projections, double del, double* bounds ) { double xi_min = ( projections[0] < projections[1] ) ? projections[0] : projections[1]; double xi_max = ( projections[0] > projections[1] ) ? projections[0] : projections[1]; @@ -247,7 +247,7 @@ TRIBOL_ENZYME_INLINE void bounds_from_projections_raw( const double* projections bounds[1] = xi_max; } -TRIBOL_ENZYME_INLINE double smooth_bound_value_raw( double bound, double del ) +TRIBOL_ENZYME_INLINE double smooth_bound( double bound, double del ) { const double xi = bound + 0.5; double xi_hat = 0.0; @@ -270,10 +270,10 @@ TRIBOL_ENZYME_INLINE double smooth_bound_value_raw( double bound, double del ) return xi_hat - 0.5; } -TRIBOL_ENZYME_INLINE void smooth_bounds_raw( const double* bounds, double del, double* smooth_bounds ) +TRIBOL_ENZYME_INLINE void smooth_bounds( const double* bounds, double del, double* smooth_bounds ) { - smooth_bounds[0] = smooth_bound_value_raw( bounds[0], del ); - smooth_bounds[1] = smooth_bound_value_raw( bounds[1], del ); + smooth_bounds[0] = smooth_bound( bounds[0], del ); + smooth_bounds[1] = smooth_bound( bounds[1], del ); } // Integrate the nodal smoothed gap and tributary area contributions over edge A. @@ -477,11 +477,11 @@ static void kernel_out_enzyme( const double* x, const void* kp_void, double* out // Recompute the integration bounds and quadrature from the current geometry. double bounds[2]; - bounds_from_projections_raw( projs, kp->del, bounds ); + bounds_from_projections( projs, kp->del, bounds ); double xi_bounds[2]; - smooth_bounds_raw( bounds, kp->del, xi_bounds ); + smooth_bounds( bounds, kp->del, xi_bounds ); QuadPoints qp; - compute_quadrature_raw( xi_bounds, kp->N, &qp ); + compute_quadrature( xi_bounds, kp->N, &qp ); Gparams gp; for ( std::size_t i = 0; i < qp.qp.size(); ++i ) { @@ -540,31 +540,19 @@ void d2_kernel( const double* x, const KernelParams* kp, double* H ) } } +// Isolate loop-local arrays to avoid a leak in Enzyme's reverse-mode tape. TRIBOL_ENZYME_INLINE double qp_penalty_kernel_qp_energy( double xiA, double w, const double* A0, const double* A1, const double* B0, const double* B1, const double* nB, double eta, double penalty, double J ) { - const double N1 = 0.5 - xiA; - const double N2 = 0.5 + xiA; - const double x1x = N1 * A0[0] + N2 * A1[0]; - const double x1y = N1 * A0[1] + N2 * A1[1]; - - const double tBx = B1[0] - B0[0]; - const double tBy = B1[1] - B0[1]; - const double dxB = x1x - B0[0]; - const double dyB = x1y - B0[1]; - const double det = tBx * nB[1] - tBy * nB[0]; + double x1[2]; + iso_map( A0, A1, xiA, x1 ); - double x2x = x1x; - double x2y = x1y; - if ( std::abs( det ) >= 1e-12 ) { - const double alpha = ( dxB * nB[1] - dyB * nB[0] ) / det; - x2x = B0[0] + alpha * tBx; - x2y = B0[1] + alpha * tBy; - } + double x2[2]; + find_intersection( B0, B1, x1, nB, x2 ); - const double dx = x1x - x2x; - const double dy = x1y - x2y; + const double dx = x1[0] - x2[0]; + const double dy = x1[1] - x2[1]; const double gn = -( dx * nB[0] + dy * nB[1] ); const double gap = gn * eta; @@ -581,11 +569,11 @@ TRIBOL_ENZYME_INLINE void qp_penalty_kernel( const double* x, const KernelParams double projs[2] = { 0.0, 0.0 }; get_projections( A0, A1, B0, B1, projs ); double bounds[2]; - bounds_from_projections_raw( projs, kp->del, bounds ); + bounds_from_projections( projs, kp->del, bounds ); double xi_bounds[2]; - smooth_bounds_raw( bounds, kp->del, xi_bounds ); + smooth_bounds( bounds, kp->del, xi_bounds ); QuadPoints qp; - compute_quadrature_raw( xi_bounds, kp->N, &qp ); + compute_quadrature( xi_bounds, kp->N, &qp ); double nB[2]; find_normal( B0, B1, nB ); @@ -599,9 +587,9 @@ TRIBOL_ENZYME_INLINE void qp_penalty_kernel( const double* x, const KernelParams const double J = line_jacobian( A0, A1 ); double value = 0.0; - value += qp_penalty_kernel_qp_energy( qp.qp[0], qp.w[0], A0, A1, B0, B1, nB, eta, kp->k, J ); - value += qp_penalty_kernel_qp_energy( qp.qp[1], qp.w[1], A0, A1, B0, B1, nB, eta, kp->k, J ); - value += qp_penalty_kernel_qp_energy( qp.qp[2], qp.w[2], A0, A1, B0, B1, nB, eta, kp->k, J ); + for ( int i = 0; i < kp->N; ++i ) { + value += qp_penalty_kernel_qp_energy( qp.qp[i], qp.w[i], A0, A1, B0, B1, nB, eta, kp->k, J ); + } *energy = value; } @@ -718,10 +706,24 @@ std::array EnergyMortarCalculator::projections( const InterfacePair& TRIBOL_ENZYME_INLINE std::array ContactSmoothing::bounds_from_projections( const std::array& proj, double del ) { - const double projections[2] = { proj[0], proj[1] }; - double bounds[2]; - bounds_from_projections_raw( projections, del, bounds ); - return { bounds[0], bounds[1] }; + double xi_min = std::min( proj[0], proj[1] ); + double xi_max = std::max( proj[0], proj[1] ); + + // Limit the integration interval to the extended range [-0.5 - del, 0.5 + del]. + if ( xi_max < -0.5 - del ) { + xi_max = -0.5 - del; + } + if ( xi_min > 0.5 + del ) { + xi_min = 0.5 + del; + } + if ( xi_min < -0.5 - del ) { + xi_min = -0.5 - del; + } + if ( xi_max > 0.5 + del ) { + xi_max = 0.5 + del; + } + + return { xi_min, xi_max }; } // Smooth the integration bounds using a C1 ramp near the ends of edge A. @@ -733,8 +735,32 @@ TRIBOL_ENZYME_INLINE std::array ContactSmoothing::smooth_bounds( cons double del ) { std::array smooth_bounds; - const double bounds_raw[2] = { bounds[0], bounds[1] }; - smooth_bounds_raw( bounds_raw, del, smooth_bounds.data() ); + for ( int i = 0; i < 2; ++i ) { + double xi = 0.0; + double xi_hat = 0.0; + + // Shift from the local coordinate interval [-0.5, 0.5] to [0, 1]. + xi = bounds[i] + 0.5; + if ( del == 0.0 ) { + xi_hat = xi; + } else { + // Apply quadratic ramps near the endpoints and leave the interior unchanged. + if ( 0.0 - del <= xi && xi <= del ) { + xi_hat = ( 1.0 / ( 4 * del ) ) * ( xi * xi ) + 0.5 * xi + del / 4.0; + } else if ( ( 1.0 - del ) <= xi && xi <= 1.0 + del ) { + double b = -1.0 / ( 4.0 * del ); + double c = 0.5 + 1.0 / ( 2.0 * del ); + double d = 1.0 - del + ( 1.0 / ( 4.0 * del ) ) * pow( 1.0 - del, 2 ) - 0.5 * ( 1.0 - del ) - + ( 1.0 - del ) / ( 2.0 * del ); + + xi_hat = b * xi * xi + c * xi + d; + } else if ( del <= xi && xi <= ( 1.0 - del ) ) { + xi_hat = xi; + } + } + // Shift the smoothed coordinate back to [-0.5, 0.5]. + smooth_bounds[i] = xi_hat - 0.5; + } return smooth_bounds; } @@ -744,8 +770,23 @@ TRIBOL_ENZYME_INLINE QuadPoints EnergyMortarCalculator::compute_quadrature( cons int N ) { QuadPoints out; - const double bounds[2] = { xi_bounds[0], xi_bounds[1] }; - compute_quadrature_raw( bounds, N, &out ); + + std::array qpoints; + std::array weights; + + determine_legendre_nodes( N, qpoints.data() ); + determine_legendre_weights( N, weights.data() ); + + const double xi_min = xi_bounds[0]; + const double xi_max = xi_bounds[1]; + // Map the reference quadrature rule to [xi_min, xi_max]. + const double J = 0.5 * ( xi_max - xi_min ); + + for ( int i = 0; i < N; ++i ) { + out.qp[i] = 0.5 * ( xi_max - xi_min ) * qpoints[i] + 0.5 * ( xi_max + xi_min ); + out.w[i] = weights[i] * J; + } + return out; } From 56c31f98d607a73eea2bed2012d17a13baaddbf0 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Fri, 21 Aug 2026 11:03:10 -0700 Subject: [PATCH 5/7] just change function signatures instead of add overloads --- .../tribol_finite_diff_energy_mortar.cpp | 16 +- src/tribol/physics/EnergyMortar.cpp | 141 +++++------------- src/tribol/physics/EnergyMortar.hpp | 18 +-- 3 files changed, 53 insertions(+), 122 deletions(-) diff --git a/src/tests/tribol_finite_diff_energy_mortar.cpp b/src/tests/tribol_finite_diff_energy_mortar.cpp index fb565b6c..3dc059e8 100644 --- a/src/tests/tribol_finite_diff_energy_mortar.cpp +++ b/src/tests/tribol_finite_diff_energy_mortar.cpp @@ -59,11 +59,13 @@ FiniteDiffResult EnergyMortarCalculator::validate_g_tilde( const InterfacePair& auto viewer2 = mesh2.getView(); auto projs0 = projections( pair, viewer1, viewer2 ); - auto bounds0 = smoother_.bounds_from_projections( projs0, p_.del ); - auto smooth_bounds0 = smoother_.smooth_bounds( bounds0, p_.del ); + double bounds0[2]; + smoother_.bounds_from_projections( projs0.data(), p_.del, bounds0 ); + double smooth_bounds0[2]; + smoother_.smooth_bounds( bounds0, p_.del, smooth_bounds0 ); QuadPoints qp0; if ( !p_.enzyme_quadrature ) { - qp0 = compute_quadrature( smooth_bounds0, p_.N ); + compute_quadrature( smooth_bounds0, p_.N, &qp0 ); } auto [g1_base, g2_base] = eval_gtilde( pair, viewer1, viewer2 ); @@ -290,9 +292,11 @@ FiniteDiffResult EnergyMortarCalculator::validate_hessian( const InterfacePair& QuadPoints qp0; if ( !p_.enzyme_quadrature ) { auto projs0 = projections( pair, viewer1, viewer2 ); - auto bounds0 = smoother_.bounds_from_projections( projs0, p_.del ); - auto smooth_bounds0 = smoother_.smooth_bounds( bounds0, p_.del ); - qp0 = compute_quadrature( smooth_bounds0, p_.N ); + double bounds0[2]; + smoother_.bounds_from_projections( projs0.data(), p_.del, bounds0 ); + double smooth_bounds0[2]; + smoother_.smooth_bounds( bounds0, p_.del, smooth_bounds0 ); + compute_quadrature( smooth_bounds0, p_.N, &qp0 ); } auto eval_from_offsets = [&]( const std::array& du ) -> std::pair { diff --git a/src/tribol/physics/EnergyMortar.cpp b/src/tribol/physics/EnergyMortar.cpp index 92e17130..e1a5378b 100644 --- a/src/tribol/physics/EnergyMortar.cpp +++ b/src/tribol/physics/EnergyMortar.cpp @@ -105,24 +105,6 @@ TRIBOL_ENZYME_INLINE void determine_legendre_weights( int N, double* W ) } } -TRIBOL_ENZYME_INLINE void compute_quadrature( const double* xi_bounds, int N, QuadPoints* out ) -{ - double qpoints[3] = { 0.0, 0.0, 0.0 }; - double weights[3] = { 0.0, 0.0, 0.0 }; - - determine_legendre_nodes( N, qpoints ); - determine_legendre_weights( N, weights ); - - const double xi_min = xi_bounds[0]; - const double xi_max = xi_bounds[1]; - const double J = 0.5 * ( xi_max - xi_min ); - - for ( int i = 0; i < N; ++i ) { - out->qp[i] = 0.5 * ( xi_max - xi_min ) * qpoints[i] + 0.5 * ( xi_max + xi_min ); - out->w[i] = weights[i] * J; - } -} - // Map a point from the 1D parent segment coordinate to physical coordinates. // Parametric space: [-0.5, 0.5] TRIBOL_ENZYME_INLINE void iso_map( const double* coord1, const double* coord2, double xi, double* mapped_coord ) @@ -225,28 +207,6 @@ TRIBOL_ENZYME_INLINE void get_projections( const double* A0, const double* A1, c projections[1] = xi_max; } -TRIBOL_ENZYME_INLINE void bounds_from_projections( const double* projections, double del, double* bounds ) -{ - double xi_min = ( projections[0] < projections[1] ) ? projections[0] : projections[1]; - double xi_max = ( projections[0] > projections[1] ) ? projections[0] : projections[1]; - - if ( xi_max < -0.5 - del ) { - xi_max = -0.5 - del; - } - if ( xi_min > 0.5 + del ) { - xi_min = 0.5 + del; - } - if ( xi_min < -0.5 - del ) { - xi_min = -0.5 - del; - } - if ( xi_max > 0.5 + del ) { - xi_max = 0.5 + del; - } - - bounds[0] = xi_min; - bounds[1] = xi_max; -} - TRIBOL_ENZYME_INLINE double smooth_bound( double bound, double del ) { const double xi = bound + 0.5; @@ -270,12 +230,6 @@ TRIBOL_ENZYME_INLINE double smooth_bound( double bound, double del ) return xi_hat - 0.5; } -TRIBOL_ENZYME_INLINE void smooth_bounds( const double* bounds, double del, double* smooth_bounds ) -{ - smooth_bounds[0] = smooth_bound( bounds[0], del ); - smooth_bounds[1] = smooth_bound( bounds[1], del ); -} - // Integrate the nodal smoothed gap and tributary area contributions over edge A. // The quadrature rule is supplied through gp, allowing this kernel to be reused // for both fixed-quadrature and geometry-dependent quadrature paths. @@ -477,11 +431,11 @@ static void kernel_out_enzyme( const double* x, const void* kp_void, double* out // Recompute the integration bounds and quadrature from the current geometry. double bounds[2]; - bounds_from_projections( projs, kp->del, bounds ); + ContactSmoothing::bounds_from_projections( projs, kp->del, bounds ); double xi_bounds[2]; - smooth_bounds( bounds, kp->del, xi_bounds ); + ContactSmoothing::smooth_bounds( bounds, kp->del, xi_bounds ); QuadPoints qp; - compute_quadrature( xi_bounds, kp->N, &qp ); + EnergyMortarCalculator::compute_quadrature( xi_bounds, kp->N, &qp ); Gparams gp; for ( std::size_t i = 0; i < qp.qp.size(); ++i ) { @@ -569,11 +523,11 @@ TRIBOL_ENZYME_INLINE void qp_penalty_kernel( const double* x, const KernelParams double projs[2] = { 0.0, 0.0 }; get_projections( A0, A1, B0, B1, projs ); double bounds[2]; - bounds_from_projections( projs, kp->del, bounds ); + ContactSmoothing::bounds_from_projections( projs, kp->del, bounds ); double xi_bounds[2]; - smooth_bounds( bounds, kp->del, xi_bounds ); + ContactSmoothing::smooth_bounds( bounds, kp->del, xi_bounds ); QuadPoints qp; - compute_quadrature( xi_bounds, kp->N, &qp ); + EnergyMortarCalculator::compute_quadrature( xi_bounds, kp->N, &qp ); double nB[2]; find_normal( B0, B1, nB ); @@ -654,10 +608,13 @@ Gparams EnergyMortarCalculator::construct_gparams( const InterfacePair& pair, co // Build the smoothed integration bounds from the projection of edge B onto edge A. auto projs = EnergyMortarCalculator::compute_projection_bounds( pair, mesh1, mesh2 ); - auto bounds = smoother_.bounds_from_projections( projs, p_.del ); - auto smooth_bounds = smoother_.smooth_bounds( bounds, p_.del ); + double bounds[2]; + smoother_.bounds_from_projections( projs.data(), p_.del, bounds ); + double smooth_bounds[2]; + smoother_.smooth_bounds( bounds, p_.del, smooth_bounds ); - auto qp = EnergyMortarCalculator::compute_quadrature( smooth_bounds, p_.N ); + QuadPoints qp; + EnergyMortarCalculator::compute_quadrature( smooth_bounds, p_.N, &qp ); const int N = static_cast( qp.qp.size() ); @@ -703,11 +660,11 @@ std::array EnergyMortarCalculator::projections( const InterfacePair& } // Clamp the projection interval to the local smoothing support around edge A. -TRIBOL_ENZYME_INLINE std::array ContactSmoothing::bounds_from_projections( const std::array& proj, - double del ) +TRIBOL_ENZYME_INLINE void ContactSmoothing::bounds_from_projections( const double* projections, double del, + double* bounds ) { - double xi_min = std::min( proj[0], proj[1] ); - double xi_max = std::max( proj[0], proj[1] ); + double xi_min = std::min( projections[0], projections[1] ); + double xi_max = std::max( projections[0], projections[1] ); // Limit the integration interval to the extended range [-0.5 - del, 0.5 + del]. if ( xi_max < -0.5 - del ) { @@ -723,7 +680,8 @@ TRIBOL_ENZYME_INLINE std::array ContactSmoothing::bounds_from_project xi_max = 0.5 + del; } - return { xi_min, xi_max }; + bounds[0] = xi_min; + bounds[1] = xi_max; } // Smooth the integration bounds using a C1 ramp near the ends of edge A. @@ -731,51 +689,21 @@ TRIBOL_ENZYME_INLINE std::array ContactSmoothing::bounds_from_project // Bounds of intergration by applying a quadratic ramping function near the ends of the paramteric // space. The smooth region/length is defined by the input del. The returned 'bounds' is the new bounds // of intergation that result after the quadratic ramping has been applied. -TRIBOL_ENZYME_INLINE std::array ContactSmoothing::smooth_bounds( const std::array& bounds, - double del ) +TRIBOL_ENZYME_INLINE void ContactSmoothing::smooth_bounds( const double* bounds, double del, double* smooth_bounds ) { - std::array smooth_bounds; - for ( int i = 0; i < 2; ++i ) { - double xi = 0.0; - double xi_hat = 0.0; - - // Shift from the local coordinate interval [-0.5, 0.5] to [0, 1]. - xi = bounds[i] + 0.5; - if ( del == 0.0 ) { - xi_hat = xi; - } else { - // Apply quadratic ramps near the endpoints and leave the interior unchanged. - if ( 0.0 - del <= xi && xi <= del ) { - xi_hat = ( 1.0 / ( 4 * del ) ) * ( xi * xi ) + 0.5 * xi + del / 4.0; - } else if ( ( 1.0 - del ) <= xi && xi <= 1.0 + del ) { - double b = -1.0 / ( 4.0 * del ); - double c = 0.5 + 1.0 / ( 2.0 * del ); - double d = 1.0 - del + ( 1.0 / ( 4.0 * del ) ) * pow( 1.0 - del, 2 ) - 0.5 * ( 1.0 - del ) - - ( 1.0 - del ) / ( 2.0 * del ); - - xi_hat = b * xi * xi + c * xi + d; - } else if ( del <= xi && xi <= ( 1.0 - del ) ) { - xi_hat = xi; - } - } - // Shift the smoothed coordinate back to [-0.5, 0.5]. - smooth_bounds[i] = xi_hat - 0.5; - } - - return smooth_bounds; + smooth_bounds[0] = smooth_bound( bounds[0], del ); + smooth_bounds[1] = smooth_bound( bounds[1], del ); } // Build a three-point Gauss-Legendre quadrature rule over the local integration bounds. -TRIBOL_ENZYME_INLINE QuadPoints EnergyMortarCalculator::compute_quadrature( const std::array& xi_bounds, - int N ) +TRIBOL_ENZYME_INLINE void EnergyMortarCalculator::compute_quadrature( const double* xi_bounds, int N, + QuadPoints* quadrature ) { - QuadPoints out; - - std::array qpoints; - std::array weights; + double qpoints[3] = { 0.0, 0.0, 0.0 }; + double weights[3] = { 0.0, 0.0, 0.0 }; - determine_legendre_nodes( N, qpoints.data() ); - determine_legendre_weights( N, weights.data() ); + determine_legendre_nodes( N, qpoints ); + determine_legendre_weights( N, weights ); const double xi_min = xi_bounds[0]; const double xi_max = xi_bounds[1]; @@ -783,11 +711,9 @@ TRIBOL_ENZYME_INLINE QuadPoints EnergyMortarCalculator::compute_quadrature( cons const double J = 0.5 * ( xi_max - xi_min ); for ( int i = 0; i < N; ++i ) { - out.qp[i] = 0.5 * ( xi_max - xi_min ) * qpoints[i] + 0.5 * ( xi_max + xi_min ); - out.w[i] = weights[i] * J; + quadrature->qp[i] = 0.5 * ( xi_max - xi_min ) * qpoints[i] + 0.5 * ( xi_max + xi_min ); + quadrature->w[i] = weights[i] * J; } - - return out; } // Evaluate the weighted normal gap at local coordinate xiA on edge A. @@ -836,10 +762,13 @@ NodalContactData EnergyMortarCalculator::compute_nodal_contact_data( const Inter auto projs = projections( pair, mesh1, mesh2 ); // Build the smoothed integration interval from the projection bounds. - auto bounds = smoother_.bounds_from_projections( projs, p_.del ); - auto smooth_bounds = smoother_.smooth_bounds( bounds, p_.del ); + double bounds[2]; + smoother_.bounds_from_projections( projs.data(), p_.del, bounds ); + double smooth_bounds[2]; + smoother_.smooth_bounds( bounds, p_.del, smooth_bounds ); - auto qp = compute_quadrature( smooth_bounds, p_.N ); + QuadPoints qp; + compute_quadrature( smooth_bounds, p_.N, &qp ); double g_tilde1 = 0.0; double g_tilde2 = 0.0; diff --git a/src/tribol/physics/EnergyMortar.hpp b/src/tribol/physics/EnergyMortar.hpp index 21316132..cf0f3217 100644 --- a/src/tribol/physics/EnergyMortar.hpp +++ b/src/tribol/physics/EnergyMortar.hpp @@ -91,17 +91,15 @@ class ContactSmoothing { public: /// Clamp the projected overlap interval to the extended smoothing support. /// - /// The input `proj` contains the local projection bounds of edge B onto edge A. - /// The returned interval is restricted to the extended local range - /// `[-0.5 - del, 0.5 + del]`. - static std::array bounds_from_projections( const std::array& proj, double del ); + /// The input `projections` contains the local projection bounds of edge B onto edge A. + /// The output `bounds` is restricted to the extended local range `[-0.5 - del, 0.5 + del]`. + static void bounds_from_projections( const double* projections, double del, double* bounds ); /// Smooth the integration bounds using the smoothing length `del`. /// - /// The returned bounds are obtained by applying the endpoint smoothing map to - /// the clamped integration interval. When `del = 0`, the bounds are returned - /// without smoothing. - static std::array smooth_bounds( const std::array& bounds, double del ); + /// The output `smooth_bounds` is obtained by applying the endpoint smoothing map to + /// the clamped integration interval. When `del = 0`, the bounds are unchanged. + static void smooth_bounds( const double* bounds, double del, double* smooth_bounds ); }; /// Evaluates Energy Mortar contact quantities for a single interface pair. @@ -132,8 +130,8 @@ class EnergyMortarCalculator { /// /// The input bounds are local coordinates on edge A. The returned quadrature /// points and weights are mapped from the reference interval to - /// `[xi_bounds[0], xi_bounds[1]]` - static QuadPoints compute_quadrature( const std::array& xi_bounds, int N ); + /// `[xi_bounds[0], xi_bounds[1]]` and written to `quadrature`. + static void compute_quadrature( const double* xi_bounds, int N, QuadPoints* quadrature ); /// Compute the nodal smoothed gap integrals and tributary areas. /// From 4ac2dc1cd1df413ce889db1ee84d21b87b18d7aa Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Fri, 21 Aug 2026 11:06:28 -0700 Subject: [PATCH 6/7] restore smooth bounds loop body --- src/tribol/physics/EnergyMortar.cpp | 51 +++++++++++++++-------------- 1 file changed, 26 insertions(+), 25 deletions(-) diff --git a/src/tribol/physics/EnergyMortar.cpp b/src/tribol/physics/EnergyMortar.cpp index e1a5378b..9ecf7b45 100644 --- a/src/tribol/physics/EnergyMortar.cpp +++ b/src/tribol/physics/EnergyMortar.cpp @@ -207,29 +207,6 @@ TRIBOL_ENZYME_INLINE void get_projections( const double* A0, const double* A1, c projections[1] = xi_max; } -TRIBOL_ENZYME_INLINE double smooth_bound( double bound, double del ) -{ - const double xi = bound + 0.5; - double xi_hat = 0.0; - - if ( del == 0.0 ) { - xi_hat = xi; - } else if ( 0.0 - del <= xi && xi <= del ) { - xi_hat = ( 1.0 / ( 4.0 * del ) ) * ( xi * xi ) + 0.5 * xi + del / 4.0; - } else if ( ( 1.0 - del ) <= xi && xi <= 1.0 + del ) { - const double b = -1.0 / ( 4.0 * del ); - const double c = 0.5 + 1.0 / ( 2.0 * del ); - const double one_minus_del = 1.0 - del; - const double d = 1.0 - del + ( 1.0 / ( 4.0 * del ) ) * one_minus_del * one_minus_del - 0.5 * one_minus_del - - one_minus_del / ( 2.0 * del ); - xi_hat = b * xi * xi + c * xi + d; - } else if ( del <= xi && xi <= ( 1.0 - del ) ) { - xi_hat = xi; - } - - return xi_hat - 0.5; -} - // Integrate the nodal smoothed gap and tributary area contributions over edge A. // The quadrature rule is supplied through gp, allowing this kernel to be reused // for both fixed-quadrature and geometry-dependent quadrature paths. @@ -691,8 +668,32 @@ TRIBOL_ENZYME_INLINE void ContactSmoothing::bounds_from_projections( const doubl // of intergation that result after the quadratic ramping has been applied. TRIBOL_ENZYME_INLINE void ContactSmoothing::smooth_bounds( const double* bounds, double del, double* smooth_bounds ) { - smooth_bounds[0] = smooth_bound( bounds[0], del ); - smooth_bounds[1] = smooth_bound( bounds[1], del ); + for ( int i = 0; i < 2; ++i ) { + double xi = 0.0; + double xi_hat = 0.0; + + // Shift from the local coordinate interval [-0.5, 0.5] to [0, 1]. + xi = bounds[i] + 0.5; + if ( del == 0.0 ) { + xi_hat = xi; + } else { + // Apply quadratic ramps near the endpoints and leave the interior unchanged. + if ( 0.0 - del <= xi && xi <= del ) { + xi_hat = ( 1.0 / ( 4 * del ) ) * ( xi * xi ) + 0.5 * xi + del / 4.0; + } else if ( ( 1.0 - del ) <= xi && xi <= 1.0 + del ) { + double b = -1.0 / ( 4.0 * del ); + double c = 0.5 + 1.0 / ( 2.0 * del ); + double d = 1.0 - del + ( 1.0 / ( 4.0 * del ) ) * pow( 1.0 - del, 2 ) - 0.5 * ( 1.0 - del ) - + ( 1.0 - del ) / ( 2.0 * del ); + + xi_hat = b * xi * xi + c * xi + d; + } else if ( del <= xi && xi <= ( 1.0 - del ) ) { + xi_hat = xi; + } + } + // Shift the smoothed coordinate back to [-0.5, 0.5]. + smooth_bounds[i] = xi_hat - 0.5; + } } // Build a three-point Gauss-Legendre quadrature rule over the local integration bounds. From f0c6b4d1baf522efcb7da14f3b09c77fcc9efea3 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Fri, 21 Aug 2026 13:44:36 -0700 Subject: [PATCH 7/7] restored smooth_bound(), needed to avoid loop-local tape reuse in enzyme leading to gradient errors on -O0 --- src/tribol/physics/EnergyMortar.cpp | 57 ++++++++++++++++------------- 1 file changed, 31 insertions(+), 26 deletions(-) diff --git a/src/tribol/physics/EnergyMortar.cpp b/src/tribol/physics/EnergyMortar.cpp index 9ecf7b45..fa8d2424 100644 --- a/src/tribol/physics/EnergyMortar.cpp +++ b/src/tribol/physics/EnergyMortar.cpp @@ -207,6 +207,35 @@ TRIBOL_ENZYME_INLINE void get_projections( const double* A0, const double* A1, c projections[1] = xi_max; } +// Isolate each endpoint to avoid incorrect loop-local tape reuse in Enzyme reverse mode. +TRIBOL_ENZYME_INLINE double smooth_bound( double bound, double del ) +{ + double xi = 0.0; + double xi_hat = 0.0; + + // Shift from the local coordinate interval [-0.5, 0.5] to [0, 1]. + xi = bound + 0.5; + if ( del == 0.0 ) { + xi_hat = xi; + } else { + // Apply quadratic ramps near the endpoints and leave the interior unchanged. + if ( 0.0 - del <= xi && xi <= del ) { + xi_hat = ( 1.0 / ( 4 * del ) ) * ( xi * xi ) + 0.5 * xi + del / 4.0; + } else if ( ( 1.0 - del ) <= xi && xi <= 1.0 + del ) { + double b = -1.0 / ( 4.0 * del ); + double c = 0.5 + 1.0 / ( 2.0 * del ); + double d = 1.0 - del + ( 1.0 / ( 4.0 * del ) ) * pow( 1.0 - del, 2 ) - 0.5 * ( 1.0 - del ) - + ( 1.0 - del ) / ( 2.0 * del ); + + xi_hat = b * xi * xi + c * xi + d; + } else if ( del <= xi && xi <= ( 1.0 - del ) ) { + xi_hat = xi; + } + } + // Shift the smoothed coordinate back to [-0.5, 0.5]. + return xi_hat - 0.5; +} + // Integrate the nodal smoothed gap and tributary area contributions over edge A. // The quadrature rule is supplied through gp, allowing this kernel to be reused // for both fixed-quadrature and geometry-dependent quadrature paths. @@ -668,32 +697,8 @@ TRIBOL_ENZYME_INLINE void ContactSmoothing::bounds_from_projections( const doubl // of intergation that result after the quadratic ramping has been applied. TRIBOL_ENZYME_INLINE void ContactSmoothing::smooth_bounds( const double* bounds, double del, double* smooth_bounds ) { - for ( int i = 0; i < 2; ++i ) { - double xi = 0.0; - double xi_hat = 0.0; - - // Shift from the local coordinate interval [-0.5, 0.5] to [0, 1]. - xi = bounds[i] + 0.5; - if ( del == 0.0 ) { - xi_hat = xi; - } else { - // Apply quadratic ramps near the endpoints and leave the interior unchanged. - if ( 0.0 - del <= xi && xi <= del ) { - xi_hat = ( 1.0 / ( 4 * del ) ) * ( xi * xi ) + 0.5 * xi + del / 4.0; - } else if ( ( 1.0 - del ) <= xi && xi <= 1.0 + del ) { - double b = -1.0 / ( 4.0 * del ); - double c = 0.5 + 1.0 / ( 2.0 * del ); - double d = 1.0 - del + ( 1.0 / ( 4.0 * del ) ) * pow( 1.0 - del, 2 ) - 0.5 * ( 1.0 - del ) - - ( 1.0 - del ) / ( 2.0 * del ); - - xi_hat = b * xi * xi + c * xi + d; - } else if ( del <= xi && xi <= ( 1.0 - del ) ) { - xi_hat = xi; - } - } - // Shift the smoothed coordinate back to [-0.5, 0.5]. - smooth_bounds[i] = xi_hat - 0.5; - } + smooth_bounds[0] = smooth_bound( bounds[0], del ); + smooth_bounds[1] = smooth_bound( bounds[1], del ); } // Build a three-point Gauss-Legendre quadrature rule over the local integration bounds.