feat: averaging of stress components using quadrature

src/coreComponents/common/DataLayouts.hpp

@@ -155,24 +155,24 @@ namespace cells
 /// Cell node map permutation when using cuda.
 using NODE_MAP_PERMUTATION = RAJA::PERM_JI;
 
-/// Cell strain permutation when using cuda
-using STRAIN_PERM = RAJA::PERM_JI;
+/// Cell tensor (i.e. average stress and strain) permutation when using cuda
+using RANK2_TENSOR_PERM = RAJA::PERM_JI;
 
 #else
 
 /// Cell node map permutation when not using cuda.
 using NODE_MAP_PERMUTATION = RAJA::PERM_IJ;
 
-/// Cell strain permutation when not using cuda
-using STRAIN_PERM = RAJA::PERM_IJ;
+/// Cell tensor (i.e. average stress and strain) permutation when not using cuda
+using RANK2_TENSOR_PERM = RAJA::PERM_IJ;
 
 #endif
 
 /// Cell node map unit stride dimension.
 static constexpr int NODE_MAP_USD = LvArray::typeManipulation::getStrideOneDimension( NODE_MAP_PERMUTATION {} );
 
 /// Cell strain unit stride dimension
-static constexpr int STRAIN_USD = LvArray::typeManipulation::getStrideOneDimension( STRAIN_PERM {} );
+static constexpr int RANK2_TENSOR_USD = LvArray::typeManipulation::getStrideOneDimension( RANK2_TENSOR_PERM {} );
 
 } // namespace cells
 

src/coreComponents/constitutive/solid/SolidBase.cpp

@@ -37,7 +37,7 @@ SolidBase::SolidBase( string const & name, Group * const parent ):
   string const voightLabels[6] = { "XX", "YY", "ZZ", "YZ", "XZ", "XY" };
 
   registerWrapper( viewKeyStruct::stressString(), &m_newStress ).
-    setPlotLevel( PlotLevel::LEVEL_0 ).
+    setPlotLevel( PlotLevel::NOPLOT ).
     setApplyDefaultValue( 0 ). // default to zero initial stress
     setDescription( "Current Material Stress" ).
     setDimLabels( 2, voightLabels );

src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsFields.hpp

@@ -42,9 +42,13 @@
 using arrayView2dLayoutIncrDisplacement = arrayView2d< real64, nodes::INCR_DISPLACEMENT_USD >;
 using arrayViewConst2dLayoutIncrDisplacement = arrayView2d< real64 const, nodes::INCR_DISPLACEMENT_USD >;
 
-using array2dLayoutStrain = array2d< real64, cells::STRAIN_PERM >;
-using arrayView2dLayoutStrain = arrayView2d< real64, cells::STRAIN_USD >;
-using arrayViewConst2dLayoutStrain = arrayView2d< real64 const, cells::STRAIN_USD >;
+using array2dLayoutStrain = array2d< real64, cells::RANK2_TENSOR_PERM >;
+using arrayView2dLayoutStrain = arrayView2d< real64, cells::RANK2_TENSOR_USD >;
+using arrayViewConst2dLayoutStrain = arrayView2d< real64 const, cells::RANK2_TENSOR_USD >;
+
+using array2dLayoutAvgStress = array2d< real64, cells::RANK2_TENSOR_PERM >;
+using arrayView2dLayoutAvgStress = arrayView2d< real64, cells::RANK2_TENSOR_USD >;
+using arrayViewConst2dLayoutAvgStress = arrayView2d< real64 const, cells::RANK2_TENSOR_USD >;
 
 using array2dLayoutVelocity = array2d< real64, nodes::VELOCITY_PERM >;
 using arrayView2dLayoutVelocity = arrayView2d< real64, nodes::VELOCITY_USD >;
@@ -79,13 +83,21 @@
                WRITE_AND_READ,
                "Incremental displacements for the current time step on the nodes" );
 
-DECLARE_FIELD( strain,
-               "strain",
+DECLARE_FIELD( averageStrain,
+               "averageStrain",
                array2dLayoutStrain,
                0,
                LEVEL_0,
                WRITE_AND_READ,
-               "Average strain in cell" );
+               "Quadrature averaged strain components in cell" );
+
+DECLARE_FIELD( averageStress,
+               "averageStress",
+               array2dLayoutAvgStress,
+               0,
+               LEVEL_0,
+               WRITE_AND_READ,
+               "Quadrature averaged stress components in cell" );
 
 DECLARE_FIELD( incrementalBubbleDisplacement,
                "incrementalBubbleDisplacement",

src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp

@@ -158,7 +158,8 @@
     {
       setConstitutiveNamesCallSuper( subRegion );
 
-      subRegion.registerField< solidMechanics::strain >( getName() ).reference().resizeDimension< 1 >( 6 );
+      subRegion.registerField< solidMechanics::averageStrain >( getName() ).reference().resizeDimension< 1 >( 6 );
+      subRegion.registerField< solidMechanics::averageStress >( getName() ).reference().resizeDimension< 1 >( 6 );
     } );
 
     NodeManager & nodes = meshLevel.getNodeManager();
@@ -944,19 +945,26 @@
       SolidBase & constitutiveRelation = getConstitutiveModel< SolidBase >( subRegion, solidMaterialName );
       constitutiveRelation.saveConvergedState();
 
-      solidMechanics::arrayView2dLayoutStrain strain = subRegion.getField< solidMechanics::strain >();
+      arrayView3d < real64 const, solid::STRESS_USD > const stress = constitutiveRelation.getStress();
+
+      solidMechanics::arrayView2dLayoutStrain avgStrain = subRegion.getField< solidMechanics::averageStrain >();
+      solidMechanics::arrayView2dLayoutAvgStress avgStress = subRegion.getField< solidMechanics::averageStress >();
 
       finiteElement::FiniteElementBase & subRegionFE = subRegion.template getReference< finiteElement::FiniteElementBase >( this->getDiscretizationName());
       finiteElement::FiniteElementDispatchHandler< BASE_FE_TYPES >::dispatch3D( subRegionFE, [&] ( auto const finiteElement )
       {
         using FE_TYPE = decltype( finiteElement );
-        AverageStrainOverQuadraturePointsKernelFactory::createAndLaunch< CellElementSubRegion, FE_TYPE, parallelDevicePolicy<> >( nodeManager,
+        AverageStressStrainOverQuadraturePointsKernelFactory::createAndLaunch< CellElementSubRegion, FE_TYPE, parallelDevicePolicy<> >( nodeManager,
                                                                                                                                   mesh.getEdgeManager(),
                                                                                                                                   mesh.getFaceManager(),
                                                                                                                                   subRegion,
                                                                                                                                   finiteElement,
                                                                                                                                   disp,
-                                                                                                                                  strain );
+                                                                                                                                  avgStrain,
+                                                                                                                                  stress,
+                                                                                                                                  avgStress );
+
+
       } );
 
 

src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp

@@ -23,7 +23,7 @@
 #include "codingUtilities/EnumStrings.hpp"
 #include "common/TimingMacros.hpp"
 #include "kernels/SolidMechanicsLagrangianFEMKernels.hpp"
-#include "kernels/StrainHelper.hpp"
+#include "kernels/StressStrainAverageKernels.hpp"
 #include "mesh/MeshForLoopInterface.hpp"
 #include "mesh/mpiCommunications/CommunicationTools.hpp"
 #include "mesh/mpiCommunications/MPI_iCommData.hpp"

src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp → src/coreComponents/physicsSolvers/solidMechanics/kernels/StressStrainAverageKernels.hpp

@@ -14,11 +14,11 @@
  */
 
 /**
- * @file StrainHelper.hpp
+ * @file StressStrainAverageKernels.hpp
  */
 
-#ifndef GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_STRAINHELPER_HPP_
-#define GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_STRAINHELPER_HPP_
+#ifndef GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_STRESSSTRAINAVERAGEKERNELS_HPP_
+#define GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_STRESSSTRAINAVERAGEKERNELS_HPP_
 
 #include "common/DataTypes.hpp"
 #include "common/GEOS_RAJA_Interface.hpp"
@@ -29,14 +29,16 @@
 
 namespace geos
 {
+
+
 /**
- * @class AverageStrainOverQuadraturePoints
+ * @class AverageStressStrainOverQuadraturePoints
  * @tparam SUBREGION_TYPE the subRegion type
  * @tparam FE_TYPE the finite element type
  */
 template< typename SUBREGION_TYPE,
           typename FE_TYPE >
-class AverageStrainOverQuadraturePoints :
+class AverageStressStrainOverQuadraturePoints :
   public AverageOverQuadraturePointsBase< SUBREGION_TYPE,
                                           FE_TYPE >
 {
@@ -59,21 +61,27 @@
    * @param finiteElementSpace the finite element space
    * @param displacement the displacement solution field
    * @param avgStrain the strain averaged over quadrature points
+   * @param stress the stress solution field
+   * @param avgStress the stress averaged over quadrature points
    */
-  AverageStrainOverQuadraturePoints( NodeManager & nodeManager,
+  AverageStressStrainOverQuadraturePoints( NodeManager & nodeManager,
                                      EdgeManager const & edgeManager,
                                      FaceManager const & faceManager,
                                      SUBREGION_TYPE const & elementSubRegion,
                                      FE_TYPE const & finiteElementSpace,
                                      fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement,
-                                     fields::solidMechanics::arrayView2dLayoutStrain const avgStrain ):
+                                     fields::solidMechanics::arrayView2dLayoutStrain const avgStrain,
+                                     arrayView3d< real64 const, solid::STRESS_USD > const stress,
+                                     fields::solidMechanics::arrayView2dLayoutAvgStress const avgStress ):
     Base( nodeManager,
           edgeManager,
           faceManager,
           elementSubRegion,
           finiteElementSpace ),
     m_displacement( displacement ),
-    m_avgStrain( avgStrain )
+    m_avgStrain( avgStrain ),
+    m_stress( stress ),
+    m_avgStress( avgStress )
   {}
 
   /**
@@ -105,6 +113,7 @@
     for( int icomp = 0; icomp < 6; ++icomp )
     {
       m_avgStrain[k][icomp] = 0.0;
+      m_avgStress[k][icomp] = 0.0;
     }
   }
 
@@ -129,6 +138,7 @@
     for( int icomp = 0; icomp < 6; ++icomp )
     {
       m_avgStrain[k][icomp] += detJxW*strain[icomp]/m_elementVolume[k];
+      m_avgStress[k][icomp] += detJxW*m_stress[k][q][icomp]/m_elementVolume[k];
     }
   }
 
@@ -166,15 +176,21 @@
   /// The average strain
   fields::solidMechanics::arrayView2dLayoutStrain const m_avgStrain;
 
+  /// The stress solution
+  arrayView3d< real64 const, solid::STRESS_USD > const m_stress;
+
+  /// The average stress
+  fields::solidMechanics::arrayView2dLayoutAvgStress const m_avgStress;
+
 };
 
 
 
 /**
- * @class AverageStrainOverQuadraturePointsKernelFactory
+ * @class AverageStressStrainOverQuadraturePointsKernelFactory
  * @brief Class to create and launch the kernel
  */
-class AverageStrainOverQuadraturePointsKernelFactory
+class AverageStressStrainOverQuadraturePointsKernelFactory
 {
 public:
 
@@ -201,13 +217,15 @@
                    SUBREGION_TYPE const & elementSubRegion,
                    FE_TYPE const & finiteElementSpace,
                    fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement,
-                   fields::solidMechanics::arrayView2dLayoutStrain const avgStrain )
+                   fields::solidMechanics::arrayView2dLayoutStrain const avgStrain,
+                   arrayView3d< real64 const, solid::STRESS_USD > const stress,
+                   fields::solidMechanics::arrayView2dLayoutAvgStress const avgStress )
   {
-    AverageStrainOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE >
+    AverageStressStrainOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE >
     kernel( nodeManager, edgeManager, faceManager, elementSubRegion, finiteElementSpace,
-            displacement, avgStrain );
+            displacement, avgStrain, stress, avgStress );
 
-    AverageStrainOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE >::template
+    AverageStressStrainOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE >::template
     kernelLaunch< POLICY >( elementSubRegion.size(), kernel );
   }
 };
@@ -217,4 +235,5 @@
 }
 
 
-#endif /* GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_STRAINHELPER_HPP_ */
+
+#endif /* GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_STRESSSTRAINAVERAGEKERNELS_HPP_ */