From 2bf033a2ac445a6e05815133ac7110887ae08ad9 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Tue, 2 Jun 2026 15:08:06 +0200 Subject: [PATCH 1/3] 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/3] 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/3] 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 */