PhaseFieldFractureSolver.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_MULTIPHYSICS_PHASEFIELDFRACTURESOLVER_HPP_
#define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_PHASEFIELDFRACTURESOLVER_HPP_
 
#include "physicsSolvers/multiphysics/CoupledSolver.hpp"
#include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp"
#include "physicsSolvers/simplePDE/PhaseFieldDamageFEM.hpp"
 
namespace geos
{
 
class PhaseFieldFractureSolver : public CoupledSolver< SolidMechanicsLagrangianFEM, PhaseFieldDamageFEM >
{
public:
 
  using Base = CoupledSolver< SolidMechanicsLagrangianFEM, PhaseFieldDamageFEM >;
  using Base::m_solvers;
  using Base::m_dofManager;
  using Base::m_localMatrix;
  using Base::m_rhs;
  using Base::m_solution;
 
  PhaseFieldFractureSolver( const string & name,
                            Group * const parent );
 
  ~PhaseFieldFractureSolver() override = default;
 
  static string catalogName()
  {
    return "PhaseFieldFracture";
  }
  string getCatalogName() const override { return catalogName(); }
 
  static string coupledSolverAttributePrefix() { return "PhaseFieldFracture"; }
 
  enum class SolverType : integer
  {
    SolidMechanics = 0,
    Damage = 1
  };
 
  virtual void postInputInitialization() override final;
 
  SolidMechanicsLagrangianFEM * solidMechanicsSolver() const
  {
    return std::get< toUnderlying( SolverType::SolidMechanics ) >( m_solvers );
  }
 
  PhaseFieldDamageFEM * damageSolver() const
  {
    return std::get< toUnderlying( SolverType::Damage ) >( m_solvers );
  }
 
  virtual void mapSolutionBetweenSolvers( DomainPartition & domain, integer const solverType ) override final;
 
protected:
 
  virtual void initializePostInitialConditionsPreSubGroups() override final {}
 
private:
 
  void mapDamageToQuadrature( DomainPartition & domain );
 
};
 
 
template< typename FE_TYPE >
struct DamageInterpolationKernel
{
  DamageInterpolationKernel( CellElementSubRegion const & subRegion ):
    m_numElems( subRegion.size() )
  {}
 
  void interpolateDamageAndGradient( arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemToNodes,
                                     arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const xNodes,
                                     arrayView1d< real64 const > const nodalDamage,
                                     arrayView2d< real64 > damageFieldOnMaterial,
                                     arrayView3d< real64 > damageGradOnMaterial )
  {
    forAll< parallelDevicePolicy<> >( m_numElems, [=] GEOS_HOST_DEVICE ( localIndex const k )
    {
      constexpr localIndex numNodesPerElement = FE_TYPE::numNodes;
      constexpr localIndex n_q_points = FE_TYPE::numQuadraturePoints;
 
      real64 xLocal[ numNodesPerElement ][ 3 ];
      real64 nodalDamageLocal[ numNodesPerElement ];
 
      for( localIndex a = 0; a < numNodesPerElement; ++a )
      {
        localIndex const localNodeIndex = elemToNodes( k, a );
 
        for( int dim=0; dim < 3; ++dim )
        {
          xLocal[a][dim] = xNodes[ localNodeIndex ][dim];
        }
 
        nodalDamageLocal[ a ] = nodalDamage[ localNodeIndex ];
      }
 
      for( localIndex q = 0; q < n_q_points; ++q )
      {
        real64 N[ numNodesPerElement ];
        FE_TYPE::calcN( q, N );
 
        real64 dNdX[ numNodesPerElement ][ 3 ];
        real64 const detJ = FE_TYPE::calcGradN( q, xLocal, dNdX );
 
        GEOS_UNUSED_VAR( detJ );
 
        real64 qDamage = 0.0;
        real64 qDamageGrad[3] = {0, 0, 0};
        finiteElement::feOps::valueAndGradient( N, dNdX, nodalDamageLocal, qDamage, qDamageGrad );
 
        damageFieldOnMaterial( k, q ) = qDamage;
 
        for( int dim=0; dim < 3; ++dim )
        {
          damageGradOnMaterial[k][q][dim] = qDamageGrad[dim];
        }
      }
 
    } );
  }
 
  localIndex m_numElems;
};
 
} /* namespace geos */
 
#endif /* GEOS_PHYSICSSOLVERS_MULTIPHYSICS_PHASEFIELDFRACTURESOLVER_HPP_ */