SolidMechanicsDisplacementJumpUpdateKernels.hpp

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2024 TotalEnergies
 * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2023-2024 Chevron
 * Copyright (c) 2019-     GEOS/GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_PHYSICSSOLVERS_CONTACT_KERNELS_SOLIDMECHANICSDISPLACEMENTJUPDATEKERNELS_HPP_
#define GEOS_PHYSICSSOLVERS_CONTACT_KERNELS_SOLIDMECHANICSDISPLACEMENTJUPDATEKERNELS_HPP_
 
#include "SolidMechanicsConformingContactKernelsBase.hpp"
#include "mesh/MeshFields.hpp"
 
namespace geos
{
 
namespace solidMechanicsConformingContactKernels
{
 
template< typename CONSTITUTIVE_TYPE,
          typename FE_TYPE >
class DispJumpUpdate :
  public ConformingContactKernelsBase< CONSTITUTIVE_TYPE,
                                       FE_TYPE >
{
public:
  using Base = ConformingContactKernelsBase< CONSTITUTIVE_TYPE,
                                             FE_TYPE >;
 
  static constexpr int numNodesPerElem = Base::maxNumTestSupportPointsPerElem;
 
  static constexpr int numUdofs = Base::numUdofs;
 
  static constexpr int numBdofs = Base::numBdofs;
 
  static constexpr int numTdofs = Base::numTdofs;
 
  using Base::m_X;
  using Base::m_finiteElementSpace;
  using Base::m_dofNumber;
  using Base::m_dofRankOffset;
  using Base::m_elemsToFaces;
  using Base::m_faceToNodes;
  using Base::m_rotationMatrix;
  using Base::m_dispJump;
 
  DispJumpUpdate( NodeManager const & nodeManager,
                  EdgeManager const & edgeManager,
                  FaceManager const & faceManager,
                  localIndex const targetRegionIndex,
                  FaceElementSubRegion & elementSubRegion,
                  FE_TYPE const & finiteElementSpace,
                  CONSTITUTIVE_TYPE & inputConstitutiveType,
                  arrayView1d< globalIndex const > const uDofNumber,
                  arrayView1d< globalIndex const > const bDofNumber,
                  globalIndex const rankOffset,
                  CRSMatrixView< real64, globalIndex const > const inputMatrix,
                  arrayView1d< real64 > const inputRhs,
                  real64 const inputDt,
                  arrayView1d< localIndex const > const & faceElementList ):
    Base( nodeManager,
          edgeManager,
          faceManager,
          targetRegionIndex,
          elementSubRegion,
          finiteElementSpace,
          inputConstitutiveType,
          uDofNumber,
          bDofNumber,
          rankOffset,
          inputMatrix,
          inputRhs,
          inputDt,
          faceElementList ),
    m_displacement( nodeManager.getField< fields::solidMechanics::totalDisplacement >()),
    m_bubbleDisp( faceManager.getField< fields::solidMechanics::totalBubbleDisplacement >() ),
    m_incrDisp( nodeManager.getField< fields::solidMechanics::incrementalDisplacement >() ),
    m_incrBubbleDisp( faceManager.getField< fields::solidMechanics::incrementalBubbleDisplacement >() ),
    m_deltaDispJump( elementSubRegion.getField< fields::contact::deltaDispJump >().toView() ),
    m_elementArea( elementSubRegion.getField< fields::elementArea >().toView() ),
    m_slip( elementSubRegion.getField< fields::contact::slip >().toView() )
  {}
 
  //***************************************************************************
 
  struct StackVariables : public Base::StackVariables
  {
 
public:
 
    GEOS_HOST_DEVICE
    StackVariables():
      Base::StackVariables(),
                                       uLocal{},
                                       bLocal{},
                                       duLocal{},
                                       dbLocal{},
                                       deltaDispJumpLocal{}
    {}
 
    real64 uLocal[numUdofs];
 
    real64 bLocal[numBdofs];
 
    real64 duLocal[numUdofs];
 
    real64 dbLocal[numBdofs];
 
    real64 deltaDispJumpLocal[numTdofs];
 
  };
 
  //***************************************************************************
 
  //START_kernelLauncher
  template< typename POLICY,
            typename KERNEL_TYPE >
  static
  real64
  kernelLaunch( localIndex const numElems,
                KERNEL_TYPE const & kernelComponent )
  {
    return Base::template kernelLaunch< POLICY, KERNEL_TYPE >( numElems, kernelComponent );
  }
  //END_kernelLauncher
 
  GEOS_HOST_DEVICE
  inline
  void setup( localIndex const k,
              StackVariables & stack ) const
  {
    constexpr int shift = numNodesPerElem * 3;
 
    int permutation[numNodesPerElem];
    m_finiteElementSpace.getPermutation( permutation );
 
    localIndex const kf0 = m_elemsToFaces[k][0];
    localIndex const kf1 = m_elemsToFaces[k][1];
    for( localIndex a=0; a<numNodesPerElem; ++a )
    {
      localIndex const kn0 = m_faceToNodes( kf0, a );
      localIndex const kn1 = m_faceToNodes( kf1, a );
 
      for( int i=0; i<3; ++i )
      {
        stack.X[ a ][ i ] = m_X[ m_faceToNodes( kf0, permutation[ a ] ) ][ i ];
        stack.uLocal[a*3+i] = m_displacement[kn0][i];
        stack.uLocal[shift + a*3+i] = m_displacement[kn1][i];
        stack.duLocal[a*3+i] = m_incrDisp[kn0][i];
        stack.duLocal[shift + a*3+i] = m_incrDisp[kn1][i];
      }
    }
 
    for( int j=0; j<3; ++j )
    {
      for( int i=0; i<3; ++i )
      {
        stack.localRotationMatrix[ i ][ j ] = m_rotationMatrix( k, i, j );
      }
    }
 
    for( int i=0; i<3; ++i )
    {
      stack.bLocal[ i ] = m_bubbleDisp[ kf0 ][i];
      stack.bLocal[ 3 + i ] = m_bubbleDisp[ kf1 ][i];
      stack.dbLocal[ i ] = m_incrBubbleDisp[ kf0 ][i];
      stack.dbLocal[ 3 + i ] = m_incrBubbleDisp[ kf1 ][i];
    }
 
  }
 
  GEOS_HOST_DEVICE
  inline
  real64 complete( localIndex const k,
                   StackVariables & stack ) const
  {
    real64 matRtAtu[3][numUdofs];
    real64 matRtAtb[3][numBdofs];
 
    // transp(R) * Atu
    LvArray::tensorOps::Rij_eq_AkiBkj< 3, numUdofs, 3 >( matRtAtu, stack.localRotationMatrix,
                                                         stack.localAtu );
    // transp(R) * Atb
    LvArray::tensorOps::Rij_eq_AkiBkj< 3, numBdofs, 3 >( matRtAtb, stack.localRotationMatrix,
                                                         stack.localAtb );
 
    // Compute the node contribute of the displacement and delta displacement jump
    LvArray::tensorOps::Ri_eq_AijBj< 3, numUdofs >( stack.dispJumpLocal, matRtAtu, stack.uLocal );
    LvArray::tensorOps::Ri_eq_AijBj< 3, numUdofs >( stack.deltaDispJumpLocal, matRtAtu, stack.duLocal );
 
    // Compute the bubble contribute of the displacement and delta displacement jump
    LvArray::tensorOps::Ri_add_AijBj< 3, numBdofs >( stack.dispJumpLocal, matRtAtb, stack.bLocal );
    LvArray::tensorOps::Ri_add_AijBj< 3, numBdofs >( stack.deltaDispJumpLocal, matRtAtb, stack.dbLocal );
 
    // Store the results
    real64 const scale = 1 / m_elementArea[k];
 
    for( int i=0; i<3; ++i )
    {
      m_dispJump[ k ][ i ] = scale * stack.dispJumpLocal[ i ];
      m_deltaDispJump[ k ][ i ] = scale * stack.deltaDispJumpLocal[ i ];
    }
 
    m_slip[k] = LvArray::math::sqrt( LvArray::math::square( m_dispJump( k, 1 ) ) + LvArray::math::square( m_dispJump( k, 2 ) ) );
 
 
    return 0.0;
  }
 
protected:
 
  arrayView2d< real64 const, nodes::TOTAL_DISPLACEMENT_USD > const m_displacement;
 
  arrayView2d< real64 const > const m_bubbleDisp;
 
  arrayView2d< real64 const, nodes::INCR_DISPLACEMENT_USD > const m_incrDisp;
 
  arrayView2d< real64 const > const m_incrBubbleDisp;
 
  arrayView2d< real64 > const m_deltaDispJump;
 
  arrayView1d< real64 const > const m_elementArea;
 
  arrayView1d< real64 > const m_slip;
 
};
 
using DispJumpUpdateFactory = finiteElement::InterfaceKernelFactory< DispJumpUpdate,
                                                                     arrayView1d< globalIndex const > const,
                                                                     arrayView1d< globalIndex const > const,
                                                                     globalIndex const,
                                                                     CRSMatrixView< real64, globalIndex const > const,
                                                                     arrayView1d< real64 > const,
                                                                     real64 const,
                                                                     arrayView1d< localIndex const > const >;
 
} // namespace SolidMechanicsALMKernels
 
} // namespace geos
 
#endif /* GEOS_PHYSICSSOLVERS_CONTACT_KERNELS_SOLIDMECHANICSDISPLACEMENTJUPDATEKERNELS_HPP_ */