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::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)
 

Private Attributes

std::vector< 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 202 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<ShapeMatricesTypeDisplacement,
ShapeMatricesTypePressure>
ProcessLib::TH2M::TH2MLocalAssembler::ShapeMatricesTypeDisplacement
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
Definition TH2MFEM.h:51
ProcessLib::TH2M::TH2MLocalAssembler::ShapeMatricesTypePressure
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
Definition TH2MFEM.h:54

Definition at line 204 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:210
ProcessLib::TH2M::TH2MLocalAssembler::_ip_data
std::vector< IpData > _ip_data
Definition TH2MFEM.h:206
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:127

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

927{
928 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
929 temperature_size + displacement_size;
930 assert(local_x.size() == matrix_size);
931
932 auto const capillary_pressure =
933 Eigen::Map<VectorType<capillary_pressure_size> const>(
934 local_x.data() + capillary_pressure_index, capillary_pressure_size);
935
936 auto const capillary_pressure_prev =
937 Eigen::Map<VectorType<capillary_pressure_size> const>(
938 local_x_prev.data() + capillary_pressure_index,
939 capillary_pressure_size);
940
941 // pointer to local_M_data vector
942 auto local_M =
943 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
944 local_M_data, matrix_size, matrix_size);
945
946 // pointer to local_K_data vector
947 auto local_K =
948 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
949 local_K_data, matrix_size, matrix_size);
950
951 // pointer to local_rhs_data vector
952 auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
953 local_rhs_data, matrix_size);
954
955 // component-formulation
956 // W - liquid phase main component
957 // C - gas phase main component
958 // pointer-matrices to the mass matrix - C component equation
959 auto MCpG = local_M.template block<C_size, gas_pressure_size>(
960 C_index, gas_pressure_index);
961 auto MCpC = local_M.template block<C_size, capillary_pressure_size>(
962 C_index, capillary_pressure_index);
963 auto MCT = local_M.template block<C_size, temperature_size>(
964 C_index, temperature_index);
965 auto MCu = local_M.template block<C_size, displacement_size>(
966 C_index, displacement_index);
967
968 // pointer-matrices to the stiffness matrix - C component equation
969 auto LCpG = local_K.template block<C_size, gas_pressure_size>(
970 C_index, gas_pressure_index);
971 auto LCpC = local_K.template block<C_size, capillary_pressure_size>(
972 C_index, capillary_pressure_index);
973 auto LCT = local_K.template block<C_size, temperature_size>(
974 C_index, temperature_index);
975
976 // pointer-matrices to the mass matrix - W component equation
977 auto MWpG = local_M.template block<W_size, gas_pressure_size>(
978 W_index, gas_pressure_index);
979 auto MWpC = local_M.template block<W_size, capillary_pressure_size>(
980 W_index, capillary_pressure_index);
981 auto MWT = local_M.template block<W_size, temperature_size>(
982 W_index, temperature_index);
983 auto MWu = local_M.template block<W_size, displacement_size>(
984 W_index, displacement_index);
985
986 // pointer-matrices to the stiffness matrix - W component equation
987 auto LWpG = local_K.template block<W_size, gas_pressure_size>(
988 W_index, gas_pressure_index);
989 auto LWpC = local_K.template block<W_size, capillary_pressure_size>(
990 W_index, capillary_pressure_index);
991 auto LWT = local_K.template block<W_size, temperature_size>(
992 W_index, temperature_index);
993
994 // pointer-matrices to the mass matrix - temperature equation
995 auto MTu = local_M.template block<temperature_size, displacement_size>(
996 temperature_index, displacement_index);
997
998 // pointer-matrices to the stiffness matrix - temperature equation
999 auto KTT = local_K.template block<temperature_size, temperature_size>(
1000 temperature_index, temperature_index);
1001
1002 // pointer-matrices to the stiffness matrix - displacement equation
1003 auto KUpG = local_K.template block<displacement_size, gas_pressure_size>(
1004 displacement_index, gas_pressure_index);
1005 auto KUpC =
1006 local_K.template block<displacement_size, capillary_pressure_size>(
1007 displacement_index, capillary_pressure_index);
1008
1009 // pointer-vectors to the right hand side terms - C-component equation
1010 auto fC = local_f.template segment<C_size>(C_index);
1011 // pointer-vectors to the right hand side terms - W-component equation
1012 auto fW = local_f.template segment<W_size>(W_index);
1013 // pointer-vectors to the right hand side terms - temperature equation
1014 auto fT = local_f.template segment<temperature_size>(temperature_index);
1015 // pointer-vectors to the right hand side terms - displacement equation
1016 auto fU = local_f.template segment<displacement_size>(displacement_index);
1017
1018 unsigned const n_integration_points =
1019 this->integration_method_.getNumberOfPoints();
1020
1021 ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const models{
1022 this->solid_material_, *this->process_data_.phase_transition_model_};
1023
1024 auto const [ip_constitutive_data, ip_constitutive_variables] =
1025 updateConstitutiveVariables(
1026 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1027 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1028 local_x_prev.size()),
1029 t, dt, models);
1030
1031 for (unsigned int_point = 0; int_point < n_integration_points; int_point++)
1032 {
1033 auto& ip = _ip_data[int_point];
1034 auto& ip_cv = ip_constitutive_variables[int_point];
1035 auto& ip_out = this->output_data_[int_point];
1036 auto& current_state = this->current_states_[int_point];
1037 auto const& prev_state = this->prev_states_[int_point];
1038
1039 auto const& Np = ip.N_p;
1040 auto const& NT = Np;
1041 auto const& Nu = ip.N_u;
1042 ParameterLib::SpatialPosition const pos{
1043 std::nullopt, this->element_.getID(), int_point,
1044 MathLib::Point3d(
1045 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
1046 ShapeMatricesTypeDisplacement>(
1047 this->element_, Nu))};
1048
1049 auto const& NpT = Np.transpose().eval();
1050 auto const& NTT = NT.transpose().eval();
1051
1052 auto const& gradNp = ip.dNdx_p;
1053 auto const& gradNT = gradNp;
1054 auto const& gradNu = ip.dNdx_u;
1055
1056 auto const& gradNpT = gradNp.transpose().eval();
1057 auto const& gradNTT = gradNT.transpose().eval();
1058
1059 auto const& w = ip.integration_weight;
1060
1061 auto const& m = Invariants::identity2;
1062
1063 auto const mT = m.transpose().eval();
1064
1065 auto const x_coord =
1066 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
1067 ShapeMatricesTypeDisplacement>(
1068 this->element_, Nu);
1069
1070 auto const Bu =
1071 LinearBMatrix::computeBMatrix<DisplacementDim,
1072 ShapeFunctionDisplacement::NPOINTS,
1073 typename BMatricesType::BMatrixType>(
1074 gradNu, Nu, x_coord, this->is_axially_symmetric_);
1075
1076 auto const BuT = Bu.transpose().eval();
1077
1078 double const pCap = Np.dot(capillary_pressure);
1079 double const pCap_prev = Np.dot(capillary_pressure_prev);
1080
1081 double const beta_T_SR = ip_cv.s_therm_exp_data.beta_T_SR;
1082
1083 auto const I =
1084 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
1085
1086 double const sD_G =
1087 ip_cv.phase_transition_data.diffusion_coefficient_vapour;
1088 double const sD_L =
1089 ip_cv.phase_transition_data.diffusion_coefficient_solute;
1090
1091 GlobalDimMatrixType const D_C_G = sD_G * I;
1092 GlobalDimMatrixType const D_W_G = sD_G * I;
1093 GlobalDimMatrixType const D_C_L = sD_L * I;
1094 GlobalDimMatrixType const D_W_L = sD_L * I;
1095
1096 auto const s_L = current_state.S_L_data.S_L;
1097 auto const s_G = 1. - s_L;
1098 auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
1099
1100 auto& alpha_B = ip_cv.biot_data();
1101 auto& beta_p_SR = ip_cv.beta_p_SR();
1102
1103 auto const& b = this->process_data_.specific_body_force;
1104
1105 // volume fraction
1106 auto const phi_G = s_G * ip_out.porosity_data.phi;
1107 auto const phi_L = s_L * ip_out.porosity_data.phi;
1108 auto const phi_S = 1. - ip_out.porosity_data.phi;
1109
1110 auto const rhoGR = ip_out.fluid_density_data.rho_GR;
1111 auto const rhoLR = ip_out.fluid_density_data.rho_LR;
1112
1113 // effective density
1114 auto const rho = phi_G * rhoGR + phi_L * rhoLR +
1115 phi_S * ip_out.solid_density_data.rho_SR;
1116
1117 // abbreviations
1118 const double rho_C_FR =
1119 s_G * current_state.constituent_density_data.rho_C_GR +
1120 s_L * current_state.constituent_density_data.rho_C_LR;
1121 const double rho_W_FR =
1122 s_G * current_state.constituent_density_data.rho_W_GR +
1123 s_L * current_state.rho_W_LR();
1124
1125 auto const rho_C_GR_dot =
1126 (current_state.constituent_density_data.rho_C_GR -
1127 prev_state.constituent_density_data->rho_C_GR) /
1128 dt;
1129 auto const rho_C_LR_dot =
1130 (current_state.constituent_density_data.rho_C_LR -
1131 prev_state.constituent_density_data->rho_C_LR) /
1132 dt;
1133 auto const rho_W_GR_dot =
1134 (current_state.constituent_density_data.rho_W_GR -
1135 prev_state.constituent_density_data->rho_W_GR) /
1136 dt;
1137 auto const rho_W_LR_dot =
1138 (current_state.rho_W_LR() - **prev_state.rho_W_LR) / dt;
1139
1140 GlobalDimMatrixType const k_over_mu_G =
1141 ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_G /
1142 ip_cv.viscosity_data.mu_GR;
1143 GlobalDimMatrixType const k_over_mu_L =
1144 ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_L /
1145 ip_cv.viscosity_data.mu_LR;
1146
1147 // ---------------------------------------------------------------------
1148 // C-component equation
1149 // ---------------------------------------------------------------------
1150
1151 MCpG.noalias() += NpT * rho_C_FR *
1152 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1153 Np * w;
1154 MCpC.noalias() -= NpT * rho_C_FR *
1155 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1156 s_L * Np * w;
1157
1158 if (this->process_data_.apply_mass_lumping)
1159 {
1160 if (pCap - pCap_prev != 0.) // avoid division by Zero
1161 {
1162 MCpC.noalias() +=
1163 NpT *
1164 (ip_out.porosity_data.phi *
1165 (current_state.constituent_density_data.rho_C_LR -
1166 current_state.constituent_density_data.rho_C_GR) -
1167 rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1168 beta_p_SR) *
1169 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1170 }
1171 }
1172
1173 MCT.noalias() -= NpT * rho_C_FR * (alpha_B - ip_out.porosity_data.phi) *
1174 beta_T_SR * Np * w;
1175 MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
1176
1177 using DisplacementDimMatrix =
1178 Eigen::Matrix<double, DisplacementDim, DisplacementDim>;
1179
1180 DisplacementDimMatrix const advection_C_G =
1181 current_state.constituent_density_data.rho_C_GR * k_over_mu_G;
1182 DisplacementDimMatrix const advection_C_L =
1183 current_state.constituent_density_data.rho_C_LR * k_over_mu_L;
1184
1185 DisplacementDimMatrix const diffusion_CGpGR =
1186 -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dpGR;
1187 DisplacementDimMatrix const diffusion_CLpGR =
1188 -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dpGR;
1189
1190 DisplacementDimMatrix const diffusion_CGpCap =
1191 -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dpCap;
1192 DisplacementDimMatrix const diffusion_CLpCap =
1193 -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dpCap;
1194
1195 DisplacementDimMatrix const diffusion_CGT =
1196 -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dT;
1197 DisplacementDimMatrix const diffusion_CLT =
1198 -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dT;
1199
1200 DisplacementDimMatrix const advection_C = advection_C_G + advection_C_L;
1201 DisplacementDimMatrix const diffusion_C_pGR =
1202 diffusion_CGpGR + diffusion_CLpGR;
1203 DisplacementDimMatrix const diffusion_C_pCap =
1204 diffusion_CGpCap + diffusion_CLpCap;
1205
1206 DisplacementDimMatrix const diffusion_C_T =
1207 diffusion_CGT + diffusion_CLT;
1208
1209 LCpG.noalias() +=
1210 gradNpT * (advection_C + diffusion_C_pGR) * gradNp * w;
1211
1212 LCpC.noalias() +=
1213 gradNpT * (diffusion_C_pCap - advection_C_L) * gradNp * w;
1214
1215 LCT.noalias() += gradNpT * (diffusion_C_T)*gradNp * w;
1216
1217 fC.noalias() +=
1218 gradNpT * (advection_C_G * rhoGR + advection_C_L * rhoLR) * b * w;
1219
1220 if (!this->process_data_.apply_mass_lumping)
1221 {
1222 fC.noalias() -=
1223 NpT *
1224 (ip_out.porosity_data.phi *
1225 (current_state.constituent_density_data.rho_C_LR -
1226 current_state.constituent_density_data.rho_C_GR) -
1227 rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1228 beta_p_SR) *
1229 s_L_dot * w;
1230 }
1231 // fC_III
1232 fC.noalias() -= NpT * ip_out.porosity_data.phi *
1233 (s_G * rho_C_GR_dot + s_L * rho_C_LR_dot) * w;
1234
1235 // ---------------------------------------------------------------------
1236 // W-component equation
1237 // ---------------------------------------------------------------------
1238
1239 MWpG.noalias() += NpT * rho_W_FR *
1240 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1241 Np * w;
1242 MWpC.noalias() -= NpT * rho_W_FR *
1243 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1244 s_L * Np * w;
1245
1246 if (this->process_data_.apply_mass_lumping)
1247 {
1248 if (pCap - pCap_prev != 0.) // avoid division by Zero
1249 {
1250 MWpC.noalias() +=
1251 NpT *
1252 (ip_out.porosity_data.phi *
1253 (current_state.rho_W_LR() -
1254 current_state.constituent_density_data.rho_W_GR) -
1255 rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1256 beta_p_SR) *
1257 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1258 }
1259 }
1260
1261 MWT.noalias() -= NpT * rho_W_FR * (alpha_B - ip_out.porosity_data.phi) *
1262 beta_T_SR * Np * w;
1263
1264 MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
1265
1266 DisplacementDimMatrix const advection_W_G =
1267 current_state.constituent_density_data.rho_W_GR * k_over_mu_G;
1268 DisplacementDimMatrix const advection_W_L =
1269 current_state.rho_W_LR() * k_over_mu_L;
1270
1271 DisplacementDimMatrix const diffusion_WGpGR =
1272 phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dpGR;
1273 DisplacementDimMatrix const diffusion_WLpGR =
1274 phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dpGR;
1275
1276 DisplacementDimMatrix const diffusion_WGpCap =
1277 phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dpCap;
1278 DisplacementDimMatrix const diffusion_WLpCap =
1279 phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dpCap;
1280
1281 DisplacementDimMatrix const diffusion_WGT =
1282 phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dT;
1283 DisplacementDimMatrix const diffusion_WLT =
1284 phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dT;
1285
1286 DisplacementDimMatrix const advection_W = advection_W_G + advection_W_L;
1287 DisplacementDimMatrix const diffusion_W_pGR =
1288 diffusion_WGpGR + diffusion_WLpGR;
1289 DisplacementDimMatrix const diffusion_W_pCap =
1290 diffusion_WGpCap + diffusion_WLpCap;
1291
1292 DisplacementDimMatrix const diffusion_W_T =
1293 diffusion_WGT + diffusion_WLT;
1294
1295 LWpG.noalias() +=
1296 gradNpT * (advection_W + diffusion_W_pGR) * gradNp * w;
1297
1298 LWpC.noalias() +=
1299 gradNpT * (diffusion_W_pCap - advection_W_L) * gradNp * w;
1300
1301 LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
1302
1303 fW.noalias() +=
1304 gradNpT * (advection_W_G * rhoGR + advection_W_L * rhoLR) * b * w;
1305
1306 if (!this->process_data_.apply_mass_lumping)
1307 {
1308 fW.noalias() -=
1309 NpT *
1310 (ip_out.porosity_data.phi *
1311 (current_state.rho_W_LR() -
1312 current_state.constituent_density_data.rho_W_GR) -
1313 rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1314 beta_p_SR) *
1315 s_L_dot * w;
1316 }
1317
1318 fW.noalias() -= NpT * ip_out.porosity_data.phi *
1319 (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
1320
1321 // ---------------------------------------------------------------------
1322 // - temperature equation
1323 // ---------------------------------------------------------------------
1324
1325 MTu.noalias() += NTT *
1326 ip_cv.effective_volumetric_enthalpy_data.rho_h_eff *
1327 mT * Bu * w;
1328
1329 KTT.noalias() +=
1330 gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w;
1331
1332 auto const rho_u_eff_dot = (current_state.internal_energy_data() -
1333 **prev_state.internal_energy_data) /
1334 dt;
1335 fT.noalias() -= NTT * rho_u_eff_dot * w;
1336
1337 fT.noalias() += gradNTT *
1338 (rhoGR * ip_out.enthalpy_data.h_G *
1339 ip_out.darcy_velocity_data.w_GS +
1340 rhoLR * ip_out.enthalpy_data.h_L *
1341 ip_out.darcy_velocity_data.w_LS) *
1342 w;
1343
1344 fT.noalias() += gradNTT *
1345 (current_state.constituent_density_data.rho_C_GR *
1346 ip_cv.phase_transition_data.hCG *
1347 ip_out.diffusion_velocity_data.d_CG +
1348 current_state.constituent_density_data.rho_W_GR *
1349 ip_cv.phase_transition_data.hWG *
1350 ip_out.diffusion_velocity_data.d_WG) *
1351 w;
1352
1353 fT.noalias() += NTT *
1354 (rhoGR * ip_out.darcy_velocity_data.w_GS.transpose() +
1355 rhoLR * ip_out.darcy_velocity_data.w_LS.transpose()) *
1356 b * w;
1357
1358 // ---------------------------------------------------------------------
1359 // - displacement equation
1360 // ---------------------------------------------------------------------
1361
1362 KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
1363
1364 KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_S_L.chi_S_L * m * Np) * w;
1365
1366 fU.noalias() -= (BuT * current_state.eff_stress_data.sigma -
1367 N_u_op(Nu).transpose() * rho * b) *
1368 w;
1369
1370 if (this->process_data_.apply_mass_lumping)
1371 {
1372 MCpG = MCpG.colwise().sum().eval().asDiagonal();
1373 MCpC = MCpC.colwise().sum().eval().asDiagonal();
1374 MWpG = MWpG.colwise().sum().eval().asDiagonal();
1375 MWpC = MWpC.colwise().sum().eval().asDiagonal();
1376 }
1377 } // int_point-loop
1378}
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:216
ProcessLib::TH2M::TH2MLocalAssembler::capillary_pressure_size
static const int capillary_pressure_size
Definition TH2MFEM.h:217
ProcessLib::TH2M::TH2MLocalAssembler::gas_pressure_index
static const int gas_pressure_index
Definition TH2MFEM.h:214
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:220
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:111
ProcessLib::TH2M::LocalAssemblerInterface::solid_material_
ConstitutiveRelations::SolidConstitutiveRelation< DisplacementDim > const & solid_material_
Definition LocalAssemblerInterface.h:131
ProcessLib::TH2M::LocalAssemblerInterface::current_states_
std::vector< typename ConstitutiveRelations::StatefulData< DisplacementDim > > current_states_
Definition LocalAssemblerInterface.h:114
ProcessLib::TH2M::LocalAssemblerInterface::prev_states_
std::vector< typename ConstitutiveRelations::StatefulDataPrev< DisplacementDim > > prev_states_
Definition LocalAssemblerInterface.h:118
ProcessLib::TH2M::LocalAssemblerInterface::output_data_
std::vector< ConstitutiveRelations::OutputData< DisplacementDim > > output_data_
Definition LocalAssemblerInterface.h:125

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

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

2247{
2248 auto const gas_pressure =
2249 local_x.template segment<gas_pressure_size>(gas_pressure_index);
2250 auto const capillary_pressure =
2251 local_x.template segment<capillary_pressure_size>(
2252 capillary_pressure_index);
2253 auto const liquid_pressure = (gas_pressure - capillary_pressure).eval();
2254
2255 NumLib::interpolateToHigherOrderNodes<
2256 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2257 DisplacementDim>(this->element_, this->is_axially_symmetric_,
2258 gas_pressure,
2259 *this->process_data_.gas_pressure_interpolated);
2260
2261 NumLib::interpolateToHigherOrderNodes<
2262 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2263 DisplacementDim>(this->element_, this->is_axially_symmetric_,
2264 capillary_pressure,
2265 *this->process_data_.capillary_pressure_interpolated);
2266
2267 NumLib::interpolateToHigherOrderNodes<
2268 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2269 DisplacementDim>(this->element_, this->is_axially_symmetric_,
2270 liquid_pressure,
2271 *this->process_data_.liquid_pressure_interpolated);
2272
2273 auto const temperature =
2274 local_x.template segment<temperature_size>(temperature_index);
2275
2276 NumLib::interpolateToHigherOrderNodes<
2277 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2278 DisplacementDim>(this->element_, this->is_axially_symmetric_,
2279 temperature,
2280 *this->process_data_.temperature_interpolated);
2281
2282 unsigned const n_integration_points =
2283 this->integration_method_.getNumberOfPoints();
2284
2285 double saturation_avg = 0;
2286
2287 ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const models{
2288 this->solid_material_, *this->process_data_.phase_transition_model_};
2289
2290 updateConstitutiveVariables(local_x, local_x_prev, t, dt, models);
2291
2292 for (unsigned ip = 0; ip < n_integration_points; ip++)
2293 {
2294 saturation_avg += this->current_states_[ip].S_L_data.S_L;
2295 }
2296 saturation_avg /= n_integration_points;
2297 (*this->process_data_.element_saturation)[this->element_.getID()] =
2298 saturation_avg;
2299}
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().

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

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

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 material_state.pushBackState();
151 }
152
153 for (unsigned ip = 0; ip < n_integration_points; ip++)
154 {
155 this->prev_states_[ip] = this->current_states_[ip];
156 }
157 }
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:123

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

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

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

792{
793 [[maybe_unused]] auto const matrix_size =
794 gas_pressure_size + capillary_pressure_size + temperature_size +
795 displacement_size;
796
797 assert(local_x.size() == matrix_size);
798
799 auto const capillary_pressure =
800 local_x.template segment<capillary_pressure_size>(
801 capillary_pressure_index);
802
803 auto const p_GR =
804 local_x.template segment<gas_pressure_size>(gas_pressure_index);
805
806 auto const temperature =
807 local_x.template segment<temperature_size>(temperature_index);
808
809 auto const displacement =
810 local_x.template segment<displacement_size>(displacement_index);
811
812 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
813 auto const& medium =
814 *this->process_data_.media_map.getMedium(this->element_.getID());
815 auto const& solid_phase = medium.phase("Solid");
816
817 ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const models{
818 this->solid_material_, *this->process_data_.phase_transition_model_};
819
820 unsigned const n_integration_points =
821 this->integration_method_.getNumberOfPoints();
822
823 for (unsigned ip = 0; ip < n_integration_points; ip++)
824 {
825 MPL::VariableArray vars;
826
827 auto& ip_data = _ip_data[ip];
828 auto& ip_out = this->output_data_[ip];
829 auto& prev_state = this->prev_states_[ip];
830 auto const& Np = ip_data.N_p;
831 auto const& NT = Np;
832 auto const& Nu = ip_data.N_u;
833 auto const& gradNu = ip_data.dNdx_u;
834 auto const x_coord =
835 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
836 ShapeMatricesTypeDisplacement>(
837 this->element_, Nu);
838 ParameterLib::SpatialPosition const pos{
839 std::nullopt, this->element_.getID(), ip,
840 MathLib::Point3d(
841 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
842 ShapeMatricesTypeDisplacement>(
843 this->element_, ip_data.N_u))};
844
845 double const pCap = Np.dot(capillary_pressure);
846 vars.capillary_pressure = pCap;
847
848 double const T = NT.dot(temperature);
849 ConstitutiveRelations::TemperatureData const T_data{
850 T, T}; // T_prev = T in initialization.
851 vars.temperature = T;
852
853 auto const Bu =
854 LinearBMatrix::computeBMatrix<DisplacementDim,
855 ShapeFunctionDisplacement::NPOINTS,
856 typename BMatricesType::BMatrixType>(
857 gradNu, Nu, x_coord, this->is_axially_symmetric_);
858
859 auto& eps = ip_out.eps_data.eps;
860 eps.noalias() = Bu * displacement;
861
862 // Set volumetric strain rate for the general case without swelling.
863 vars.volumetric_strain = Invariants::trace(eps);
864
865 double const S_L =
866 medium.property(MPL::PropertyType::saturation)
867 .template value<double>(
868 vars, pos, t, std::numeric_limits<double>::quiet_NaN());
869 this->prev_states_[ip].S_L_data->S_L = S_L;
870
871 // TODO (naumov) Double computation of C_el might be avoided if
872 // updateConstitutiveVariables is called before. But it might interfere
873 // with eps_m initialization.
874 ConstitutiveRelations::ElasticTangentStiffnessData<DisplacementDim>
875 C_el_data;
876 models.elastic_tangent_stiffness_model.eval({pos, t, dt}, T_data,
877 C_el_data);
878 auto const& C_el = C_el_data.stiffness_tensor;
879
880 // Set eps_m_prev from potentially non-zero eps and sigma_sw from
881 // restart.
882 auto const& sigma_sw = this->current_states_[ip].swelling_data.sigma_sw;
883 prev_state.mechanical_strain_data->eps_m.noalias() =
884 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
885 ? eps + C_el.inverse() * sigma_sw
886 : eps;
887
888 if (this->process_data_.initial_stress.isTotalStress())
889 {
890 auto const alpha_b =
891 medium.property(MPL::PropertyType::biot_coefficient)
892 .template value<double>(vars, pos, t, 0.0 /*dt*/);
893
894 vars.liquid_saturation = S_L;
895 double const bishop =
896 medium.property(MPL::PropertyType::bishops_effective_stress)
897 .template value<double>(vars, pos, t, 0.0 /*dt*/);
898
899 this->current_states_[ip].eff_stress_data.sigma.noalias() +=
900 alpha_b * Np.dot(p_GR - bishop * capillary_pressure) *
901 Invariants::identity2;
902 this->prev_states_[ip].eff_stress_data =
903 this->current_states_[ip].eff_stress_data;
904 }
905 }
906
907 // local_x_prev equal to local_x s.t. the local_x_dot is zero.
908 updateConstitutiveVariables(local_x, local_x, t, 0, models);
909
910 for (unsigned ip = 0; ip < n_integration_points; ip++)
911 {
912 this->material_states_[ip].pushBackState();
913 this->prev_states_[ip] = this->current_states_[ip];
914 }
915}

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

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

References ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::beta_p_SR_model, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::biot_model, 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::ViscosityModel::eval(), ProcessLib::TH2M::ConstitutiveRelations::SolidDensityModel::eval(), ProcessLib::TH2M::ConstitutiveRelations::SolidHeatCapacityModel::eval(), NumLib::interpolateCoordinates(), NumLib::interpolateXCoordinate(), ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::mechanical_strain_model, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::permeability_model, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::phase_transition_model, 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 >::solid_heat_capacity_model, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::swelling_model, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::thermal_conductivity_model, ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::total_stress_model, and ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels< DisplacementDim >::viscosity_model.

Member Data Documentation

◆ _ip_data

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

Definition at line 206 of file TH2MFEM.h.

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

◆ _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 210 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 224 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 225 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 216 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 217 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 220 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 221 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 214 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 215 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 218 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 219 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 226 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 227 of file TH2MFEM.h.

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