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::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();
 
@@ -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(), [=] ( 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,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 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];
+                }
+
+
+              }
+
+              // 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;
+
+                //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 );
               }
 
               averageYoungModulus += 0.5*E;

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

@@ -302,6 +302,8 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase
 
     constexpr static char const * tolTauLimitString() { return "tolTauLimit"; }
 
+    constexpr static char const * anisotropicString() { return "anisotropic"; }
+
   };
 
   /// Tolerance for the sliding check: the tangential traction must exceed (1 + m_slidingCheckTolerance) * t_lim to activate the sliding
@@ -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 */