feat: averaging of stress components using quadrature

src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsFields.hpp

@@ -46,6 +46,10 @@
 using arrayView2dLayoutStrain = arrayView2d< real64, cells::STRAIN_USD >;
 using arrayViewConst2dLayoutStrain = arrayView2d< real64 const, cells::STRAIN_USD >;
 
+using array2dLayoutAvgStress = array2d< real64, cells::STRAIN_PERM >;
+using arrayView2dLayoutAvgStress = arrayView2d< real64, cells::STRAIN_USD >;
+using arrayViewConst2dLayoutAvgStress = arrayView2d< real64 const, cells::STRAIN_USD >;
+
 using array2dLayoutVelocity = array2d< real64, nodes::VELOCITY_PERM >;
 using arrayView2dLayoutVelocity = arrayView2d< real64, nodes::VELOCITY_USD >;
 using arrayViewConst2dLayoutVelocity = arrayView2d< real64 const, nodes::VELOCITY_USD >;
@@ -87,6 +91,14 @@
                WRITE_AND_READ,
                "Average strain in cell" );
 
+DECLARE_FIELD( averageStress,
+               "averageStress",
+               array2dLayoutAvgStress,
+               0,
+               LEVEL_0,
+               WRITE_AND_READ,
+               "Quadrature averaged stress components in cell" );
+
 DECLARE_FIELD( incrementalBubbleDisplacement,
                "incrementalBubbleDisplacement",
                array2d< real64 >,

src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp

@@ -159,6 +159,7 @@
       setConstitutiveNamesCallSuper( subRegion );
 
       subRegion.registerField< solidMechanics::strain >( getName() ).reference().resizeDimension< 1 >( 6 );
+      subRegion.registerField< solidMechanics::averageStress >( getName() ).reference().resizeDimension< 1 >( 6 );
     } );
 
     NodeManager & nodes = meshLevel.getNodeManager();
@@ -944,7 +945,10 @@
       SolidBase & constitutiveRelation = getConstitutiveModel< SolidBase >( subRegion, solidMaterialName );
       constitutiveRelation.saveConvergedState();
 
+      arrayView3d < real64 const, solid::STRESS_USD > const stress = constitutiveRelation.getStress();
+
       solidMechanics::arrayView2dLayoutStrain strain = subRegion.getField< solidMechanics::strain >();
+      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 )
@@ -957,6 +961,16 @@
                                                                                                                                   finiteElement,
                                                                                                                                   disp,
                                                                                                                                   strain );
+
+
+        AverageStressOverQuadraturePointsKernelFactory::createAndLaunch< CellElementSubRegion, FE_TYPE, parallelDevicePolicy<> >( nodeManager,
+                                                                                                                                  mesh.getEdgeManager(),
+                                                                                                                                  mesh.getFaceManager(),
+                                                                                                                                  subRegion,
+                                                                                                                                  finiteElement,
+                                                                                                                                  stress,
+                                                                                                                                  avgStress );
+
       } );
 
 

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

@@ -214,6 +214,182 @@
 
 
 
+
+
+
+/**
+ * @class AverageStressOverQuadraturePoints
+ * @tparam SUBREGION_TYPE the subRegion type
+ * @tparam FE_TYPE the finite element type
+ */
+template< typename SUBREGION_TYPE,
+          typename FE_TYPE >
+class AverageStressOverQuadraturePoints :
+  public AverageOverQuadraturePointsBase< SUBREGION_TYPE,
+                                          FE_TYPE >
+{
+public:
+
+  /// Alias for the base class;
+  using Base = AverageOverQuadraturePointsBase< SUBREGION_TYPE,
+                                                FE_TYPE >;
+
+  using Base::m_elementVolume;
+  using Base::m_elemsToNodes;
+  using Base::m_finiteElementSpace;
+
+  /**
+   * @brief Constructor for the class
+   * @param nodeManager the node manager
+   * @param edgeManager the edge manager
+   * @param faceManager the face manager
+   * @param elementSubRegion the element subRegion
+   * @param finiteElementSpace the finite element space
+   * @param stress the stress solution field
+   * @param avgStress the stress averaged over quadrature points
+   */
+  AverageStressOverQuadraturePoints( NodeManager & nodeManager,
+                                     EdgeManager const & edgeManager,
+                                     FaceManager const & faceManager,
+                                     SUBREGION_TYPE const & elementSubRegion,
+                                     FE_TYPE const & finiteElementSpace,
+                                     arrayView3d< real64 const, solid::STRESS_USD > const stress,
+                                     fields::solidMechanics::arrayView2dLayoutAvgStress const avgStress ):
+    Base( nodeManager,
+          edgeManager,
+          faceManager,
+          elementSubRegion,
+          finiteElementSpace ),
+    m_stress( stress ),
+    m_avgStress( avgStress )
+  {}
+
+  /**
+   * @copydoc finiteElement::KernelBase::StackVariables
+   */
+  struct StackVariables : Base::StackVariables
+  {  };
+
+  /**
+   * @brief Performs the setup phase for the kernel.
+   * @param k The element index.
+   * @param stack The StackVariable object that hold the stack variables.
+   */
+  GEOS_HOST_DEVICE
+  void setup( localIndex const k,
+              StackVariables & stack ) const
+  {
+    Base::setup( k, stack );
+
+    for( int icomp = 0; icomp < 6; ++icomp )
+    {
+      m_avgStress[k][icomp] = 0.0;
+    }
+  }
+
+  /**
+   * @brief Increment the average property with the contribution of the property at this quadrature point
+   * @param k The element index
+   * @param q The quadrature point index
+   * @param stack The StackVariables object that hold the stack variables.
+   */
+  GEOS_HOST_DEVICE
+  void quadraturePointKernel( localIndex const k,
+                              localIndex const q,
+                              StackVariables & stack ) const
+  {
+    //real64 const weight = FE_TYPE::transformedQuadratureWeight( q, stack.xLocal, stack.feStack ) / m_elementVolume[k];
+
+    real64 dNdX[ FE_TYPE::maxSupportPoints ][3];
+    real64 const detJxW = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal, stack.feStack, dNdX );
+
+    for( int icomp = 0; icomp < 6; ++icomp )
+    {
+      m_avgStress[k][icomp] += detJxW*m_stress[k][q][icomp]/m_elementVolume[k];
+    }
+  }
+
+  /**
+   * @brief Launch the kernel over the elements in the subRegion
+   * @tparam POLICY the kernel policy
+   * @tparam KERNEL_TYPE the type of kernel
+   * @param numElems the number of elements in the subRegion
+   * @param kernelComponent the kernel component
+   */
+  template< typename POLICY,
+            typename KERNEL_TYPE >
+  static void
+  kernelLaunch( localIndex const numElems,
+                KERNEL_TYPE const & kernelComponent )
+  {
+    forAll< POLICY >( numElems,
+                      [=] GEOS_HOST_DEVICE ( localIndex const k )
+    {
+      typename KERNEL_TYPE::StackVariables stack;
+
+      kernelComponent.setup( k, stack );
+      for( integer q = 0; q < FE_TYPE::numQuadraturePoints; ++q )
+      {
+        kernelComponent.quadraturePointKernel( k, q, stack );
+      }
+    } );
+  }
+
+protected:
+
+  /// The stress solution
+  arrayView3d< real64 const, solid::STRESS_USD > const m_stress;
+
+  /// The average stress
+  fields::solidMechanics::arrayView2dLayoutAvgStress const m_avgStress;
+
+};
+
+
+
+/**
+ * @class AverageStressOverQuadraturePointsKernelFactory
+ * @brief Class to create and launch the kernel
+ */
+class AverageStressOverQuadraturePointsKernelFactory
+{
+public:
+
+  /**
+   * @brief Create a new kernel and launch
+   * @tparam SUBREGION_TYPE the subRegion type
+   * @tparam FE_TYPE the finite element type
+   * @tparam POLICY the kernel policy
+   * @param nodeManager the node manager
+   * @param edgeManager the edge manager
+   * @param faceManager the face manager
+   * @param elementSubRegion the element subRegion
+   * @param finiteElementSpace the finite element space
+   * @param property the property at quadrature points
+   * @param averageProperty the property averaged over quadrature points
+   */
+  template< typename SUBREGION_TYPE,
+            typename FE_TYPE,
+            typename POLICY >
+  static void
+  createAndLaunch( NodeManager & nodeManager,
+                   EdgeManager const & edgeManager,
+                   FaceManager const & faceManager,
+                   SUBREGION_TYPE const & elementSubRegion,
+                   FE_TYPE const & finiteElementSpace,
+                   arrayView3d< real64 const, solid::STRESS_USD > const  stress,
+                   fields::solidMechanics::arrayView2dLayoutAvgStress const avgStress )
+  {
+    AverageStressOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE >
+    kernel( nodeManager, edgeManager, faceManager, elementSubRegion, finiteElementSpace,
+            stress, avgStress );
+
+    AverageStressOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE >::template
+    kernelLaunch< POLICY >( elementSubRegion.size(), kernel );
+  }
+};
+
+
 }