From 2bf033a2ac445a6e05815133ac7110887ae08ad9 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Tue, 2 Jun 2026 15:08:06 +0200 Subject: [PATCH 1/9] right tetra scaling --- .../contact/SolidMechanicsAugmentedLagrangianContact.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp index ff850c0d744..1c25458fc00 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -41,6 +41,7 @@ #include "finiteElement/FiniteElementDiscretization.hpp" #include "mesh/DomainPartition.hpp" +#include #include #if defined( GEOS_USE_CUDA ) @@ -2006,12 +2007,14 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio real64 const nu = ( 3.0 * K - 2.0 * G ) / ( 2.0 * ( 3.0 * K + G ) ); real64 const M = K + 4.0 / 3.0 * G; - real64 const charLength = pow( volume, 1.0 / 3.0 ); + // real64 const charLength = pow( volume, 1.0 / 3.0 ); + real64 const charLength = pow( 6*std::sqrt(2)*volume, 1.0 / 3.0 ); // Combine E and nu to obtain a stiffness approximation (like it was an hexahedron) for( localIndex j = 0; j < 3; ++j ) { - stiffDiagApprox[ i ][ j ] = E / ( ( 1.0 + nu )*( 1.0 - 2.0*nu ) ) * 4.0 / 9.0 * ( 2.0 - 3.0 * nu ) * charLength; + // stiffDiagApprox[ i ][ j ] = E / ( ( 1.0 + nu )*( 1.0 - 2.0*nu ) ) * 4.0 / 9.0 * ( 2.0 - 3.0 * nu ) * charLength; + stiffDiagApprox[ i ][ j ] = E / ( ( 1.0 + nu )*( 1.0 - 2.0*nu ) ) * ( 2.0 - 3.0 * nu ) * charLength; } averageYoungModulus += 0.5*E; From 7a3b5c0c8fbf9021089fc4d084fd77e6573f6628 Mon Sep 17 00:00:00 2001 From: jafranc Date: Wed, 3 Jun 2026 15:11:05 +0200 Subject: [PATCH 2/9] rough dispatch --- ...SolidMechanicsAugmentedLagrangianContact.cpp | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp index 1c25458fc00..8e3231a23e4 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -1976,6 +1976,10 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + // Triangle face elements bound tetrahedral bulk elements; quadrilateral faces bound hexahedral ones. + // The charLength formula differs between the two element types. + bool const isTriangle = subRegion.size() > 0 && subRegion.getElementType( 0 ) == ElementType::Triangle; + forAll< parallelHostPolicy >( subRegion.size(), [=] ( localIndex const kfe ) { @@ -2007,14 +2011,19 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio real64 const nu = ( 3.0 * K - 2.0 * G ) / ( 2.0 * ( 3.0 * K + G ) ); real64 const M = K + 4.0 / 3.0 * G; - // real64 const charLength = pow( volume, 1.0 / 3.0 ); - real64 const charLength = pow( 6*std::sqrt(2)*volume, 1.0 / 3.0 ); + // For tetrahedra (triangle faces): charLength = edge = (6*sqrt(2)*V)^(1/3) + // For hexahedra (quadrilateral faces): charLength = (V)^(1/3) + real64 const charLength = isTriangle + ? pow( 6 * std::sqrt( 2 ) * volume, 1.0 / 3.0 ) + : pow( volume, 1.0 / 3.0 ); // Combine E and nu to obtain a stiffness approximation (like it was an hexahedron) for( localIndex j = 0; j < 3; ++j ) { - // stiffDiagApprox[ i ][ j ] = E / ( ( 1.0 + nu )*( 1.0 - 2.0*nu ) ) * 4.0 / 9.0 * ( 2.0 - 3.0 * nu ) * charLength; - stiffDiagApprox[ i ][ j ] = E / ( ( 1.0 + nu )*( 1.0 - 2.0*nu ) ) * ( 2.0 - 3.0 * nu ) * charLength; + + stiffDiagApprox[ i ][ j ] = isTriangle + ? E / ( ( 1.0 + nu )*( 1.0 - 2.0*nu ) ) * ( 2.0 - 3.0 * nu ) * charLength + : E / ( ( 1.0 + nu )*( 1.0 - 2.0*nu ) ) * 4.0 / 9.0 * ( 2.0 - 3.0 * nu ) * charLength; } averageYoungModulus += 0.5*E; From 85b5f1fcbe6c3459b6aacebb58b1d7339c0628c8 Mon Sep 17 00:00:00 2001 From: jafranc Date: Fri, 19 Jun 2026 08:34:24 -0500 Subject: [PATCH 3/9] restore anisotropic box scaling --- ...lidMechanicsAugmentedLagrangianContact.cpp | 55 +++++++++++++++++-- ...lidMechanicsAugmentedLagrangianContact.hpp | 5 ++ 2 files changed, 56 insertions(+), 4 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp index 8e3231a23e4..99e3eca70a6 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -106,6 +106,11 @@ SolidMechanicsAugmentedLagrangianContact::SolidMechanicsAugmentedLagrangianConta setApplyDefaultValue( 5.e-02 ). setDescription( "Tolerance for the sliding check" ); + registerWrapper( viewKeyStruct::symmetricString(), &m_isAnisotropic ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 1 ). + setDescription( "Flag to use anisotropic scaling in tolerances and penalties computations" ); + // Set the default linear solver parameters LinearSolverParameters & linSolParams = m_linearSolverParameters.get(); @@ -1933,6 +1938,7 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio string_array const & ) { FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); ElementRegionManager & elemManager = mesh.getElemManager(); // Get the "face to element" map (valid for the entire mesh) @@ -1955,6 +1961,11 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio using NodeMapViewType = arrayView2d< localIndex const, cells::NODE_MAP_USD >; ElementRegionManager::ElementViewAccessor< NodeMapViewType > const elemToNode = elemManager.constructViewAccessor< CellElementSubRegion::NodeMapType, NodeMapViewType >( ElementSubRegionBase::viewKeyStruct::nodeListString() ); + + ElementRegionManager::ElementViewConst< NodeMapViewType > const elemToNodeView = elemToNode.toNestedViewConst(); + + // Get the coordinates for all nodes + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodePosition = nodeManager.referencePosition(); elemManager.forElementSubRegions< FaceElementSubRegion >( [&]( FaceElementSubRegion & subRegion ) { @@ -1997,6 +2008,7 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio for( localIndex i = 0; i < 2; ++i ) { + localIndex const faceIndex = elemsToFaces[kfe][i]; localIndex const er = faceToElemRegion[faceIndex][0]; localIndex const esr = faceToElemSubRegion[faceIndex][0]; @@ -2011,19 +2023,54 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio real64 const nu = ( 3.0 * K - 2.0 * G ) / ( 2.0 * ( 3.0 * K + G ) ); real64 const M = K + 4.0 / 3.0 * G; + real64 bbox[3]{}; + + if(m_isAnisotropic) + { + NodeMapViewType const & cellElemsToNodes = elemToNodeView[er][esr]; + localIndex const numNodesPerElem = cellElemsToNodes.size( 1 ); + + real64 maxSize[3]; + real64 minSize[3]; + for( localIndex j = 0; j < 3; ++j ) + { + maxSize[j] = nodePosition[cellElemsToNodes[ei][0]][j]; + minSize[j] = nodePosition[cellElemsToNodes[ei][0]][j]; + } + + for( localIndex a = 1; a < numNodesPerElem; ++a ) + { + for( localIndex j = 0; j < 3; ++j ) + { + maxSize[j] = fmax( maxSize[j], nodePosition[cellElemsToNodes[ei][a]][j] ); + minSize[j] = fmin( minSize[j], nodePosition[cellElemsToNodes[ei][a]][j] ); + } + } + + for( localIndex j = 0; j < 3; ++j ) + { + bbox[j] = maxSize[j] - minSize[j]; + } + + + } + + // For anisotropic factor: XYZ-aligned bbox length // For tetrahedra (triangle faces): charLength = edge = (6*sqrt(2)*V)^(1/3) // For hexahedra (quadrilateral faces): charLength = (V)^(1/3) - real64 const charLength = isTriangle + real64 const charLength = m_isAnisotropic ? bbox[0] : ( isTriangle ? pow( 6 * std::sqrt( 2 ) * volume, 1.0 / 3.0 ) - : pow( volume, 1.0 / 3.0 ); + : pow( volume, 1.0 / 3.0 ) ); // Combine E and nu to obtain a stiffness approximation (like it was an hexahedron) for( localIndex j = 0; j < 3; ++j ) { - stiffDiagApprox[ i ][ j ] = isTriangle + //TODO (jafranc) once stabilized, get rid of this ugly ternary + stiffDiagApprox[ i ][ j ] = m_isAnisotropic ? E / ( ( 1.0 + nu )*( 1.0 - 2.0*nu ) ) * 4.0 / 9.0 * ( 2.0 - 3.0 * nu ) * volume / bbox[j] / bbox[j] + : ( isTriangle ? E / ( ( 1.0 + nu )*( 1.0 - 2.0*nu ) ) * ( 2.0 - 3.0 * nu ) * charLength - : E / ( ( 1.0 + nu )*( 1.0 - 2.0*nu ) ) * 4.0 / 9.0 * ( 2.0 - 3.0 * nu ) * charLength; + : E / ( ( 1.0 + nu )*( 1.0 - 2.0*nu ) ) * 4.0 / 9.0 * ( 2.0 - 3.0 * nu ) * charLength ); } averageYoungModulus += 0.5*E; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp index 02af0550280..6da88bb436c 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp @@ -301,6 +301,8 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase constexpr static char const * tolNormalTracFacString() { return "tolNormalTrac"; } constexpr static char const * tolTauLimitString() { return "tolTauLimit"; } + + constexpr static char const * anisotropicString() { return "anisotropic"; } }; @@ -330,6 +332,9 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase /// Factor to adjust the tolerance for normal traction real64 m_tolNormalTracFac = 0.5; + /// Flag for anisotropic scaling in Tolerances and Penalties + int m_isAnisotropic = 1; + }; } /* namespace geos */ From e0c45f4ba90a6d66c70335beffa5100b0a1443a6 Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Thu, 27 Aug 2026 11:30:00 -0500 Subject: [PATCH 4/9] uncrustify --- ...lidMechanicsAugmentedLagrangianContact.cpp | 44 +++++++++---------- ...lidMechanicsAugmentedLagrangianContact.hpp | 2 +- 2 files changed, 23 insertions(+), 23 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp index 70d5528d5e8..2e4668620ee 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -110,7 +110,7 @@ SolidMechanicsAugmentedLagrangianContact::SolidMechanicsAugmentedLagrangianConta setInputFlag( InputFlags::OPTIONAL ). setApplyDefaultValue( 1 ). setDescription( "Flag to use anisotropic scaling in tolerances and penalties computations" ); - + // Set the default linear solver parameters LinearSolverParameters & linSolParams = m_linearSolverParameters.get(); @@ -1961,10 +1961,10 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio using NodeMapViewType = arrayView2d< localIndex const, cells::NODE_MAP_USD >; ElementRegionManager::ElementViewAccessor< NodeMapViewType > const elemToNode = elemManager.constructViewAccessor< CellElementSubRegion::NodeMapType, NodeMapViewType >( ElementSubRegionBase::viewKeyStruct::nodeListString() ); - + ElementRegionManager::ElementViewConst< NodeMapViewType > const elemToNodeView = elemToNode.toNestedViewConst(); - // Get the coordinates for all nodes + // Get the coordinates for all nodes arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodePosition = nodeManager.referencePosition(); elemManager.forElementSubRegions< FaceElementSubRegion >( [&]( FaceElementSubRegion & subRegion ) @@ -2024,33 +2024,33 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio real64 const M = K + 4.0 / 3.0 * G; real64 bbox[3]{}; - - if(m_isAnisotropic) + + if( m_isAnisotropic ) { NodeMapViewType const & cellElemsToNodes = elemToNodeView[er][esr]; localIndex const numNodesPerElem = cellElemsToNodes.size( 1 ); real64 maxSize[3]; real64 minSize[3]; - for( localIndex j = 0; j < 3; ++j ) - { - maxSize[j] = nodePosition[cellElemsToNodes[ei][0]][j]; - minSize[j] = nodePosition[cellElemsToNodes[ei][0]][j]; - } - - for( localIndex a = 1; a < numNodesPerElem; ++a ) - { for( localIndex j = 0; j < 3; ++j ) { - maxSize[j] = fmax( maxSize[j], nodePosition[cellElemsToNodes[ei][a]][j] ); - minSize[j] = fmin( minSize[j], nodePosition[cellElemsToNodes[ei][a]][j] ); + maxSize[j] = nodePosition[cellElemsToNodes[ei][0]][j]; + minSize[j] = nodePosition[cellElemsToNodes[ei][0]][j]; } - } - for( localIndex j = 0; j < 3; ++j ) - { - bbox[j] = maxSize[j] - minSize[j]; - } + for( localIndex a = 1; a < numNodesPerElem; ++a ) + { + for( localIndex j = 0; j < 3; ++j ) + { + maxSize[j] = fmax( maxSize[j], nodePosition[cellElemsToNodes[ei][a]][j] ); + minSize[j] = fmin( minSize[j], nodePosition[cellElemsToNodes[ei][a]][j] ); + } + } + + for( localIndex j = 0; j < 3; ++j ) + { + bbox[j] = maxSize[j] - minSize[j]; + } } @@ -2058,14 +2058,14 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio // For anisotropic factor: XYZ-aligned bbox length // For tetrahedra (triangle faces): charLength = edge = (6*sqrt(2)*V)^(1/3) // For hexahedra (quadrilateral faces): charLength = (V)^(1/3) - real64 const charLength = m_isAnisotropic ? bbox[0] : ( isTriangle + real64 const charLength = m_isAnisotropic ? bbox[0] : ( isTriangle ? pow( 6 * std::sqrt( 2 ) * volume, 1.0 / 3.0 ) : pow( volume, 1.0 / 3.0 ) ); // Combine E and nu to obtain a stiffness approximation (like it was an hexahedron) for( localIndex j = 0; j < 3; ++j ) { - + //TODO (jafranc) once stabilized, get rid of this ugly ternary stiffDiagApprox[ i ][ j ] = m_isAnisotropic ? E / ( ( 1.0 + nu )*( 1.0 - 2.0*nu ) ) * 4.0 / 9.0 * ( 2.0 - 3.0 * nu ) * volume / bbox[j] / bbox[j] : ( isTriangle diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp index 6da88bb436c..d6deeb42028 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp @@ -301,7 +301,7 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase constexpr static char const * tolNormalTracFacString() { return "tolNormalTrac"; } constexpr static char const * tolTauLimitString() { return "tolTauLimit"; } - + constexpr static char const * anisotropicString() { return "anisotropic"; } }; From 93128ee49e828d10223c58a5c99926821c02086a Mon Sep 17 00:00:00 2001 From: jafranc Date: Fri, 28 Aug 2026 09:18:07 +0200 Subject: [PATCH 5/9] adress jian's comments --- ...lidMechanicsAugmentedLagrangianContact.cpp | 8 ++-- ...lidMechanicsAugmentedLagrangianContact.hpp | 40 +++++++------------ 2 files changed, 20 insertions(+), 28 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp index 2e4668620ee..163ab4cd5d5 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -106,7 +106,7 @@ SolidMechanicsAugmentedLagrangianContact::SolidMechanicsAugmentedLagrangianConta setApplyDefaultValue( 5.e-02 ). setDescription( "Tolerance for the sliding check" ); - registerWrapper( viewKeyStruct::symmetricString(), &m_isAnisotropic ). + registerWrapper( viewKeyStruct::isAnisotropicString(), &m_isAnisotropic ). setInputFlag( InputFlags::OPTIONAL ). setApplyDefaultValue( 1 ). setDescription( "Flag to use anisotropic scaling in tolerances and penalties computations" ); @@ -765,7 +765,7 @@ void SolidMechanicsAugmentedLagrangianContact::assembleForceResidualPressureCont ElementRegionManager const & elementRegionManager = mesh.getElemManager(); elementRegionManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&] - ( localIndex const regionIndex, auto & elementSubRegion ) + ( localIndex const regionIndex, auto & elementSubRegion ) { if( elementSubRegion.template hasWrapper< string >( FlowSolverBase::viewKeyStruct::solidNamesString() ) ) { @@ -2050,6 +2050,9 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio for( localIndex j = 0; j < 3; ++j ) { bbox[j] = maxSize[j] - minSize[j]; + // Avoid division by zero in case of degenerate elements + if( bbox[j] < 1e-12 ) + GEOS_ERROR( GEOS_FMT( "SolidMechanicsAugmentedLagrangianContact::computeTolerances: degenerate element detected with zero size in direction {}", j ), getDataContext() ); } @@ -2066,7 +2069,6 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio for( localIndex j = 0; j < 3; ++j ) { - //TODO (jafranc) once stabilized, get rid of this ugly ternary stiffDiagApprox[ i ][ j ] = m_isAnisotropic ? E / ( ( 1.0 + nu )*( 1.0 - 2.0*nu ) ) * 4.0 / 9.0 * ( 2.0 - 3.0 * nu ) * volume / bbox[j] / bbox[j] : ( isTriangle ? E / ( ( 1.0 + nu )*( 1.0 - 2.0*nu ) ) * ( 2.0 - 3.0 * nu ) * charLength diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp index d6deeb42028..44d7c589aef 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp @@ -13,7 +13,7 @@ * ------------------------------------------------------------------------------------------------------------ */ -/* + /* * SolidMechanicsAugmentedLagrangianContact.hpp * */ @@ -118,8 +118,7 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase * @brief Loop over the finite element type on the fracture subregions of meshName and apply callback. * @tparam LAMBDA The callback function type * @param meshName The mesh name. - * @param lambda The callback function. Take the finite element type name and - * the list of face element of the same type. + * @param lambda The callback function. Take the finite element type name and the list of face element of the same type. */ template< typename LAMBDA > void forFiniteElementOnFractureSubRegions( string const & meshName, LAMBDA && lambda ) const @@ -143,8 +142,7 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase * @brief Loop over the finite element type on the stick fracture subregions of meshName and apply callback. * @tparam LAMBDA The callback function type * @param meshName The mesh name. - * @param lambda The callback function. Take the finite element type name and - * the list of face element of the same type. + * @param lambda The callback function. Take the finite element type name and the list of face element of the same type. */ template< typename LAMBDA > void forFiniteElementOnStickFractureSubRegions( string const & meshName, LAMBDA && lambda ) const @@ -170,8 +168,7 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase * @brief Loop over the finite element type on the slip fracture subregions of meshName and apply callback. * @tparam LAMBDA The callback function type * @param meshName The mesh name. - * @param lambda The callback function. Take the finite element type name and - * the list of face element of the same type. + * @param lambda The callback function. Take the finite element type name and the list of face element of the same type. */ template< typename LAMBDA > void forFiniteElementOnSlipFractureSubRegions( string const & meshName, LAMBDA && lambda ) const @@ -194,23 +191,20 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase } /** - * @brief Create the list of finite elements of the same type - * for each FaceElementSubRegion (Triangle or Quadrilateral) - * and of the same fracture state (Stick or Slip). + * @brief Create the list of finite elements of the same type for each FaceElementSubRegion (Triangle or Quadrilateral) and of the same + *fracture state (Stick or Slip). * @param domain The physical domain object */ void updateStickSlipList( DomainPartition const & domain ); /** - * @brief Create the list of finite elements of the same type - * for each FaceElementSubRegion (Triangle or Quadrilateral). + * @brief Create the list of finite elements of the same type for each FaceElementSubRegion (Triangle or Quadrilateral). * @param domain The physical domain object */ void createFaceTypeList( DomainPartition const & domain ); /** - * @brief Create the list of elements belonging to CellElementSubRegion - * that are enriched with the bubble basis functions + * @brief Create the list of elements belonging to CellElementSubRegion that are enriched with the bubble basis functions * @param domain The physical domain object */ void createBubbleCellList( DomainPartition & domain ) const; @@ -224,8 +218,7 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase void validateTetrahedralQuadrature( Group & meshBodies ); /** - * @brief add the number of non-zero elements induced by the coupling between - * nodal and bubble displacement. + * @brief add the number of non-zero elements induced by the coupling between nodal and bubble displacement. * @param domain the physical domain object * @param dofManager degree-of-freedom manager associated with the linear system * @param rowLengths the array containing the number of non-zero elements for each row @@ -251,15 +244,12 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase * @param domain The physical domain object * * This function computes the initial traction on each fracture element by: - * 1. Getting the stress tensor from both adjacent volume elements (one on each side of the fracture) - * 2. Computing the traction vector as t = sigma * n (where n is the face normal) for each side - * 3. Averaging the tractions from both sides - * 4. Rotating the averaged traction to the local coordinate system of the fracture - * 5. Validating the traction against the Coulomb friction law and warning if inconsistent + * 1. Getting the stress tensor from both adjacent volume elements (one on each side of the fracture) 2. Computing the traction vector as + *t = sigma * n (where n is the face normal) for each side 3. Averaging the tractions from both sides 4. Rotating the averaged traction to + *the local coordinate system of the fracture 5. Validating the traction against the Coulomb friction law and warning if inconsistent * - * This initialization ensures that the ALM traction field starts with a physically - * consistent value rather than zero, which is important for proper convergence - * when the domain is under stress. + * This initialization ensures that the ALM traction field starts with a physically consistent value rather than zero, which is important + *for proper convergence when the domain is under stress. */ void initializeTractionFromAdjacentCellStress( DomainPartition & domain ) const; @@ -302,7 +292,7 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase constexpr static char const * tolTauLimitString() { return "tolTauLimit"; } - constexpr static char const * anisotropicString() { return "anisotropic"; } + constexpr static char const * isAnisotropicString() { return "anisotropic"; } }; From 26c07780f62d2e0e8bb10e8a73a54ce7c98a0bb6 Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Fri, 28 Aug 2026 10:59:24 -0500 Subject: [PATCH 6/9] uncrustify --- .../SolidMechanicsAugmentedLagrangianContact.cpp | 2 +- .../SolidMechanicsAugmentedLagrangianContact.hpp | 11 ++++++----- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp index 48f330d681c..0baa6ca60a6 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -765,7 +765,7 @@ void SolidMechanicsAugmentedLagrangianContact::assembleForceResidualPressureCont ElementRegionManager const & elementRegionManager = mesh.getElemManager(); elementRegionManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&] - ( localIndex const regionIndex, auto & elementSubRegion ) + ( localIndex const regionIndex, auto & elementSubRegion ) { if( elementSubRegion.template hasWrapper< string >( FlowSolverBase::viewKeyStruct::solidNamesString() ) ) { diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp index 44d7c589aef..873b3b347c1 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp @@ -13,7 +13,7 @@ * ------------------------------------------------------------------------------------------------------------ */ - /* +/* * SolidMechanicsAugmentedLagrangianContact.hpp * */ @@ -192,7 +192,7 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase /** * @brief Create the list of finite elements of the same type for each FaceElementSubRegion (Triangle or Quadrilateral) and of the same - *fracture state (Stick or Slip). + * fracture state (Stick or Slip). * @param domain The physical domain object */ void updateStickSlipList( DomainPartition const & domain ); @@ -245,11 +245,12 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase * * This function computes the initial traction on each fracture element by: * 1. Getting the stress tensor from both adjacent volume elements (one on each side of the fracture) 2. Computing the traction vector as - *t = sigma * n (where n is the face normal) for each side 3. Averaging the tractions from both sides 4. Rotating the averaged traction to - *the local coordinate system of the fracture 5. Validating the traction against the Coulomb friction law and warning if inconsistent + * t = sigma * n (where n is the face normal) for each side 3. Averaging the tractions from both sides 4. Rotating the averaged traction + *to + * the local coordinate system of the fracture 5. Validating the traction against the Coulomb friction law and warning if inconsistent * * This initialization ensures that the ALM traction field starts with a physically consistent value rather than zero, which is important - *for proper convergence when the domain is under stress. + * for proper convergence when the domain is under stress. */ void initializeTractionFromAdjacentCellStress( DomainPartition & domain ) const; From a4cf0b90409b8fa779f074b2097db3a19bfa8446 Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Fri, 28 Aug 2026 11:06:12 -0500 Subject: [PATCH 7/9] wording style --- .../contact/SolidMechanicsAugmentedLagrangianContact.hpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp index 873b3b347c1..63e01cccabc 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp @@ -244,10 +244,11 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase * @param domain The physical domain object * * This function computes the initial traction on each fracture element by: - * 1. Getting the stress tensor from both adjacent volume elements (one on each side of the fracture) 2. Computing the traction vector as - * t = sigma * n (where n is the face normal) for each side 3. Averaging the tractions from both sides 4. Rotating the averaged traction - *to - * the local coordinate system of the fracture 5. Validating the traction against the Coulomb friction law and warning if inconsistent + * 1. Getting the stress tensor from both adjacent volume elements (one on each side of the fracture) + * 2. Computing the traction vector as t = sigma * n (where n is the face normal) for each side + * 3. Averaging the tractions from both sides + * 4. Rotating the averaged traction to the local coordinate system of the fracture + * 5. Validating the traction against the Coulomb friction law and warning if inconsistent * * This initialization ensures that the ALM traction field starts with a physically consistent value rather than zero, which is important * for proper convergence when the domain is under stress. From 8db0815a0d269342f330af8b40d351f9094065b9 Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Fri, 28 Aug 2026 11:11:53 -0500 Subject: [PATCH 8/9] uncrustify again --- .../contact/SolidMechanicsAugmentedLagrangianContact.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp index 63e01cccabc..405b2004d7d 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp @@ -244,10 +244,10 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase * @param domain The physical domain object * * This function computes the initial traction on each fracture element by: - * 1. Getting the stress tensor from both adjacent volume elements (one on each side of the fracture) - * 2. Computing the traction vector as t = sigma * n (where n is the face normal) for each side - * 3. Averaging the tractions from both sides - * 4. Rotating the averaged traction to the local coordinate system of the fracture + * 1. Getting the stress tensor from both adjacent volume elements (one on each side of the fracture) + * 2. Computing the traction vector as t = sigma * n (where n is the face normal) for each side + * 3. Averaging the tractions from both sides + * 4. Rotating the averaged traction to the local coordinate system of the fracture * 5. Validating the traction against the Coulomb friction law and warning if inconsistent * * This initialization ensures that the ALM traction field starts with a physically consistent value rather than zero, which is important From ed488ae3706d3fe3c3a7b4d2198a7460d98123e5 Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Tue, 1 Sep 2026 10:12:21 -0500 Subject: [PATCH 9/9] rebaseline --- .integrated_tests.yaml | 2 +- BASELINE_NOTES.md | 5 +++++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index e59c7f27847..dbbaf8edc6b 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr4096-17439-18734f3 + baseline: integratedTests/baseline_integratedTests-pr4083-17503-96f5eac allow_fail: all: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index b813fa7e566..d64a5aa9594 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -5,7 +5,12 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #4096 (2026-09-01) +===================== +Dispatch scaling wrt to fracture element type for ALM contact solver + PR #4096 (2026-08-27) +===================== Refactor the phase-field module and add new benchmarks for baseline PR #4127 (2026-08-23)