From ec65f7e89d9b99451afa587f0269f9fe2937d6e3 Mon Sep 17 00:00:00 2001 From: Chris White Date: Mon, 6 Jul 2026 16:46:10 -0700 Subject: [PATCH 1/3] use project specific option, add needed mpi dependency --- src/tests/CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt index 20164e57..40c719dc 100644 --- a/src/tests/CMakeLists.txt +++ b/src/tests/CMakeLists.txt @@ -252,7 +252,8 @@ if( MFEM_FOUND ) set(mfem_smoke_depends mfem gtest) # MFEM's openmp dependency isn't being picked up by FindMFEM.cmake. Force it to be added for now. - blt_list_append(TO mfem_smoke_depends ELEMENTS blt::openmp IF ENABLE_OPENMP ) + blt_list_append(TO mfem_smoke_depends ELEMENTS blt::openmp IF TRIBOL_USE_OPENMP ) + blt_list_append(TO mfem_smoke_depends ELEMENTS blt::mpi IF TRIBOL_USE_MPI ) blt_add_executable( NAME mfem_smoke_test SOURCES mfem_smoke.cpp From d4566e8c2ddeedb00deff22ee5e2665441b2bbf8 Mon Sep 17 00:00:00 2001 From: Chris White Date: Mon, 6 Jul 2026 16:46:33 -0700 Subject: [PATCH 2/3] fix warnings that were promoted to errors in the smith build --- src/tests/tribol_mfem_jacobian.cpp | 6 +- src/tribol/geom/CompGeom.hpp | 156 +++++++++++++------------ src/tribol/geom/GeomUtilities.hpp | 34 +++--- src/tribol/mesh/CouplingScheme.hpp | 24 ++-- src/tribol/mesh/MeshData.hpp | 2 +- src/tribol/mesh/MethodCouplingData.hpp | 2 +- src/tribol/mesh/MfemData.hpp | 20 ++-- 7 files changed, 127 insertions(+), 117 deletions(-) diff --git a/src/tests/tribol_mfem_jacobian.cpp b/src/tests/tribol_mfem_jacobian.cpp index 2fba02b3..065da083 100644 --- a/src/tests/tribol_mfem_jacobian.cpp +++ b/src/tests/tribol_mfem_jacobian.cpp @@ -695,7 +695,7 @@ TEST_F( MfemJacobianTest, jacobian_primary_dual_assembles_in_ho_and_lor ) EXPECT_EQ( mat->Height(), ctx.mesh_data->GetParentCoords().ParFESpace()->GetTrueVSize() ); EXPECT_EQ( mat->Width(), ctx.submesh_data->GetSubmeshFESpace().GetTrueVSize() ); - const int local_pairs = contributions[0].row_elem_ids.size(); + int local_pairs = contributions[0].row_elem_ids.size(); int global_pairs = 0; MPI_Allreduce( &local_pairs, &global_pairs, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD ); if ( global_pairs > 0 ) { @@ -724,7 +724,7 @@ TEST_F( MfemJacobianTest, jacobian_dual_primary_assembles_in_ho_and_lor ) EXPECT_EQ( mat->Height(), ctx.submesh_data->GetSubmeshFESpace().GetTrueVSize() ); EXPECT_EQ( mat->Width(), ctx.mesh_data->GetParentCoords().ParFESpace()->GetTrueVSize() ); - const int local_pairs = contributions[0].row_elem_ids.size(); + int local_pairs = contributions[0].row_elem_ids.size(); int global_pairs = 0; MPI_Allreduce( &local_pairs, &global_pairs, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD ); if ( global_pairs > 0 ) { @@ -795,7 +795,7 @@ TEST_F( MfemJacobianTest, jacobian_primary_primary_aggregates_mortar_and_nonmort expected += AssembleSolverBlockJacobian( ctx, { nm }, SolverBlock::Primary, SolverBlock::Primary ); expected += AssembleSolverBlockJacobian( ctx, { nn }, SolverBlock::Primary, SolverBlock::Primary ); - const int local_pairs = + int local_pairs = mm.row_elem_ids.size() + mn.row_elem_ids.size() + nm.row_elem_ids.size() + nn.row_elem_ids.size(); int global_pairs = 0; MPI_Allreduce( &local_pairs, &global_pairs, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD ); diff --git a/src/tribol/geom/CompGeom.hpp b/src/tribol/geom/CompGeom.hpp index f0ef6036..b9a3e2c1 100644 --- a/src/tribol/geom/CompGeom.hpp +++ b/src/tribol/geom/CompGeom.hpp @@ -190,7 +190,6 @@ class ContactPlanePair : public CompGeomPair { TRIBOL_HOST_DEVICE inline void getFace1Coords( RealT* x1, int num_coords ) const; TRIBOL_HOST_DEVICE inline void getFace2Coords( RealT* x2, int num_coords ) const; - ; TRIBOL_HOST_DEVICE inline void getFace1ProjectedCoords( RealT* x1_proj, int num_coords ) const; @@ -297,14 +296,14 @@ class ContactPlanePair : public CompGeomPair { * * \return Face id */ - TRIBOL_HOST_DEVICE inline int getCpElementId1() const { return m_pair->m_element_id1; } + TRIBOL_HOST_DEVICE inline IndexT getCpElementId1() const { return m_pair->m_element_id1; } /*! * \brief Get the id of the second element that forms the contact plane * * \return Face id */ - TRIBOL_HOST_DEVICE inline int getCpElementId2() const { return m_pair->m_element_id2; } + TRIBOL_HOST_DEVICE inline IndexT getCpElementId2() const { return m_pair->m_element_id2; } /*! * \brief Set the first contact plane element id @@ -512,8 +511,8 @@ class CommonPlanePair : public ContactPlanePair { * */ TRIBOL_HOST_DEVICE inline bool exceedsMaxAutoInterpen( const MeshData::Viewer& mesh1, const MeshData::Viewer& mesh2, - const int faceId1, const int faceId2, const Parameters& params, - const RealT gap ); + const IndexT faceId1, const IndexT faceId2, + const Parameters& params, const RealT gap ); }; // end class CommonPlanePair @@ -740,27 +739,27 @@ class CompGeom { * * @return common plane object */ - TRIBOL_HOST_DEVICE inline CommonPlanePair& getCommonPlane( int id ) const { return m_common_plane_pairs[id]; } + TRIBOL_HOST_DEVICE inline CommonPlanePair& getCommonPlane( IndexT id ) const { return m_common_plane_pairs[id]; } /** * @brief Get a single mortar plane from the array view * * @return mortar plane object */ - TRIBOL_HOST_DEVICE inline MortarPlanePair& getMortarPlane( int id ) const { return m_mortar_plane_pairs[id]; } + TRIBOL_HOST_DEVICE inline MortarPlanePair& getMortarPlane( IndexT id ) const { return m_mortar_plane_pairs[id]; } /** * @brief Get a single aligned mortar plane from the array view * * @return algined mortar plane object */ - TRIBOL_HOST_DEVICE inline AlignedMortarPlanePair& getAlignedMortarPlane( int id ) const + TRIBOL_HOST_DEVICE inline AlignedMortarPlanePair& getAlignedMortarPlane( IndexT id ) const { return m_aligned_mortar_plane_pairs[id]; } template - TRIBOL_HOST_DEVICE inline T& getPlane( int id ); + TRIBOL_HOST_DEVICE inline T& getPlane( IndexT id ); private: ArrayViewT m_common_plane_pairs; @@ -807,7 +806,7 @@ class CompGeom { * * @return common plane object */ - const CommonPlanePair& getCommonPlane( int id ) const { return m_common_plane_pairs[id]; } + const CommonPlanePair& getCommonPlane( IndexT id ) const { return m_common_plane_pairs[id]; } /** * @brief Get the list of mortar plane pairs @@ -821,7 +820,7 @@ class CompGeom { * * @return mortar plane object */ - const MortarPlanePair& getMortarPlane( int id ) const { return m_mortar_plane_pairs[id]; } + const MortarPlanePair& getMortarPlane( IndexT id ) const { return m_mortar_plane_pairs[id]; } /** * @brief Get the list of aligned mortar plane pairs @@ -835,7 +834,7 @@ class CompGeom { * * @return aligned mortar plane object */ - const AlignedMortarPlanePair& getAlignedMortarPlane( int id ) const { return m_aligned_mortar_plane_pairs[id]; } + const AlignedMortarPlanePair& getAlignedMortarPlane( IndexT id ) const { return m_aligned_mortar_plane_pairs[id]; } /** * @brief Allocate contact plane pairs arrays based on contact method @@ -865,7 +864,7 @@ class CompGeom { } // end switch } - int getNumActivePairs( const ContactMethod method ) const + IndexT getNumActivePairs( const ContactMethod method ) const { switch ( method ) { case COMMON_PLANE: { @@ -961,7 +960,7 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CheckInterfacePairByMethod( InterfacePair& pair, const MeshData::Viewer& mesh1, const MeshData::Viewer& mesh2, const Parameters& params, ContactCase TRIBOL_UNUSED_PARAM( cCase ), bool& isInteracting, CompGeom::Viewer& cg, IndexT* plane_ct ) { - auto dim = mesh1.spatialDimension(); + auto dim = static_cast( mesh1.spatialDimension() ); T my_plane( &pair, params, dim ); FaceGeomException face_err = NO_FACE_GEOM_EXCEPTION; if ( dim == 3 ) { @@ -1028,7 +1027,7 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CheckInterfacePair( InterfacePair& p //------------------------------------------------------------------------------ TRIBOL_HOST_DEVICE inline bool CommonPlanePair::exceedsMaxAutoInterpen( const MeshData::Viewer& mesh1, const MeshData::Viewer& mesh2, - const int faceId1, const int faceId2, + const IndexT faceId1, const IndexT faceId2, const Parameters& params, const RealT gap ) { if ( params.auto_contact_check ) { @@ -1152,17 +1151,19 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::checkFacePair( cons // mark the convexity of each face constexpr int max_nodes = 4; + const int num_nodes_1 = static_cast( mesh1.numberOfNodesPerElement() ); + const int num_nodes_2 = static_cast( mesh2.numberOfNodesPerElement() ); RealT x1_loc[max_nodes]; RealT y1_loc[max_nodes]; RealT x2_loc[max_nodes]; RealT y2_loc[max_nodes]; GlobalTo2DLocalCoords( &m_x1_prime[0], &m_y1_prime[0], &m_z1_prime[0], m_e1X, m_e1Y, m_e1Z, m_e2X, m_e2Y, m_e2Z, m_cX, - m_cY, m_cZ, &x1_loc[0], &y1_loc[0], mesh1.numberOfNodesPerElement() ); + m_cY, m_cZ, &x1_loc[0], &y1_loc[0], num_nodes_1 ); GlobalTo2DLocalCoords( &m_x2_prime[0], &m_y2_prime[0], &m_z2_prime[0], m_e1X, m_e1Y, m_e1Z, m_e2X, m_e2Y, m_e2Z, m_cX, - m_cY, m_cZ, &x2_loc[0], &y2_loc[0], mesh2.numberOfNodesPerElement() ); + m_cY, m_cZ, &x2_loc[0], &y2_loc[0], num_nodes_2 ); - m_face1_convex = IsConvex( x1_loc, y1_loc, mesh1.numberOfNodesPerElement() ); - m_face2_convex = IsConvex( x2_loc, y2_loc, mesh2.numberOfNodesPerElement() ); + m_face1_convex = IsConvex( x1_loc, y1_loc, num_nodes_1 ); + m_face2_convex = IsConvex( x2_loc, y2_loc, num_nodes_2 ); // explicitly call compute overlap routine for common plane so CUDA can determine stack size FaceGeomException interpen_err = CommonPlanePair::computeOverlap3D( @@ -1190,13 +1191,13 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::checkEdgePair( cons // set face "prime" coordinates to current configuration coordinates so we can pull these // off the plane object later for ( int i = 0; i < mesh1.numberOfNodesPerElement(); ++i ) { - const int nodeId = mesh1.getGlobalNodeId( element_id1, i ); + const IndexT nodeId = mesh1.getGlobalNodeId( element_id1, i ); m_x1_prime[i] = mesh1.getPosition()[0][nodeId]; m_y1_prime[i] = mesh1.getPosition()[1][nodeId]; } for ( int i = 0; i < mesh2.numberOfNodesPerElement(); ++i ) { - const int nodeId = mesh2.getGlobalNodeId( element_id2, i ); + const IndexT nodeId = mesh2.getGlobalNodeId( element_id2, i ); m_x2_prime[i] = mesh2.getPosition()[0][nodeId]; m_y2_prime[i] = mesh2.getPosition()[1][nodeId]; } @@ -1247,14 +1248,14 @@ TRIBOL_HOST_DEVICE inline FaceGeomException MortarPlanePair::checkFacePair( cons // set face "prime" coordinates to current configuration coordinates for full overlap methods for ( int i = 0; i < mesh1.numberOfNodesPerElement(); ++i ) { - const int nodeId = mesh1.getGlobalNodeId( element_id1, i ); + const IndexT nodeId = mesh1.getGlobalNodeId( element_id1, i ); m_x1_prime[i] = mesh1.getPosition()[0][nodeId]; m_y1_prime[i] = mesh1.getPosition()[1][nodeId]; m_z1_prime[i] = mesh1.getPosition()[2][nodeId]; } for ( int i = 0; i < mesh2.numberOfNodesPerElement(); ++i ) { - const int nodeId = mesh2.getGlobalNodeId( element_id2, i ); + const IndexT nodeId = mesh2.getGlobalNodeId( element_id2, i ); m_x2_prime[i] = mesh2.getPosition()[0][nodeId]; m_y2_prime[i] = mesh2.getPosition()[1][nodeId]; m_z2_prime[i] = mesh2.getPosition()[2][nodeId]; @@ -1290,12 +1291,14 @@ TRIBOL_HOST_DEVICE inline FaceGeomException MortarPlanePair::computeOverlap3D( c { IndexT element_id1 = this->getCpElementId1(); IndexT element_id2 = this->getCpElementId2(); + const int num_nodes_1 = static_cast( m1.numberOfNodesPerElement() ); + const int num_nodes_2 = static_cast( m2.numberOfNodesPerElement() ); // project face vertex coordinates to contact plane ProjectPointsToPlane( x1, y1, z1, this->m_nX, this->m_nY, this->m_nZ, this->m_cX, this->m_cY, this->m_cZ, - &m_x1_bar[0], &m_y1_bar[0], &m_z1_bar[0], m1.numberOfNodesPerElement() ); + &m_x1_bar[0], &m_y1_bar[0], &m_z1_bar[0], num_nodes_1 ); ProjectPointsToPlane( x2, y2, z2, this->m_nX, this->m_nY, this->m_nZ, this->m_cX, this->m_cY, this->m_cZ, - &m_x2_bar[0], &m_y2_bar[0], &m_z2_bar[0], m2.numberOfNodesPerElement() ); + &m_x2_bar[0], &m_y2_bar[0], &m_z2_bar[0], num_nodes_2 ); // project the projected global nodal coordinates onto local // contact plane 2D coordinate system. @@ -1306,9 +1309,9 @@ TRIBOL_HOST_DEVICE inline FaceGeomException MortarPlanePair::computeOverlap3D( c RealT y2_bar_local[max_nodes_per_elem]; this->globalTo2DLocalCoords( &m_x1_bar[0], &m_y1_bar[0], &m_z1_bar[0], &x1_bar_local[0], &y1_bar_local[0], - m1.numberOfNodesPerElement() ); + num_nodes_1 ); this->globalTo2DLocalCoords( &m_x2_bar[0], &m_y2_bar[0], &m_z2_bar[0], &x2_bar_local[0], &y2_bar_local[0], - m2.numberOfNodesPerElement() ); + num_nodes_2 ); // compute the full intersection polygon vertex coordinates RealT* X1 = &x1_bar_local[0]; @@ -1319,15 +1322,15 @@ TRIBOL_HOST_DEVICE inline FaceGeomException MortarPlanePair::computeOverlap3D( c // assuming each face's vertices are ordered WRT that face's outward unit normal, // reorder face 2 vertices to be consistent with face 1. DO NOT CALL POLYREORDER() // to do this. - ElemReverse( X2, Y2, m2.numberOfNodesPerElement() ); + ElemReverse( X2, Y2, num_nodes_2 ); // compute intersection polygon and area. RealT pos_tol = this->m_params.len_collapse_ratio * axom::utilities::max( m1.getFaceRadius()[element_id1], m2.getFaceRadius()[element_id2] ); RealT len_tol = pos_tol; FaceGeomException inter_err = - Intersection2DPolygon( X1, Y1, m1.numberOfNodesPerElement(), X2, Y2, m2.numberOfNodesPerElement(), pos_tol, - len_tol, this->m_polyLocX, this->m_polyLocY, this->m_numPolyVert, this->m_area, false ); + Intersection2DPolygon( X1, Y1, num_nodes_1, X2, Y2, num_nodes_2, pos_tol, len_tol, this->m_polyLocX, + this->m_polyLocY, this->m_numPolyVert, this->m_area, false ); if ( inter_err != NO_FACE_GEOM_EXCEPTION ) { return inter_err; @@ -1409,14 +1412,14 @@ TRIBOL_HOST_DEVICE inline FaceGeomException AlignedMortarPlanePair::checkFacePai // set face "prime" coordinates to current configuration coordinates for full overlap methods for ( int i = 0; i < mesh1.numberOfNodesPerElement(); ++i ) { - const int nodeId = mesh1.getGlobalNodeId( element_id1, i ); + const IndexT nodeId = mesh1.getGlobalNodeId( element_id1, i ); m_x1_prime[i] = mesh1.getPosition()[0][nodeId]; m_y1_prime[i] = mesh1.getPosition()[1][nodeId]; m_z1_prime[i] = mesh1.getPosition()[2][nodeId]; } for ( int i = 0; i < mesh2.numberOfNodesPerElement(); ++i ) { - const int nodeId = mesh2.getGlobalNodeId( element_id2, i ); + const IndexT nodeId = mesh2.getGlobalNodeId( element_id2, i ); m_x2_prime[i] = mesh2.getPosition()[0][nodeId]; m_y2_prime[i] = mesh2.getPosition()[1][nodeId]; m_z2_prime[i] = mesh2.getPosition()[2][nodeId]; @@ -1455,18 +1458,20 @@ TRIBOL_HOST_DEVICE inline FaceGeomException AlignedMortarPlanePair::computeOverl const MeshData::Viewer& m2 ) { IndexT element_id2 = this->getCpElementId2(); + const int num_nodes_1 = static_cast( m1.numberOfNodesPerElement() ); + const int num_nodes_2 = static_cast( m2.numberOfNodesPerElement() ); // project face vertex coordinates to contact plane ProjectPointsToPlane( x1, y1, z1, this->m_nX, this->m_nY, this->m_nZ, this->m_cX, this->m_cY, this->m_cZ, - &m_x1_bar[0], &m_y1_bar[0], &m_z1_bar[0], m1.numberOfNodesPerElement() ); + &m_x1_bar[0], &m_y1_bar[0], &m_z1_bar[0], num_nodes_1 ); ProjectPointsToPlane( x2, y2, z2, this->m_nX, this->m_nY, this->m_nZ, this->m_cX, this->m_cY, this->m_cZ, - &m_x2_bar[0], &m_y2_bar[0], &m_z2_bar[0], m2.numberOfNodesPerElement() ); + &m_x2_bar[0], &m_y2_bar[0], &m_z2_bar[0], num_nodes_2 ); // Compute face centroids using projected coordinates passed in RealT cx1, cy1, cz1; RealT cx2, cy2, cz2; - VertexAvgCentroid( x1, y1, z1, m1.numberOfNodesPerElement(), cx1, cy1, cz1 ); - VertexAvgCentroid( x2, y2, z2, m2.numberOfNodesPerElement(), cx2, cy2, cz2 ); + VertexAvgCentroid( x1, y1, z1, num_nodes_1, cx1, cy1, cz1 ); + VertexAvgCentroid( x2, y2, z2, num_nodes_2, cx2, cy2, cz2 ); // compute the gap vector between face centroids. Then, project the gap vector // on the contact plane unit normal. The magnitude of the gap vector should be very @@ -1486,8 +1491,8 @@ TRIBOL_HOST_DEVICE inline FaceGeomException AlignedMortarPlanePair::computeOverl // if we are here we have contact between two aligned faces; per mortar method take // the non-mortar (mesh2) face as the contact plane. For Aligned mortar we can use // face 2 as the overlap plane and overlap - this->m_numPolyVert = m2.numberOfNodesPerElement(); - for ( int a = 0; a < m2.numberOfNodesPerElement(); ++a ) { + this->m_numPolyVert = num_nodes_2; + for ( int a = 0; a < num_nodes_2; ++a ) { this->m_polyX[a] = x2[a]; this->m_polyY[a] = y2[a]; this->m_polyZ[a] = z2[a]; @@ -2007,6 +2012,8 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::computeOverlap3D( c // set up vertex id arrays to indicate which face vertices pass through // contact plane (i.e. lie on the other side) StackArrayT mesh( { &m1, &m2 } ); + StackArrayT num_nodes_per_elem( + { static_cast( m1.numberOfNodesPerElement() ), static_cast( m2.numberOfNodesPerElement() ) } ); int interpenVertex1[max_nodes_per_elem]; int interpenVertex2[max_nodes_per_elem]; @@ -2016,9 +2023,9 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::computeOverlap3D( c } StackArrayT element_id( { getCpElementId1(), getCpElementId2() } ); - StackArrayT num_intersections( { 0, 0 } ); - StackArrayT num_intersections_inside( { 0, 0 } ); - StackArrayT num_nodes_otherside( { 0, 0 } ); + StackArrayT num_intersections( { 0, 0 } ); + StackArrayT num_intersections_inside( { 0, 0 } ); + StackArrayT num_nodes_otherside( { 0, 0 } ); // compute interpen data for convex face-pairs only if ( m_face1_convex && m_face2_convex ) { @@ -2050,27 +2057,27 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::computeOverlap3D( c // get the other face normal and centroid for line-face-plane intersections RealT fn[max_dim], cx[max_dim]; - RealT num_nodes_other; + int num_nodes_other; if ( i == 0 ) { mesh[1]->getFaceNormal( element_id[1], fn ); mesh[1]->getFaceCentroid( element_id[1], cx ); - num_nodes_other = mesh[1]->numberOfNodesPerElement(); + num_nodes_other = static_cast( mesh[1]->numberOfNodesPerElement() ); } else { mesh[0]->getFaceNormal( element_id[0], fn ); mesh[0]->getFaceCentroid( element_id[0], cx ); - num_nodes_other = mesh[0]->numberOfNodesPerElement(); + num_nodes_other = static_cast( mesh[0]->numberOfNodesPerElement() ); } int k = 0; int k_inside = 0; int k_otherside = 0; - for ( int j = 0; j < mesh[i]->numberOfNodesPerElement(); ++j ) // loop over face segments + for ( int j = 0; j < num_nodes_per_elem[i]; ++j ) // loop over face segments { bool intersection_is_node = false; // determine local segment vertex ids int ja = j; - int jb = ( j == ( mesh[i]->numberOfNodesPerElement() - 1 ) ) ? 0 : ( j + 1 ); + int jb = ( j == ( num_nodes_per_elem[i] - 1 ) ) ? 0 : ( j + 1 ); // initialize current entry in the vertex id list interpenVertex[ja] = -1; @@ -2092,12 +2099,12 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::computeOverlap3D( c // Debug print faces to screen to catch unforeseen degenerate face configurations etc. SLIC_WARNING( "Degenerate face configuration detected with number of line-plane intersections > 2." ); SLIC_INFO( "Planar coordinates for face 1 in CommonPlanePair::computeOverlap3D(): " ); - for ( int a = 0; a < mesh[0]->numberOfNodesPerElement(); ++a ) { + for ( int a = 0; a < num_nodes_per_elem[0]; ++a ) { std::cout << x1[a] << ", " << y1[a] << ", " << z1[a] << std::endl; } SLIC_INFO( "Planar coordinates for face 2 in CommonPlanePair::computeOverlap3D(): " ); - for ( int b = 0; b < mesh[1]->numberOfNodesPerElement(); ++b ) { + for ( int b = 0; b < num_nodes_per_elem[1]; ++b ) { std::cout << x2[b] << ", " << y2[b] << ", " << z2[b] << std::endl; } @@ -2145,7 +2152,7 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::computeOverlap3D( c // get the local coordinates of the current intersection point RealT xInter_local, yInter_local; Points3DTo2D( &xInter[2 * i + k], &yInter[2 * i + k], &zInter[2 * i + k], fn[0], fn[1], fn[2], cx[0], cx[1], - cx[2], 1.0, &xInter_local, &yInter_local ); + cx[2], 1, &xInter_local, &yInter_local ); // get the local coordinates of the other face's centroid RealT cx_other_local, cy_other_local; @@ -2194,7 +2201,7 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::computeOverlap3D( c // count the number of vertices (intersection points and interpen points) for the clipped // portion of the i^th face that interpenetrates the opposing face. numV[i] = k; // could be zero intersection points - for ( int vid = 0; vid < mesh[i]->numberOfNodesPerElement(); ++vid ) { + for ( int vid = 0; vid < num_nodes_per_elem[i]; ++vid ) { // increment total vertex counter if ids match if ( interpenVertex[vid] == vid ) ++numV[i]; @@ -2281,7 +2288,7 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::computeOverlap3D( c // populate the face 1 vertices that cross the contact plane int k = num_intersections[0]; - for ( int m = 0; m < mesh[0]->numberOfNodesPerElement(); ++m ) { + for ( int m = 0; m < num_nodes_per_elem[0]; ++m ) { if ( interpenVertex1[m] != -1 ) { // set nodal coordinates to the "prime" coords sent into this routine // associated with the average face plane @@ -2294,7 +2301,7 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::computeOverlap3D( c // populate the face 2 vertices that cross the contact plane k = num_intersections[1]; - for ( int n = 0; n < mesh[1]->numberOfNodesPerElement(); ++n ) { + for ( int n = 0; n < num_nodes_per_elem[1]; ++n ) { if ( interpenVertex2[n] != -1 ) { // set nodal coordinates to the "prime" coords sent into this routine // associated with the average face plane @@ -2325,10 +2332,10 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::computeOverlap3D( c } // end if (!m_fullOverlap) if ( m_fullOverlap ) { // populate the face vertex array with the face coordinates themselves - numV[0] = mesh[0]->numberOfNodesPerElement(); - numV[1] = mesh[1]->numberOfNodesPerElement(); + numV[0] = num_nodes_per_elem[0]; + numV[1] = num_nodes_per_elem[1]; // face 1 - for ( int m = 0; m < mesh[0]->numberOfNodesPerElement(); ++m ) { + for ( int m = 0; m < num_nodes_per_elem[0]; ++m ) { // use face averaged "prime" coordinates for consistency cfx1[m] = x1[m]; cfy1[m] = y1[m]; @@ -2336,7 +2343,7 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::computeOverlap3D( c } // face 2 - for ( int n = 0; n < mesh[1]->numberOfNodesPerElement(); ++n ) { + for ( int n = 0; n < num_nodes_per_elem[1]; ++n ) { // use face averaged "prime" coordinates for consistency cfx2[n] = x2[n]; cfy2[n] = y2[n]; @@ -2414,12 +2421,12 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::computeOverlap3D( c } // Project averaged face coordinates sent into this routine to the common plane and store - for ( int i = 0; i < mesh[0]->numberOfNodesPerElement(); ++i ) { + for ( int i = 0; i < num_nodes_per_elem[0]; ++i ) { ProjectPointToPlane( x1[i], y1[i], z1[i], m_nX, m_nY, m_nZ, m_cX, m_cY, m_cZ, m_x1_bar[i], m_y1_bar[i], m_z1_bar[i] ); } - for ( int i = 0; i < mesh[1]->numberOfNodesPerElement(); ++i ) { + for ( int i = 0; i < num_nodes_per_elem[1]; ++i ) { ProjectPointToPlane( x2[i], y2[i], z2[i], m_nX, m_nY, m_nZ, m_cX, m_cY, m_cZ, m_x2_bar[i], m_y2_bar[i], m_z2_bar[i] ); } @@ -2456,7 +2463,8 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::projectPointsAndCom // sanity check #ifdef TRIBOL_USE_HOST if ( m_fullOverlap ) { - if ( num_vert_1 != m1.numberOfNodesPerElement() || num_vert_2 != m2.numberOfNodesPerElement() ) { + if ( num_vert_1 != static_cast( m1.numberOfNodesPerElement() ) || + num_vert_2 != static_cast( m2.numberOfNodesPerElement() ) ) { SLIC_ERROR( "CommonPlanePair::projectPointsAndComputeOverlap(): full overlap requires " << "input number of vertices to match number of nodes per element." ); } @@ -2549,8 +2557,10 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::computeOverlap2D( c { // all edge-edge interactions suitable for an interpenetration overlap // calculation are edges that intersect at a single point - int edgeId1 = getCpElementId1(); - int edgeId2 = getCpElementId2(); + IndexT edgeId1 = getCpElementId1(); + IndexT edgeId2 = getCpElementId2(); + const int num_nodes_1 = static_cast( m1.numberOfNodesPerElement() ); + const int num_nodes_2 = static_cast( m2.numberOfNodesPerElement() ); constexpr int max_dim = 2; RealT cx1[max_dim], cx2[max_dim]; @@ -2600,9 +2610,9 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::computeOverlap2D( c int interId2 = -1; int k1 = 0; int k2 = 0; - for ( int i = 0; i < m1.numberOfNodesPerElement(); ++i ) { - int nodeId1 = m1.getGlobalNodeId( edgeId1, i ); - int nodeId2 = m2.getGlobalNodeId( edgeId2, i ); + for ( int i = 0; i < num_nodes_1; ++i ) { + IndexT nodeId1 = m1.getGlobalNodeId( edgeId1, i ); + IndexT nodeId2 = m2.getGlobalNodeId( edgeId2, i ); RealT lvx1 = m1.getPosition()[0][nodeId1] - cx2[0]; RealT lvy1 = m1.getPosition()[1][nodeId1] - cx2[1]; RealT lvx2 = m2.getPosition()[0][nodeId2] - cx1[0]; @@ -2644,13 +2654,13 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::computeOverlap2D( c // populate arrays holding the interpenetrating edge portions to be // used in computing the overlap - int nodeInter1 = m1.getGlobalNodeId( edgeId1, interId1 ); + IndexT nodeInter1 = m1.getGlobalNodeId( edgeId1, interId1 ); x1_to_project[0] = m1.getPosition()[0][nodeInter1]; x1_to_project[1] = xInter; y1_to_project[0] = m1.getPosition()[1][nodeInter1]; y1_to_project[1] = yInter; - int nodeInter2 = m2.getGlobalNodeId( edgeId2, interId2 ); + IndexT nodeInter2 = m2.getGlobalNodeId( edgeId2, interId2 ); x2_to_project[0] = m2.getPosition()[0][nodeInter2]; x2_to_project[1] = xInter; y2_to_project[0] = m2.getPosition()[1][nodeInter2]; @@ -2728,8 +2738,8 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::computeOverlap2D( c ProjectEdgeNodesToSegment( m1, edgeId1, this->m_nX, this->m_nY, this->m_cX, this->m_cY, &projX1[0], &projY1[0] ); ProjectEdgeNodesToSegment( m2, edgeId2, this->m_nX, this->m_nY, this->m_cX, this->m_cY, &projX2[0], &projY2[0] ); FaceGeomException segError = - CheckSegOverlap( &projX1[0], &projY1[0], &projX2[0], &projY2[0], m1.numberOfNodesPerElement(), - m2.numberOfNodesPerElement(), &m_polyX[0], &m_polyY[0], m_area ); + CheckSegOverlap( &projX1[0], &projY1[0], &projX2[0], &projY2[0], num_nodes_1, num_nodes_2, &m_polyX[0], + &m_polyY[0], m_area ); if ( segError != NO_FACE_GEOM_EXCEPTION ) { return segError; @@ -2772,11 +2782,11 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::computeOverlap2D( c } // end if (m_fullOverlap) // Project edge coordinates to the common plane and store - for ( int i = 0; i < m1.numberOfNodesPerElement(); ++i ) { + for ( int i = 0; i < num_nodes_1; ++i ) { ProjectPointToSegment( m_x1_prime[i], m_y1_prime[i], m_nX, m_nY, m_cX, m_cY, m_x1_bar[i], m_y1_bar[i] ); } - for ( int i = 0; i < m2.numberOfNodesPerElement(); ++i ) { + for ( int i = 0; i < num_nodes_2; ++i ) { ProjectPointToSegment( m_x2_prime[i], m_y2_prime[i], m_nX, m_nY, m_cX, m_cY, m_x2_bar[i], m_y2_bar[i] ); } @@ -2786,19 +2796,19 @@ TRIBOL_HOST_DEVICE inline FaceGeomException CommonPlanePair::computeOverlap2D( c //------------------------------------------------------------------------------ template <> -TRIBOL_HOST_DEVICE inline CommonPlanePair& CompGeom::Viewer::getPlane( int id ) +TRIBOL_HOST_DEVICE inline CommonPlanePair& CompGeom::Viewer::getPlane( IndexT id ) { return m_common_plane_pairs[id]; }; template <> -TRIBOL_HOST_DEVICE inline MortarPlanePair& CompGeom::Viewer::getPlane( int id ) +TRIBOL_HOST_DEVICE inline MortarPlanePair& CompGeom::Viewer::getPlane( IndexT id ) { return m_mortar_plane_pairs[id]; }; template <> -TRIBOL_HOST_DEVICE inline AlignedMortarPlanePair& CompGeom::Viewer::getPlane( int id ) +TRIBOL_HOST_DEVICE inline AlignedMortarPlanePair& CompGeom::Viewer::getPlane( IndexT id ) { return m_aligned_mortar_plane_pairs[id]; }; diff --git a/src/tribol/geom/GeomUtilities.hpp b/src/tribol/geom/GeomUtilities.hpp index 4d72714e..a1332076 100644 --- a/src/tribol/geom/GeomUtilities.hpp +++ b/src/tribol/geom/GeomUtilities.hpp @@ -54,7 +54,7 @@ TRIBOL_HOST_DEVICE inline void ComputeLocalBasis( RealT nx, RealT ny, RealT nz, * * \pre length(pX), length(pY), length(pZ) >= number of nodes on face */ -TRIBOL_HOST_DEVICE inline void ProjectFaceNodesToPlane( const MeshData::Viewer& mesh, int faceId, RealT nrmlX, +TRIBOL_HOST_DEVICE inline void ProjectFaceNodesToPlane( const MeshData::Viewer& mesh, IndexT faceId, RealT nrmlX, RealT nrmlY, RealT nrmlZ, RealT cX, RealT cY, RealT cZ, RealT* pX, RealT* pY, RealT* pZ ); @@ -72,7 +72,7 @@ TRIBOL_HOST_DEVICE inline void ProjectFaceNodesToPlane( const MeshData::Viewer& * \param [out] pY pointer to array of projected nodal y-coordinates * */ -TRIBOL_HOST_DEVICE inline void ProjectEdgeNodesToSegment( const MeshData::Viewer& mesh, int edgeId, RealT nrmlX, +TRIBOL_HOST_DEVICE inline void ProjectEdgeNodesToSegment( const MeshData::Viewer& mesh, IndexT edgeId, RealT nrmlX, RealT nrmlY, RealT cX, RealT cY, RealT* pX, RealT* pY ); /*! @@ -1869,14 +1869,14 @@ TRIBOL_HOST_DEVICE inline void ComputeLocalBasis( RealT nx, RealT ny, RealT nz, } //------------------------------------------------------------------------------ -TRIBOL_HOST_DEVICE inline void ProjectFaceNodesToPlane( const MeshData::Viewer& mesh, int faceId, RealT nrmlX, +TRIBOL_HOST_DEVICE inline void ProjectFaceNodesToPlane( const MeshData::Viewer& mesh, IndexT faceId, RealT nrmlX, RealT nrmlY, RealT nrmlZ, RealT cX, RealT cY, RealT cZ, RealT* pX, RealT* pY, RealT* pZ ) { // loop over nodes and project onto the plane defined by the point-normal // input arguments for ( int i = 0; i < mesh.numberOfNodesPerElement(); ++i ) { - const int nodeId = mesh.getGlobalNodeId( faceId, i ); + const IndexT nodeId = mesh.getGlobalNodeId( faceId, i ); ProjectPointToPlane( mesh.getPosition()[0][nodeId], mesh.getPosition()[1][nodeId], mesh.getPosition()[2][nodeId], nrmlX, nrmlY, nrmlZ, cX, cY, cZ, pX[i], pY[i], pZ[i] ); } @@ -1886,11 +1886,11 @@ TRIBOL_HOST_DEVICE inline void ProjectFaceNodesToPlane( const MeshData::Viewer& } // end ProjectFaceNodesToPlane() //------------------------------------------------------------------------------ -TRIBOL_HOST_DEVICE inline void ProjectEdgeNodesToSegment( const MeshData::Viewer& mesh, int edgeId, RealT nrmlX, +TRIBOL_HOST_DEVICE inline void ProjectEdgeNodesToSegment( const MeshData::Viewer& mesh, IndexT edgeId, RealT nrmlX, RealT nrmlY, RealT cX, RealT cY, RealT* pX, RealT* pY ) { for ( int i = 0; i < mesh.numberOfNodesPerElement(); ++i ) { - const int nodeId = mesh.getGlobalNodeId( edgeId, i ); + const IndexT nodeId = mesh.getGlobalNodeId( edgeId, i ); ProjectPointToSegment( mesh.getPosition()[0][nodeId], mesh.getPosition()[1][nodeId], nrmlX, nrmlY, cX, cY, pX[i], pY[i] ); } @@ -1960,17 +1960,17 @@ TRIBOL_HOST_DEVICE inline void PolyInterYCentroid( const int namax, const RealT* RealT vol; // calculate origin shift to avoid roundoff errors - // TODO figure out numeric limits associated with RealT that also works on device - RealT xorg = FLT_MAX; - RealT yorg = FLT_MAX; - RealT xa_min = FLT_MAX; - RealT xa_max = -FLT_MAX; - RealT ya_min = FLT_MAX; - RealT ya_max = -FLT_MAX; - RealT xb_min = FLT_MAX; - RealT xb_max = -FLT_MAX; - RealT yb_min = FLT_MAX; - RealT yb_max = -FLT_MAX; + constexpr RealT max_real = DBL_MAX; + RealT xorg = max_real; + RealT yorg = max_real; + RealT xa_min = max_real; + RealT xa_max = -max_real; + RealT ya_min = max_real; + RealT ya_max = -max_real; + RealT xb_min = max_real; + RealT xb_max = -max_real; + RealT yb_min = max_real; + RealT yb_max = -max_real; RealT qy = 0.0; diff --git a/src/tribol/mesh/CouplingScheme.hpp b/src/tribol/mesh/CouplingScheme.hpp index e7dbe61a..8b62ca86 100644 --- a/src/tribol/mesh/CouplingScheme.hpp +++ b/src/tribol/mesh/CouplingScheme.hpp @@ -128,7 +128,7 @@ class CouplingScheme { * * @return spatial dimension */ - TRIBOL_HOST_DEVICE int spatialDimension() const { return m_mesh1.spatialDimension(); } + TRIBOL_HOST_DEVICE int spatialDimension() const { return static_cast( m_mesh1.spatialDimension() ); } /** * @brief Return a view of the first mesh in the coupling scheme @@ -192,7 +192,7 @@ class CouplingScheme { * * @return the gap tolerance for the common plane method */ - TRIBOL_HOST_DEVICE inline RealT getGapTol( int fid1, int fid2 ) const; + TRIBOL_HOST_DEVICE inline RealT getGapTol( IndexT fid1, IndexT fid2 ) const; /** * @brief Get a view of the computational geometry container @@ -284,21 +284,21 @@ class CouplingScheme { * * @return unique coupling scheme ID */ - int getId() const { return m_id; } + IndexT getId() const { return m_id; } /** * @brief Get the integer ID of the first mesh * * @return unique ID of the first mesh */ - int getMeshId1() const { return m_mesh_id1; } + IndexT getMeshId1() const { return m_mesh_id1; } /** * @brief Get the integer ID of the second mesh * * @return unique ID of the second mesh */ - int getMeshId2() const { return m_mesh_id2; } + IndexT getMeshId2() const { return m_mesh_id2; } /** * @brief Get the Parameters struct @@ -513,7 +513,7 @@ class CouplingScheme { * * @return number of active interface pairs */ - int getNumActivePairs() const { return m_cg_pairs.getNumActivePairs( getContactMethod() ); } + IndexT getNumActivePairs() const { return m_cg_pairs.getNumActivePairs( getContactMethod() ); } /** * @brief Return the contact plane given by id @@ -1009,7 +1009,7 @@ TRIBOL_HOST_DEVICE inline ContactPlanePair& CouplingScheme::Viewer::getContactPl } //------------------------------------------------------------------------------ -TRIBOL_HOST_DEVICE inline RealT CouplingScheme::Viewer::getGapTol( int fid1, int fid2 ) const +TRIBOL_HOST_DEVICE inline RealT CouplingScheme::Viewer::getGapTol( IndexT fid1, IndexT fid2 ) const { RealT gap_tol = 0.; // add debug warning if this routine is called for interface methods @@ -1072,9 +1072,9 @@ TRIBOL_HOST_DEVICE inline bool CouplingScheme::Viewer::pruneMethodFacePair( cons auto& mesh1 = this->getMesh1View(); auto& mesh2 = this->getMesh2View(); - int dim = mesh1.spatialDimension(); - int num_nodes_face_1 = mesh1.numberOfNodesPerElement(); - int num_nodes_face_2 = mesh2.numberOfNodesPerElement(); + int dim = static_cast( mesh1.spatialDimension() ); + int num_nodes_face_1 = static_cast( mesh1.numberOfNodesPerElement() ); + int num_nodes_face_2 = static_cast( mesh2.numberOfNodesPerElement() ); RealT fn1[max_dim], cx1[max_dim]; mesh1.getFaceNormal( fid1, fn1 ); @@ -1094,7 +1094,7 @@ TRIBOL_HOST_DEVICE inline bool CouplingScheme::Viewer::pruneMethodFacePair( cons RealT z2[max_nodes_per_face]; for ( int i = 0; i < mesh1.numberOfNodesPerElement(); ++i ) { - const int nodeId_1 = mesh1.getGlobalNodeId( fid1, i ); + const IndexT nodeId_1 = mesh1.getGlobalNodeId( fid1, i ); x1[i] = mesh1.getPosition()[0][nodeId_1]; y1[i] = mesh1.getPosition()[1][nodeId_1]; if ( dim == 3 ) { @@ -1103,7 +1103,7 @@ TRIBOL_HOST_DEVICE inline bool CouplingScheme::Viewer::pruneMethodFacePair( cons } for ( int i = 0; i < mesh2.numberOfNodesPerElement(); ++i ) { - const int nodeId_2 = mesh2.getGlobalNodeId( fid2, i ); + const IndexT nodeId_2 = mesh2.getGlobalNodeId( fid2, i ); x2[i] = mesh2.getPosition()[0][nodeId_2]; y2[i] = mesh2.getPosition()[1][nodeId_2]; if ( dim == 3 ) { diff --git a/src/tribol/mesh/MeshData.hpp b/src/tribol/mesh/MeshData.hpp index cd8a4bd1..75d958b3 100644 --- a/src/tribol/mesh/MeshData.hpp +++ b/src/tribol/mesh/MeshData.hpp @@ -159,7 +159,7 @@ class MeshData { * * @return spatial dimension */ - TRIBOL_HOST_DEVICE int spatialDimension() const { return m_position.size(); } + TRIBOL_HOST_DEVICE IndexT spatialDimension() const { return m_position.size(); } /** * @brief Number of nodes in the mesh diff --git a/src/tribol/mesh/MethodCouplingData.hpp b/src/tribol/mesh/MethodCouplingData.hpp index 32e7e045..5aa77c07 100644 --- a/src/tribol/mesh/MethodCouplingData.hpp +++ b/src/tribol/mesh/MethodCouplingData.hpp @@ -231,7 +231,7 @@ class MethodData { * * See @ref getElementBlockJacobians for a definition of the blocks. */ - int getNSpaces() const { return m_blockJSpaces.size(); } + IndexT getNSpaces() const { return m_blockJSpaces.size(); } /*! * \brief Get the element ids for each entry of the getBlockJ 2D ArrayT diff --git a/src/tribol/mesh/MfemData.hpp b/src/tribol/mesh/MfemData.hpp index fb652657..e7c03830 100644 --- a/src/tribol/mesh/MfemData.hpp +++ b/src/tribol/mesh/MfemData.hpp @@ -672,16 +672,16 @@ class MfemMeshData { /** * @brief Get the number of elements in the first Tribol registered mesh * - * @return int + * @return IndexT */ - int GetMesh1NE() const { return GetUpdateData().conn_1_.size() / GetUpdateData().num_verts_per_elem_; } + IndexT GetMesh1NE() const { return GetUpdateData().conn_1_.size() / GetUpdateData().num_verts_per_elem_; } /** * @brief Get the number of elements in the second Tribol registered mesh * - * @return int + * @return IndexT */ - int GetMesh2NE() const { return GetUpdateData().conn_2_.size() / GetUpdateData().num_verts_per_elem_; } + IndexT GetMesh2NE() const { return GetUpdateData().conn_2_.size() / GetUpdateData().num_verts_per_elem_; } /** * @brief Get the total number of vertices in both Tribol registered meshes @@ -1790,10 +1790,10 @@ struct PackedPairJacobianContribs { const Array1D* row_elem_map{ nullptr }; /// Tribol element-id -> redecomp element-id map for column element ids const Array1D* col_elem_map{ nullptr }; - Array1D row_elem_ids; ///< Tribol element IDs for rows - Array1D col_elem_ids; ///< Tribol element IDs for columns + Array1D row_elem_ids; ///< Tribol element IDs for rows + Array1D col_elem_ids; ///< Tribol element IDs for columns Array1D jacobian_data; ///< Flattened Jacobian data - Array1D value_offsets; ///< Offsets into jacobian_data for each element + Array1D value_offsets; ///< Offsets into jacobian_data for each element PackedPairJacobianContribs() = default; @@ -1817,7 +1817,7 @@ struct PackedPairJacobianContribs { * @param n_jacobian_scalar_values Total number of scalar Jacobian values that will be appended across all pairs. * This is used to reserve capacity in @ref jacobian_data. */ - void reserve( int n_pairs, int n_jacobian_scalar_values ) + void reserve( IndexT n_pairs, IndexT n_jacobian_scalar_values ) { row_elem_ids.reserve( n_pairs ); col_elem_ids.reserve( n_pairs ); @@ -1833,7 +1833,7 @@ struct PackedPairJacobianContribs { * @param data Pointer to a contiguous, column-major dense block of Jacobian values. * @param size Number of scalar entries in @p data (typically `n_row_dofs * n_col_dofs`). */ - void append( int row_elem_id, int col_elem_id, const double* data, int size ) + void append( int row_elem_id, int col_elem_id, const double* data, IndexT size ) { row_elem_ids.push_back( row_elem_id ); col_elem_ids.push_back( col_elem_id ); @@ -1846,7 +1846,7 @@ struct PackedPairJacobianContribs { /** * @brief Number of element contributions stored in this packed block */ - int numEntries() const { return row_elem_ids.size(); } + IndexT numEntries() const { return row_elem_ids.size(); } }; /** From a5c589d8f61f90a8cf4726d17c9d8e822629cd95 Mon Sep 17 00:00:00 2001 From: Chris White Date: Tue, 7 Jul 2026 13:19:45 -0700 Subject: [PATCH 3/3] use numeric limits --- src/tribol/geom/GeomUtilities.hpp | 26 +++++++++++++++----------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/src/tribol/geom/GeomUtilities.hpp b/src/tribol/geom/GeomUtilities.hpp index a1332076..54b0990d 100644 --- a/src/tribol/geom/GeomUtilities.hpp +++ b/src/tribol/geom/GeomUtilities.hpp @@ -8,6 +8,7 @@ #include #include +#include #include "axom/slic.hpp" @@ -1960,17 +1961,20 @@ TRIBOL_HOST_DEVICE inline void PolyInterYCentroid( const int namax, const RealT* RealT vol; // calculate origin shift to avoid roundoff errors - constexpr RealT max_real = DBL_MAX; - RealT xorg = max_real; - RealT yorg = max_real; - RealT xa_min = max_real; - RealT xa_max = -max_real; - RealT ya_min = max_real; - RealT ya_max = -max_real; - RealT xb_min = max_real; - RealT xb_max = -max_real; - RealT yb_min = max_real; - RealT yb_max = -max_real; + RealT realt_max = std::numeric_limits::max(); + + // clang-format off + RealT xorg = realt_max; + RealT yorg = realt_max; + RealT xa_min = realt_max; + RealT xa_max = -realt_max; + RealT ya_min = realt_max; + RealT ya_max = -realt_max; + RealT xb_min = realt_max; + RealT xb_max = -realt_max; + RealT yb_min = realt_max; + RealT yb_max = -realt_max; + // clang-format on RealT qy = 0.0;