OGS
List of all members | Public Types | Public Member Functions | Static Public Attributes | Private Types | Private Member Functions | Private Attributes | Static Private Attributes
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
 
using ShapeMatricesTypePressure
 
template<int N>
using VectorType
 
template<int M, int N>
using MatrixType
 
using GlobalDimMatrixType
 
using 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)
 
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
 
virtual void setInitialConditions (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, 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 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 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, int const process_id)
 
void postNonLinearSolver (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, 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 ~ExtrapolatableElement ()=default
 

Static Public Attributes

static int const KelvinVectorSize
 
static constexpr auto & N_u_op
 

Private Types

using BMatricesType
 
using IpData
 

Private Member Functions

std::size_t setIPDataInitialConditions (std::string_view const name, double const *values, int const integration_order) override
 
void setInitialConditionsConcrete (Eigen::VectorXd const local_x, double const t, int const process_id) 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, 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 & 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
private
Initial value:
BMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>

Definition at line 249 of file TH2MFEM.h.

◆ GlobalDimMatrixType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::GlobalDimMatrixType
Initial value:
typename ShapeMatricesTypePressure::GlobalDimMatrixType
EigenFixedShapeMatrixPolicy::GlobalDimMatrixType
MatrixType< GlobalDim, GlobalDim > GlobalDimMatrixType
Definition ShapeMatrixPolicy.h:84

Definition at line 65 of file TH2MFEM.h.

◆ GlobalDimVectorType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::GlobalDimVectorType
Initial value:
typename ShapeMatricesTypePressure::GlobalDimVectorType
EigenFixedShapeMatrixPolicy::GlobalDimVectorType
VectorType< GlobalDim > GlobalDimVectorType
Definition ShapeMatrixPolicy.h:85

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
private
Initial value:
IntegrationPointData<BMatricesType, ShapeMatricesTypeDisplacement,
ShapeMatricesTypePressure, DisplacementDim,
ShapeFunctionDisplacement::NPOINTS>
ProcessLib::TH2M::TH2MLocalAssembler::ShapeMatricesTypeDisplacement
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
Definition TH2MFEM.h:51
ProcessLib::TH2M::TH2MLocalAssembler::BMatricesType
BMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > BMatricesType
Definition TH2MFEM.h:249
ProcessLib::TH2M::TH2MLocalAssembler::ShapeMatricesTypePressure
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
Definition TH2MFEM.h:54

Definition at line 251 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
Initial value:
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
Initial value:
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
Initial value:
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
Initial value:
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 =
47 this->integration_method_.getNumberOfPoints();
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,
54 ShapeMatricesTypeDisplacement,
55 DisplacementDim>(e, is_axially_symmetric,
56 this->integration_method_);
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 =
68 this->integration_method_.getWeightedPoint(ip).getWeight() *
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}
MathLib::WeightedPoint::getWeight
double getWeight() const
Definition WeightedPoint.h:80
NumLib::GenericIntegrationMethod::getWeightedPoint
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
Definition GenericIntegrationMethod.h:36
ProcessLib::TH2M::TH2MLocalAssembler::_secondary_data
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
Definition TH2MFEM.h:259
ProcessLib::TH2M::TH2MLocalAssembler::_ip_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
Definition TH2MFEM.h:255
NumLib::initShapeMatrices
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)
Definition InitShapeMatrices.h:53
ProcessLib::HydroMechanics::SecondaryData::N_u
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u
Definition HydroMechanicsFEM.h:161
ProcessLib::TH2M::LocalAssemblerInterface::integration_method_
NumLib::GenericIntegrationMethod const & integration_method_
Definition LocalAssemblerInterface.h:169

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 983 of file TH2MFEM-impl.h.

989{
990 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
991 temperature_size + displacement_size;
992 assert(local_x.size() == matrix_size);
993
994 auto const capillary_pressure =
995 Eigen::Map<VectorType<capillary_pressure_size> const>(
996 local_x.data() + capillary_pressure_index, capillary_pressure_size);
997
998 auto const capillary_pressure_prev =
999 Eigen::Map<VectorType<capillary_pressure_size> const>(
1000 local_x_prev.data() + capillary_pressure_index,
1001 capillary_pressure_size);
1002
1003 // pointer to local_M_data vector
1004 auto local_M =
1005 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1006 local_M_data, matrix_size, matrix_size);
1007
1008 // pointer to local_K_data vector
1009 auto local_K =
1010 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1011 local_K_data, matrix_size, matrix_size);
1012
1013 // pointer to local_rhs_data vector
1014 auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
1015 local_rhs_data, matrix_size);
1016
1017 // component-formulation
1018 // W - liquid phase main component
1019 // C - gas phase main component
1020 // pointer-matrices to the mass matrix - C component equation
1021 auto MCpG = local_M.template block<C_size, gas_pressure_size>(
1022 C_index, gas_pressure_index);
1023 auto MCpC = local_M.template block<C_size, capillary_pressure_size>(
1024 C_index, capillary_pressure_index);
1025 auto MCT = local_M.template block<C_size, temperature_size>(
1026 C_index, temperature_index);
1027 auto MCu = local_M.template block<C_size, displacement_size>(
1028 C_index, displacement_index);
1029
1030 // pointer-matrices to the stiffness matrix - C component equation
1031 auto LCpG = local_K.template block<C_size, gas_pressure_size>(
1032 C_index, gas_pressure_index);
1033 auto LCpC = local_K.template block<C_size, capillary_pressure_size>(
1034 C_index, capillary_pressure_index);
1035 auto LCT = local_K.template block<C_size, temperature_size>(
1036 C_index, temperature_index);
1037
1038 // pointer-matrices to the mass matrix - W component equation
1039 auto MWpG = local_M.template block<W_size, gas_pressure_size>(
1040 W_index, gas_pressure_index);
1041 auto MWpC = local_M.template block<W_size, capillary_pressure_size>(
1042 W_index, capillary_pressure_index);
1043 auto MWT = local_M.template block<W_size, temperature_size>(
1044 W_index, temperature_index);
1045 auto MWu = local_M.template block<W_size, displacement_size>(
1046 W_index, displacement_index);
1047
1048 // pointer-matrices to the stiffness matrix - W component equation
1049 auto LWpG = local_K.template block<W_size, gas_pressure_size>(
1050 W_index, gas_pressure_index);
1051 auto LWpC = local_K.template block<W_size, capillary_pressure_size>(
1052 W_index, capillary_pressure_index);
1053 auto LWT = local_K.template block<W_size, temperature_size>(
1054 W_index, temperature_index);
1055
1056 // pointer-matrices to the mass matrix - temperature equation
1057 auto MTu = local_M.template block<temperature_size, displacement_size>(
1058 temperature_index, displacement_index);
1059
1060 // pointer-matrices to the stiffness matrix - temperature equation
1061 auto KTT = local_K.template block<temperature_size, temperature_size>(
1062 temperature_index, temperature_index);
1063
1064 // pointer-matrices to the stiffness matrix - displacement equation
1065 auto KUpG = local_K.template block<displacement_size, gas_pressure_size>(
1066 displacement_index, gas_pressure_index);
1067 auto KUpC =
1068 local_K.template block<displacement_size, capillary_pressure_size>(
1069 displacement_index, capillary_pressure_index);
1070
1071 // pointer-vectors to the right hand side terms - C-component equation
1072 auto fC = local_f.template segment<C_size>(C_index);
1073 // pointer-vectors to the right hand side terms - W-component equation
1074 auto fW = local_f.template segment<W_size>(W_index);
1075 // pointer-vectors to the right hand side terms - temperature equation
1076 auto fT = local_f.template segment<temperature_size>(temperature_index);
1077 // pointer-vectors to the right hand side terms - displacement equation
1078 auto fU = local_f.template segment<displacement_size>(displacement_index);
1079
1080 unsigned const n_integration_points =
1081 this->integration_method_.getNumberOfPoints();
1082
1083 ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const models{
1084 this->solid_material_, *this->process_data_.phase_transition_model_};
1085
1086 auto const [ip_constitutive_data, ip_constitutive_variables] =
1087 updateConstitutiveVariables(
1088 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1089 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1090 local_x_prev.size()),
1091 t, dt, models);
1092
1093 for (unsigned int_point = 0; int_point < n_integration_points; int_point++)
1094 {
1095 auto& ip = _ip_data[int_point];
1096 auto& ip_cv = ip_constitutive_variables[int_point];
1097 auto& ip_out = this->output_data_[int_point];
1098 auto& current_state = this->current_states_[int_point];
1099 auto const& prev_state = this->prev_states_[int_point];
1100
1101 auto const& Np = ip.N_p;
1102 auto const& NT = Np;
1103 auto const& Nu = ip.N_u;
1104 ParameterLib::SpatialPosition const pos{
1105 std::nullopt, this->element_.getID(), int_point,
1106 MathLib::Point3d(
1107 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
1108 ShapeMatricesTypeDisplacement>(
1109 this->element_, Nu))};
1110
1111 auto const& NpT = Np.transpose().eval();
1112 auto const& NTT = NT.transpose().eval();
1113
1114 auto const& gradNp = ip.dNdx_p;
1115 auto const& gradNT = gradNp;
1116 auto const& gradNu = ip.dNdx_u;
1117
1118 auto const& gradNpT = gradNp.transpose().eval();
1119 auto const& gradNTT = gradNT.transpose().eval();
1120
1121 auto const& w = ip.integration_weight;
1122
1123 auto const& m = Invariants::identity2;
1124
1125 auto const mT = m.transpose().eval();
1126
1127 auto const x_coord =
1128 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
1129 ShapeMatricesTypeDisplacement>(
1130 this->element_, Nu);
1131
1132 auto const Bu =
1133 LinearBMatrix::computeBMatrix<DisplacementDim,
1134 ShapeFunctionDisplacement::NPOINTS,
1135 typename BMatricesType::BMatrixType>(
1136 gradNu, Nu, x_coord, this->is_axially_symmetric_);
1137
1138 auto const BuT = Bu.transpose().eval();
1139
1140 double const pCap = Np.dot(capillary_pressure);
1141 double const pCap_prev = Np.dot(capillary_pressure_prev);
1142
1143 double const beta_T_SR = ip_cv.s_therm_exp_data.beta_T_SR;
1144
1145 auto const I =
1146 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
1147
1148 double const sD_G =
1149 ip_cv.phase_transition_data.diffusion_coefficient_vapour;
1150 double const sD_L =
1151 ip_cv.phase_transition_data.diffusion_coefficient_solute;
1152
1153 GlobalDimMatrixType const D_C_G = sD_G * I;
1154 GlobalDimMatrixType const D_W_G = sD_G * I;
1155 GlobalDimMatrixType const D_C_L = sD_L * I;
1156 GlobalDimMatrixType const D_W_L = sD_L * I;
1157
1158 auto const s_L = current_state.S_L_data.S_L;
1159 auto const s_G = 1. - s_L;
1160 auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
1161
1162 auto& alpha_B = ip_cv.biot_data();
1163 auto& beta_p_SR = ip_cv.beta_p_SR();
1164
1165 auto const& b = this->process_data_.specific_body_force;
1166
1167 // volume fraction
1168 auto const phi_G = s_G * ip_out.porosity_data.phi;
1169 auto const phi_L = s_L * ip_out.porosity_data.phi;
1170 auto const phi_S = 1. - ip_out.porosity_data.phi;
1171
1172 auto const rhoGR = ip_out.fluid_density_data.rho_GR;
1173 auto const rhoLR = ip_out.fluid_density_data.rho_LR;
1174
1175 // effective density
1176 auto const rho = phi_G * rhoGR + phi_L * rhoLR +
1177 phi_S * ip_out.solid_density_data.rho_SR;
1178
1179 // abbreviations
1180 const double rho_C_FR =
1181 s_G * current_state.constituent_density_data.rho_C_GR +
1182 s_L * current_state.constituent_density_data.rho_C_LR;
1183 const double rho_W_FR =
1184 s_G * current_state.constituent_density_data.rho_W_GR +
1185 s_L * current_state.rho_W_LR();
1186
1187 auto const rho_C_GR_dot =
1188 (current_state.constituent_density_data.rho_C_GR -
1189 prev_state.constituent_density_data->rho_C_GR) /
1190 dt;
1191 auto const rho_C_LR_dot =
1192 (current_state.constituent_density_data.rho_C_LR -
1193 prev_state.constituent_density_data->rho_C_LR) /
1194 dt;
1195 auto const rho_W_GR_dot =
1196 (current_state.constituent_density_data.rho_W_GR -
1197 prev_state.constituent_density_data->rho_W_GR) /
1198 dt;
1199 auto const rho_W_LR_dot =
1200 (current_state.rho_W_LR() - **prev_state.rho_W_LR) / dt;
1201
1202 auto const rho_h_eff = ip.rho_G_h_G + ip.rho_L_h_L + ip.rho_S_h_S;
1203
1204 auto const rho_u_eff_dot = (ip.rho_u_eff - ip.rho_u_eff_prev) / dt;
1205
1206 GlobalDimMatrixType const k_over_mu_G =
1207 ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_G /
1208 ip_cv.viscosity_data.mu_GR;
1209 GlobalDimMatrixType const k_over_mu_L =
1210 ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_L /
1211 ip_cv.viscosity_data.mu_LR;
1212
1213 // ---------------------------------------------------------------------
1214 // C-component equation
1215 // ---------------------------------------------------------------------
1216
1217 MCpG.noalias() += NpT * rho_C_FR *
1218 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1219 Np * w;
1220 MCpC.noalias() -= NpT * rho_C_FR *
1221 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1222 s_L * Np * w;
1223
1224 if (this->process_data_.apply_mass_lumping)
1225 {
1226 if (pCap - pCap_prev != 0.) // avoid division by Zero
1227 {
1228 MCpC.noalias() +=
1229 NpT *
1230 (ip_out.porosity_data.phi *
1231 (current_state.constituent_density_data.rho_C_LR -
1232 current_state.constituent_density_data.rho_C_GR) -
1233 rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1234 beta_p_SR) *
1235 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1236 }
1237 }
1238
1239 MCT.noalias() -= NpT * rho_C_FR * (alpha_B - ip_out.porosity_data.phi) *
1240 beta_T_SR * Np * w;
1241 MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
1242
1243 using DisplacementDimMatrix =
1244 Eigen::Matrix<double, DisplacementDim, DisplacementDim>;
1245
1246 DisplacementDimMatrix const advection_C_G =
1247 current_state.constituent_density_data.rho_C_GR * k_over_mu_G;
1248 DisplacementDimMatrix const advection_C_L =
1249 current_state.constituent_density_data.rho_C_LR * k_over_mu_L;
1250
1251 DisplacementDimMatrix const diffusion_CGpGR =
1252 -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dpGR;
1253 DisplacementDimMatrix const diffusion_CLpGR =
1254 -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dpGR;
1255
1256 DisplacementDimMatrix const diffusion_CGpCap =
1257 -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dpCap;
1258 DisplacementDimMatrix const diffusion_CLpCap =
1259 -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dpCap;
1260
1261 DisplacementDimMatrix const diffusion_CGT =
1262 -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dT;
1263 DisplacementDimMatrix const diffusion_CLT =
1264 -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dT;
1265
1266 DisplacementDimMatrix const advection_C = advection_C_G + advection_C_L;
1267 DisplacementDimMatrix const diffusion_C_pGR =
1268 diffusion_CGpGR + diffusion_CLpGR;
1269 DisplacementDimMatrix const diffusion_C_pCap =
1270 diffusion_CGpCap + diffusion_CLpCap;
1271
1272 DisplacementDimMatrix const diffusion_C_T =
1273 diffusion_CGT + diffusion_CLT;
1274
1275 LCpG.noalias() +=
1276 gradNpT * (advection_C + diffusion_C_pGR) * gradNp * w;
1277
1278 LCpC.noalias() +=
1279 gradNpT * (diffusion_C_pCap - advection_C_L) * gradNp * w;
1280
1281 LCT.noalias() += gradNpT * (diffusion_C_T)*gradNp * w;
1282
1283 fC.noalias() +=
1284 gradNpT * (advection_C_G * rhoGR + advection_C_L * rhoLR) * b * w;
1285
1286 if (!this->process_data_.apply_mass_lumping)
1287 {
1288 fC.noalias() -=
1289 NpT *
1290 (ip_out.porosity_data.phi *
1291 (current_state.constituent_density_data.rho_C_LR -
1292 current_state.constituent_density_data.rho_C_GR) -
1293 rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1294 beta_p_SR) *
1295 s_L_dot * w;
1296 }
1297 // fC_III
1298 fC.noalias() -= NpT * ip_out.porosity_data.phi *
1299 (s_G * rho_C_GR_dot + s_L * rho_C_LR_dot) * w;
1300
1301 // ---------------------------------------------------------------------
1302 // W-component equation
1303 // ---------------------------------------------------------------------
1304
1305 MWpG.noalias() += NpT * rho_W_FR *
1306 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1307 Np * w;
1308 MWpC.noalias() -= NpT * rho_W_FR *
1309 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1310 s_L * Np * w;
1311
1312 if (this->process_data_.apply_mass_lumping)
1313 {
1314 if (pCap - pCap_prev != 0.) // avoid division by Zero
1315 {
1316 MWpC.noalias() +=
1317 NpT *
1318 (ip_out.porosity_data.phi *
1319 (current_state.rho_W_LR() -
1320 current_state.constituent_density_data.rho_W_GR) -
1321 rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1322 beta_p_SR) *
1323 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1324 }
1325 }
1326
1327 MWT.noalias() -= NpT * rho_W_FR * (alpha_B - ip_out.porosity_data.phi) *
1328 beta_T_SR * Np * w;
1329
1330 MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
1331
1332 DisplacementDimMatrix const advection_W_G =
1333 current_state.constituent_density_data.rho_W_GR * k_over_mu_G;
1334 DisplacementDimMatrix const advection_W_L =
1335 current_state.rho_W_LR() * k_over_mu_L;
1336
1337 DisplacementDimMatrix const diffusion_WGpGR =
1338 phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dpGR;
1339 DisplacementDimMatrix const diffusion_WLpGR =
1340 phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dpGR;
1341
1342 DisplacementDimMatrix const diffusion_WGpCap =
1343 phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dpCap;
1344 DisplacementDimMatrix const diffusion_WLpCap =
1345 phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dpCap;
1346
1347 DisplacementDimMatrix const diffusion_WGT =
1348 phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dT;
1349 DisplacementDimMatrix const diffusion_WLT =
1350 phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dT;
1351
1352 DisplacementDimMatrix const advection_W = advection_W_G + advection_W_L;
1353 DisplacementDimMatrix const diffusion_W_pGR =
1354 diffusion_WGpGR + diffusion_WLpGR;
1355 DisplacementDimMatrix const diffusion_W_pCap =
1356 diffusion_WGpCap + diffusion_WLpCap;
1357
1358 DisplacementDimMatrix const diffusion_W_T =
1359 diffusion_WGT + diffusion_WLT;
1360
1361 LWpG.noalias() +=
1362 gradNpT * (advection_W + diffusion_W_pGR) * gradNp * w;
1363
1364 LWpC.noalias() +=
1365 gradNpT * (diffusion_W_pCap - advection_W_L) * gradNp * w;
1366
1367 LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
1368
1369 fW.noalias() +=
1370 gradNpT * (advection_W_G * rhoGR + advection_W_L * rhoLR) * b * w;
1371
1372 if (!this->process_data_.apply_mass_lumping)
1373 {
1374 fW.noalias() -=
1375 NpT *
1376 (ip_out.porosity_data.phi *
1377 (current_state.rho_W_LR() -
1378 current_state.constituent_density_data.rho_W_GR) -
1379 rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1380 beta_p_SR) *
1381 s_L_dot * w;
1382 }
1383
1384 fW.noalias() -= NpT * ip_out.porosity_data.phi *
1385 (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
1386
1387 // ---------------------------------------------------------------------
1388 // - temperature equation
1389 // ---------------------------------------------------------------------
1390
1391 MTu.noalias() += NTT * rho_h_eff * mT * Bu * w;
1392
1393 KTT.noalias() += gradNTT * ip.lambda * gradNT * w;
1394
1395 fT.noalias() -= NTT * rho_u_eff_dot * w;
1396
1397 fT.noalias() += gradNTT *
1398 (rhoGR * ip_out.enthalpy_data.h_G * ip.w_GS +
1399 rhoLR * ip_out.enthalpy_data.h_L * ip.w_LS) *
1400 w;
1401
1402 fT.noalias() += gradNTT *
1403 (current_state.constituent_density_data.rho_C_GR *
1404 ip_cv.phase_transition_data.hCG * ip.d_CG +
1405 current_state.constituent_density_data.rho_W_GR *
1406 ip_cv.phase_transition_data.hWG * ip.d_WG) *
1407 w;
1408
1409 fT.noalias() +=
1410 NTT * (rhoGR * ip.w_GS.transpose() + rhoLR * ip.w_LS.transpose()) *
1411 b * w;
1412
1413 // ---------------------------------------------------------------------
1414 // - displacement equation
1415 // ---------------------------------------------------------------------
1416
1417 KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
1418
1419 KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_S_L.chi_S_L * m * Np) * w;
1420
1421 fU.noalias() -= (BuT * current_state.eff_stress_data.sigma -
1422 N_u_op(Nu).transpose() * rho * b) *
1423 w;
1424
1425 if (this->process_data_.apply_mass_lumping)
1426 {
1427 MCpG = MCpG.colwise().sum().eval().asDiagonal();
1428 MCpC = MCpC.colwise().sum().eval().asDiagonal();
1429 MWpG = MWpG.colwise().sum().eval().asDiagonal();
1430 MWpC = MWpC.colwise().sum().eval().asDiagonal();
1431 }
1432 } // int_point-loop
1433}
MeshLib::Element::getID
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:89
ProcessLib::BMatrixPolicyType::BMatrixType
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition BMatrixPolicy.h:55
ProcessLib::TH2M::TH2MLocalAssembler::N_u_op
static constexpr auto & N_u_op
Definition TH2MFEM.h:75
ProcessLib::TH2M::TH2MLocalAssembler::capillary_pressure_index
static const int capillary_pressure_index
Definition TH2MFEM.h:265
ProcessLib::TH2M::TH2MLocalAssembler::capillary_pressure_size
static const int capillary_pressure_size
Definition TH2MFEM.h:266
ProcessLib::TH2M::TH2MLocalAssembler::gas_pressure_index
static const int gas_pressure_index
Definition TH2MFEM.h:263
ProcessLib::TH2M::TH2MLocalAssembler::updateConstitutiveVariables
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
ProcessLib::TH2M::TH2MLocalAssembler::GlobalDimMatrixType
typename ShapeMatricesTypePressure::GlobalDimMatrixType GlobalDimMatrixType
Definition TH2MFEM.h:65
ProcessLib::TH2M::TH2MLocalAssembler::displacement_index
static const int displacement_index
Definition TH2MFEM.h:269
MaterialLib::Fluid::PropertyVariableType::rho
@ rho
density. For some models, rho substitutes p
NumLib::interpolateXCoordinate
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
Definition InitShapeMatrices.h:71
NumLib::interpolateCoordinates
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
Definition InitShapeMatrices.h:82
ProcessLib::LinearBMatrix::computeBMatrix
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
MathLib::KelvinVector::Invariants::identity2
static Eigen::Matrix< double, KelvinVectorSize, 1 > const identity2
Kelvin mapping of 2nd order identity tensor.
Definition KelvinVector.h:107
ProcessLib::TH2M::LocalAssemblerInterface::process_data_
TH2MProcessData< DisplacementDim > & process_data_
Definition LocalAssemblerInterface.h:153
ProcessLib::TH2M::LocalAssemblerInterface::solid_material_
ConstitutiveRelations::SolidConstitutiveRelation< DisplacementDim > const & solid_material_
Definition LocalAssemblerInterface.h:173
ProcessLib::TH2M::LocalAssemblerInterface::current_states_
std::vector< typename ConstitutiveRelations::StatefulData< DisplacementDim > > current_states_
Definition LocalAssemblerInterface.h:156
ProcessLib::TH2M::LocalAssemblerInterface::prev_states_
std::vector< typename ConstitutiveRelations::StatefulDataPrev< DisplacementDim > > prev_states_
Definition LocalAssemblerInterface.h:160
ProcessLib::TH2M::LocalAssemblerInterface::output_data_
std::vector< ConstitutiveRelations::OutputData< DisplacementDim > > output_data_
Definition LocalAssemblerInterface.h:167

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 1440 of file TH2MFEM-impl.h.

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

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 2429 of file TH2MFEM-impl.h.

2433{
2434 auto const gas_pressure =
2435 local_x.template segment<gas_pressure_size>(gas_pressure_index);
2436 auto const capillary_pressure =
2437 local_x.template segment<capillary_pressure_size>(
2438 capillary_pressure_index);
2439 auto const liquid_pressure = (gas_pressure - capillary_pressure).eval();
2440
2441 NumLib::interpolateToHigherOrderNodes<
2442 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2443 DisplacementDim>(this->element_, this->is_axially_symmetric_,
2444 gas_pressure,
2445 *this->process_data_.gas_pressure_interpolated);
2446
2447 NumLib::interpolateToHigherOrderNodes<
2448 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2449 DisplacementDim>(this->element_, this->is_axially_symmetric_,
2450 capillary_pressure,
2451 *this->process_data_.capillary_pressure_interpolated);
2452
2453 NumLib::interpolateToHigherOrderNodes<
2454 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2455 DisplacementDim>(this->element_, this->is_axially_symmetric_,
2456 liquid_pressure,
2457 *this->process_data_.liquid_pressure_interpolated);
2458
2459 auto const temperature =
2460 local_x.template segment<temperature_size>(temperature_index);
2461
2462 NumLib::interpolateToHigherOrderNodes<
2463 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2464 DisplacementDim>(this->element_, this->is_axially_symmetric_,
2465 temperature,
2466 *this->process_data_.temperature_interpolated);
2467
2468 unsigned const n_integration_points =
2469 this->integration_method_.getNumberOfPoints();
2470
2471 double saturation_avg = 0;
2472
2473 ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const models{
2474 this->solid_material_, *this->process_data_.phase_transition_model_};
2475
2476 updateConstitutiveVariables(local_x, local_x_prev, t, dt, models);
2477
2478 for (unsigned ip = 0; ip < n_integration_points; ip++)
2479 {
2480 saturation_avg += this->current_states_[ip].S_L_data.S_L;
2481 }
2482 saturation_avg /= n_integration_points;
2483 (*this->process_data_.element_saturation)[this->element_.getID()] =
2484 saturation_avg;
2485}
NumLib::interpolateToHigherOrderNodes
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
Definition Interpolation.h:107
ProcessLib::Graph::eval
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 2273 of file TH2MFEM-impl.h.

2279{
2280 unsigned const n_integration_points =
2281 this->integration_method_.getNumberOfPoints();
2282
2283 cache.clear();
2284 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
2285 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2286 cache, DisplacementDim, n_integration_points);
2287
2288 for (unsigned ip = 0; ip < n_integration_points; ip++)
2289 {
2290 cache_matrix.col(ip).noalias() = _ip_data[ip].w_GS;
2291 }
2292
2293 return cache;
2294}
MathLib::createZeroedMatrix
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 2299 of file TH2MFEM-impl.h.

2305{
2306 unsigned const n_integration_points =
2307 this->integration_method_.getNumberOfPoints();
2308
2309 cache.clear();
2310 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
2311 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2312 cache, DisplacementDim, n_integration_points);
2313
2314 for (unsigned ip = 0; ip < n_integration_points; ip++)
2315 {
2316 cache_matrix.col(ip).noalias() = _ip_data[ip].w_LS;
2317 }
2318
2319 return cache;
2320}

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 2351 of file TH2MFEM-impl.h.

2357{
2358 unsigned const n_integration_points =
2359 this->integration_method_.getNumberOfPoints();
2360
2361 cache.clear();
2362 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
2363 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2364 cache, DisplacementDim, n_integration_points);
2365
2366 for (unsigned ip = 0; ip < n_integration_points; ip++)
2367 {
2368 cache_matrix.col(ip).noalias() = _ip_data[ip].d_CG;
2369 }
2370
2371 return cache;
2372}

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 2403 of file TH2MFEM-impl.h.

2409{
2410 unsigned const n_integration_points =
2411 this->integration_method_.getNumberOfPoints();
2412
2413 cache.clear();
2414 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
2415 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2416 cache, DisplacementDim, n_integration_points);
2417
2418 for (unsigned ip = 0; ip < n_integration_points; ip++)
2419 {
2420 cache_matrix.col(ip).noalias() = _ip_data[ip].d_WL;
2421 }
2422
2423 return cache;
2424}

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 2377 of file TH2MFEM-impl.h.

2383{
2384 unsigned const n_integration_points =
2385 this->integration_method_.getNumberOfPoints();
2386
2387 cache.clear();
2388 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
2389 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2390 cache, DisplacementDim, n_integration_points);
2391
2392 for (unsigned ip = 0; ip < n_integration_points; ip++)
2393 {
2394 cache_matrix.col(ip).noalias() = _ip_data[ip].d_CL;
2395 }
2396
2397 return cache;
2398}

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 2325 of file TH2MFEM-impl.h.

2331{
2332 unsigned const n_integration_points =
2333 this->integration_method_.getNumberOfPoints();
2334
2335 cache.clear();
2336 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
2337 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2338 cache, DisplacementDim, n_integration_points);
2339
2340 for (unsigned ip = 0; ip < n_integration_points; ip++)
2341 {
2342 cache_matrix.col(ip).noalias() = _ip_data[ip].d_WG;
2343 }
2344
2345 return cache;
2346}

References MathLib::createZeroedMatrix().

◆ 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

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

Definition at line 238 of file TH2MFEM.h.

243 {
244 return ProcessLib::getIntegrationPointScalarData(_ip_data, &IpData::h_S,
245 cache);
246 }
ProcessLib::getIntegrationPointScalarData
std::vector< double > const & getIntegrationPointScalarData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
Definition SetOrGetIntegrationPointData.h:176

References ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data, ProcessLib::getIntegrationPointScalarData(), and ProcessLib::TH2M::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::h_S.

◆ 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 184 of file TH2MFEM.h.

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

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 117 of file TH2MFEM.h.

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

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 ,
int const  )
inlineoverrideprivatevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 160 of file TH2MFEM.h.

164 {
165 unsigned const n_integration_points =
166 this->integration_method_.getNumberOfPoints();
167
168 for (unsigned ip = 0; ip < n_integration_points; ip++)
169 {
170 _ip_data[ip].pushBackState();
171 this->material_states_[ip].pushBackState();
172 }
173
174 for (unsigned ip = 0; ip < n_integration_points; ip++)
175 {
176 this->prev_states_[ip] = this->current_states_[ip];
177 }
178 }

References ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data, ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >::current_states_, NumLib::GenericIntegrationMethod::getNumberOfPoints(), ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >::integration_method_, ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >::material_states_, and ProcessLib::TH2M::LocalAssemblerInterface< DisplacementDim >::prev_states_.

◆ setInitialConditionsConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::setInitialConditionsConcrete ( Eigen::VectorXd const local_x,
double const t,
int const process_id )
overrideprivatevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

851{
852 [[maybe_unused]] auto const matrix_size =
853 gas_pressure_size + capillary_pressure_size + temperature_size +
854 displacement_size;
855
856 assert(local_x.size() == matrix_size);
857
858 auto const capillary_pressure =
859 local_x.template segment<capillary_pressure_size>(
860 capillary_pressure_index);
861
862 auto const p_GR =
863 local_x.template segment<gas_pressure_size>(gas_pressure_index);
864
865 auto const temperature =
866 local_x.template segment<temperature_size>(temperature_index);
867
868 auto const displacement =
869 local_x.template segment<displacement_size>(displacement_index);
870
871 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
872 auto const& medium =
873 *this->process_data_.media_map.getMedium(this->element_.getID());
874 auto const& solid_phase = medium.phase("Solid");
875
876 ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const models{
877 this->solid_material_, *this->process_data_.phase_transition_model_};
878
879 unsigned const n_integration_points =
880 this->integration_method_.getNumberOfPoints();
881
882 for (unsigned ip = 0; ip < n_integration_points; ip++)
883 {
884 MPL::VariableArray vars;
885
886 auto& ip_data = _ip_data[ip];
887 auto& ip_out = this->output_data_[ip];
888 auto& prev_state = this->prev_states_[ip];
889 auto const& Np = ip_data.N_p;
890 auto const& NT = Np;
891 auto const& Nu = ip_data.N_u;
892 auto const& gradNu = ip_data.dNdx_u;
893 auto const x_coord =
894 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
895 ShapeMatricesTypeDisplacement>(
896 this->element_, Nu);
897 ParameterLib::SpatialPosition const pos{
898 std::nullopt, this->element_.getID(), ip,
899 MathLib::Point3d(
900 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
901 ShapeMatricesTypeDisplacement>(
902 this->element_, ip_data.N_u))};
903
904 double const pCap = Np.dot(capillary_pressure);
905 vars.capillary_pressure = pCap;
906
907 double const T = NT.dot(temperature);
908 ConstitutiveRelations::TemperatureData const T_data{
909 T, T}; // T_prev = T in initialization.
910 vars.temperature = T;
911
912 auto const Bu =
913 LinearBMatrix::computeBMatrix<DisplacementDim,
914 ShapeFunctionDisplacement::NPOINTS,
915 typename BMatricesType::BMatrixType>(
916 gradNu, Nu, x_coord, this->is_axially_symmetric_);
917
918 auto& eps = ip_out.eps_data.eps;
919 eps.noalias() = Bu * displacement;
920
921 // Set volumetric strain rate for the general case without swelling.
922 vars.volumetric_strain = Invariants::trace(eps);
923
924 double const S_L =
925 medium.property(MPL::PropertyType::saturation)
926 .template value<double>(
927 vars, pos, t, std::numeric_limits<double>::quiet_NaN());
928 this->prev_states_[ip].S_L_data->S_L = S_L;
929
930 // TODO (naumov) Double computation of C_el might be avoided if
931 // updateConstitutiveVariables is called before. But it might interfere
932 // with eps_m initialization.
933 ConstitutiveRelations::ElasticTangentStiffnessData<DisplacementDim>
934 C_el_data;
935 models.elastic_tangent_stiffness_model.eval({pos, t, dt}, T_data,
936 C_el_data);
937 auto const& C_el = C_el_data.stiffness_tensor;
938
939 // Set eps_m_prev from potentially non-zero eps and sigma_sw from
940 // restart.
941 auto const& sigma_sw = this->current_states_[ip].swelling_data.sigma_sw;
942 prev_state.mechanical_strain_data->eps_m.noalias() =
943 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
944 ? eps + C_el.inverse() * sigma_sw
945 : eps;
946
947 if (this->process_data_.initial_stress.isTotalStress())
948 {
949 auto const alpha_b =
950 medium.property(MPL::PropertyType::biot_coefficient)
951 .template value<double>(vars, pos, t, 0.0 /*dt*/);
952
953 vars.liquid_saturation = S_L;
954 double const bishop =
955 medium.property(MPL::PropertyType::bishops_effective_stress)
956 .template value<double>(vars, pos, t, 0.0 /*dt*/);
957
958 this->current_states_[ip].eff_stress_data.sigma.noalias() +=
959 alpha_b * Np.dot(p_GR - bishop * capillary_pressure) *
960 Invariants::identity2;
961 this->prev_states_[ip].eff_stress_data =
962 this->current_states_[ip].eff_stress_data;
963 }
964 }
965
966 // local_x_prev equal to local_x s.t. the local_x_dot is zero.
967 updateConstitutiveVariables(local_x, local_x, t, 0, models);
968
969 for (unsigned ip = 0; ip < n_integration_points; ip++)
970 {
971 auto& ip_data = _ip_data[ip];
972 ip_data.pushBackState();
973
974 this->material_states_[ip].pushBackState();
975 this->prev_states_[ip] = this->current_states_[ip];
976 }
977}

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 787 of file TH2MFEM-impl.h.

790{
791 if (integration_order !=
792 static_cast<int>(this->integration_method_.getIntegrationOrder()))
793 {
794 OGS_FATAL(
795 "Setting integration point initial conditions; The integration "
796 "order of the local assembler for element {:d} is different "
797 "from the integration order in the initial condition.",
798 this->element_.getID());
799 }
800
801 if (name == "sigma" && this->process_data_.initial_stress.value)
802 {
803 OGS_FATAL(
804 "Setting initial conditions for stress from integration "
805 "point data and from a parameter '{:s}' is not possible "
806 "simultaneously.",
807 this->process_data_.initial_stress.value->name);
808 }
809
810 if (name.starts_with("material_state_variable_"))
811 {
812 name.remove_prefix(24);
813 DBUG("Setting material state variable '{:s}'", name);
814
815 auto const& internal_variables =
816 this->solid_material_.getInternalVariables();
817 if (auto const iv = std::find_if(
818 begin(internal_variables), end(internal_variables),
819 [&name](auto const& iv) { return iv.name == name; });
820 iv != end(internal_variables))
821 {
822 DBUG("Setting material state variable '{:s}'", name);
823 return ProcessLib::setIntegrationPointDataMaterialStateVariables(
824 values, this->material_states_,
825 &ConstitutiveRelations::MaterialStateData<
826 DisplacementDim>::material_state_variables,
827 iv->reference);
828 }
829
830 WARN(
831 "Could not find variable {:s} in solid material model's "
832 "internal variables.",
833 name);
834 return 0;
835 }
836
837 // TODO this logic could be pulled out of the local assembler into the
838 // process. That might lead to a slightly better performance due to less
839 // string comparisons.
840 return ProcessLib::Reflection::reflectSetIPData<DisplacementDim>(
841 name, values, this->current_states_);
842}
OGS_FATAL
#define OGS_FATAL(...)
Definition Error.h:26
DBUG
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
WARN
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40
ProcessLib::setIntegrationPointDataMaterialStateVariables
std::size_t setIntegrationPointDataMaterialStateVariables(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType member, std::function< std::span< double >(MaterialStateVariables &)> get_values_span)
Definition SetOrGetIntegrationPointData.h:270

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

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, ProcessLib::TH2M::ConstitutiveRelations::BiotModel::eval(), ProcessLib::TH2M::ConstitutiveRelations::SaturationModel::eval(), ProcessLib::TH2M::ConstitutiveRelations::PureLiquidDensityModel::eval(), ProcessLib::TH2M::ConstitutiveRelations::PhaseTransitionModel::eval(), ProcessLib::TH2M::ConstitutiveRelations::PorosityModel::eval(), ProcessLib::TH2M::ConstitutiveRelations::BishopsModel::eval(), ProcessLib::TH2M::ConstitutiveRelations::SolidDensityModel::eval(), MaterialPropertyLib::VariableArray::gas_phase_pressure, NumLib::interpolateCoordinates(), NumLib::interpolateXCoordinate(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::mechanical_strain_model, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::permeability_model, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::phase_transition_model, MaterialPropertyLib::VariableArray::porosity, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::porosity_model, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::pure_liquid_density_model, 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 >::solid_density_model, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::swelling_model, MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::thermal_conductivity, 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 255 of file TH2MFEM.h.

Referenced by ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::TH2MLocalAssembler(), ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtEnthalpySolid(), 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

Definition at line 259 of file TH2MFEM.h.

Referenced by ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::TH2MLocalAssembler(), and ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getShapeMatrix().

◆ C_index

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

Definition at line 273 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 274 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 265 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 266 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 269 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 270 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 263 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 264 of file TH2MFEM.h.

◆ KelvinVectorSize

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
int const ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::KelvinVectorSize
static
Initial value:
=
MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim)
MathLib::KelvinVector::kelvin_vector_dimensions
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:
= MathLib::eigenBlockMatrixView<
DisplacementDim,
typename ShapeMatricesTypeDisplacement::NodalRowVectorType>
MathLib::eigenBlockMatrixView
constexpr Eigen::CwiseNullaryOp< EigenBlockMatrixViewFunctor< D, M >, typename EigenBlockMatrixViewFunctor< D, M >::Matrix > eigenBlockMatrixView(const Eigen::MatrixBase< M > &matrix)
Definition EigenBlockMatrixView.h:49
EigenFixedShapeMatrixPolicy::NodalRowVectorType
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
Definition ShapeMatrixPolicy.h:78

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 267 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 268 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 275 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 276 of file TH2MFEM.h.

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