OGS
ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > Class Template Reference

Detailed Description

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
class ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >

Definition at line 48 of file TH2MFEM.h.

#include <TH2MFEM.h>

Inheritance diagram for ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >:
[legend]
Collaboration diagram for ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >:
[legend]

Public Types

using ShapeMatricesTypeDisplacement = ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim >
 
using ShapeMatricesTypePressure = ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim >
 
template<int N>
using VectorType = typename ShapeMatricesTypePressure::template VectorType< N >
 
template<int M, int N>
using MatrixType = typename ShapeMatricesTypePressure::template MatrixType< M, N >
 
using GlobalDimMatrixType = typename ShapeMatricesTypePressure::GlobalDimMatrixType
 
using GlobalDimVectorType = typename ShapeMatricesTypePressure::GlobalDimVectorType
 
using SymmetricTensor = Eigen::Matrix< double, KelvinVectorSize, 1 >
 
using Invariants = MathLib::KelvinVector::Invariants< KelvinVectorSize >
 

Public Member Functions

 TH2MLocalAssembler (TH2MLocalAssembler const &)=delete
 
 TH2MLocalAssembler (TH2MLocalAssembler &&)=delete
 
 TH2MLocalAssembler (MeshLib::Element const &e, std::size_t const, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, TH2MProcessData< DisplacementDim > &process_data)
 
- Public Member Functions inherited from ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >
 LocalAssemblerInterface (MeshLib::Element const &e, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, TH2MProcessData< DisplacementDim > &process_data)
 
virtual std::size_t setIPDataInitialConditions (std::string_view name, double const *values, int const integration_order)=0
 
virtual std::vector< double > const & getIntPtDarcyVelocityGas (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtDarcyVelocityLiquid (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtDiffusionVelocityVapourGas (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtDiffusionVelocityGasGas (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtDiffusionVelocitySoluteLiquid (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtDiffusionVelocityLiquidLiquid (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtLiquidDensity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtGasDensity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtSolidDensity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtVapourPressure (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtPorosity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtMoleFractionGas (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtMassFractionGas (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtMassFractionLiquid (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtRelativePermeabilityGas (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtRelativePermeabilityLiquid (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtIntrinsicPermeability (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtEnthalpyGas (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtEnthalpyLiquid (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtEnthalpySolid (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
unsigned getNumberOfIntegrationPoints () const
 
int getMaterialID () const
 
std::vector< double > getMaterialStateVariableInternalState (std::function< std::span< double >(typename ConstitutiveRelations::SolidConstitutiveRelation< DisplacementDim >::MaterialStateVariables &)> const &get_values_span, int const &n_components) const
 
ConstitutiveRelations::SolidConstitutiveRelation< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt (unsigned integration_point) const
 
- Public Member Functions inherited from ProcessLib::LocalAssemblerInterface
virtual ~LocalAssemblerInterface ()=default
 
void setInitialConditions (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, bool const use_monolithic_scheme, int const process_id)
 
virtual void initialize (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
 
virtual void preAssemble (double const, double const, std::vector< double > const &)
 
virtual void assemble (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
 
virtual void assembleForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
 
virtual void assembleWithJacobian (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
 
virtual void assembleWithJacobianForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
 
virtual void computeSecondaryVariable (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, int const process_id)
 
virtual void preTimestep (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
 
virtual void postTimestep (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, bool const use_monolithic_scheme, int const process_id)
 
void postNonLinearSolver (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, bool const use_monolithic_scheme, int const process_id)
 
virtual Eigen::Vector3d getFlux (MathLib::Point3d const &, double const, std::vector< double > const &) const
 
virtual Eigen::Vector3d getFlux (MathLib::Point3d const &, double const, std::vector< std::vector< double > > const &) const
 Fits to staggered scheme.
 
- Public Member Functions inherited from NumLib::ExtrapolatableElement
virtual Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix (const unsigned integration_point) const =0
 Provides the shape matrix at the given integration point.
 
virtual ~ExtrapolatableElement ()=default
 

Static Public Attributes

static int const KelvinVectorSize
 
static constexpr auto & N_u_op
 

Private Types

using BMatricesType = BMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim >
 
using IpData = IntegrationPointData< BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, ShapeFunctionDisplacement::NPOINTS >
 

Private Member Functions

std::size_t setIPDataInitialConditions (std::string_view const name, double const *values, int const integration_order) override
 
void setInitialConditionsConcrete (std::vector< double > const &local_x_data, double const t, bool const, int const) override
 
void assemble (double const, double const, std::vector< double > const &, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
 
void assembleWithJacobian (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &, std::vector< double > &, std::vector< double > &local_rhs_data, std::vector< double > &local_Jac_data) override
 
void initializeConcrete () override
 
void postTimestepConcrete (Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, bool const, int const) override
 
void computeSecondaryVariableConcrete (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev) override
 
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix (const unsigned integration_point) const override
 Provides the shape matrix at the given integration point.
 
std::vector< double > const & getIntPtDarcyVelocityGas (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::vector< double > const & getIntPtDarcyVelocityLiquid (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::vector< double > const & getIntPtDiffusionVelocityVapourGas (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::vector< double > const & getIntPtDiffusionVelocityGasGas (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::vector< double > const & getIntPtDiffusionVelocitySoluteLiquid (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::vector< double > const & getIntPtDiffusionVelocityLiquidLiquid (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::tuple< std::vector< ConstitutiveRelations::ConstitutiveData< DisplacementDim > >, std::vector< ConstitutiveRelations::ConstitutiveTempData< DisplacementDim > > > updateConstitutiveVariables (Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, ConstitutiveRelations::ConstitutiveModels< DisplacementDim > const &models)
 
virtual std::vector< double > const & getIntPtLiquidDensity (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtGasDensity (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtSolidDensity (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtVapourPressure (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtPorosity (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtMoleFractionGas (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtMassFractionGas (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtMassFractionLiquid (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtRelativePermeabilityGas (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtRelativePermeabilityLiquid (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtIntrinsicPermeability (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtEnthalpyGas (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtEnthalpyLiquid (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtEnthalpySolid (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 

Private Attributes

std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
 
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
 

Static Private Attributes

static const int gas_pressure_index = 0
 
static const int gas_pressure_size = ShapeFunctionPressure::NPOINTS
 
static const int capillary_pressure_index = ShapeFunctionPressure::NPOINTS
 
static const int capillary_pressure_size = ShapeFunctionPressure::NPOINTS
 
static const int temperature_index = 2 * ShapeFunctionPressure::NPOINTS
 
static const int temperature_size = ShapeFunctionPressure::NPOINTS
 
static const int displacement_index = ShapeFunctionPressure::NPOINTS * 3
 
static const int displacement_size
 
static const int C_index = 0
 
static const int C_size = ShapeFunctionPressure::NPOINTS
 
static const int W_index = ShapeFunctionPressure::NPOINTS
 
static const int W_size = ShapeFunctionPressure::NPOINTS
 

Additional Inherited Members

- Static Public Member Functions inherited from ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >
static auto getReflectionDataForOutput ()
 
- Public Attributes inherited from ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >
TH2MProcessData< DisplacementDim > & process_data_
 
std::vector< typename ConstitutiveRelations::StatefulData< DisplacementDim > > current_states_
 
std::vector< typename ConstitutiveRelations::StatefulDataPrev< DisplacementDim > > prev_states_
 
std::vector< ConstitutiveRelations::MaterialStateData< DisplacementDim > > material_states_
 
std::vector< ConstitutiveRelations::OutputData< DisplacementDim > > output_data_
 
NumLib::GenericIntegrationMethod const & integration_method_
 
MeshLib::Element const & element_
 
bool const is_axially_symmetric_
 
ConstitutiveRelations::SolidConstitutiveRelation< DisplacementDim > const & solid_material_
 

Member Typedef Documentation

◆ BMatricesType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::BMatricesType = BMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>
private

Definition at line 376 of file TH2MFEM.h.

◆ GlobalDimMatrixType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::GlobalDimMatrixType = typename ShapeMatricesTypePressure::GlobalDimMatrixType

Definition at line 65 of file TH2MFEM.h.

◆ GlobalDimVectorType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::GlobalDimVectorType = typename ShapeMatricesTypePressure::GlobalDimVectorType

Definition at line 68 of file TH2MFEM.h.

◆ Invariants

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>

Definition at line 79 of file TH2MFEM.h.

◆ IpData

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::IpData = IntegrationPointData<BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, ShapeFunctionDisplacement::NPOINTS>
private

Definition at line 378 of file TH2MFEM.h.

◆ MatrixType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
template<int M, int N>
using ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::MatrixType = typename ShapeMatricesTypePressure::template MatrixType<M, N>

Definition at line 62 of file TH2MFEM.h.

◆ ShapeMatricesTypeDisplacement

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ShapeMatricesTypeDisplacement = ShapeMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>

Definition at line 51 of file TH2MFEM.h.

◆ ShapeMatricesTypePressure

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ShapeMatricesTypePressure = ShapeMatrixPolicyType<ShapeFunctionPressure, DisplacementDim>

Definition at line 54 of file TH2MFEM.h.

◆ SymmetricTensor

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>

Definition at line 73 of file TH2MFEM.h.

◆ VectorType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
template<int N>
using ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::VectorType = typename ShapeMatricesTypePressure::template VectorType<N>

Definition at line 58 of file TH2MFEM.h.

Constructor & Destructor Documentation

◆ TH2MLocalAssembler() [1/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::TH2MLocalAssembler ( TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > const &  )
delete

◆ TH2MLocalAssembler() [2/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::TH2MLocalAssembler ( TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > &&  )
delete

◆ TH2MLocalAssembler() [3/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::TH2MLocalAssembler ( MeshLib::Element const &  e,
std::size_t const  ,
NumLib::GenericIntegrationMethod const &  integration_method,
bool const  is_axially_symmetric,
TH2MProcessData< DisplacementDim > &  process_data 
)

Definition at line 36 of file TH2MFEM-impl.h.

43 : LocalAssemblerInterface<DisplacementDim>(
44 e, integration_method, is_axially_symmetric, process_data)
45{
46 unsigned const n_integration_points =
48
49 _ip_data.resize(n_integration_points);
50 _secondary_data.N_u.resize(n_integration_points);
51
52 auto const shape_matrices_u =
53 NumLib::initShapeMatrices<ShapeFunctionDisplacement,
55 DisplacementDim>(e, is_axially_symmetric,
57
58 auto const shape_matrices_p =
59 NumLib::initShapeMatrices<ShapeFunctionPressure,
60 ShapeMatricesTypePressure, DisplacementDim>(
61 e, is_axially_symmetric, this->integration_method_);
62
63 for (unsigned ip = 0; ip < n_integration_points; ip++)
64 {
65 auto& ip_data = _ip_data[ip];
66 auto const& sm_u = shape_matrices_u[ip];
67 ip_data.integration_weight =
69 sm_u.integralMeasure * sm_u.detJ;
70
71 ip_data.N_u = sm_u.N;
72 ip_data.dNdx_u = sm_u.dNdx;
73
74 ip_data.N_p = shape_matrices_p[ip].N;
75 ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
76
77 _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
78 }
79}
double getWeight() const
Definition: WeightedPoint.h:80
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
Definition: TH2MFEM.h:55
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
Definition: TH2MFEM.h:386
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
Definition: TH2MFEM.h:52
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
Definition: TH2MFEM.h:382
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u
NumLib::GenericIntegrationMethod const & integration_method_

References ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data, ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_secondary_data, NumLib::GenericIntegrationMethod::getNumberOfPoints(), MathLib::WeightedPoint::getWeight(), NumLib::GenericIntegrationMethod::getWeightedPoint(), NumLib::initShapeMatrices(), ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >::integration_method_, and ProcessLib::HydroMechanics::SecondaryData< ShapeMatrixType >::N_u.

Member Function Documentation

◆ assemble()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assemble ( double const  t,
double const  dt,
std::vector< double > const &  local_x,
std::vector< double > const &  local_x_prev,
std::vector< double > &  local_M_data,
std::vector< double > &  local_K_data,
std::vector< double > &  local_rhs_data 
)
overrideprivatevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 1022 of file TH2MFEM-impl.h.

1028{
1029 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
1031 assert(local_x.size() == matrix_size);
1032
1033 auto const capillary_pressure =
1034 Eigen::Map<VectorType<capillary_pressure_size> const>(
1036
1037 auto const capillary_pressure_prev =
1038 Eigen::Map<VectorType<capillary_pressure_size> const>(
1039 local_x_prev.data() + capillary_pressure_index,
1041
1042 // pointer to local_M_data vector
1043 auto local_M =
1044 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1045 local_M_data, matrix_size, matrix_size);
1046
1047 // pointer to local_K_data vector
1048 auto local_K =
1049 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1050 local_K_data, matrix_size, matrix_size);
1051
1052 // pointer to local_rhs_data vector
1053 auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
1054 local_rhs_data, matrix_size);
1055
1056 // component-formulation
1057 // W - liquid phase main component
1058 // C - gas phase main component
1059 // pointer-matrices to the mass matrix - C component equation
1060 auto MCpG = local_M.template block<C_size, gas_pressure_size>(
1062 auto MCpC = local_M.template block<C_size, capillary_pressure_size>(
1064 auto MCT = local_M.template block<C_size, temperature_size>(
1066 auto MCu = local_M.template block<C_size, displacement_size>(
1068
1069 // pointer-matrices to the stiffness matrix - C component equation
1070 auto LCpG = local_K.template block<C_size, gas_pressure_size>(
1072 auto LCpC = local_K.template block<C_size, capillary_pressure_size>(
1074 auto LCT = local_K.template block<C_size, temperature_size>(
1076
1077 // pointer-matrices to the mass matrix - W component equation
1078 auto MWpG = local_M.template block<W_size, gas_pressure_size>(
1080 auto MWpC = local_M.template block<W_size, capillary_pressure_size>(
1082 auto MWT = local_M.template block<W_size, temperature_size>(
1084 auto MWu = local_M.template block<W_size, displacement_size>(
1086
1087 // pointer-matrices to the stiffness matrix - W component equation
1088 auto LWpG = local_K.template block<W_size, gas_pressure_size>(
1090 auto LWpC = local_K.template block<W_size, capillary_pressure_size>(
1092 auto LWT = local_K.template block<W_size, temperature_size>(
1094
1095 // pointer-matrices to the mass matrix - temperature equation
1096 auto MTu = local_M.template block<temperature_size, displacement_size>(
1098
1099 // pointer-matrices to the stiffness matrix - temperature equation
1100 auto KTT = local_K.template block<temperature_size, temperature_size>(
1102
1103 // pointer-matrices to the stiffness matrix - displacement equation
1104 auto KUpG = local_K.template block<displacement_size, gas_pressure_size>(
1106 auto KUpC =
1107 local_K.template block<displacement_size, capillary_pressure_size>(
1109
1110 // pointer-vectors to the right hand side terms - C-component equation
1111 auto fC = local_f.template segment<C_size>(C_index);
1112 // pointer-vectors to the right hand side terms - W-component equation
1113 auto fW = local_f.template segment<W_size>(W_index);
1114 // pointer-vectors to the right hand side terms - temperature equation
1115 auto fT = local_f.template segment<temperature_size>(temperature_index);
1116 // pointer-vectors to the right hand side terms - displacement equation
1117 auto fU = local_f.template segment<displacement_size>(displacement_index);
1118
1119 unsigned const n_integration_points =
1121
1122 ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const models{
1123 this->solid_material_};
1124
1125 auto const [ip_constitutive_data, ip_constitutive_variables] =
1127 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1128 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1129 local_x_prev.size()),
1130 t, dt, models);
1131
1132 for (unsigned int_point = 0; int_point < n_integration_points; int_point++)
1133 {
1134 auto& ip = _ip_data[int_point];
1135 auto& ip_cv = ip_constitutive_variables[int_point];
1136 auto& current_state = this->current_states_[int_point];
1137 auto const& prev_state = this->prev_states_[int_point];
1138
1139 auto const& Np = ip.N_p;
1140 auto const& NT = Np;
1141 auto const& Nu = ip.N_u;
1143 std::nullopt, this->element_.getID(), int_point,
1145 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
1147 this->element_, Nu))};
1148
1149 auto const& NpT = Np.transpose().eval();
1150 auto const& NTT = NT.transpose().eval();
1151
1152 auto const& gradNp = ip.dNdx_p;
1153 auto const& gradNT = gradNp;
1154 auto const& gradNu = ip.dNdx_u;
1155
1156 auto const& gradNpT = gradNp.transpose().eval();
1157 auto const& gradNTT = gradNT.transpose().eval();
1158
1159 auto const& w = ip.integration_weight;
1160
1161 auto const& m = Invariants::identity2;
1162
1163 auto const mT = m.transpose().eval();
1164
1165 auto const x_coord =
1166 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
1168 this->element_, Nu);
1169
1170 auto const Bu =
1171 LinearBMatrix::computeBMatrix<DisplacementDim,
1172 ShapeFunctionDisplacement::NPOINTS,
1174 gradNu, Nu, x_coord, this->is_axially_symmetric_);
1175
1176 auto const BuT = Bu.transpose().eval();
1177
1178 double const pCap = Np.dot(capillary_pressure);
1179 double const pCap_prev = Np.dot(capillary_pressure_prev);
1180
1181 double const beta_T_SR = ip_cv.s_therm_exp_data.beta_T_SR;
1182
1183 auto const I =
1184 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
1185
1186 const double sD_G = ip.diffusion_coefficient_vapour;
1187 const double sD_L = ip.diffusion_coefficient_solute;
1188
1189 GlobalDimMatrixType const D_C_G = sD_G * I;
1190 GlobalDimMatrixType const D_W_G = sD_G * I;
1191 GlobalDimMatrixType const D_C_L = sD_L * I;
1192 GlobalDimMatrixType const D_W_L = sD_L * I;
1193
1194 auto& k_S = ip.k_S;
1195
1196 auto const s_L = current_state.S_L_data.S_L;
1197 auto const s_G = 1. - s_L;
1198 auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
1199
1200 auto& alpha_B = ip_cv.biot_data();
1201 auto& beta_p_SR = ip_cv.beta_p_SR();
1202
1203 auto const& b = this->process_data_.specific_body_force;
1204
1205 // porosity
1206 auto& phi = ip.phi;
1207
1208 // volume fraction
1209 auto const phi_G = s_G * phi;
1210 auto const phi_L = s_L * phi;
1211 auto const phi_S = 1. - phi;
1212
1213 // solid phase density
1214 auto& rho_SR = ip.rhoSR;
1215 // effective density
1216 auto const rho = phi_G * ip.rhoGR + phi_L * ip.rhoLR + phi_S * rho_SR;
1217
1218 // abbreviations
1219 const double rho_C_FR = s_G * ip.rhoCGR + s_L * ip.rhoCLR;
1220 const double rho_W_FR = s_G * ip.rhoWGR + s_L * ip.rhoWLR;
1221
1222 // phase specific enthalpies
1223 auto& h_G = ip.h_G;
1224 auto& h_L = ip.h_L;
1225
1226 auto const rho_C_GR_dot = (ip.rhoCGR - ip.rhoCGR_prev) / dt;
1227 auto const rho_C_LR_dot = (ip.rhoCLR - ip.rhoCLR_prev) / dt;
1228 auto const rho_W_GR_dot = (ip.rhoWGR - ip.rhoWGR_prev) / dt;
1229 auto const rho_W_LR_dot = (ip.rhoWLR - ip.rhoWLR_prev) / dt;
1230
1231 auto const rho_h_eff = ip.rho_G_h_G + ip.rho_L_h_L + ip.rho_S_h_S;
1232
1233 auto const rho_u_eff_dot = (ip.rho_u_eff - ip.rho_u_eff_prev) / dt;
1234
1235 GlobalDimMatrixType const k_over_mu_G = k_S * ip.k_rel_G / ip.muGR;
1236 GlobalDimMatrixType const k_over_mu_L = k_S * ip.k_rel_L / ip.muLR;
1237
1238 // ---------------------------------------------------------------------
1239 // C-component equation
1240 // ---------------------------------------------------------------------
1241
1242 MCpG.noalias() += NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * Np * w;
1243 MCpC.noalias() -=
1244 NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1245
1246 if (this->process_data_.apply_mass_lumping)
1247 {
1248 if (pCap - pCap_prev != 0.) // avoid division by Zero
1249 {
1250 MCpC.noalias() +=
1251 NpT *
1252 (phi * (ip.rhoCLR - ip.rhoCGR) -
1253 rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1254 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1255 }
1256 }
1257
1258 MCT.noalias() -= NpT * rho_C_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1259 MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
1260
1261 using DisplacementDimMatrix =
1262 Eigen::Matrix<double, DisplacementDim, DisplacementDim>;
1263
1264 DisplacementDimMatrix const advection_C_G = ip.rhoCGR * k_over_mu_G;
1265 DisplacementDimMatrix const advection_C_L = ip.rhoCLR * k_over_mu_L;
1266
1267 DisplacementDimMatrix const diffusion_CGpGR =
1268 -phi_G * ip.rhoGR * D_C_G * ip.dxmWG_dpGR;
1269 DisplacementDimMatrix const diffusion_CLpGR =
1270 -phi_L * ip.rhoLR * D_C_L * ip.dxmWL_dpGR;
1271
1272 DisplacementDimMatrix const diffusion_CGpCap =
1273 -phi_G * ip.rhoGR * D_C_G * ip.dxmWG_dpCap;
1274 DisplacementDimMatrix const diffusion_CLpCap =
1275 -phi_L * ip.rhoLR * D_C_L * ip.dxmWL_dpCap;
1276
1277 DisplacementDimMatrix const diffusion_CGT =
1278 -phi_G * ip.rhoGR * D_C_G * ip.dxmWG_dT;
1279 DisplacementDimMatrix const diffusion_CLT =
1280 -phi_L * ip.rhoLR * D_C_L * ip.dxmWL_dT;
1281
1282 DisplacementDimMatrix const advection_C = advection_C_G + advection_C_L;
1283 DisplacementDimMatrix const diffusion_C_pGR =
1284 diffusion_CGpGR + diffusion_CLpGR;
1285 DisplacementDimMatrix const diffusion_C_pCap =
1286 diffusion_CGpCap + diffusion_CLpCap;
1287
1288 DisplacementDimMatrix const diffusion_C_T =
1289 diffusion_CGT + diffusion_CLT;
1290
1291 LCpG.noalias() +=
1292 gradNpT * (advection_C + diffusion_C_pGR) * gradNp * w;
1293
1294 LCpC.noalias() +=
1295 gradNpT * (diffusion_C_pCap - advection_C_L) * gradNp * w;
1296
1297 LCT.noalias() += gradNpT * (diffusion_C_T)*gradNp * w;
1298
1299 fC.noalias() += gradNpT *
1300 (advection_C_G * ip.rhoGR + advection_C_L * ip.rhoLR) *
1301 b * w;
1302
1303 if (!this->process_data_.apply_mass_lumping)
1304 {
1305 fC.noalias() -= NpT *
1306 (phi * (ip.rhoCLR - ip.rhoCGR) -
1307 rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1308 s_L_dot * w;
1309 }
1310 // fC_III
1311 fC.noalias() -=
1312 NpT * phi * (s_G * rho_C_GR_dot + s_L * rho_C_LR_dot) * w;
1313
1314 // ---------------------------------------------------------------------
1315 // W-component equation
1316 // ---------------------------------------------------------------------
1317
1318 MWpG.noalias() += NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * Np * w;
1319 MWpC.noalias() -=
1320 NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1321
1322 if (this->process_data_.apply_mass_lumping)
1323 {
1324 if (pCap - pCap_prev != 0.) // avoid division by Zero
1325 {
1326 MWpC.noalias() +=
1327 NpT *
1328 (phi * (ip.rhoWLR - ip.rhoWGR) -
1329 rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1330 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1331 }
1332 }
1333
1334 MWT.noalias() -= NpT * rho_W_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1335
1336 MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
1337
1338 DisplacementDimMatrix const advection_W_G = ip.rhoWGR * k_over_mu_G;
1339 DisplacementDimMatrix const advection_W_L = ip.rhoWLR * k_over_mu_L;
1340
1341 DisplacementDimMatrix const diffusion_WGpGR =
1342 phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dpGR;
1343 DisplacementDimMatrix const diffusion_WLpGR =
1344 phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dpGR;
1345
1346 DisplacementDimMatrix const diffusion_WGpCap =
1347 phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dpCap;
1348 DisplacementDimMatrix const diffusion_WLpCap =
1349 phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dpCap;
1350
1351 DisplacementDimMatrix const diffusion_WGT =
1352 phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dT;
1353 DisplacementDimMatrix const diffusion_WLT =
1354 phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dT;
1355
1356 DisplacementDimMatrix const advection_W = advection_W_G + advection_W_L;
1357 DisplacementDimMatrix const diffusion_W_pGR =
1358 diffusion_WGpGR + diffusion_WLpGR;
1359 DisplacementDimMatrix const diffusion_W_pCap =
1360 diffusion_WGpCap + diffusion_WLpCap;
1361
1362 DisplacementDimMatrix const diffusion_W_T =
1363 diffusion_WGT + diffusion_WLT;
1364
1365 LWpG.noalias() +=
1366 gradNpT * (advection_W + diffusion_W_pGR) * gradNp * w;
1367
1368 LWpC.noalias() +=
1369 gradNpT * (diffusion_W_pCap - advection_W_L) * gradNp * w;
1370
1371 LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
1372
1373 fW.noalias() += gradNpT *
1374 (advection_W_G * ip.rhoGR + advection_W_L * ip.rhoLR) *
1375 b * w;
1376
1377 if (!this->process_data_.apply_mass_lumping)
1378 {
1379 fW.noalias() -= NpT *
1380 (phi * (ip.rhoWLR - ip.rhoWGR) -
1381 rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1382 s_L_dot * w;
1383 }
1384
1385 fW.noalias() -=
1386 NpT * phi * (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
1387
1388 // ---------------------------------------------------------------------
1389 // - temperature equation
1390 // ---------------------------------------------------------------------
1391
1392 MTu.noalias() += NTT * rho_h_eff * mT * Bu * w;
1393
1394 KTT.noalias() += gradNTT * ip.lambda * gradNT * w;
1395
1396 fT.noalias() -= NTT * rho_u_eff_dot * w;
1397
1398 fT.noalias() +=
1399 gradNTT * (ip.rhoGR * h_G * ip.w_GS + ip.rhoLR * h_L * ip.w_LS) * w;
1400
1401 fT.noalias() +=
1402 gradNTT *
1403 (ip.rhoCGR * ip.h_CG * ip.d_CG + ip.rhoWGR * ip.h_WG * ip.d_WG) * w;
1404
1405 fT.noalias() +=
1406 NTT *
1407 (ip.rhoGR * ip.w_GS.transpose() + ip.rhoLR * ip.w_LS.transpose()) *
1408 b * w;
1409
1410 // ---------------------------------------------------------------------
1411 // - displacement equation
1412 // ---------------------------------------------------------------------
1413
1414 KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
1415
1416 KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_S_L.chi_S_L * m * Np) * w;
1417
1418 fU.noalias() -= (BuT * current_state.eff_stress_data.sigma -
1419 N_u_op(Nu).transpose() * rho * b) *
1420 w;
1421
1422 if (this->process_data_.apply_mass_lumping)
1423 {
1424 MCpG = MCpG.colwise().sum().eval().asDiagonal();
1425 MCpC = MCpC.colwise().sum().eval().asDiagonal();
1426 MWpG = MWpG.colwise().sum().eval().asDiagonal();
1427 MWpC = MWpC.colwise().sum().eval().asDiagonal();
1428 }
1429 } // int_point-loop
1430}
std::size_t getID() const
Returns the ID of the element.
Definition: Element.h:89
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition: BMatrixPolicy.h:55
static const int gas_pressure_size
Definition: TH2MFEM.h:391
static constexpr auto & N_u_op
Definition: TH2MFEM.h:75
static const int capillary_pressure_index
Definition: TH2MFEM.h:392
typename ShapeMatricesTypePressure::GlobalDimMatrixType GlobalDimMatrixType
Definition: TH2MFEM.h:66
static const int capillary_pressure_size
Definition: TH2MFEM.h:393
static const int temperature_size
Definition: TH2MFEM.h:395
static const int gas_pressure_index
Definition: TH2MFEM.h:390
static const int displacement_size
Definition: TH2MFEM.h:397
std::tuple< std::vector< ConstitutiveRelations::ConstitutiveData< DisplacementDim > >, std::vector< ConstitutiveRelations::ConstitutiveTempData< DisplacementDim > > > updateConstitutiveVariables(Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, ConstitutiveRelations::ConstitutiveModels< DisplacementDim > const &models)
Definition: TH2MFEM-impl.h:88
static const int temperature_index
Definition: TH2MFEM.h:394
static const int displacement_index
Definition: TH2MFEM.h:396
@ rho
density. For some models, rho substitutes p
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
Definition: LinearBMatrix.h:42
static Eigen::Matrix< double, KelvinVectorSize, 1 > const identity2
Kelvin mapping of 2nd order identity tensor.
Definition: KelvinVector.h:107
TH2MProcessData< DisplacementDim > & process_data_
ConstitutiveRelations::SolidConstitutiveRelation< DisplacementDim > const & solid_material_
std::vector< typename ConstitutiveRelations::StatefulData< DisplacementDim > > current_states_
std::vector< typename ConstitutiveRelations::StatefulDataPrev< DisplacementDim > > prev_states_

References ProcessLib::LinearBMatrix::computeBMatrix(), NumLib::interpolateCoordinates(), and NumLib::interpolateXCoordinate().

◆ assembleWithJacobian()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobian ( double const  t,
double const  dt,
std::vector< double > const &  local_x,
std::vector< double > const &  local_x_prev,
std::vector< double > &  ,
std::vector< double > &  ,
std::vector< double > &  local_rhs_data,
std::vector< double > &  local_Jac_data 
)
overrideprivatevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 1437 of file TH2MFEM-impl.h.

1445{
1446 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
1448 assert(local_x.size() == matrix_size);
1449
1450 auto const temperature = Eigen::Map<VectorType<temperature_size> const>(
1451 local_x.data() + temperature_index, temperature_size);
1452
1453 auto const gas_pressure = Eigen::Map<VectorType<gas_pressure_size> const>(
1454 local_x.data() + gas_pressure_index, gas_pressure_size);
1455
1456 auto const capillary_pressure =
1457 Eigen::Map<VectorType<capillary_pressure_size> const>(
1459
1460 auto const displacement = Eigen::Map<VectorType<displacement_size> const>(
1461 local_x.data() + displacement_index, displacement_size);
1462
1463 auto const gas_pressure_prev =
1464 Eigen::Map<VectorType<gas_pressure_size> const>(
1465 local_x_prev.data() + gas_pressure_index, gas_pressure_size);
1466
1467 auto const capillary_pressure_prev =
1468 Eigen::Map<VectorType<capillary_pressure_size> const>(
1469 local_x_prev.data() + capillary_pressure_index,
1471
1472 auto const temperature_prev =
1473 Eigen::Map<VectorType<temperature_size> const>(
1474 local_x_prev.data() + temperature_index, temperature_size);
1475
1476 auto const displacement_prev =
1477 Eigen::Map<VectorType<displacement_size> const>(
1478 local_x_prev.data() + displacement_index, displacement_size);
1479
1480 auto local_Jac =
1481 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1482 local_Jac_data, matrix_size, matrix_size);
1483
1484 auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
1485 local_rhs_data, matrix_size);
1486
1487 // component-formulation
1488 // W - liquid phase main component
1489 // C - gas phase main component
1490
1491 // C component equation matrices
1492 MatrixType<C_size, gas_pressure_size> MCpG =
1493 MatrixType<C_size, gas_pressure_size>::Zero(C_size, gas_pressure_size);
1494 MatrixType<C_size, capillary_pressure_size> MCpC =
1495 MatrixType<C_size, capillary_pressure_size>::Zero(
1497 MatrixType<C_size, temperature_size> MCT =
1498 MatrixType<C_size, temperature_size>::Zero(C_size, temperature_size);
1499 MatrixType<C_size, displacement_size> MCu =
1500 MatrixType<C_size, displacement_size>::Zero(C_size, displacement_size);
1501
1502 MatrixType<C_size, gas_pressure_size> LCpG =
1503 MatrixType<C_size, gas_pressure_size>::Zero(C_size, gas_pressure_size);
1504 MatrixType<C_size, capillary_pressure_size> LCpC =
1505 MatrixType<C_size, capillary_pressure_size>::Zero(
1507 MatrixType<C_size, temperature_size> LCT =
1508 MatrixType<C_size, temperature_size>::Zero(C_size, temperature_size);
1509
1510 // mass matrix - W component equation
1511 MatrixType<W_size, gas_pressure_size> MWpG =
1512 MatrixType<W_size, gas_pressure_size>::Zero(W_size, gas_pressure_size);
1513 MatrixType<W_size, capillary_pressure_size> MWpC =
1514 MatrixType<W_size, capillary_pressure_size>::Zero(
1516 MatrixType<W_size, temperature_size> MWT =
1517 MatrixType<W_size, temperature_size>::Zero(W_size, temperature_size);
1518 MatrixType<W_size, displacement_size> MWu =
1519 MatrixType<W_size, displacement_size>::Zero(W_size, displacement_size);
1520
1521 // stiffness matrix - W component equation
1522 MatrixType<W_size, gas_pressure_size> LWpG =
1523 MatrixType<W_size, gas_pressure_size>::Zero(W_size, gas_pressure_size);
1524 MatrixType<W_size, capillary_pressure_size> LWpC =
1525 MatrixType<W_size, capillary_pressure_size>::Zero(
1527 MatrixType<W_size, temperature_size> LWT =
1528 MatrixType<W_size, temperature_size>::Zero(W_size, temperature_size);
1529
1530 // mass matrix - temperature equation
1531 MatrixType<temperature_size, displacement_size> MTu =
1532 MatrixType<temperature_size, displacement_size>::Zero(
1534
1535 // stiffness matrix - temperature equation
1536 MatrixType<temperature_size, temperature_size> KTT =
1537 MatrixType<temperature_size, temperature_size>::Zero(temperature_size,
1539
1540 // stiffness matrices - displacement equation coupling into pressures
1541 MatrixType<displacement_size, gas_pressure_size> KUpG =
1542 MatrixType<displacement_size, gas_pressure_size>::Zero(
1544 MatrixType<displacement_size, capillary_pressure_size> KUpC =
1545 MatrixType<displacement_size, capillary_pressure_size>::Zero(
1547
1548 // pointer-vectors to the right hand side terms - C-component equation
1549 auto fC = local_f.template segment<C_size>(C_index);
1550 // pointer-vectors to the right hand side terms - W-component equation
1551 auto fW = local_f.template segment<W_size>(W_index);
1552 // pointer-vectors to the right hand side terms - temperature equation
1553 auto fT = local_f.template segment<temperature_size>(temperature_index);
1554 // pointer-vectors to the right hand side terms - displacement equation
1555 auto fU = local_f.template segment<displacement_size>(displacement_index);
1556
1557 unsigned const n_integration_points =
1559
1560 ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const models{
1561 this->solid_material_};
1562
1563 auto const [ip_constitutive_data, ip_constitutive_variables] =
1565 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1566 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1567 local_x_prev.size()),
1568 t, dt, models);
1569
1570 for (unsigned int_point = 0; int_point < n_integration_points; int_point++)
1571 {
1572 auto& ip = _ip_data[int_point];
1573 auto& ip_cd = ip_constitutive_data[int_point];
1574 auto& ip_cv = ip_constitutive_variables[int_point];
1575 auto& current_state = this->current_states_[int_point];
1576 auto& prev_state = this->prev_states_[int_point];
1577
1578 auto const& Np = ip.N_p;
1579 auto const& NT = Np;
1580 auto const& Nu = ip.N_u;
1582 std::nullopt, this->element_.getID(), int_point,
1584 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
1586 this->element_, Nu))};
1587
1588 auto const& NpT = Np.transpose().eval();
1589 auto const& NTT = NT.transpose().eval();
1590
1591 auto const& gradNp = ip.dNdx_p;
1592 auto const& gradNT = gradNp;
1593 auto const& gradNu = ip.dNdx_u;
1594
1595 auto const& gradNpT = gradNp.transpose().eval();
1596 auto const& gradNTT = gradNT.transpose().eval();
1597
1598 auto const& w = ip.integration_weight;
1599
1600 auto const& m = Invariants::identity2;
1601 auto const mT = m.transpose().eval();
1602
1603 auto const x_coord =
1604 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
1606 this->element_, Nu);
1607
1608 auto const Bu =
1609 LinearBMatrix::computeBMatrix<DisplacementDim,
1610 ShapeFunctionDisplacement::NPOINTS,
1612 gradNu, Nu, x_coord, this->is_axially_symmetric_);
1613
1614 auto const BuT = Bu.transpose().eval();
1615
1616 double const div_u_dot =
1617 Invariants::trace(Bu * (displacement - displacement_prev) / dt);
1618
1619 double const pGR = Np.dot(gas_pressure);
1620 double const pCap = Np.dot(capillary_pressure);
1621 double const T = NT.dot(temperature);
1622
1623 GlobalDimVectorType const gradpGR = gradNp * gas_pressure;
1624 GlobalDimVectorType const gradpCap = gradNp * capillary_pressure;
1625 GlobalDimVectorType const gradT = gradNT * temperature;
1626
1627 double const pGR_prev = Np.dot(gas_pressure_prev);
1628 double const pCap_prev = Np.dot(capillary_pressure_prev);
1629 double const T_prev = NT.dot(temperature_prev);
1630 double const beta_T_SR = ip_cv.s_therm_exp_data.beta_T_SR;
1631
1632 auto const I =
1633 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
1634
1635 const double sD_G = ip.diffusion_coefficient_vapour;
1636 const double sD_L = ip.diffusion_coefficient_solute;
1637
1638 GlobalDimMatrixType const D_C_G = sD_G * I;
1639 GlobalDimMatrixType const D_W_G = sD_G * I;
1640 GlobalDimMatrixType const D_C_L = sD_L * I;
1641 GlobalDimMatrixType const D_W_L = sD_L * I;
1642
1643 auto& k_S = ip.k_S;
1644
1645 auto& s_L = current_state.S_L_data.S_L;
1646 auto const s_G = 1. - s_L;
1647 auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
1648
1649 auto const alpha_B = ip_cv.biot_data();
1650 auto const beta_p_SR = ip_cv.beta_p_SR();
1651
1652 auto const& b = this->process_data_.specific_body_force;
1653
1654 // porosity
1655 auto& phi = ip.phi;
1656
1657 // volume fraction
1658 auto const phi_G = s_G * phi;
1659 auto const phi_L = s_L * phi;
1660 auto const phi_S = 1. - phi;
1661
1662 // solid phase density
1663 auto& rho_SR = ip.rhoSR;
1664 // effective density
1665 auto const rho = phi_G * ip.rhoGR + phi_L * ip.rhoLR + phi_S * rho_SR;
1666
1667 // abbreviations
1668 const double rho_C_FR = s_G * ip.rhoCGR + s_L * ip.rhoCLR;
1669 const double rho_W_FR = s_G * ip.rhoWGR + s_L * ip.rhoWLR;
1670
1671 // phase specific enthalpies
1672 auto& h_G = ip.h_G;
1673 auto& h_L = ip.h_L;
1674
1675 auto const rho_C_GR_dot = (ip.rhoCGR - ip.rhoCGR_prev) / dt;
1676 auto const rho_C_LR_dot = (ip.rhoCLR - ip.rhoCLR_prev) / dt;
1677 auto const rho_W_GR_dot = (ip.rhoWGR - ip.rhoWGR_prev) / dt;
1678 auto const rho_W_LR_dot = (ip.rhoWLR - ip.rhoWLR_prev) / dt;
1679
1680 auto const rho_h_eff = ip.rho_G_h_G + ip.rho_L_h_L + ip.rho_S_h_S;
1681
1682 auto const rho_u_eff_dot = (ip.rho_u_eff - ip.rho_u_eff_prev) / dt;
1683
1684 GlobalDimMatrixType const k_over_mu_G = k_S * ip.k_rel_G / ip.muGR;
1685 GlobalDimMatrixType const k_over_mu_L = k_S * ip.k_rel_L / ip.muLR;
1686
1687 // ---------------------------------------------------------------------
1688 // C-component equation
1689 // ---------------------------------------------------------------------
1690
1691 MCpG.noalias() += NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * Np * w;
1692 MCpC.noalias() -=
1693 NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1694
1695 if (this->process_data_.apply_mass_lumping)
1696 {
1697 if (pCap - pCap_prev != 0.) // avoid division by Zero
1698 {
1699 MCpC.noalias() +=
1700 NpT *
1701 (phi * (ip.rhoCLR - ip.rhoCGR) -
1702 rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1703 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1704 }
1705 }
1706
1707 MCT.noalias() -= NpT * rho_C_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1708 // d (fC_4_MCT * T_dot)/d T
1709 local_Jac
1710 .template block<C_size, temperature_size>(C_index,
1712 .noalias() += NpT * ip_cv.dfC_4_MCT_dT * (T - T_prev) / dt * NT * w;
1713
1714 MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
1715 // d (fC_4_MCu * u_dot)/d T
1716 local_Jac
1717 .template block<C_size, temperature_size>(C_index,
1719 .noalias() += NpT * ip_cv.dfC_4_MCu_dT * div_u_dot * NT * w;
1720
1721 GlobalDimMatrixType const advection_C_G = ip.rhoCGR * k_over_mu_G;
1722 GlobalDimMatrixType const advection_C_L = ip.rhoCLR * k_over_mu_L;
1723 GlobalDimMatrixType const diffusion_C_G_p =
1724 -phi_G * ip.rhoGR * D_C_G * ip.dxmWG_dpGR;
1725 GlobalDimMatrixType const diffusion_C_L_p =
1726 -phi_L * ip.rhoLR * D_C_L * ip.dxmWL_dpLR;
1727 GlobalDimMatrixType const diffusion_C_G_T =
1728 -phi_G * ip.rhoGR * D_C_G * ip.dxmWG_dT;
1729 GlobalDimMatrixType const diffusion_C_L_T =
1730 -phi_L * ip.rhoLR * D_C_L * ip.dxmWL_dT;
1731
1732 GlobalDimMatrixType const advection_C = advection_C_G + advection_C_L;
1733 GlobalDimMatrixType const diffusion_C_p =
1734 diffusion_C_G_p + diffusion_C_L_p;
1735 GlobalDimMatrixType const diffusion_C_T =
1736 diffusion_C_G_T + diffusion_C_L_T;
1737
1738 LCpG.noalias() += gradNpT * (advection_C + diffusion_C_p) * gradNp * w;
1739
1740 // d (fC_4_LCpG * grad p_GR)/d p_GR
1741 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1742 gradNpT *
1743 (ip_cv.dadvection_C_dp_GR
1744 // + ip_cv.ddiffusion_C_p_dp_GR TODO (naumov)
1745 ) *
1746 gradpGR * Np * w;
1747
1748 // d (fC_4_LCpG * grad p_GR)/d p_cap
1749 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
1750 gradNpT *
1751 (ip_cv.dadvection_C_dp_cap
1752 // + ip_cv.ddiffusion_C_p_dp_GR TODO (naumov)
1753 ) *
1754 gradpGR * Np * w;
1755
1756 // d (fC_4_LCpG * grad p_GR)/d T
1757 local_Jac
1758 .template block<C_size, temperature_size>(C_index,
1760 .noalias() += gradNpT * ip_cv.dfC_4_LCpG_dT * gradpGR * NT * w;
1761
1762 // d (fC_4_MCpG * p_GR_dot)/d p_GR
1763 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1764 NpT * ip_cv.dfC_4_MCpG_dp_GR * (pGR - pGR_prev) / dt * Np * w;
1765
1766 // d (fC_4_MCpG * p_GR_dot)/d T
1767 local_Jac
1768 .template block<C_size, temperature_size>(C_index,
1770 .noalias() +=
1771 NpT * ip_cv.dfC_4_MCpG_dT * (pGR - pGR_prev) / dt * NT * w;
1772
1773 LCpC.noalias() -=
1774 gradNpT * (advection_C_L + diffusion_C_L_p) * gradNp * w;
1775
1776 /* TODO (naumov) This part is not tested by any of the current ctests.
1777 // d (fC_4_LCpC * grad p_cap)/d p_GR
1778 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1779 gradNpT *
1780 (ip_cv.dfC_4_LCpC_a_dp_GR
1781 // + ip_cv.dfC_4_LCpC_d_dp_GR TODO (naumov)
1782 ) *
1783 gradpCap * Np * w;
1784 // d (fC_4_LCpC * grad p_cap)/d p_cap
1785 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
1786 gradNpT *
1787 (ip_cv.dfC_4_LCpC_a_dp_cap
1788 // + ip_cv.dfC_4_LCpC_d_dp_cap TODO (naumov)
1789 ) *
1790 gradpCap * Np * w;
1791
1792 local_Jac
1793 .template block<C_size, temperature_size>(C_index,
1794 temperature_index)
1795 .noalias() += gradNpT *
1796 (ip_cv.dfC_4_LCpC_a_dT
1797 // + ip_cv.dfC_4_LCpC_d_dT TODO (naumov)
1798 ) *
1799 gradpCap * Np * w;
1800 */
1801
1802 LCT.noalias() += gradNpT * diffusion_C_T * gradNp * w;
1803
1804 // fC_1
1805 fC.noalias() += gradNpT *
1806 (advection_C_G * ip.rhoGR + advection_C_L * ip.rhoLR) *
1807 b * w;
1808
1809 if (!this->process_data_.apply_mass_lumping)
1810 {
1811 // fC_2 = \int a * s_L_dot
1812 auto const a = phi * (ip.rhoCLR - ip.rhoCGR) -
1813 rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR;
1814 fC.noalias() -= NpT * a * s_L_dot * w;
1815
1816 local_Jac.template block<C_size, C_size>(C_index, C_index)
1817 .noalias() +=
1818 NpT *
1819 (ip_cv.dfC_2a_dp_GR * s_L_dot /*- a * (ds_L_dp_GR = 0) / dt*/) *
1820 Np * w;
1821
1822 local_Jac.template block<C_size, W_size>(C_index, W_index)
1823 .noalias() +=
1824 NpT *
1825 (ip_cv.dfC_2a_dp_cap * s_L_dot + a * ip_cv.dS_L_dp_cap() / dt) *
1826 Np * w;
1827
1828 local_Jac
1829 .template block<C_size, temperature_size>(C_index,
1831 .noalias() += NpT * ip_cv.dfC_2a_dT * s_L_dot * NT * w;
1832 }
1833 {
1834 // fC_3 = \int phi * a
1835 double const a = s_G * rho_C_GR_dot + s_L * rho_C_LR_dot;
1836 fC.noalias() -= NpT * phi * a * w;
1837
1838 local_Jac.template block<C_size, C_size>(C_index, C_index)
1839 .noalias() += NpT * phi * ip_cv.dfC_3a_dp_GR * Np * w;
1840
1841 local_Jac.template block<C_size, W_size>(C_index, W_index)
1842 .noalias() += NpT * phi * ip_cv.dfC_3a_dp_cap * Np * w;
1843
1844 local_Jac
1845 .template block<C_size, temperature_size>(C_index,
1847 .noalias() += NpT *
1848 (
1849#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
1850 ip.dphi_dT * a +
1851#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
1852 phi * ip_cv.dfC_3a_dT) *
1853 NT * w;
1854 }
1855 // ---------------------------------------------------------------------
1856 // W-component equation
1857 // ---------------------------------------------------------------------
1858
1859 MWpG.noalias() += NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * Np * w;
1860 MWpC.noalias() -=
1861 NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1862
1863 if (this->process_data_.apply_mass_lumping)
1864 {
1865 if (pCap - pCap_prev != 0.) // avoid division by Zero
1866 {
1867 MWpC.noalias() +=
1868 NpT *
1869 (phi * (ip.rhoWLR - ip.rhoWGR) -
1870 rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1871 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1872 }
1873 }
1874
1875 MWT.noalias() -= NpT * rho_W_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1876
1877 MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
1878
1879 GlobalDimMatrixType const advection_W_G = ip.rhoWGR * k_over_mu_G;
1880 GlobalDimMatrixType const advection_W_L = ip.rhoWLR * k_over_mu_L;
1881 GlobalDimMatrixType const diffusion_W_G_p =
1882 phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dpGR;
1883 GlobalDimMatrixType const diffusion_W_L_p =
1884 phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dpLR;
1885 GlobalDimMatrixType const diffusion_W_G_T =
1886 phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dT;
1887 GlobalDimMatrixType const diffusion_W_L_T =
1888 phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dT;
1889
1890 GlobalDimMatrixType const advection_W = advection_W_G + advection_W_L;
1891 GlobalDimMatrixType const diffusion_W_p =
1892 diffusion_W_G_p + diffusion_W_L_p;
1893 GlobalDimMatrixType const diffusion_W_T =
1894 diffusion_W_G_T + diffusion_W_L_T;
1895
1896 LWpG.noalias() += gradNpT * (advection_W + diffusion_W_p) * gradNp * w;
1897
1898 // fW_4 LWpG' parts; LWpG = \int grad (a + d) grad
1899 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1900 gradNpT * (ip_cv.dfW_4_LWpG_a_dp_GR + ip_cv.dfW_4_LWpG_d_dp_GR) *
1901 gradpGR * Np * w;
1902
1903 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1904 gradNpT * (ip_cv.dfW_4_LWpG_a_dp_cap + ip_cv.dfW_4_LWpG_d_dp_cap) *
1905 gradpGR * Np * w;
1906
1907 local_Jac
1908 .template block<W_size, temperature_size>(W_index,
1910 .noalias() += gradNpT *
1911 (ip_cv.dfW_4_LWpG_a_dT + ip_cv.dfW_4_LWpG_d_dT) *
1912 gradpGR * NT * w;
1913
1914 LWpC.noalias() -=
1915 gradNpT * (advection_W_L + diffusion_W_L_p) * gradNp * w;
1916
1917 // fW_4 LWp_cap' parts; LWpC = \int grad (a + d) grad
1918 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() -=
1919 gradNpT * (ip_cv.dfW_4_LWpC_a_dp_GR + ip_cv.dfW_4_LWpC_d_dp_GR) *
1920 gradpCap * Np * w;
1921
1922 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() -=
1923 gradNpT * (ip_cv.dfW_4_LWpC_a_dp_cap + ip_cv.dfW_4_LWpC_d_dp_cap) *
1924 gradpCap * Np * w;
1925
1926 local_Jac
1927 .template block<W_size, temperature_size>(W_index,
1929 .noalias() -= gradNpT *
1930 (ip_cv.dfW_4_LWpC_a_dT + ip_cv.dfW_4_LWpC_d_dT) *
1931 gradpCap * NT * w;
1932
1933 LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
1934
1935 // fW_1
1936 fW.noalias() += gradNpT *
1937 (advection_W_G * ip.rhoGR + advection_W_L * ip.rhoLR) *
1938 b * w;
1939
1940 // fW_2 = \int (f - g) * s_L_dot
1941 if (!this->process_data_.apply_mass_lumping)
1942 {
1943 double const f = phi * (ip.rhoWLR - ip.rhoWGR);
1944 double const g = rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR;
1945
1946 fW.noalias() -= NpT * (f - g) * s_L_dot * w;
1947
1948 local_Jac.template block<W_size, C_size>(W_index, C_index)
1949 .noalias() += NpT * (ip_cv.dfW_2a_dp_GR - ip_cv.dfW_2b_dp_GR) *
1950 s_L_dot * Np * w;
1951
1952 // sign negated because of dp_cap = -dp_LR
1953 // TODO (naumov) Had to change the sign to get equal Jacobian WW
1954 // blocks in A2 Test. Where is the error?
1955 local_Jac.template block<W_size, W_size>(W_index, W_index)
1956 .noalias() +=
1957 NpT *
1958 ((ip_cv.dfW_2a_dp_cap - ip_cv.dfW_2b_dp_cap) * s_L_dot +
1959 (f - g) * ip_cv.dS_L_dp_cap() / dt) *
1960 Np * w;
1961
1962 local_Jac
1963 .template block<W_size, temperature_size>(W_index,
1965 .noalias() +=
1966 NpT * (ip_cv.dfW_2a_dT - ip_cv.dfW_2b_dT) * s_L_dot * Np * w;
1967 }
1968
1969 // fW_3 = \int phi * a
1970 fW.noalias() -=
1971 NpT * phi * (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
1972
1973 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1974 NpT * phi * ip_cv.dfW_3a_dp_GR * Np * w;
1975
1976 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1977 NpT * phi * ip_cv.dfW_3a_dp_cap * Np * w;
1978
1979 local_Jac
1980 .template block<W_size, temperature_size>(W_index,
1982 .noalias() +=
1983 NpT *
1984 (
1985#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
1986 ip.dphi_dT * (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) +
1987#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
1988 phi * ip_cv.dfW_3a_dT) *
1989 NT * w;
1990
1991 // ---------------------------------------------------------------------
1992 // - temperature equation
1993 // ---------------------------------------------------------------------
1994
1995 MTu.noalias() += NTT * rho_h_eff * mT * Bu * w;
1996
1997 // dfT_4/dp_GR
1998 // d (MTu * u_dot)/dp_GR
1999 local_Jac
2000 .template block<temperature_size, C_size>(temperature_index,
2001 C_index)
2002 .noalias() += NTT * ip_cv.drho_h_eff_dp_GR * div_u_dot * NT * w;
2003
2004 // dfT_4/dp_cap
2005 // d (MTu * u_dot)/dp_cap
2006 local_Jac
2007 .template block<temperature_size, W_size>(temperature_index,
2008 W_index)
2009 .noalias() -= NTT * ip_cv.drho_h_eff_dp_cap * div_u_dot * NT * w;
2010
2011 // dfT_4/dT
2012 // d (MTu * u_dot)/dT
2013 local_Jac
2014 .template block<temperature_size, temperature_size>(
2016 .noalias() += NTT * ip_cv.drho_h_eff_dT * div_u_dot * NT * w;
2017
2018 KTT.noalias() += gradNTT * ip.lambda * gradNT * w;
2019
2020 // d KTT/dp_GR * T
2021 // TODO (naumov) always zero if lambda_xR have no derivatives wrt. p_GR.
2022 // dlambda_dp_GR =
2023 // (dphi_G_dp_GR = 0) * lambdaGR + phi_G * dlambda_GR_dp_GR +
2024 // (dphi_L_dp_GR = 0) * lambdaLR + phi_L * dlambda_LR_dp_GR +
2025 // (dphi_S_dp_GR = 0) * lambdaSR + phi_S * dlambda_SR_dp_GR +
2026 // = 0
2027 //
2028 // Since dlambda/dp_GR is 0 the derivative is omitted:
2029 // local_Jac
2030 // .template block<temperature_size, C_size>(temperature_index,
2031 // C_index)
2032 // .noalias() += gradNTT * dlambda_dp_GR * gradT * Np * w;
2033
2034 // d KTT/dp_cap * T
2035 local_Jac
2036 .template block<temperature_size, W_size>(temperature_index,
2037 W_index)
2038 .noalias() += gradNTT * ip_cv.dlambda_dp_cap * gradT * Np * w;
2039
2040 // d KTT/dT * T
2041 local_Jac
2042 .template block<temperature_size, temperature_size>(
2044 .noalias() += gradNTT * ip_cv.dlambda_dT * gradT * NT * w;
2045
2046 // fT_1
2047 fT.noalias() -= NTT * rho_u_eff_dot * w;
2048
2049 // dfT_1/dp_GR
2050 local_Jac
2051 .template block<temperature_size, C_size>(temperature_index,
2052 C_index)
2053 .noalias() += NpT / dt * ip_cv.drho_u_eff_dp_GR * Np * w;
2054
2055 // dfT_1/dp_cap
2056 local_Jac
2057 .template block<temperature_size, W_size>(temperature_index,
2058 W_index)
2059 .noalias() += NpT / dt * ip_cv.drho_u_eff_dp_cap * Np * w;
2060
2061 // dfT_1/dT
2062 // MTT
2063 local_Jac
2064 .template block<temperature_size, temperature_size>(
2066 .noalias() += NTT * ip_cv.drho_u_eff_dT / dt * NT * w;
2067
2068 // fT_2
2069 fT.noalias() +=
2070 gradNTT * (ip.rhoGR * h_G * ip.w_GS + ip.rhoLR * h_L * ip.w_LS) * w;
2071
2072 // dfT_2/dp_GR
2073 local_Jac
2074 .template block<temperature_size, C_size>(temperature_index,
2075 C_index)
2076 .noalias() -=
2077 // dfT_2/dp_GR first part
2078 gradNTT * ip_cv.drho_GR_h_w_eff_dp_GR_Npart * Np * w +
2079 // dfT_2/dp_GR second part
2080 gradNTT * ip_cv.drho_GR_h_w_eff_dp_GR_gradNpart * gradNp * w;
2081
2082 // dfT_2/dp_cap
2083 local_Jac
2084 .template block<temperature_size, W_size>(temperature_index,
2085 W_index)
2086 .noalias() -=
2087 // first part of dfT_2/dp_cap
2088 gradNTT * (-ip_cv.drho_LR_h_w_eff_dp_cap_Npart) * Np * w +
2089 // second part of dfT_2/dp_cap
2090 gradNTT * (-ip_cv.drho_LR_h_w_eff_dp_cap_gradNpart) * gradNp * w;
2091
2092 // dfT_2/dT
2093 local_Jac
2094 .template block<temperature_size, temperature_size>(
2096 .noalias() -= gradNTT * ip_cv.drho_GR_h_w_eff_dT * NT * w;
2097
2098 // fT_3
2099 fT.noalias() +=
2100 NTT *
2101 (ip.rhoGR * ip.w_GS.transpose() + ip.rhoLR * ip.w_LS.transpose()) *
2102 b * w;
2103
2104 fT.noalias() +=
2105 gradNTT *
2106 (ip.rhoCGR * ip.h_CG * ip.d_CG + ip.rhoWGR * ip.h_WG * ip.d_WG) * w;
2107
2108 // ---------------------------------------------------------------------
2109 // - displacement equation
2110 // ---------------------------------------------------------------------
2111
2112 KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
2113
2114 // dfU_2/dp_GR = dKUpG/dp_GR * p_GR + KUpG. The former is zero, the
2115 // latter is handled below.
2116
2117 KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_S_L.chi_S_L * m * Np) * w;
2118
2119 // dfU_2/dp_cap = dKUpC/dp_cap * p_cap + KUpC. The former is handled
2120 // here, the latter below.
2121 local_Jac
2122 .template block<displacement_size, W_size>(displacement_index,
2123 W_index)
2124 .noalias() += BuT * alpha_B * ip_cv.chi_S_L.dchi_dS_L *
2125 ip_cv.dS_L_dp_cap() * pCap * m * Np * w;
2126
2127 local_Jac
2128 .template block<displacement_size, displacement_size>(
2130 .noalias() += BuT * ip_cd.s_mech_data.stiffness_tensor * Bu * w;
2131
2132 // fU_1
2133 fU.noalias() -= (BuT * current_state.eff_stress_data.sigma -
2134 N_u_op(Nu).transpose() * rho * b) *
2135 w;
2136
2137 // KuT
2138 local_Jac
2139 .template block<displacement_size, temperature_size>(
2141 .noalias() -=
2142 BuT *
2143 (ip_cd.s_mech_data.stiffness_tensor *
2144 ip_cv.s_therm_exp_data.solid_linear_thermal_expansivity_vector) *
2145 NT * w;
2146
2147 /* TODO (naumov) Test with gravity needed to check this Jacobian part.
2148 local_Jac
2149 .template block<displacement_size, temperature_size>(
2150 displacement_index, temperature_index)
2151 .noalias() += N_u_op(Nu).transpose() * ip_cv.drho_dT * b *
2152 N_u_op(Nu).transpose() * w;
2153 */
2154
2155 if (this->process_data_.apply_mass_lumping)
2156 {
2157 MCpG = MCpG.colwise().sum().eval().asDiagonal();
2158 MCpC = MCpC.colwise().sum().eval().asDiagonal();
2159 MWpG = MWpG.colwise().sum().eval().asDiagonal();
2160 MWpC = MWpC.colwise().sum().eval().asDiagonal();
2161 }
2162 } // int_point-loop
2163
2164 // --- Gas ---
2165 // fC_4
2166 fC.noalias() -= LCpG * gas_pressure + LCpC * capillary_pressure +
2167 LCT * temperature +
2168 MCpG * (gas_pressure - gas_pressure_prev) / dt +
2169 MCpC * (capillary_pressure - capillary_pressure_prev) / dt +
2170 MCT * (temperature - temperature_prev) / dt +
2171 MCu * (displacement - displacement_prev) / dt;
2172
2173 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
2174 LCpG + MCpG / dt;
2175 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
2176 LCpC + MCpC / dt;
2177 local_Jac
2178 .template block<C_size, temperature_size>(C_index, temperature_index)
2179 .noalias() += LCT + MCT / dt;
2180 local_Jac
2181 .template block<C_size, displacement_size>(C_index, displacement_index)
2182 .noalias() += MCu / dt;
2183
2184 // --- Capillary pressure ---
2185 // fW_4
2186 fW.noalias() -= LWpG * gas_pressure + LWpC * capillary_pressure +
2187 LWT * temperature +
2188 MWpG * (gas_pressure - gas_pressure_prev) / dt +
2189 MWpC * (capillary_pressure - capillary_pressure_prev) / dt +
2190 MWT * (temperature - temperature_prev) / dt +
2191 MWu * (displacement - displacement_prev) / dt;
2192
2193 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
2194 LWpC + MWpC / dt;
2195 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
2196 LWpG + MWpG / dt;
2197 local_Jac
2198 .template block<W_size, temperature_size>(W_index, temperature_index)
2199 .noalias() += LWT + MWT / dt;
2200 local_Jac
2201 .template block<W_size, displacement_size>(W_index, displacement_index)
2202 .noalias() += MWu / dt;
2203
2204 // --- Temperature ---
2205 // fT_4
2206 fT.noalias() -=
2207 KTT * temperature + MTu * (displacement - displacement_prev) / dt;
2208
2209 local_Jac
2210 .template block<temperature_size, temperature_size>(temperature_index,
2212 .noalias() += KTT;
2213 local_Jac
2214 .template block<temperature_size, displacement_size>(temperature_index,
2216 .noalias() += MTu / dt;
2217
2218 // --- Displacement ---
2219 // fU_2
2220 fU.noalias() -= KUpG * gas_pressure + KUpC * capillary_pressure;
2221
2222 local_Jac
2223 .template block<displacement_size, C_size>(displacement_index, C_index)
2224 .noalias() += KUpG;
2225 local_Jac
2226 .template block<displacement_size, W_size>(displacement_index, W_index)
2227 .noalias() += KUpC;
2228}
typename ShapeMatricesTypePressure::GlobalDimVectorType GlobalDimVectorType
Definition: TH2MFEM.h:69
static double trace(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Trace of the corresponding tensor.

References ProcessLib::LinearBMatrix::computeBMatrix(), NumLib::interpolateCoordinates(), and NumLib::interpolateXCoordinate().

◆ computeSecondaryVariableConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::computeSecondaryVariableConcrete ( double const  t,
double const  dt,
Eigen::VectorXd const &  local_x,
Eigen::VectorXd const &  local_x_prev 
)
overrideprivatevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 2389 of file TH2MFEM-impl.h.

2393{
2394 auto const gas_pressure =
2395 local_x.template segment<gas_pressure_size>(gas_pressure_index);
2396 auto const capillary_pressure =
2397 local_x.template segment<capillary_pressure_size>(
2399 auto const liquid_pressure = (gas_pressure - capillary_pressure).eval();
2400
2402 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2403 DisplacementDim>(this->element_, this->is_axially_symmetric_,
2404 gas_pressure,
2405 *this->process_data_.gas_pressure_interpolated);
2406
2408 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2409 DisplacementDim>(this->element_, this->is_axially_symmetric_,
2411 *this->process_data_.capillary_pressure_interpolated);
2412
2414 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2415 DisplacementDim>(this->element_, this->is_axially_symmetric_,
2416 liquid_pressure,
2417 *this->process_data_.liquid_pressure_interpolated);
2418
2419 auto const temperature =
2420 local_x.template segment<temperature_size>(temperature_index);
2421
2423 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2424 DisplacementDim>(this->element_, this->is_axially_symmetric_,
2425 temperature,
2426 *this->process_data_.temperature_interpolated);
2427
2428 unsigned const n_integration_points =
2430
2431 double saturation_avg = 0;
2432
2433 ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const models{
2434 this->solid_material_};
2435
2436 updateConstitutiveVariables(local_x, local_x_prev, t, dt, models);
2437
2438 for (unsigned ip = 0; ip < n_integration_points; ip++)
2439 {
2440 saturation_avg += this->current_states_[ip].S_L_data.S_L;
2441 }
2442 saturation_avg /= n_integration_points;
2443 (*this->process_data_.element_saturation)[this->element_.getID()] =
2444 saturation_avg;
2445}
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
auto eval(Function &f, Tuples &... ts) -> typename detail::GetFunctionReturnType< decltype(&Function::eval)>::type
Definition: Apply.h:267

References NumLib::interpolateToHigherOrderNodes().

◆ getIntPtDarcyVelocityGas()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtDarcyVelocityGas ( const double  t,
std::vector< GlobalVector * > const &  x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  dof_table,
std::vector< double > &  cache 
) const
overrideprivatevirtual

Implements ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >.

Definition at line 2233 of file TH2MFEM-impl.h.

2239{
2240 unsigned const n_integration_points =
2242
2243 cache.clear();
2244 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
2245 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2246 cache, DisplacementDim, n_integration_points);
2247
2248 for (unsigned ip = 0; ip < n_integration_points; ip++)
2249 {
2250 cache_matrix.col(ip).noalias() = _ip_data[ip].w_GS;
2251 }
2252
2253 return cache;
2254}
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:32

References MathLib::createZeroedMatrix().

◆ getIntPtDarcyVelocityLiquid()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtDarcyVelocityLiquid ( const double  t,
std::vector< GlobalVector * > const &  x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  dof_table,
std::vector< double > &  cache 
) const
overrideprivatevirtual

Implements ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >.

Definition at line 2259 of file TH2MFEM-impl.h.

2265{
2266 unsigned const n_integration_points =
2268
2269 cache.clear();
2270 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
2271 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2272 cache, DisplacementDim, n_integration_points);
2273
2274 for (unsigned ip = 0; ip < n_integration_points; ip++)
2275 {
2276 cache_matrix.col(ip).noalias() = _ip_data[ip].w_LS;
2277 }
2278
2279 return cache;
2280}

References MathLib::createZeroedMatrix().

◆ getIntPtDiffusionVelocityGasGas()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtDiffusionVelocityGasGas ( const double  t,
std::vector< GlobalVector * > const &  x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  dof_table,
std::vector< double > &  cache 
) const
overrideprivatevirtual

Implements ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >.

Definition at line 2311 of file TH2MFEM-impl.h.

2317{
2318 unsigned const n_integration_points =
2320
2321 cache.clear();
2322 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
2323 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2324 cache, DisplacementDim, n_integration_points);
2325
2326 for (unsigned ip = 0; ip < n_integration_points; ip++)
2327 {
2328 cache_matrix.col(ip).noalias() = _ip_data[ip].d_CG;
2329 }
2330
2331 return cache;
2332}

References MathLib::createZeroedMatrix().

◆ getIntPtDiffusionVelocityLiquidLiquid()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtDiffusionVelocityLiquidLiquid ( const double  t,
std::vector< GlobalVector * > const &  x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  dof_table,
std::vector< double > &  cache 
) const
overrideprivatevirtual

Implements ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >.

Definition at line 2363 of file TH2MFEM-impl.h.

2369{
2370 unsigned const n_integration_points =
2372
2373 cache.clear();
2374 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
2375 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2376 cache, DisplacementDim, n_integration_points);
2377
2378 for (unsigned ip = 0; ip < n_integration_points; ip++)
2379 {
2380 cache_matrix.col(ip).noalias() = _ip_data[ip].d_WL;
2381 }
2382
2383 return cache;
2384}

References MathLib::createZeroedMatrix().

◆ getIntPtDiffusionVelocitySoluteLiquid()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtDiffusionVelocitySoluteLiquid ( const double  t,
std::vector< GlobalVector * > const &  x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  dof_table,
std::vector< double > &  cache 
) const
overrideprivatevirtual

Implements ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >.

Definition at line 2337 of file TH2MFEM-impl.h.

2343{
2344 unsigned const n_integration_points =
2346
2347 cache.clear();
2348 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
2349 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2350 cache, DisplacementDim, n_integration_points);
2351
2352 for (unsigned ip = 0; ip < n_integration_points; ip++)
2353 {
2354 cache_matrix.col(ip).noalias() = _ip_data[ip].d_CL;
2355 }
2356
2357 return cache;
2358}

References MathLib::createZeroedMatrix().

◆ getIntPtDiffusionVelocityVapourGas()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtDiffusionVelocityVapourGas ( const double  t,
std::vector< GlobalVector * > const &  x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  dof_table,
std::vector< double > &  cache 
) const
overrideprivatevirtual

Implements ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >.

Definition at line 2285 of file TH2MFEM-impl.h.

2291{
2292 unsigned const n_integration_points =
2294
2295 cache.clear();
2296 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
2297 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2298 cache, DisplacementDim, n_integration_points);
2299
2300 for (unsigned ip = 0; ip < n_integration_points; ip++)
2301 {
2302 cache_matrix.col(ip).noalias() = _ip_data[ip].d_WG;
2303 }
2304
2305 return cache;
2306}

References MathLib::createZeroedMatrix().

◆ getIntPtEnthalpyGas()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
virtual std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtEnthalpyGas ( const double  ,
std::vector< GlobalVector * > const &  ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  ,
std::vector< double > &  cache 
) const
inlineoverrideprivatevirtual

◆ getIntPtEnthalpyLiquid()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
virtual std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtEnthalpyLiquid ( const double  ,
std::vector< GlobalVector * > const &  ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  ,
std::vector< double > &  cache 
) const
inlineoverrideprivatevirtual

◆ getIntPtEnthalpySolid()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
virtual std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtEnthalpySolid ( const double  ,
std::vector< GlobalVector * > const &  ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  ,
std::vector< double > &  cache 
) const
inlineoverrideprivatevirtual

◆ getIntPtGasDensity()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
virtual std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtGasDensity ( const double  ,
std::vector< GlobalVector * > const &  ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  ,
std::vector< double > &  cache 
) const
inlineoverrideprivatevirtual

◆ getIntPtIntrinsicPermeability()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
virtual std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtIntrinsicPermeability ( const double  ,
std::vector< GlobalVector * > const &  ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  ,
std::vector< double > &  cache 
) const
inlineoverrideprivatevirtual

◆ getIntPtLiquidDensity()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
virtual std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtLiquidDensity ( const double  ,
std::vector< GlobalVector * > const &  ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  ,
std::vector< double > &  cache 
) const
inlineoverrideprivatevirtual

◆ getIntPtMassFractionGas()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
virtual std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtMassFractionGas ( const double  ,
std::vector< GlobalVector * > const &  ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  ,
std::vector< double > &  cache 
) const
inlineoverrideprivatevirtual

◆ getIntPtMassFractionLiquid()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
virtual std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtMassFractionLiquid ( const double  ,
std::vector< GlobalVector * > const &  ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  ,
std::vector< double > &  cache 
) const
inlineoverrideprivatevirtual

◆ getIntPtMoleFractionGas()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
virtual std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtMoleFractionGas ( const double  ,
std::vector< GlobalVector * > const &  ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  ,
std::vector< double > &  cache 
) const
inlineoverrideprivatevirtual

◆ getIntPtPorosity()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
virtual std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtPorosity ( const double  ,
std::vector< GlobalVector * > const &  ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  ,
std::vector< double > &  cache 
) const
inlineoverrideprivatevirtual

◆ getIntPtRelativePermeabilityGas()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
virtual std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtRelativePermeabilityGas ( const double  ,
std::vector< GlobalVector * > const &  ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  ,
std::vector< double > &  cache 
) const
inlineoverrideprivatevirtual

◆ getIntPtRelativePermeabilityLiquid()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
virtual std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtRelativePermeabilityLiquid ( const double  ,
std::vector< GlobalVector * > const &  ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  ,
std::vector< double > &  cache 
) const
inlineoverrideprivatevirtual

◆ getIntPtSolidDensity()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
virtual std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtSolidDensity ( const double  ,
std::vector< GlobalVector * > const &  ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  ,
std::vector< double > &  cache 
) const
inlineoverrideprivatevirtual

◆ getIntPtVapourPressure()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
virtual std::vector< double > const & ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtVapourPressure ( const double  ,
std::vector< GlobalVector * > const &  ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  ,
std::vector< double > &  cache 
) const
inlineoverrideprivatevirtual

◆ getShapeMatrix()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
Eigen::Map< const Eigen::RowVectorXd > ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getShapeMatrix ( const unsigned  integration_point) const
inlineoverrideprivatevirtual

Provides the shape matrix at the given integration point.

Implements NumLib::ExtrapolatableElement.

Definition at line 186 of file TH2MFEM.h.

188 {
189 auto const& N_u = _secondary_data.N_u[integration_point];
190
191 // assumes N is stored contiguously in memory
192 return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
193 }

References ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_secondary_data, and ProcessLib::HydroMechanics::SecondaryData< ShapeMatrixType >::N_u.

◆ initializeConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::initializeConcrete ( )
inlineoverrideprivatevirtual

Set initial stress from parameter.

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 118 of file TH2MFEM.h.

119 {
120 unsigned const n_integration_points =
122
123 for (unsigned ip = 0; ip < n_integration_points; ip++)
124 {
125 auto& ip_data = _ip_data[ip];
126
127 ParameterLib::SpatialPosition const x_position{
128 std::nullopt, this->element_.getID(), ip,
130 ShapeFunctionDisplacement,
132 ip_data.N_u))};
133
135 if (this->process_data_.initial_stress.value)
136 {
137 this->current_states_[ip].eff_stress_data.sigma.noalias() =
139 DisplacementDim>(
140 (*this->process_data_.initial_stress.value)(
141 std::numeric_limits<
142 double>::quiet_NaN() /* time independent */,
143 x_position));
144 }
145
146 double const t = 0; // TODO (naumov) pass t from top
147 auto& material_state = this->material_states_[ip];
148 this->solid_material_.initializeInternalStateVariables(
149 t, x_position, *material_state.material_state_variables);
150
151 ip_data.pushBackState();
152 material_state.pushBackState();
153 }
154
155 for (unsigned ip = 0; ip < n_integration_points; ip++)
156 {
157 this->prev_states_[ip] = this->current_states_[ip];
158 }
159 }
Eigen::Matrix< double, Eigen::MatrixBase< Derived >::RowsAtCompileTime, 1 > symmetricTensorToKelvinVector(Eigen::MatrixBase< Derived > const &v)
Definition: KelvinVector.h:201
std::vector< ConstitutiveRelations::MaterialStateData< DisplacementDim > > material_states_

References ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data, ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >::current_states_, ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >::element_, MeshLib::Element::getID(), NumLib::GenericIntegrationMethod::getNumberOfPoints(), ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >::integration_method_, NumLib::interpolateCoordinates(), ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >::material_states_, ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >::prev_states_, ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >::process_data_, ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >::solid_material_, and MathLib::KelvinVector::symmetricTensorToKelvinVector().

◆ postTimestepConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::postTimestepConcrete ( Eigen::VectorXd const &  ,
Eigen::VectorXd const &  ,
double const  ,
double const  ,
bool const  ,
int const   
)
inlineoverrideprivatevirtual

◆ setInitialConditionsConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::setInitialConditionsConcrete ( std::vector< double > const &  local_x_data,
double const  t,
bool const  ,
int const   
)
overrideprivatevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 886 of file TH2MFEM-impl.h.

891{
892 [[maybe_unused]] auto const matrix_size =
895
896 assert(local_x_data.size() == matrix_size);
897
898 auto const local_x = Eigen::Map<Eigen::VectorXd const>(local_x_data.data(),
899 local_x_data.size());
900 auto const capillary_pressure =
901 local_x.template segment<capillary_pressure_size>(
903
904 auto const p_GR =
905 local_x.template segment<gas_pressure_size>(gas_pressure_index);
906
907 auto const temperature =
908 local_x.template segment<temperature_size>(temperature_index);
909
910 auto const displacement =
911 local_x.template segment<displacement_size>(displacement_index);
912
913 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
914 auto const& medium =
915 *this->process_data_.media_map.getMedium(this->element_.getID());
916 auto const& solid_phase = medium.phase("Solid");
917
918 ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const models{
919 this->solid_material_};
920
921 unsigned const n_integration_points =
923
924 for (unsigned ip = 0; ip < n_integration_points; ip++)
925 {
927
928 auto& ip_data = _ip_data[ip];
929 auto& prev_state = this->prev_states_[ip];
930 auto const& Np = ip_data.N_p;
931 auto const& NT = Np;
932 auto const& Nu = ip_data.N_u;
933 auto const& gradNu = ip_data.dNdx_u;
934 auto const x_coord =
935 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
937 this->element_, Nu);
939 std::nullopt, this->element_.getID(), ip,
941 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
943 this->element_, ip_data.N_u))};
944
945 double const pCap = Np.dot(capillary_pressure);
946 vars.capillary_pressure = pCap;
947
948 double const T = NT.dot(temperature);
949 TemperatureData const T_data{T, T}; // T_prev = T in initialization.
950 vars.temperature = T;
951
952 auto const Bu =
953 LinearBMatrix::computeBMatrix<DisplacementDim,
954 ShapeFunctionDisplacement::NPOINTS,
956 gradNu, Nu, x_coord, this->is_axially_symmetric_);
957
958 auto& eps = this->output_data_[ip].eps_data.eps;
959 eps.noalias() = Bu * displacement;
960
961 // Set volumetric strain rate for the general case without swelling.
963
964 double const S_L =
965 medium.property(MPL::PropertyType::saturation)
966 .template value<double>(
967 vars, pos, t, std::numeric_limits<double>::quiet_NaN());
968 this->prev_states_[ip].S_L_data->S_L = S_L;
969
970 // TODO (naumov) Double computation of C_el might be avoided if
971 // updateConstitutiveVariables is called before. But it might interfere
972 // with eps_m initialization.
973 ConstitutiveRelations::ElasticTangentStiffnessData<DisplacementDim>
974 C_el_data;
975 models.elastic_tangent_stiffness_model.eval({pos, t, dt}, T_data,
976 C_el_data);
977 auto const& C_el = C_el_data.stiffness_tensor;
978
979 // Set eps_m_prev from potentially non-zero eps and sigma_sw from
980 // restart.
981 auto const& sigma_sw = this->current_states_[ip].swelling_data.sigma_sw;
982 prev_state.mechanical_strain_data->eps_m.noalias() =
983 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
984 ? eps + C_el.inverse() * sigma_sw
985 : eps;
986
987 if (this->process_data_.initial_stress.isTotalStress())
988 {
989 auto const alpha_b =
990 medium.property(MPL::PropertyType::biot_coefficient)
991 .template value<double>(vars, pos, t, 0.0 /*dt*/);
992
993 vars.liquid_saturation = S_L;
994 double const bishop =
995 medium.property(MPL::PropertyType::bishops_effective_stress)
996 .template value<double>(vars, pos, t, 0.0 /*dt*/);
997
998 this->current_states_[ip].eff_stress_data.sigma.noalias() +=
999 alpha_b * Np.dot(p_GR - bishop * capillary_pressure) *
1001 this->prev_states_[ip].eff_stress_data =
1002 this->current_states_[ip].eff_stress_data;
1003 }
1004 }
1005
1006 // local_x_prev equal to local_x s.t. the local_x_dot is zero.
1007 updateConstitutiveVariables(local_x, local_x, t, 0, models);
1008
1009 for (unsigned ip = 0; ip < n_integration_points; ip++)
1010 {
1011 auto& ip_data = _ip_data[ip];
1012 ip_data.pushBackState();
1013
1014 this->material_states_[ip].pushBackState();
1015 }
1016}
std::vector< ConstitutiveRelations::OutputData< DisplacementDim > > output_data_

References MaterialPropertyLib::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), NumLib::interpolateCoordinates(), NumLib::interpolateXCoordinate(), MaterialPropertyLib::VariableArray::liquid_saturation, ProcessLib::TH2M::ConstitutiveRelations::ElasticTangentStiffnessData< DisplacementDim >::stiffness_tensor, MaterialPropertyLib::VariableArray::temperature, and MaterialPropertyLib::VariableArray::volumetric_strain.

◆ setIPDataInitialConditions()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::size_t ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::setIPDataInitialConditions ( std::string_view const  name,
double const *  values,
int const  integration_order 
)
overrideprivatevirtual
Returns
the number of read integration points.

Implements ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >.

Definition at line 826 of file TH2MFEM-impl.h.

829{
830 if (integration_order !=
831 static_cast<int>(this->integration_method_.getIntegrationOrder()))
832 {
833 OGS_FATAL(
834 "Setting integration point initial conditions; The integration "
835 "order of the local assembler for element {:d} is different "
836 "from the integration order in the initial condition.",
837 this->element_.getID());
838 }
839
840 if (name == "sigma" && this->process_data_.initial_stress.value)
841 {
842 OGS_FATAL(
843 "Setting initial conditions for stress from integration "
844 "point data and from a parameter '{:s}' is not possible "
845 "simultaneously.",
846 this->process_data_.initial_stress.value->name);
847 }
848
849 if (name.starts_with("material_state_variable_"))
850 {
851 name.remove_prefix(24);
852 DBUG("Setting material state variable '{:s}'", name);
853
854 auto const& internal_variables =
855 this->solid_material_.getInternalVariables();
856 if (auto const iv = std::find_if(
857 begin(internal_variables), end(internal_variables),
858 [&name](auto const& iv) { return iv.name == name; });
859 iv != end(internal_variables))
860 {
861 DBUG("Setting material state variable '{:s}'", name);
863 values, this->material_states_,
864 &ConstitutiveRelations::MaterialStateData<
865 DisplacementDim>::material_state_variables,
866 iv->reference);
867 }
868
869 WARN(
870 "Could not find variable {:s} in solid material model's "
871 "internal variables.",
872 name);
873 return 0;
874 }
875
876 // TODO this logic could be pulled out of the local assembler into the
877 // process. That might lead to a slightly better performance due to less
878 // string comparisons.
879 return ProcessLib::Reflection::reflectSetIPData<DisplacementDim>(
880 name, values, this->current_states_);
881}
#define OGS_FATAL(...)
Definition: Error.h:26
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:30
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:40
std::size_t setIntegrationPointDataMaterialStateVariables(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType member, std::function< std::span< double >(MaterialStateVariables &)> get_values_span)

References DBUG(), OGS_FATAL, ProcessLib::setIntegrationPointDataMaterialStateVariables(), and WARN().

◆ updateConstitutiveVariables()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::tuple< std::vector< ConstitutiveRelations::ConstitutiveData< DisplacementDim > >, std::vector< ConstitutiveRelations::ConstitutiveTempData< DisplacementDim > > > ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::updateConstitutiveVariables ( Eigen::VectorXd const &  local_x,
Eigen::VectorXd const &  local_x_prev,
double const  t,
double const  dt,
ConstitutiveRelations::ConstitutiveModels< DisplacementDim > const &  models 
)
private

Definition at line 87 of file TH2MFEM-impl.h.

93{
94 [[maybe_unused]] auto const matrix_size =
97
98 assert(local_x.size() == matrix_size);
99
100 auto const gas_pressure =
101 local_x.template segment<gas_pressure_size>(gas_pressure_index);
102 auto const capillary_pressure =
103 local_x.template segment<capillary_pressure_size>(
105
106 auto const temperature =
107 local_x.template segment<temperature_size>(temperature_index);
108 auto const temperature_prev =
109 local_x_prev.template segment<temperature_size>(temperature_index);
110
111 auto const displacement =
112 local_x.template segment<displacement_size>(displacement_index);
113 auto const displacement_prev =
114 local_x_prev.template segment<displacement_size>(displacement_index);
115
116 auto const& medium =
117 *this->process_data_.media_map.getMedium(this->element_.getID());
118 auto const& gas_phase = medium.phase("Gas");
119 auto const& liquid_phase = medium.phase("AqueousLiquid");
120 auto const& solid_phase = medium.phase("Solid");
121 MediaData media_data{medium};
122
123 unsigned const n_integration_points =
125
126 std::vector<ConstitutiveRelations::ConstitutiveData<DisplacementDim>>
127 ip_constitutive_data(n_integration_points);
128 std::vector<ConstitutiveRelations::ConstitutiveTempData<DisplacementDim>>
129 ip_constitutive_variables(n_integration_points);
130
131 PhaseTransitionModelVariables ptmv;
132
133 for (unsigned ip = 0; ip < n_integration_points; ip++)
134 {
135 auto& ip_data = _ip_data[ip];
136 auto& ip_cv = ip_constitutive_variables[ip];
137 auto& ip_cd = ip_constitutive_data[ip];
138 auto& current_state = this->current_states_[ip];
139 auto& prev_state = this->prev_states_[ip];
140
141 auto const& Np = ip_data.N_p;
142 auto const& NT = Np;
143 auto const& Nu = ip_data.N_u;
144 auto const& gradNu = ip_data.dNdx_u;
145 auto const& gradNp = ip_data.dNdx_p;
147 std::nullopt, this->element_.getID(), ip,
149 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
151 this->element_, Nu))};
152 auto const x_coord =
153 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
155 this->element_, Nu);
156
157 double const T = NT.dot(temperature);
158 double const T_prev = NT.dot(temperature_prev);
159 TemperatureData const T_data{T, T_prev};
160 double const pGR = Np.dot(gas_pressure);
161 double const pCap = Np.dot(capillary_pressure);
162 double const pLR = pGR - pCap;
163 GlobalDimVectorType const gradpGR = gradNp * gas_pressure;
164 GlobalDimVectorType const gradpCap = gradNp * capillary_pressure;
165 GlobalDimVectorType const gradT = gradNp * temperature;
166
168 MPL::VariableArray vars_prev;
169 vars.temperature = T;
170 vars.gas_phase_pressure = pGR;
171 vars.capillary_pressure = pCap;
172 vars.liquid_phase_pressure = pLR;
173
174 // medium properties
175 models.elastic_tangent_stiffness_model.eval({pos, t, dt}, T_data,
176 ip_cv.C_el_data);
177
178 models.biot_model.eval({pos, t, dt}, media_data, ip_cv.biot_data);
179
180 auto const Bu =
181 LinearBMatrix::computeBMatrix<DisplacementDim,
182 ShapeFunctionDisplacement::NPOINTS,
184 gradNu, Nu, x_coord, this->is_axially_symmetric_);
185
186 auto& eps = this->output_data_[ip].eps_data.eps;
187 eps.noalias() = Bu * displacement;
188
189 // Set volumetric strain for the general case without swelling.
191
192 models.S_L_model.eval({pos, t, dt}, media_data,
194 current_state.S_L_data, ip_cv.dS_L_dp_cap);
195
196 vars.liquid_saturation = current_state.S_L_data.S_L;
197 vars_prev.liquid_saturation = prev_state.S_L_data->S_L;
198
199 models.chi_S_L_model.eval({pos, t, dt}, media_data,
200 current_state.S_L_data, ip_cv.chi_S_L);
201
202 // solid phase compressibility
203 models.beta_p_SR_model.eval({pos, t, dt}, ip_cv.biot_data,
204 ip_cv.C_el_data, ip_cv.beta_p_SR);
205
206 // If there is swelling stress rate, compute swelling stress.
207 models.swelling_model.eval(
208 {pos, t, dt}, media_data, ip_cv.C_el_data, current_state.S_L_data,
209 prev_state.S_L_data, prev_state.swelling_data,
210 current_state.swelling_data, ip_cv.swelling_data);
211
212 // solid phase linear thermal expansion coefficient
213 models.s_therm_exp_model.eval({pos, t, dt}, media_data,
214 ip_cv.s_therm_exp_data);
215
216 models.mechanical_strain_model.eval(
217 T_data, ip_cv.s_therm_exp_data, this->output_data_[ip].eps_data,
218 Bu * displacement_prev, prev_state.mechanical_strain_data,
219 ip_cv.swelling_data, current_state.mechanical_strain_data);
220
221 models.s_mech_model.eval(
222 {pos, t, dt}, T_data, current_state.mechanical_strain_data,
223 prev_state.mechanical_strain_data, prev_state.eff_stress_data,
224 current_state.eff_stress_data, this->material_states_[ip],
225 ip_cd.s_mech_data, ip_cv.equivalent_plastic_strain_data);
226
227 models.total_stress_model.eval(
228 current_state.eff_stress_data, ip_cv.biot_data, ip_cv.chi_S_L,
230 ip_cv.total_stress_data);
231
232 // relative permeability
233 // Set mechanical variables for the intrinsic permeability model
234 // For stress dependent permeability.
235 vars.total_stress.emplace<SymmetricTensor>(
237 ip_cv.total_stress_data.sigma_total));
238
239 vars.equivalent_plastic_strain = *ip_cv.equivalent_plastic_strain_data;
240
243 eps);
244
245 // intrinsic permeability
246 ip_data.k_S = MPL::formEigenTensor<DisplacementDim>(
247 medium.property(MPL::PropertyType::permeability)
248 .value(vars, pos, t, dt));
249
250 ip_data.k_rel_G =
251 medium
252 .property(
253 MPL::PropertyType::relative_permeability_nonwetting_phase)
254 .template value<double>(vars, pos, t, dt);
255
256 auto const dk_rel_G_ds_L =
257 medium[MPL::PropertyType::relative_permeability_nonwetting_phase]
258 .template dValue<double>(vars, MPL::Variable::liquid_saturation,
259 pos, t, dt);
260
261 ip_data.k_rel_L =
262 medium.property(MPL::PropertyType::relative_permeability)
263 .template value<double>(vars, pos, t, dt);
264
265 auto const dk_rel_L_ds_L =
266 medium[MPL::PropertyType::relative_permeability]
267 .template dValue<double>(vars, MPL::Variable::liquid_saturation,
268 pos, t, dt);
269
272 current_state.mechanical_strain_data.eps_m);
273
274 auto const rho_ref_SR =
275 solid_phase.property(MPL::PropertyType::density)
276 .template value<double>(
277 vars, pos, t, std::numeric_limits<double>::quiet_NaN());
278
279 double const T0 = this->process_data_.reference_temperature(t, pos)[0];
280 double const delta_T(T - T0);
281 ip_data.thermal_volume_strain =
282 ip_cv.s_therm_exp_data.beta_T_SR * delta_T;
283
284 // initial porosity
285 auto const phi_0 = medium.property(MPL::PropertyType::porosity)
286 .template value<double>(vars, pos, t, dt);
287
288 auto const phi_S_0 = 1. - phi_0;
289
290#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
291 auto const& m = Invariants::identity2;
292 double const div_u = m.transpose() * eps;
293
294 const double phi_S = phi_S_0 * (1. + ip_data.thermal_volume_strain -
295 ip_cv.biot_data() * div_u);
296#else // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
297 const double phi_S = phi_S_0;
298#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
299
300 // porosity
301 ip_data.phi = 1. - phi_S;
302 vars.porosity = ip_data.phi;
303
304 // solid phase density
305#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
306 auto const rhoSR = rho_ref_SR * (1. - ip_data.thermal_volume_strain +
307 (ip_cv.biot_data() - 1.) * div_u);
308#else // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
309 auto const rhoSR = rho_ref_SR;
310#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
311
312 // constitutive model object as specified in process creation
313 auto& ptm = *this->process_data_.phase_transition_model_;
314 ptmv = ptm.updateConstitutiveVariables(ptmv, &medium, vars, pos, t, dt);
315 auto const& c = ptmv;
316
317 auto const phi_L = current_state.S_L_data.S_L * ip_data.phi;
318 auto const phi_G = (1. - current_state.S_L_data.S_L) * ip_data.phi;
319
320 // thermal conductivity
321 ip_data.lambda = MaterialPropertyLib::formEigenTensor<DisplacementDim>(
322 medium
323 .property(
325 .value(vars, pos, t, dt));
326
327 auto const cpS =
328 solid_phase.property(MPL::PropertyType::specific_heat_capacity)
329 .template value<double>(vars, pos, t, dt);
330 ip_data.h_S = cpS * T;
331 auto const u_S = ip_data.h_S;
332
333 ip_data.rho_u_eff = phi_G * c.rhoGR * c.uG + phi_L * c.rhoLR * c.uL +
334 phi_S * rhoSR * u_S;
335
336 ip_data.rho_G_h_G = phi_G * c.rhoGR * c.hG;
337 ip_data.rho_L_h_L = phi_L * c.rhoLR * c.hL;
338 ip_data.rho_S_h_S = phi_S * rhoSR * ip_data.h_S;
339
340 ip_data.muGR = c.muGR;
341 ip_data.muLR = c.muLR;
342
343 ip_data.rhoGR = c.rhoGR;
344 ip_data.rhoLR = c.rhoLR;
345 ip_data.rhoSR = rhoSR;
346
347 ip_data.rhoCGR = c.rhoCGR;
348 ip_data.rhoCLR = c.rhoCLR;
349 ip_data.rhoWGR = c.rhoWGR;
350 ip_data.rhoWLR = c.rhoWLR;
351
352 ip_data.dxmWG_dpGR = c.dxmWG_dpGR;
353 ip_data.dxmWG_dpCap = c.dxmWG_dpCap;
354 ip_data.dxmWG_dT = c.dxmWG_dT;
355
356 ip_data.dxmWL_dpGR = c.dxmWL_dpGR;
357 ip_data.dxmWL_dpCap = c.dxmWL_dpCap;
358 ip_data.dxmWL_dT = c.dxmWL_dT;
359
360 ip_data.dxmWL_dpLR = c.dxmWL_dpLR;
361
362 // for variable output
363 ip_data.xnCG = 1. - c.xnWG;
364 ip_data.xmCG = 1. - c.xmWG;
365 ip_data.xmWG = c.xmWG;
366 ip_data.xmWL = c.xmWL;
367 auto const xmCL = 1. - c.xmWL;
368
369 ip_data.diffusion_coefficient_vapour = c.diffusion_coefficient_vapour;
370 ip_data.diffusion_coefficient_solute = c.diffusion_coefficient_solute;
371
372 ip_data.h_G = c.hG;
373 ip_data.h_CG = c.hCG;
374 ip_data.h_WG = c.hWG;
375 ip_data.h_L = c.hL;
376 ip_data.pWGR = c.pWGR;
377
378 const GlobalDimVectorType gradxmWG = ip_data.dxmWG_dpGR * gradpGR +
379 ip_data.dxmWG_dpCap * gradpCap +
380 ip_data.dxmWG_dT * gradT;
381 const GlobalDimVectorType gradxmCG = -gradxmWG;
382
383 const GlobalDimVectorType gradxmWL = ip_data.dxmWL_dpGR * gradpGR +
384 ip_data.dxmWL_dpCap * gradpCap +
385 ip_data.dxmWL_dT * gradT;
386 const GlobalDimVectorType gradxmCL = -gradxmWL;
387
388 // Todo: factor -phiAlpha / xmZetaAlpha * DZetaAlpha can be evaluated in
389 // the respective phase transition model, here only the multiplication
390 // with the gradient of the mass fractions should take place.
391
392 ip_data.d_CG = ip_data.xmCG == 0.
393 ? 0. * gradxmCG // Keep d_CG's dimension and prevent
394 // division by zero
395 : -phi_G / ip_data.xmCG *
396 ip_data.diffusion_coefficient_vapour *
397 gradxmCG;
398
399 ip_data.d_WG = ip_data.xmWG == 0.
400 ? 0. * gradxmWG // Keep d_WG's dimension and prevent
401 // division by zero
402 : -phi_G / ip_data.xmWG *
403 ip_data.diffusion_coefficient_vapour *
404 gradxmWG;
405
406 ip_data.d_CL = xmCL == 0. ? 0. * gradxmCL // Keep d_CL's dimension and
407 // prevent division by zero
408 : -phi_L / xmCL *
409 ip_data.diffusion_coefficient_solute *
410 gradxmCL;
411
412 ip_data.d_WL = ip_data.xmWL == 0.
413 ? 0. * gradxmWL // Keep d_WG's dimension and prevent
414 // division by zero
415 : -phi_L / ip_data.xmWL *
416 ip_data.diffusion_coefficient_solute *
417 gradxmWL;
418
419 // ---------------------------------------------------------------------
420 // Derivatives for Jacobian
421 // ---------------------------------------------------------------------
422 auto const drho_LR_dT =
423 liquid_phase.property(MPL::PropertyType::density)
424 .template dValue<double>(vars, MPL::Variable::temperature, pos,
425 t, dt);
426 auto const drho_SR_dT =
427 solid_phase.property(MPL::PropertyType::density)
428 .template dValue<double>(vars, MPL::Variable::temperature,
429 pos, t, dt)
430#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
431 * (1. - ip_data.thermal_volume_strain +
432 (ip_cv.biot_data() - 1.) * div_u) -
433 rho_ref_SR * ip_cv.s_therm_exp_data.beta_T_SR
434#endif
435 ;
436
437 // porosity
438 auto const dphi_0_dT =
439 medium[MPL::PropertyType::porosity].template dValue<double>(
440 vars, MPL::Variable::temperature, pos, t, dt);
441
442 auto const dphi_S_0_dT = -dphi_0_dT;
443 const double dphi_S_dT = dphi_S_0_dT
444#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
445 * (1. + ip_data.thermal_volume_strain -
446 ip_cv.biot_data() * div_u) +
447 phi_S_0 * ip_cv.s_therm_exp_data.beta_T_SR
448#endif
449 ;
450
451 ip_cv.drho_u_eff_dT =
452 phi_G * c.drho_GR_dT * c.uG + phi_G * c.rhoGR * c.du_G_dT +
453 phi_L * drho_LR_dT * c.uL + phi_L * c.rhoLR * c.du_L_dT +
454 phi_S * drho_SR_dT * u_S + phi_S * rhoSR * cpS +
455 dphi_S_dT * rhoSR * u_S;
456
457 // ds_L_dp_GR = 0;
458 // ds_G_dp_GR = -ds_L_dp_GR;
459 double const ds_G_dp_cap = -ip_cv.dS_L_dp_cap();
460
461 // dphi_G_dp_GR = -ds_L_dp_GR * ip_data.phi = 0;
462 double const dphi_G_dp_cap = -ip_cv.dS_L_dp_cap() * ip_data.phi;
463 // dphi_L_dp_GR = ds_L_dp_GR * ip_data.phi = 0;
464 double const dphi_L_dp_cap = ip_cv.dS_L_dp_cap() * ip_data.phi;
465
466 auto const lambdaGR =
467 gas_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
468 ? MPL::formEigenTensor<DisplacementDim>(
469 gas_phase
470 .property(MPL::PropertyType::thermal_conductivity)
471 .value(vars, pos, t, dt))
472 : MPL::formEigenTensor<DisplacementDim>(0.);
473
474 auto const dlambda_GR_dT =
475 gas_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
476 ? MPL::formEigenTensor<DisplacementDim>(
477 gas_phase[MPL::PropertyType::thermal_conductivity].dValue(
478 vars, MPL::Variable::temperature, pos, t, dt))
479 : MPL::formEigenTensor<DisplacementDim>(0.);
480
481 auto const lambdaLR =
482 liquid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
483 ? MPL::formEigenTensor<DisplacementDim>(
484 liquid_phase
485 .property(MPL::PropertyType::thermal_conductivity)
486 .value(vars, pos, t, dt))
487 : MPL::formEigenTensor<DisplacementDim>(0.);
488
489 auto const dlambda_LR_dT =
490 liquid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
491 ? MPL::formEigenTensor<DisplacementDim>(
492 liquid_phase[MPL::PropertyType::thermal_conductivity]
493 .dValue(vars, MPL::Variable::temperature, pos, t, dt))
494 : MPL::formEigenTensor<DisplacementDim>(0.);
495
496 auto const lambdaSR =
497 solid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
498 ? MPL::formEigenTensor<DisplacementDim>(
499 solid_phase
500 .property(MPL::PropertyType::thermal_conductivity)
501 .value(vars, pos, t, dt))
502 : MPL::formEigenTensor<DisplacementDim>(0.);
503
504 auto const dlambda_SR_dT =
505 solid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
506 ? MPL::formEigenTensor<DisplacementDim>(
507 solid_phase[MPL::PropertyType::thermal_conductivity]
508 .dValue(vars, MPL::Variable::temperature, pos, t, dt))
509 : MPL::formEigenTensor<DisplacementDim>(0.);
510
511 ip_cv.dlambda_dp_cap =
512 dphi_G_dp_cap * lambdaGR + dphi_L_dp_cap * lambdaLR;
513
514 ip_cv.dlambda_dT = phi_G * dlambda_GR_dT + phi_L * dlambda_LR_dT +
515 phi_S * dlambda_SR_dT + dphi_S_dT * lambdaSR;
516
517 // From p_LR = p_GR - p_cap it follows for
518 // drho_LR/dp_GR = drho_LR/dp_LR * dp_LR/dp_GR
519 // = drho_LR/dp_LR * (dp_GR/dp_GR - dp_cap/dp_GR)
520 // = drho_LR/dp_LR * (1 - 0)
521 double const drho_LR_dp_GR = c.drho_LR_dp_LR;
522 double const drho_LR_dp_cap = -c.drho_LR_dp_LR;
523 // drho_GR_dp_cap = 0;
524
525 ip_cv.drho_h_eff_dp_GR =
526 /*(dphi_G_dp_GR = 0) * c.rhoGR * c.hG +*/ phi_G * c.drho_GR_dp_GR *
527 c.hG +
528 /*(dphi_L_dp_GR = 0) * c.rhoLR * c.hL +*/ phi_L * drho_LR_dp_GR *
529 c.hL;
530 ip_cv.drho_h_eff_dp_cap = dphi_G_dp_cap * c.rhoGR * c.hG +
531 /*phi_G * (drho_GR_dp_cap = 0) * c.hG +*/
532 dphi_L_dp_cap * c.rhoLR * c.hL +
533 phi_L * drho_LR_dp_cap * c.hL;
534
535 // TODO (naumov) Extend for temperature dependent porosities.
536 constexpr double dphi_G_dT = 0;
537 constexpr double dphi_L_dT = 0;
538 ip_cv.drho_h_eff_dT =
539 dphi_G_dT * c.rhoGR * c.hG + phi_G * c.drho_GR_dT * c.hG +
540 phi_G * c.rhoGR * c.dh_G_dT + dphi_L_dT * c.rhoLR * c.hL +
541 phi_L * drho_LR_dT * c.hL + phi_L * c.rhoLR * c.dh_L_dT +
542 dphi_S_dT * rhoSR * ip_data.h_S + phi_S * drho_SR_dT * ip_data.h_S +
543 phi_S * rhoSR * cpS;
544
545 ip_cv.drho_u_eff_dp_GR =
546 /*(dphi_G_dp_GR = 0) * c.rhoGR * c.uG +*/
547 phi_G * c.drho_GR_dp_GR * c.uG + phi_G * c.rhoGR * c.du_G_dp_GR +
548 /*(dphi_L_dp_GR = 0) * c.rhoLR * c.uL +*/
549 phi_L * drho_LR_dp_GR * c.uL + phi_L * c.rhoLR * c.du_L_dp_GR;
550
551 ip_cv.drho_u_eff_dp_cap = dphi_G_dp_cap * c.rhoGR * c.uG +
552 /*phi_G * (drho_GR_dp_cap = 0) * c.uG +*/
553 dphi_L_dp_cap * c.rhoLR * c.uL +
554 phi_L * drho_LR_dp_cap * c.uL +
555 phi_L * c.rhoLR * c.du_L_dp_cap;
556
557 auto const& b = this->process_data_.specific_body_force;
558 GlobalDimMatrixType const k_over_mu_G =
559 ip_data.k_S * ip_data.k_rel_G / ip_data.muGR;
560 GlobalDimMatrixType const k_over_mu_L =
561 ip_data.k_S * ip_data.k_rel_L / ip_data.muLR;
562
563 // dk_over_mu_G_dp_GR =
564 // ip_data.k_S * dk_rel_G_ds_L * (ds_L_dp_GR = 0) / ip_data.muGR =
565 // 0;
566 // dk_over_mu_L_dp_GR =
567 // ip_data.k_S * dk_rel_L_ds_L * (ds_L_dp_GR = 0) / ip_data.muLR =
568 // 0;
569 ip_cv.dk_over_mu_G_dp_cap =
570 ip_data.k_S * dk_rel_G_ds_L * ip_cv.dS_L_dp_cap() / ip_data.muGR;
571 ip_cv.dk_over_mu_L_dp_cap =
572 ip_data.k_S * dk_rel_L_ds_L * ip_cv.dS_L_dp_cap() / ip_data.muLR;
573
574 ip_data.w_GS = k_over_mu_G * c.rhoGR * b - k_over_mu_G * gradpGR;
575 ip_data.w_LS = k_over_mu_L * gradpCap + k_over_mu_L * c.rhoLR * b -
576 k_over_mu_L * gradpGR;
577
578 ip_cv.drho_GR_h_w_eff_dp_GR_Npart =
579 c.drho_GR_dp_GR * c.hG * ip_data.w_GS +
580 c.rhoGR * c.hG * k_over_mu_G * c.drho_GR_dp_GR * b;
581 ip_cv.drho_GR_h_w_eff_dp_GR_gradNpart =
582 -c.rhoGR * c.hG * k_over_mu_G - c.rhoLR * c.hL * k_over_mu_L;
583
584 ip_cv.drho_LR_h_w_eff_dp_cap_Npart =
585 -drho_LR_dp_cap * c.hL * ip_data.w_LS -
586 c.rhoLR * c.hL * k_over_mu_L * drho_LR_dp_cap * b;
587 ip_cv.drho_LR_h_w_eff_dp_cap_gradNpart =
588 // TODO (naumov) why the minus sign??????
589 -c.rhoLR * c.hL * k_over_mu_L;
590
591 ip_cv.drho_GR_h_w_eff_dT = c.drho_GR_dT * c.hG * ip_data.w_GS +
592 c.rhoGR * c.dh_G_dT * ip_data.w_GS +
593 drho_LR_dT * c.hL * ip_data.w_LS +
594 c.rhoLR * c.dh_L_dT * ip_data.w_LS;
595 // TODO (naumov) + k_over_mu_G * drho_GR_dT * b + k_over_mu_L *
596 // drho_LR_dT * b
597
598 // Derivatives of s_G * rho_C_GR_dot + s_L * rho_C_LR_dot abbreviated
599 // here with S_rho_C_eff.
600 double const s_L = current_state.S_L_data.S_L;
601 double const s_G = 1. - s_L;
602 double const rho_C_FR = s_G * ip_data.rhoCGR + s_L * ip_data.rhoCLR;
603 double const rho_W_FR = s_G * ip_data.rhoWGR + s_L * ip_data.rhoWLR;
604 // TODO (naumov) Extend for partially saturated media.
605 constexpr double drho_C_GR_dp_cap = 0;
606 if (dt == 0.)
607 {
608 ip_cv.dfC_3a_dp_GR = 0.;
609 ip_cv.dfC_3a_dp_cap = 0.;
610 ip_cv.dfC_3a_dT = 0.;
611 }
612 else
613 {
614 double const rho_C_GR_dot =
615 (ip_data.rhoCGR - ip_data.rhoCGR_prev) / dt;
616 double const rho_C_LR_dot =
617 (ip_data.rhoCLR - ip_data.rhoCLR_prev) / dt;
618 ip_cv.dfC_3a_dp_GR =
619 /*(ds_G_dp_GR = 0) * rho_C_GR_dot +*/ s_G * c.drho_C_GR_dp_GR /
620 dt +
621 /*(ds_L_dp_GR = 0) * rho_C_LR_dot +*/ s_L * c.drho_C_LR_dp_GR /
622 dt;
623 ip_cv.dfC_3a_dp_cap = ds_G_dp_cap * rho_C_GR_dot +
624 s_G * drho_C_GR_dp_cap / dt +
625 ip_cv.dS_L_dp_cap() * rho_C_LR_dot -
626 s_L * c.drho_C_LR_dp_LR / dt;
627 ip_cv.dfC_3a_dT =
628 s_G * c.drho_C_GR_dT / dt + s_L * c.drho_C_LR_dT / dt;
629 }
630
631 double const drho_C_FR_dp_GR =
632 /*(ds_G_dp_GR = 0) * ip_data.rhoCGR +*/ s_G * c.drho_C_GR_dp_GR +
633 /*(ds_L_dp_GR = 0) * ip_data.rhoCLR +*/ s_L * c.drho_C_LR_dp_GR;
634 ip_cv.dfC_4_MCpG_dp_GR = drho_C_FR_dp_GR *
635 (ip_cv.biot_data() - ip_data.phi) *
636 ip_cv.beta_p_SR();
637
638 double const drho_C_FR_dT = s_G * c.drho_C_GR_dT + s_L * c.drho_C_LR_dT;
639 ip_cv.dfC_4_MCpG_dT =
640 drho_C_FR_dT * (ip_cv.biot_data() - ip_data.phi) * ip_cv.beta_p_SR()
641#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
642 - rho_C_FR * ip_data.dphi_dT * ip_cv.beta_p_SR()
643#endif
644 ;
645
646 ip_cv.dfC_4_MCT_dT = drho_C_FR_dT * (ip_cv.biot_data() - ip_data.phi) *
647 ip_cv.s_therm_exp_data.beta_T_SR
648#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
649 +
650 rho_C_FR * (ip_cv.biot_data() - ip_data.dphi_dT) *
651 ip_cv.s_therm_exp_data.beta_T_SR
652#endif
653 ;
654
655 ip_cv.dfC_4_MCu_dT = drho_C_FR_dT * ip_cv.biot_data();
656
657 ip_cv.dfC_2a_dp_GR = -ip_data.phi * c.drho_C_GR_dp_GR -
658 drho_C_FR_dp_GR * pCap *
659 (ip_cv.biot_data() - ip_data.phi) *
660 ip_cv.beta_p_SR();
661
662 double const drho_C_FR_dp_cap =
663 ds_G_dp_cap * ip_data.rhoCGR + s_G * drho_C_GR_dp_cap +
664 ip_cv.dS_L_dp_cap() * ip_data.rhoCLR - s_L * c.drho_C_LR_dp_LR;
665
666 ip_cv.dfC_2a_dp_cap =
667 ip_data.phi * (-c.drho_C_LR_dp_LR - drho_C_GR_dp_cap) -
668 drho_C_FR_dp_cap * pCap * (ip_cv.biot_data() - ip_data.phi) *
669 ip_cv.beta_p_SR() +
670 rho_C_FR * (ip_cv.biot_data() - ip_data.phi) * ip_cv.beta_p_SR();
671
672 ip_cv.dfC_2a_dT =
673#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
674 ip_data.dphi_dT * (ip_data.rhoCLR - ip_data.rhoCGR) +
675#endif
676 ip_data.phi * (c.drho_C_LR_dT - c.drho_C_GR_dT) -
677 drho_C_FR_dT * pCap * (ip_cv.biot_data() - ip_data.phi) *
678 ip_cv.beta_p_SR()
679#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
680 + rho_C_FR * pCap * ip_data.dphi_dT * ip_cv.beta_p_SR()
681#endif
682 ;
683
684 ip_cv.dadvection_C_dp_GR = c.drho_C_GR_dp_GR * k_over_mu_G
685 // + rhoCGR * (dk_over_mu_G_dp_GR = 0)
686 // + rhoCLR * (dk_over_mu_L_dp_GR = 0)
687 + c.drho_C_LR_dp_GR * k_over_mu_L;
688
689 ip_cv.dadvection_C_dp_cap =
690 //(drho_C_GR_dp_cap = 0) * k_over_mu_G
691 ip_data.rhoCGR * ip_cv.dk_over_mu_G_dp_cap +
692 (-c.drho_C_LR_dp_LR) * k_over_mu_L +
693 ip_data.rhoCLR * ip_cv.dk_over_mu_L_dp_cap;
694
695 ip_cv.dfC_4_LCpG_dT =
696 c.drho_C_GR_dT * k_over_mu_G + c.drho_C_LR_dT * k_over_mu_L
697 // + ip_cv.ddiffusion_C_p_dT TODO (naumov)
698 ;
699
700 double const drho_W_FR_dp_GR =
701 /*(ds_G_dp_GR = 0) * ip_data.rhoWGR +*/ s_G * c.drho_W_GR_dp_GR +
702 /*(ds_L_dp_GR = 0) * ip_data.rhoWLR +*/ s_L * c.drho_W_LR_dp_GR;
703 double const drho_W_FR_dp_cap =
704 ds_G_dp_cap * ip_data.rhoWGR + s_G * c.drho_W_GR_dp_cap +
705 ip_cv.dS_L_dp_cap() * ip_data.rhoWLR - s_L * c.drho_W_LR_dp_LR;
706 double const drho_W_FR_dT = s_G * c.drho_W_GR_dT + s_L * c.drho_W_LR_dT;
707
708 ip_cv.dfW_2a_dp_GR =
709 ip_data.phi * (c.drho_W_LR_dp_GR - c.drho_W_GR_dp_GR);
710 ip_cv.dfW_2b_dp_GR = drho_W_FR_dp_GR * pCap *
711 (ip_cv.biot_data() - ip_data.phi) *
712 ip_cv.beta_p_SR();
713 ip_cv.dfW_2a_dp_cap =
714 ip_data.phi * (-c.drho_W_LR_dp_LR - c.drho_W_GR_dp_cap);
715 ip_cv.dfW_2b_dp_cap =
716 drho_W_FR_dp_cap * pCap * (ip_cv.biot_data() - ip_data.phi) *
717 ip_cv.beta_p_SR() +
718 rho_W_FR * (ip_cv.biot_data() - ip_data.phi) * ip_cv.beta_p_SR();
719
720 ip_cv.dfW_2a_dT =
721#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
722 ip_data.dphi_dT * (ip_data.rhoWLR - ip_data.rhoWGR) +
723#endif
724 ip_data.phi * (c.drho_W_LR_dT - c.drho_W_GR_dT);
725 ip_cv.dfW_2b_dT =
726 drho_W_FR_dT * pCap * (ip_cv.biot_data() - ip_data.phi) *
727 ip_cv.beta_p_SR()
728#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
729 - rho_W_FR * pCap * ip_data.dphi_dT * ip_cv.beta_p_SR()
730#endif
731 ;
732
733 if (dt == 0.)
734 {
735 ip_cv.dfW_3a_dp_GR = 0.;
736 ip_cv.dfW_3a_dp_cap = 0.;
737 ip_cv.dfW_3a_dT = 0.;
738 }
739 else
740 {
741 double const rho_W_GR_dot =
742 (ip_data.rhoWGR - ip_data.rhoWGR_prev) / dt;
743 double const rho_W_LR_dot =
744 (ip_data.rhoWLR - ip_data.rhoWLR_prev) / dt;
745
746 ip_cv.dfW_3a_dp_GR =
747 /*(ds_G_dp_GR = 0) * rho_W_GR_dot +*/ s_G * c.drho_W_GR_dp_GR /
748 dt +
749 /*(ds_L_dp_GR = 0) * rho_W_LR_dot +*/ s_L * c.drho_W_LR_dp_GR /
750 dt;
751 ip_cv.dfW_3a_dp_cap = ds_G_dp_cap * rho_W_GR_dot +
752 s_G * c.drho_W_GR_dp_cap / dt +
753 ip_cv.dS_L_dp_cap() * rho_W_LR_dot -
754 s_L * c.drho_W_LR_dp_LR / dt;
755 ip_cv.dfW_3a_dT =
756 s_G * c.drho_W_GR_dT / dt + s_L * c.drho_W_LR_dT / dt;
757 }
758
759 ip_cv.dfW_4_LWpG_a_dp_GR = c.drho_W_GR_dp_GR * k_over_mu_G
760 // + rhoWGR * (dk_over_mu_G_dp_GR = 0)
761 + c.drho_W_LR_dp_GR * k_over_mu_L
762 // + rhoWLR * (dk_over_mu_L_dp_GR = 0)
763 ;
764 ip_cv.dfW_4_LWpG_a_dp_cap = c.drho_W_GR_dp_cap * k_over_mu_G +
765 ip_data.rhoWGR * ip_cv.dk_over_mu_G_dp_cap +
766 -c.drho_W_LR_dp_LR * k_over_mu_L +
767 ip_data.rhoWLR * ip_cv.dk_over_mu_L_dp_cap;
768
769 ip_cv.dfW_4_LWpG_a_dT =
770 c.drho_W_GR_dT * k_over_mu_G
771 //+ rhoWGR * (dk_over_mu_G_dT != 0 TODO for mu_G(T))
772 + c.drho_W_LR_dT * k_over_mu_L
773 //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_G(T))
774 ;
775
776 // TODO (naumov) for dxmW*/d* != 0
777 ip_cv.dfW_4_LWpG_d_dp_GR =
778 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
779 ip_cv.dfW_4_LWpG_d_dp_cap =
780 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
781 ip_cv.dfW_4_LWpG_d_dT =
782 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
783
784 ip_cv.dfW_4_LWpC_a_dp_GR = c.drho_W_LR_dp_GR * k_over_mu_L
785 //+ rhoWLR * (dk_over_mu_L_dp_GR = 0)
786 ;
787 ip_cv.dfW_4_LWpC_a_dp_cap = -c.drho_W_LR_dp_LR * k_over_mu_L +
788 ip_data.rhoWLR * ip_cv.dk_over_mu_L_dp_cap;
789 ip_cv.dfW_4_LWpC_a_dT = c.drho_W_LR_dT * k_over_mu_L
790 //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_L(T))
791 ;
792
793 // TODO (naumov) for dxmW*/d* != 0
794 ip_cv.dfW_4_LWpC_d_dp_GR =
795 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
796 ip_cv.dfW_4_LWpC_d_dp_cap =
797 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
798 ip_cv.dfW_4_LWpC_d_dT =
799 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
800
801 ip_cv.dfC_4_LCpC_a_dp_GR = c.drho_C_LR_dp_GR * k_over_mu_L
802 //+ rhoCLR * (dk_over_mu_L_dp_GR = 0)
803 ;
804 ip_cv.dfC_4_LCpC_a_dp_cap = -c.drho_C_LR_dp_LR * k_over_mu_L +
805 ip_data.rhoCLR * ip_cv.dk_over_mu_L_dp_cap;
806 ip_cv.dfC_4_LCpC_a_dT = c.drho_W_LR_dT * k_over_mu_L
807 //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_L(T))
808 ;
809
810 // TODO (naumov) for dxmW*/d* != 0
811 ip_cv.dfC_4_LCpC_d_dp_GR =
812 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
813 ip_cv.dfC_4_LCpC_d_dp_cap =
814 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
815 ip_cv.dfC_4_LCpC_d_dT =
816 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
817 }
818
819 return {ip_constitutive_data, ip_constitutive_variables};
820}
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > mechanical_strain
Definition: VariableType.h:182
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > total_stress
Definition: VariableType.h:201
Eigen::Matrix< double, KelvinVectorSize, 1 > SymmetricTensor
Definition: TH2MFEM.h:73
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:49
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
BaseLib::StrongType< double, struct CapillaryPressureTag > CapillaryPressureData
Definition: Base.h:52
BaseLib::StrongType< double, struct GasPressureTag > GasPressureData
Definition: Base.h:50

References ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::beta_p_SR_model, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::biot_model, MaterialPropertyLib::VariableArray::capillary_pressure, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::chi_S_L_model, ProcessLib::LinearBMatrix::computeBMatrix(), ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::elastic_tangent_stiffness_model, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, ProcessLib::TH2M::ConstitutiveRelations::BiotModel::eval(), ProcessLib::TH2M::ConstitutiveRelations::SaturationModel::eval(), ProcessLib::TH2M::ConstitutiveRelations::BishopsModel::eval(), MaterialPropertyLib::VariableArray::gas_phase_pressure, NumLib::interpolateCoordinates(), NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, MaterialPropertyLib::VariableArray::mechanical_strain, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::mechanical_strain_model, MaterialPropertyLib::VariableArray::porosity, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::S_L_model, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::s_mech_model, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::s_therm_exp_model, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::swelling_model, MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::thermal_conductivity, MaterialPropertyLib::VariableArray::total_stress, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::total_stress_model, and MaterialPropertyLib::VariableArray::volumetric_strain.

Member Data Documentation

◆ _ip_data

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector<IpData, Eigen::aligned_allocator<IpData> > ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data
private

Definition at line 382 of file TH2MFEM.h.

Referenced by ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::TH2MLocalAssembler(), ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtEnthalpyGas(), ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtEnthalpyLiquid(), ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtEnthalpySolid(), ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtGasDensity(), ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtIntrinsicPermeability(), ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtLiquidDensity(), ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtMassFractionGas(), ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtMassFractionLiquid(), ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtMoleFractionGas(), ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtPorosity(), ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtRelativePermeabilityGas(), ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtRelativePermeabilityLiquid(), ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtSolidDensity(), ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtVapourPressure(), ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::initializeConcrete(), and ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::postTimestepConcrete().

◆ _secondary_data

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType> ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_secondary_data
private

◆ C_index

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
const int ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::C_index = 0
staticprivate

Definition at line 400 of file TH2MFEM.h.

◆ C_size

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
const int ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::C_size = ShapeFunctionPressure::NPOINTS
staticprivate

Definition at line 401 of file TH2MFEM.h.

◆ capillary_pressure_index

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
const int ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::capillary_pressure_index = ShapeFunctionPressure::NPOINTS
staticprivate

Definition at line 392 of file TH2MFEM.h.

◆ capillary_pressure_size

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
const int ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::capillary_pressure_size = ShapeFunctionPressure::NPOINTS
staticprivate

Definition at line 393 of file TH2MFEM.h.

◆ displacement_index

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
const int ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::displacement_index = ShapeFunctionPressure::NPOINTS * 3
staticprivate

Definition at line 396 of file TH2MFEM.h.

◆ displacement_size

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
const int ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::displacement_size
staticprivate
Initial value:
=
ShapeFunctionDisplacement::NPOINTS * DisplacementDim

Definition at line 397 of file TH2MFEM.h.

◆ gas_pressure_index

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
const int ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::gas_pressure_index = 0
staticprivate

Definition at line 390 of file TH2MFEM.h.

◆ gas_pressure_size

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
const int ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::gas_pressure_size = ShapeFunctionPressure::NPOINTS
staticprivate

Definition at line 391 of file TH2MFEM.h.

◆ KelvinVectorSize

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
int const ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::KelvinVectorSize
static
Initial value:
=
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:24

Definition at line 71 of file TH2MFEM.h.

◆ N_u_op

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
constexpr auto& ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::N_u_op
staticconstexpr
Initial value:
DisplacementDim,
constexpr Eigen::CwiseNullaryOp< EigenBlockMatrixViewFunctor< D, M >, typename EigenBlockMatrixViewFunctor< D, M >::Matrix > eigenBlockMatrixView(const Eigen::MatrixBase< M > &matrix)
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType

Definition at line 75 of file TH2MFEM.h.

◆ temperature_index

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
const int ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::temperature_index = 2 * ShapeFunctionPressure::NPOINTS
staticprivate

Definition at line 394 of file TH2MFEM.h.

◆ temperature_size

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
const int ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::temperature_size = ShapeFunctionPressure::NPOINTS
staticprivate

Definition at line 395 of file TH2MFEM.h.

◆ W_index

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
const int ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::W_index = ShapeFunctionPressure::NPOINTS
staticprivate

Definition at line 402 of file TH2MFEM.h.

◆ W_size

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
const int ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::W_size = ShapeFunctionPressure::NPOINTS
staticprivate

Definition at line 403 of file TH2MFEM.h.


The documentation for this class was generated from the following files: