feat: plastic strain output

src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsFields.hpp

@@ -87,6 +87,14 @@ DECLARE_FIELD( strain,
                WRITE_AND_READ,
                "Average strain in cell" );
 
+DECLARE_FIELD( plasticStrain,
+               "plasticStrain",
+               array2dLayoutStrain,
+               0,
+               LEVEL_0,
+               WRITE_AND_READ,
+               "Average plastic strain in cell" );
+
 DECLARE_FIELD( incrementalBubbleDisplacement,
                "incrementalBubbleDisplacement",
                array2d< real64 >,

src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp

@@ -160,7 +160,9 @@ void SolidMechanicsLagrangianFEM::registerDataOnMesh( Group & meshBodies )
     {
       setConstitutiveNamesCallSuper( subRegion );
 
+      subRegion.registerField< solidMechanics::strain >( getName() ).reference().resizeDimension< 1 >( 6 );
       subRegion.registerField< solidMechanics::strain >( getName() ).setDimLabels( 1, voightLabels ).reference().resizeDimension< 1 >( 6 );
+
     } );
 
     NodeManager & nodes = meshLevel.getNodeManager();
@@ -947,18 +949,29 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS
       constitutiveRelation.saveConvergedState();
 
       solidMechanics::arrayView2dLayoutStrain strain = subRegion.getField< solidMechanics::strain >();
+      solidMechanics::arrayView2dLayoutStrain plasticStrain = subRegion.getField< solidMechanics::plasticStrain >();
+
+      constitutive::ConstitutivePassThru< SolidBase >::execute(constitutiveRelation, [&] (auto & solidModel)
+      {
+
+        using SOLID_TYPE = TYPEOFREF( solidModel );
 
       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,
+        AverageStrainOverQuadraturePointsKernelFactory::createAndLaunch< FE_TYPE, SOLID_TYPE, parallelDevicePolicy<> >( nodeManager,
                                                                                                                                   mesh.getEdgeManager(),
                                                                                                                                   mesh.getFaceManager(),
                                                                                                                                   subRegion,
                                                                                                                                   finiteElement,
+                                                                                                                                  solidModel,
                                                                                                                                   disp,
-                                                                                                                                  strain );
+                                                                                                                                  strain,
+                                                                                                                                  plasticStrain );
+      } );
+
+
       } );
 
 

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

@@ -23,6 +23,7 @@
 #include "common/DataTypes.hpp"
 #include "common/GEOS_RAJA_Interface.hpp"
 #include "finiteElement/FiniteElementDispatch.hpp"
+#include "constitutive/ConstitutivePassThru.hpp"
 #include "mesh/CellElementSubRegion.hpp"
 #include "mesh/utilities/AverageOverQuadraturePointsKernel.hpp"
 #include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp"
@@ -33,17 +34,18 @@ namespace geos
  * @class AverageStrainOverQuadraturePoints
  * @tparam SUBREGION_TYPE the subRegion type
  * @tparam FE_TYPE the finite element type
+ * @tparam SOLID_TYPE the solid mechanics constitutuve type
  */
-template< typename SUBREGION_TYPE,
-          typename FE_TYPE >
+template< typename FE_TYPE,
+          typename SOLID_TYPE >
 class AverageStrainOverQuadraturePoints :
-  public AverageOverQuadraturePointsBase< SUBREGION_TYPE,
+  public AverageOverQuadraturePointsBase< CellElementSubRegion,
                                           FE_TYPE >
 {
 public:
 
   /// Alias for the base class;
-  using Base = AverageOverQuadraturePointsBase< SUBREGION_TYPE,
+  using Base = AverageOverQuadraturePointsBase< CellElementSubRegion,
                                                 FE_TYPE >;
 
   using Base::m_elementVolume;
@@ -63,17 +65,21 @@ class AverageStrainOverQuadraturePoints :
   AverageStrainOverQuadraturePoints( NodeManager & nodeManager,
                                      EdgeManager const & edgeManager,
                                      FaceManager const & faceManager,
-                                     SUBREGION_TYPE const & elementSubRegion,
+                                     CellElementSubRegion const & elementSubRegion,
                                      FE_TYPE const & finiteElementSpace,
+                                     SOLID_TYPE const & solidModel,
                                      fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement,
-                                     fields::solidMechanics::arrayView2dLayoutStrain const avgStrain ):
+                                     fields::solidMechanics::arrayView2dLayoutStrain const avgStrain,
+                                     fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain):
     Base( nodeManager,
           edgeManager,
           faceManager,
           elementSubRegion,
           finiteElementSpace ),
+    m_solidUpdate(solidModel.createKernelUpdates()),
     m_displacement( displacement ),
-    m_avgStrain( avgStrain )
+    m_avgStrain( avgStrain ),
+    m_avgPlasticStrain( avgPlasticStrain )
   {}
 
   /**
@@ -105,6 +111,7 @@ class AverageStrainOverQuadraturePoints :
     for( int icomp = 0; icomp < 6; ++icomp )
     {
       m_avgStrain[k][icomp] = 0.0;
+      m_avgPlasticStrain[k][icomp] = 0.0;
     }
   }
 
@@ -126,9 +133,13 @@ class AverageStrainOverQuadraturePoints :
     real64 strain[6] = {0.0};
     FE_TYPE::symmetricGradient( dNdX, stack.uLocal, strain );
 
+    real64 elasticStrain[6] = {0.0};
+    m_solidUpdate.getElasticStrain(k, q, elasticStrain);
+
     for( int icomp = 0; icomp < 6; ++icomp )
     {
       m_avgStrain[k][icomp] += detJxW*strain[icomp]/m_elementVolume[k];
+      m_avgPlasticStrain[k][icomp] += detJxW*(strain[icomp] - elasticStrain[icomp])/m_elementVolume[k];
     }
   }
 
@@ -160,12 +171,18 @@ class AverageStrainOverQuadraturePoints :
 
 protected:
 
+  /// The material
+  typename SOLID_TYPE::KernelWrapper const m_solidUpdate;
+
   /// The displacement solution
   fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const m_displacement;
 
   /// The average strain
   fields::solidMechanics::arrayView2dLayoutStrain const m_avgStrain;
 
+  /// The average plastic strain
+  fields::solidMechanics::arrayView2dLayoutStrain const m_avgPlasticStrain;
+
 };
 
 
@@ -182,6 +199,7 @@ class AverageStrainOverQuadraturePointsKernelFactory
    * @brief Create a new kernel and launch
    * @tparam SUBREGION_TYPE the subRegion type
    * @tparam FE_TYPE the finite element type
+   * @tparam SOLID_TYPE the constitutive type
    * @tparam POLICY the kernel policy
    * @param nodeManager the node manager
    * @param edgeManager the edge manager
@@ -191,23 +209,25 @@ class AverageStrainOverQuadraturePointsKernelFactory
    * @param property the property at quadrature points
    * @param averageProperty the property averaged over quadrature points
    */
-  template< typename SUBREGION_TYPE,
-            typename FE_TYPE,
+  template< typename FE_TYPE,
+            typename SOLID_TYPE,
             typename POLICY >
   static void
   createAndLaunch( NodeManager & nodeManager,
                    EdgeManager const & edgeManager,
                    FaceManager const & faceManager,
-                   SUBREGION_TYPE const & elementSubRegion,
+                   CellElementSubRegion const & elementSubRegion,
                    FE_TYPE const & finiteElementSpace,
+                   SOLID_TYPE const & solidModel,
                    fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement,
-                   fields::solidMechanics::arrayView2dLayoutStrain const avgStrain )
+                   fields::solidMechanics::arrayView2dLayoutStrain const avgStrain,
+                   fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain)
   {
-    AverageStrainOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE >
+    AverageStrainOverQuadraturePoints< FE_TYPE, SOLID_TYPE >
     kernel( nodeManager, edgeManager, faceManager, elementSubRegion, finiteElementSpace,
-            displacement, avgStrain );
+            solidModel, displacement, avgStrain, avgPlasticStrain );
 
-    AverageStrainOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE >::template
+    AverageStrainOverQuadraturePoints< FE_TYPE, SOLID_TYPE >::template
     kernelLaunch< POLICY >( elementSubRegion.size(), kernel );
   }
 };