SolidMechanicsLagrangianFEM.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * 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_SOLIDMECHANICS_SOLIDMECHANICSLAGRANGIANFEM_HPP_
#define GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_SOLIDMECHANICSLAGRANGIANFEM_HPP_
 
#include "common/format/EnumStrings.hpp"
#include "common/TimingMacros.hpp"
#include "kernels/SolidMechanicsLagrangianFEMKernels.hpp"
#include "kernels/StrainHelper.hpp"
#include "mesh/mpiCommunications/CommunicationTools.hpp"
#include "mesh/mpiCommunications/MPI_iCommData.hpp"
#include "physicsSolvers/PhysicsSolverBase.hpp"
 
#include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp"
 
namespace geos
{
 
class SolidMechanicsLagrangianFEM : public PhysicsSolverBase
{
public:
 
  static string coupledSolverAttributePrefix() { return "solid"; }
 
  enum class TimeIntegrationOption : integer
  {
    QuasiStatic,      
    ImplicitDynamic,  
    ExplicitDynamic   
  };
 
  SolidMechanicsLagrangianFEM( const string & name,
                               Group * const parent );
 
  static string catalogName() { return "SolidMechanicsLagrangianFEM"; }
  string getCatalogName() const override { return catalogName(); }
 
  virtual void initializePreSubGroups() override;
 
  virtual void registerDataOnMesh( Group & meshBodies ) override;
 
  virtual
  real64 solverStep( real64 const & time_n,
                     real64 const & dt,
                     integer const cycleNumber,
                     DomainPartition & domain ) override;
 
  virtual
  real64 explicitStep( real64 const & time_n,
                       real64 const & dt,
                       integer const cycleNumber,
                       DomainPartition & domain ) override;
 
  virtual void
  implicitStepSetup( real64 const & time_n,
                     real64 const & dt,
                     DomainPartition & domain ) override;
 
  virtual void
  setupDofs( DomainPartition const & domain,
             DofManager & dofManager ) const override;
 
  virtual void
  setupSystem( DomainPartition & domain,
               DofManager & dofManager,
               CRSMatrix< real64, globalIndex > & localMatrix,
               ParallelVector & rhs,
               ParallelVector & solution,
               bool const setSparsity = false ) override;
 
  virtual std::unique_ptr< PreconditionerBase< LAInterface > >
  createPreconditioner( DomainPartition & domain ) const override;
 
  virtual void
  assembleSystem( real64 const time,
                  real64 const dt,
                  DomainPartition & domain,
                  DofManager const & dofManager,
                  CRSMatrixView< real64, globalIndex const > const & localMatrix,
                  arrayView1d< real64 > const & localRhs ) override;
 
  virtual void solveLinearSystem( DofManager const & dofManager,
                                  ParallelMatrix & matrix,
                                  ParallelVector & rhs,
                                  ParallelVector & solution ) override;
 
  virtual void
  applySystemSolution( DofManager const & dofManager,
                       arrayView1d< real64 const > const & localSolution,
                       real64 const scalingFactor,
                       real64 const dt,
                       DomainPartition & domain ) override;
 
  virtual void updateState( DomainPartition & domain ) override
  {
    // There should be nothing to update
    GEOS_UNUSED_VAR( domain );
  };
 
  virtual void applyBoundaryConditions( real64 const time,
                                        real64 const dt,
                                        DomainPartition & domain,
                                        DofManager const & dofManager,
                                        CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                        arrayView1d< real64 > const & localRhs ) override;
 
  virtual real64
  calculateResidualNorm( real64 const & time_n,
                         real64 const & dt,
                         DomainPartition const & domain,
                         DofManager const & dofManager,
                         arrayView1d< real64 const > const & localRhs ) override;
 
  virtual void resetStateToBeginningOfStep( DomainPartition & domain ) override;
 
  virtual void implicitStepComplete( real64 const & time,
                                     real64 const & dt,
                                     DomainPartition & domain ) override;
 
  template< typename TYPE_LIST,
            typename KERNEL_WRAPPER,
            typename ... PARAMS >
  real64 assemblyLaunch( MeshLevel & mesh,
                         DofManager const & dofManager,
                         string_array const & regionNames,
                         string const & materialNamesString,
                         CRSMatrixView< real64, globalIndex const > const & localMatrix,
                         arrayView1d< real64 > const & localRhs,
                         real64 const dt,
                         PARAMS && ... params );
 
  real64 explicitKernelDispatch( MeshLevel & mesh,
                                 string_array const & targetRegions,
                                 string const & finiteElementName,
                                 real64 const dt,
                                 std::string const & elementListName );
 
  void applyDisplacementBCImplicit( real64 const time,
                                    DofManager const & dofManager,
                                    DomainPartition & domain,
                                    CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                    arrayView1d< real64 > const & localRhs );
 
  void applyTractionBC( real64 const time,
                        DofManager const & dofManager,
                        DomainPartition & domain,
                        arrayView1d< real64 > const & localRhs );
 
  void applyChomboPressure( DofManager const & dofManager,
                            DomainPartition & domain,
                            arrayView1d< real64 > const & localRhs );
 
 
  void applyContactConstraint( DofManager const & dofManager,
                               DomainPartition & domain,
                               CRSMatrixView< real64, globalIndex const > const & localMatrix,
                               arrayView1d< real64 > const & localRhs );
 
  virtual real64
  scalingForSystemSolution( DomainPartition & domain,
                            DofManager const & dofManager,
                            arrayView1d< real64 const > const & localSolution ) override;
 
  void enableFixedStressPoromechanicsUpdate();
 
  virtual void saveSequentialIterationState( DomainPartition & domain ) override;
 
  struct viewKeyStruct : PhysicsSolverBase::viewKeyStruct
  {
    static constexpr char const * newmarkGammaString() { return "newmarkGamma"; }
    static constexpr char const * newmarkBetaString() { return "newmarkBeta"; }
    static constexpr char const * massDampingString() { return "massDamping"; }
    static constexpr char const * stiffnessDampingString() { return "stiffnessDamping"; }
    static constexpr char const * timeIntegrationOptionString() { return "timeIntegrationOption"; }
    static constexpr char const * maxNumResolvesString() { return "maxNumResolves"; }
    static constexpr char const * strainTheoryString() { return "strainTheory"; }
    static constexpr char const * solidMaterialNamesString() { return "solidMaterialNames"; }
    static constexpr char const * contactRelationNameString() { return "contactRelationName"; }
    static constexpr char const * noContactRelationNameString() { return "NOCONTACT"; }
    static constexpr char const * maxForceString() { return "maxForce"; }
    static constexpr char const * elemsAttachedToSendOrReceiveNodesString() { return "elemsAttachedToSendOrReceiveNodes"; }
    static constexpr char const * elemsNotAttachedToSendOrReceiveNodesString() { return "elemsNotAttachedToSendOrReceiveNodes"; }
    static constexpr char const * surfaceGeneratorNameString() { return "surfaceGeneratorName"; }
 
    static constexpr char const * sendOrReceiveNodesString() { return "sendOrReceiveNodes";}
    static constexpr char const * nonSendOrReceiveNodesString() { return "nonSendOrReceiveNodes";}
    static constexpr char const * targetNodesString() { return "targetNodes";}
    static constexpr char const * forceString() { return "Force";}
 
    static constexpr char const * contactPenaltyStiffnessString() { return "contactPenaltyStiffness"; }
 
  };
 
  SortedArray< localIndex > & getElemsAttachedToSendOrReceiveNodes( ElementSubRegionBase & subRegion )
  {
    return subRegion.getReference< SortedArray< localIndex > >( viewKeyStruct::elemsAttachedToSendOrReceiveNodesString() );
  }
 
  SortedArray< localIndex > & getElemsNotAttachedToSendOrReceiveNodes( ElementSubRegionBase & subRegion )
  {
    return subRegion.getReference< SortedArray< localIndex > >( viewKeyStruct::elemsNotAttachedToSendOrReceiveNodesString() );
  }
 
  real64 & getMaxForce() { return m_maxForce; }
  real64 const & getMaxForce() const { return m_maxForce; }
 
  void computeRigidBodyModes( DomainPartition & domain ) const;
 
  arrayView1d< ParallelVector > const & getRigidBodyModes( DomainPartition & domain ) const
  {
    computeRigidBodyModes( domain );
    return m_rigidBodyModes;
  }
 
  /*
   * @brief Utility function to set the stress initialization flag
   * @param[in] performStressInitialization true if the solver has to initialize stress, false otherwise
   */
  void setStressInitialization( bool const performStressInitialization )
  {
    m_performStressInitialization = performStressInitialization;
  }
 
  TimeIntegrationOption timeIntegrationOption() const { return m_timeIntegrationOption; }
 
protected:
  virtual void postInputInitialization() override;
 
  virtual void initializePostInitialConditionsPreSubGroups() override;
 
  virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override;
 
  real64 m_newmarkGamma;
  real64 m_newmarkBeta;
  real64 m_massDamping;
  real64 m_stiffnessDamping;
  TimeIntegrationOption m_timeIntegrationOption;
  real64 m_maxForce = 0.0;
  integer m_maxNumResolves;
  integer m_strainTheory;
 
  bool m_isFixedStressPoromechanicsUpdate;
  bool m_performStressInitialization;
 
  mutable array1d< ParallelVector > m_rigidBodyModes;
 
  real64 m_contactPenaltyStiffness;
 
private:
 
  string m_contactRelationName;
 
  PhysicsSolverBase *m_surfaceGenerator;
  string m_surfaceGeneratorName;
};
 
ENUM_STRINGS( SolidMechanicsLagrangianFEM::TimeIntegrationOption,
              "QuasiStatic",
              "ImplicitDynamic",
              "ExplicitDynamic" );
 
//**********************************************************************************************************************
//**********************************************************************************************************************
//**********************************************************************************************************************
 
 
template< typename TYPE_LIST,
          typename KERNEL_WRAPPER,
          typename ... PARAMS >
real64 SolidMechanicsLagrangianFEM::assemblyLaunch( MeshLevel & mesh,
                                                    DofManager const & dofManager,
                                                    string_array const & regionNames,
                                                    string const & materialNamesString,
                                                    CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                                    arrayView1d< real64 > const & localRhs,
                                                    real64 const dt,
                                                    PARAMS && ... params )
{
  GEOS_MARK_FUNCTION;
 
  NodeManager const & nodeManager = mesh.getNodeManager();
 
  string const dofKey = dofManager.getKey( fields::solidMechanics::totalDisplacement::key() );
  arrayView1d< globalIndex const > const & dofNumber = nodeManager.getReference< globalIndex_array >( dofKey );
 
  real64 const gravityVectorData[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( gravityVector() );
 
  KERNEL_WRAPPER kernelWrapper( dofNumber,
                                dofManager.rankOffset(),
                                localMatrix,
                                localRhs,
                                dt,
                                gravityVectorData,
                                std::forward< PARAMS >( params )... );
 
  return finiteElement::
           regionBasedKernelApplication< parallelDevicePolicy< >,
                                         TYPE_LIST >( mesh,
                                                      regionNames,
                                                      this->getDiscretizationName(),
                                                      materialNamesString,
                                                      kernelWrapper );
 
}
 
} /* namespace geos */
 
#endif /* GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_SOLIDMECHANICSLAGRANGIANFEM_HPP_ */