From 7c556606c3574bd3e737ef700fe51b105a028976 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Sun, 15 Feb 2026 23:28:13 +0100 Subject: [PATCH 1/3] try local newell --- .../mesh/utilities/ComputationalGeometry.hpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp index 61de4181531..b04c3b89fed 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp @@ -245,14 +245,17 @@ real64 centroid_3DPolygon( arraySlice1d< localIndex const > const pointsIndices, GEOS_ERROR_IF_LT( numberOfPoints, 2 ); - real64 current[ 3 ], next[ 3 ], crossProduct[ 3 ]; + real64 current[ 3 ], next[ 3 ], origin[ 3 ], crossProduct[ 3 ]; LvArray::tensorOps::copy< 3 >( next, points[ pointsIndices[ numberOfPoints - 1 ] ] ); + LvArray::tensorOps::copy< 3>( origin, points[ pointsIndices[ 0 ]] ); - for( localIndex a=0; a( current, next ); + LvArray::tensorOps::scaledAdd<3>(current, origin, -1.); LvArray::tensorOps::copy< 3 >( next, points[ pointsIndices[ a ] ] ); + LvArray::tensorOps::scaledAdd<3>(next, origin, -1.); LvArray::tensorOps::crossProduct( crossProduct, current, next ); From e1996b1a4db9961003485a3bf7af471be29d7820 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Mon, 16 Feb 2026 14:28:19 +0100 Subject: [PATCH 2/3] Local Newell / debug --- .../unitTests/testComputationalGeometry.cpp | 67 ++++++++++++++++ .../mesh/utilities/ComputationalGeometry.hpp | 76 ++++++++++++++++++- 2 files changed, 139 insertions(+), 4 deletions(-) diff --git a/src/coreComponents/mesh/unitTests/testComputationalGeometry.cpp b/src/coreComponents/mesh/unitTests/testComputationalGeometry.cpp index 403d5f25483..8b2f6998cfd 100644 --- a/src/coreComponents/mesh/unitTests/testComputationalGeometry.cpp +++ b/src/coreComponents/mesh/unitTests/testComputationalGeometry.cpp @@ -17,6 +17,7 @@ * @file testComputationalGeometry.cpp */ #include "../utilities/ComputationalGeometry.hpp" +#include #include #include @@ -70,4 +71,70 @@ TEST( testComputationalGeometry, checkCentroid3DPolygon ) } } +TEST( testComputationalGeometry, checkHighAspectRatio) +{ + constexpr localIndex numNodes = 3; + + array2d< real64, nodes::REFERENCE_POSITION_PERM > points; + points.resize(numNodes, 3); + array1d< localIndex > indices(3); + indices(0) = 0; + indices(1) = 1; + indices(2) = 2; + + //mimic a perfect triangle but with huge offset + constexpr real64 xa = 1.e4, ya = 1.e-2, za = 1.e0; + constexpr real64 xb = 1.e-1, yb = 1.e2, zb = -1.e-2; + constexpr real64 TOL = 1e-6; + + points(0,0) = 1.e6; + points(0,1) = 1.e6; + points(0,2) = 1.e6; + + points(1,0) = 1.e6 + xa; + points(1,1) = 1.e6 + ya; + points(1,2) = 1.e6 + za; + + points(2,0) = 1.e6 + xb; + points(2,1) = 1.e6 + yb; + points(2,2) = 1.e6 + zb; + + // std::cout << "Check length " << a << " " << std::sqrt(c*c + dz*dz) << std::endl; + + real64 faceCenter[ 3 ], faceNormal[ 3 ]; + real64 faceCenterOld[ 3 ], faceNormalOld[ 3 ]; + auto faceArea = computationalGeometry::centroid_3DPolygon(indices.toSliceConst(), points.toViewConst(), faceCenter, faceNormal, TOL); + auto faceAreaOld = computationalGeometry::centroid_3DPolygon_Old(indices.toSliceConst(), points.toViewConst(), faceCenterOld, faceNormalOld, TOL); + + constexpr real64 EXPECTED_AREA = 0.5*std::sqrt(LvArray::math::square(xa*yb-ya*xb) + + LvArray::math::square(ya*zb-za*yb) + + LvArray::math::square(za*xb-xa*zb) ); + + constexpr real64 EXPECTED_CENTER[3] = { 1.e6 + (xa+xb)/3 , 1.e6 + (ya+yb)/3. , 1.e6 + (za+zb)/3}; + + constexpr real64 EXPECTED_NORMAL[3] = { (ya*zb-za*yb)/2/EXPECTED_AREA, (za*xb-xa*zb)/2/EXPECTED_AREA, (xa*yb-ya*xb)/2/EXPECTED_AREA}; + std::cout << "AREA " << EXPECTED_AREA << " -- " << faceArea << "-- " << faceAreaOld << std::endl; + // + std::cout << "CENTER " << EXPECTED_CENTER[0] << " -- " << faceCenter[0] << " -- " << faceCenterOld[0]<< std::endl; + std::cout << "CENTER " << EXPECTED_CENTER[1] << " -- " << faceCenter[1] << " -- " << faceCenterOld[1]<< std::endl; + std::cout << "CENTER " << EXPECTED_CENTER[2] << " -- " << faceCenter[2] << " -- " << faceCenterOld[2]<< std::endl; + // + std::cout << "NORMAL " << EXPECTED_NORMAL[0] << " -- " << faceNormal[0] << " -- " << faceNormalOld[0] <( diameter ); } +template< typename CENTER_TYPE, typename NORMAL_TYPE > +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +real64 centroid_3DPolygon_Old( arraySlice1d< localIndex const > const pointsIndices, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & points, + CENTER_TYPE && center, + NORMAL_TYPE && normal, + real64 const areaTolerance = 0.0 ) +{ + real64 area = 0.0; + LvArray::tensorOps::fill< 3 >( center, 0 ); + LvArray::tensorOps::fill< 3 >( normal, 0 ); + + localIndex const numberOfPoints = pointsIndices.size(); + + std::cout << " Entering old impl \n"; + GEOS_ERROR_IF_LT( numberOfPoints, 2 ); + + real64 current[ 3 ], next[ 3 ], crossProduct[ 3 ]; + + LvArray::tensorOps::copy< 3 >( next, points[ pointsIndices[ numberOfPoints - 1 ] ] ); + + for( localIndex a=0; a( current, next ); + std::cout << "next[init-00] :: " << next[0] << " " << next[1] << " " << next[2] << std::endl; + LvArray::tensorOps::copy< 3 >( next, points[ pointsIndices[ a ] ] ); + std::cout << "next[init-00] :: " << next[0] << " " << next[1] << " " << next[2] << std::endl; + + LvArray::tensorOps::crossProduct( crossProduct, current, next ); + + LvArray::tensorOps::add< 3 >( normal, crossProduct ); + LvArray::tensorOps::add< 3 >( center, next ); + } + + area = LvArray::tensorOps::l2Norm< 3 >( normal ); + LvArray::tensorOps::scale< 3 >( center, 1.0 / numberOfPoints ); + + if( area > areaTolerance ) + { + LvArray::tensorOps::normalize< 3 >( normal ); + area *= 0.5; + } + else if( area < -areaTolerance ) + { + for( localIndex a=0; a const pointsIndices, real64 current[ 3 ], next[ 3 ], origin[ 3 ], crossProduct[ 3 ]; LvArray::tensorOps::copy< 3 >( next, points[ pointsIndices[ numberOfPoints - 1 ] ] ); + std::cout << "next[init] :: " << next[0] << " " << next[1] << " " << next[2] << std::endl; LvArray::tensorOps::copy< 3>( origin, points[ pointsIndices[ 0 ]] ); + std::cout << "origin[init] :: " << origin[0] << " " << origin[1] << " " << origin[2] << std::endl; - for( localIndex a=1; a( current, next ); + std::cout << "-----\n\t iter : " << a << std::endl; + LvArray::tensorOps::copy< 3 >( current, points[ pointsIndices[ a++ ]] ); LvArray::tensorOps::scaledAdd<3>(current, origin, -1.); - LvArray::tensorOps::copy< 3 >( next, points[ pointsIndices[ a ] ] ); + std::cout << "current :: " << current[0] << " " << current[1] << " " << current[2] << std::endl; + LvArray::tensorOps::copy< 3 >( next, points[ pointsIndices[ a % numberOfPoints ] ] ); + // GEOS_LOG_RANK(GEOS_FMT("next: {}, {}, {}", next[0], next[1], next[2])); LvArray::tensorOps::scaledAdd<3>(next, origin, -1.); + std::cout << "next :: " << next[0] << " " << next[1] << " " << next[2] << std::endl; LvArray::tensorOps::crossProduct( crossProduct, current, next ); + // GEOS_LOG_RANK(GEOS_FMT("crossProduct: {}", crossProduct)); + std::cout << "CrossProduct :: " << crossProduct[0] << " " << crossProduct[1] << " " << crossProduct[2] << std::endl; LvArray::tensorOps::add< 3 >( normal, crossProduct ); LvArray::tensorOps::add< 3 >( center, next ); @@ -265,6 +332,7 @@ real64 centroid_3DPolygon( arraySlice1d< localIndex const > const pointsIndices, area = LvArray::tensorOps::l2Norm< 3 >( normal ); LvArray::tensorOps::scale< 3 >( center, 1.0 / numberOfPoints ); + LvArray::tensorOps::scaledAdd< 3 >( center, origin, 1. ); if( area > areaTolerance ) { From 2129e1f83da944b55f4643cb742e2411f16e5dd4 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Mon, 16 Feb 2026 14:32:01 +0100 Subject: [PATCH 3/3] clean up --- .../unitTests/testComputationalGeometry.cpp | 12 ---- .../mesh/utilities/ComputationalGeometry.hpp | 67 ------------------- 2 files changed, 79 deletions(-) diff --git a/src/coreComponents/mesh/unitTests/testComputationalGeometry.cpp b/src/coreComponents/mesh/unitTests/testComputationalGeometry.cpp index 8b2f6998cfd..f8134983b95 100644 --- a/src/coreComponents/mesh/unitTests/testComputationalGeometry.cpp +++ b/src/coreComponents/mesh/unitTests/testComputationalGeometry.cpp @@ -99,12 +99,9 @@ TEST( testComputationalGeometry, checkHighAspectRatio) points(2,1) = 1.e6 + yb; points(2,2) = 1.e6 + zb; - // std::cout << "Check length " << a << " " << std::sqrt(c*c + dz*dz) << std::endl; - real64 faceCenter[ 3 ], faceNormal[ 3 ]; real64 faceCenterOld[ 3 ], faceNormalOld[ 3 ]; auto faceArea = computationalGeometry::centroid_3DPolygon(indices.toSliceConst(), points.toViewConst(), faceCenter, faceNormal, TOL); - auto faceAreaOld = computationalGeometry::centroid_3DPolygon_Old(indices.toSliceConst(), points.toViewConst(), faceCenterOld, faceNormalOld, TOL); constexpr real64 EXPECTED_AREA = 0.5*std::sqrt(LvArray::math::square(xa*yb-ya*xb) + LvArray::math::square(ya*zb-za*yb) + @@ -113,15 +110,6 @@ TEST( testComputationalGeometry, checkHighAspectRatio) constexpr real64 EXPECTED_CENTER[3] = { 1.e6 + (xa+xb)/3 , 1.e6 + (ya+yb)/3. , 1.e6 + (za+zb)/3}; constexpr real64 EXPECTED_NORMAL[3] = { (ya*zb-za*yb)/2/EXPECTED_AREA, (za*xb-xa*zb)/2/EXPECTED_AREA, (xa*yb-ya*xb)/2/EXPECTED_AREA}; - std::cout << "AREA " << EXPECTED_AREA << " -- " << faceArea << "-- " << faceAreaOld << std::endl; - // - std::cout << "CENTER " << EXPECTED_CENTER[0] << " -- " << faceCenter[0] << " -- " << faceCenterOld[0]<< std::endl; - std::cout << "CENTER " << EXPECTED_CENTER[1] << " -- " << faceCenter[1] << " -- " << faceCenterOld[1]<< std::endl; - std::cout << "CENTER " << EXPECTED_CENTER[2] << " -- " << faceCenter[2] << " -- " << faceCenterOld[2]<< std::endl; - // - std::cout << "NORMAL " << EXPECTED_NORMAL[0] << " -- " << faceNormal[0] << " -- " << faceNormalOld[0] <( diameter ); } -template< typename CENTER_TYPE, typename NORMAL_TYPE > -GEOS_HOST_DEVICE -GEOS_FORCE_INLINE -real64 centroid_3DPolygon_Old( arraySlice1d< localIndex const > const pointsIndices, - arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & points, - CENTER_TYPE && center, - NORMAL_TYPE && normal, - real64 const areaTolerance = 0.0 ) -{ - real64 area = 0.0; - LvArray::tensorOps::fill< 3 >( center, 0 ); - LvArray::tensorOps::fill< 3 >( normal, 0 ); - - localIndex const numberOfPoints = pointsIndices.size(); - - std::cout << " Entering old impl \n"; - GEOS_ERROR_IF_LT( numberOfPoints, 2 ); - - real64 current[ 3 ], next[ 3 ], crossProduct[ 3 ]; - - LvArray::tensorOps::copy< 3 >( next, points[ pointsIndices[ numberOfPoints - 1 ] ] ); - - for( localIndex a=0; a( current, next ); - std::cout << "next[init-00] :: " << next[0] << " " << next[1] << " " << next[2] << std::endl; - LvArray::tensorOps::copy< 3 >( next, points[ pointsIndices[ a ] ] ); - std::cout << "next[init-00] :: " << next[0] << " " << next[1] << " " << next[2] << std::endl; - - LvArray::tensorOps::crossProduct( crossProduct, current, next ); - - LvArray::tensorOps::add< 3 >( normal, crossProduct ); - LvArray::tensorOps::add< 3 >( center, next ); - } - - area = LvArray::tensorOps::l2Norm< 3 >( normal ); - LvArray::tensorOps::scale< 3 >( center, 1.0 / numberOfPoints ); - - if( area > areaTolerance ) - { - LvArray::tensorOps::normalize< 3 >( normal ); - area *= 0.5; - } - else if( area < -areaTolerance ) - { - for( localIndex a=0; a const pointsIndices, real64 current[ 3 ], next[ 3 ], origin[ 3 ], crossProduct[ 3 ]; LvArray::tensorOps::copy< 3 >( next, points[ pointsIndices[ numberOfPoints - 1 ] ] ); - std::cout << "next[init] :: " << next[0] << " " << next[1] << " " << next[2] << std::endl; LvArray::tensorOps::copy< 3>( origin, points[ pointsIndices[ 0 ]] ); - std::cout << "origin[init] :: " << origin[0] << " " << origin[1] << " " << origin[2] << std::endl; for( localIndex a=0; a( current, points[ pointsIndices[ a++ ]] ); LvArray::tensorOps::scaledAdd<3>(current, origin, -1.); - std::cout << "current :: " << current[0] << " " << current[1] << " " << current[2] << std::endl; LvArray::tensorOps::copy< 3 >( next, points[ pointsIndices[ a % numberOfPoints ] ] ); - // GEOS_LOG_RANK(GEOS_FMT("next: {}, {}, {}", next[0], next[1], next[2])); LvArray::tensorOps::scaledAdd<3>(next, origin, -1.); - std::cout << "next :: " << next[0] << " " << next[1] << " " << next[2] << std::endl; LvArray::tensorOps::crossProduct( crossProduct, current, next ); - // GEOS_LOG_RANK(GEOS_FMT("crossProduct: {}", crossProduct)); - std::cout << "CrossProduct :: " << crossProduct[0] << " " << crossProduct[1] << " " << crossProduct[2] << std::endl; LvArray::tensorOps::add< 3 >( normal, crossProduct ); LvArray::tensorOps::add< 3 >( center, next );