refact: Contact Mech - dispatch scaling wrt to fracture element type

src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp

@@ -41,6 +41,7 @@
 #include "finiteElement/FiniteElementDiscretization.hpp"
 #include "mesh/DomainPartition.hpp"
 
+#include <cmath>
 #include <stdio.h>
 
 #if defined( GEOS_USE_CUDA )
@@ -105,6 +106,11 @@ SolidMechanicsAugmentedLagrangianContact::SolidMechanicsAugmentedLagrangianConta
     setApplyDefaultValue( 5.e-02 ).
     setDescription( "Tolerance for the sliding check" );
 
+  registerWrapper( viewKeyStruct::isAnisotropicString(), &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();
 
@@ -1932,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 +1962,11 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio
     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 )
     {
       if( subRegion.hasField< contact::traction >() )
@@ -1975,6 +1987,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(), [=, this] ( localIndex const kfe )
         {
 
@@ -1992,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];
@@ -2006,12 +2023,56 @@ 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 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];
+                  // 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() );
+                }
+
+
+              }
+
+              // 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
+                  ? 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 ] = 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 );
               }
 
               averageYoungModulus += 0.5*E;

src/coreComponents/physicsSolvers/solidMechanics/contact/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
@@ -257,9 +250,8 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase
    * 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,6 +294,8 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase
 
     constexpr static char const * tolTauLimitString() { return "tolTauLimit"; }
 
+    constexpr static char const * isAnisotropicString() { return "anisotropic"; }
+
   };
 
   /// Tolerance for the sliding check: the tangential traction must exceed (1 + m_slidingCheckTolerance) * t_lim to activate the sliding
@@ -330,6 +324,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 */