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

Detailed Description

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

Definition at line 48 of file TH2MFEM.h.

#include <TH2MFEM.h>

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

Public Types

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

Public Member Functions

 TH2MLocalAssembler (TH2MLocalAssembler const &)=delete
 
 TH2MLocalAssembler (TH2MLocalAssembler &&)=delete
 
 TH2MLocalAssembler (MeshLib::Element const &e, std::size_t const, bool const is_axially_symmetric, unsigned const integration_order, TH2MProcessData< DisplacementDim > &process_data)
 
- Public Member Functions inherited from ProcessLib::LocalAssemblerInterface
virtual ~LocalAssemblerInterface ()=default
 
void setInitialConditions (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, bool const use_monolithic_scheme, int const process_id)
 
virtual void initialize (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
 
virtual void preAssemble (double const, double const, std::vector< double > const &)
 
virtual void assembleForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, 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_xdot, const double dxdot_dx, const double dx_dx, 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_dot, 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, double const t, double const dt)
 
void postNonLinearSolver (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, GlobalVector const &xdot, double const t, double const dt, bool const use_monolithic_scheme, int const process_id)
 
virtual std::vector< double > interpolateNodalValuesToIntegrationPoints (std::vector< double > const &)
 
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. More...
 
- Public Member Functions inherited from NumLib::ExtrapolatableElement
virtual ~ExtrapolatableElement ()=default
 

Static Public Attributes

static int const KelvinVectorSize
 

Private Types

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

Private Member Functions

std::size_t setIPDataInitialConditions (std::string const &name, double const *values, int const integration_order) override
 
void setInitialConditionsConcrete (std::vector< double > const &local_x, double const t, bool const, int const) override
 
void assemble (double const, double const, std::vector< double > const &, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
 
void assembleWithJacobian (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_xdot, const double, const double, 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 &, double const, double const) override
 
void computeSecondaryVariableConcrete (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_dot) override
 
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix (const unsigned integration_point) const override
 Provides the shape matrix at the given integration point. More...
 
std::vector< double > const & getIntPtDarcyVelocityGas (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::vector< double > const & getIntPtDarcyVelocityLiquid (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::vector< ConstitutiveVariables< DisplacementDim > > updateConstitutiveVariables (Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_dot, double const t, double const dt)
 
std::size_t setSigma (double const *values)
 
std::vector< double > getSigma () const override
 
std::vector< double > const & getIntPtSigma (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
std::vector< double > getEpsilon () const override
 
virtual std::vector< double > const & getIntPtEpsilon (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtLiquidDensity (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtGasDensity (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtSolidDensity (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtVapourPressure (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtPorosity (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
std::vector< double > getSaturation () const override
 
virtual std::vector< double > const & getIntPtSaturation (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtMoleFractionGas (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtMassFractionGas (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtMassFractionLiquid (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtRelativePermeabilityGas (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtRelativePermeabilityLiquid (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtEnthalpyGas (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtEnthalpyLiquid (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
virtual std::vector< double > const & getIntPtEnthalpySolid (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 

Private Attributes

TH2MProcessData< DisplacementDim > & _process_data
 
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
 
IntegrationMethod _integration_method
 
MeshLib::Element const & _element
 
bool const _is_axially_symmetric
 
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
 

Member Typedef Documentation

◆ BMatricesType

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

Definition at line 388 of file TH2MFEM.h.

◆ GlobalDimMatrixType

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

Definition at line 65 of file TH2MFEM.h.

◆ GlobalDimVectorType

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

Definition at line 68 of file TH2MFEM.h.

◆ Invariants

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

Definition at line 75 of file TH2MFEM.h.

◆ IpData

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

Definition at line 390 of file TH2MFEM.h.

◆ MatrixType

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

Definition at line 62 of file TH2MFEM.h.

◆ ShapeMatricesTypeDisplacement

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

Definition at line 51 of file TH2MFEM.h.

◆ ShapeMatricesTypePressure

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

Definition at line 54 of file TH2MFEM.h.

◆ SymmetricTensor

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

Definition at line 73 of file TH2MFEM.h.

◆ VectorType

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

Definition at line 58 of file TH2MFEM.h.

Constructor & Destructor Documentation

◆ TH2MLocalAssembler() [1/3]

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

◆ TH2MLocalAssembler() [2/3]

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

◆ TH2MLocalAssembler() [3/3]

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

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

38  : _process_data(process_data),
39  _integration_method(integration_order),
40  _element(e),
41  _is_axially_symmetric(is_axially_symmetric)
42 {
43  unsigned const n_integration_points =
44  _integration_method.getNumberOfPoints();
45 
46  _ip_data.reserve(n_integration_points);
47  _secondary_data.N_u.resize(n_integration_points);
48 
49  auto const shape_matrices_u =
50  NumLib::initShapeMatrices<ShapeFunctionDisplacement,
52  DisplacementDim>(e, is_axially_symmetric,
54 
55  auto const shape_matrices_p =
56  NumLib::initShapeMatrices<ShapeFunctionPressure,
57  ShapeMatricesTypePressure, DisplacementDim>(
58  e, is_axially_symmetric, _integration_method);
59 
60  auto const& solid_material =
62  _process_data.solid_materials, _process_data.material_ids,
63  e.getID());
64 
65  for (unsigned ip = 0; ip < n_integration_points; ip++)
66  {
67  _ip_data.emplace_back(solid_material);
68  auto& ip_data = _ip_data[ip];
69  auto const& sm_u = shape_matrices_u[ip];
70  ip_data.integration_weight =
71  _integration_method.getWeightedPoint(ip).getWeight() *
72  sm_u.integralMeasure * sm_u.detJ;
73 
74  ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
75  DisplacementDim, displacement_size>::Zero(DisplacementDim,
77  for (int i = 0; i < DisplacementDim; ++i)
78  {
79  ip_data.N_u_op
80  .template block<1, displacement_size / DisplacementDim>(
81  i, i * displacement_size / DisplacementDim)
82  .noalias() = sm_u.N;
83  }
84 
85  ip_data.N_u = sm_u.N;
86  ip_data.dNdx_u = sm_u.dNdx;
87 
88  ip_data.N_p = shape_matrices_p[ip].N;
89  ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
90 
91  _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
92  }
93 }
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
Definition: TH2MFEM.h:401
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
Definition: TH2MFEM.h:55
typename ShapeMatricesTypePressure::template MatrixType< M, N > MatrixType
Definition: TH2MFEM.h:63
MeshLib::Element const & _element
Definition: TH2MFEM.h:397
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
Definition: TH2MFEM.h:52
static const int displacement_size
Definition: TH2MFEM.h:412
TH2MProcessData< DisplacementDim > & _process_data
Definition: TH2MFEM.h:386
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
Definition: TH2MFEM.h:394
IntegrationMethod _integration_method
Definition: TH2MFEM.h:396
MechanicsBase< DisplacementDim > & selectSolidConstitutiveRelation(std::map< int, std::unique_ptr< MechanicsBase< DisplacementDim >>> const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u

References ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::_integration_method, ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::_ip_data, ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::_process_data, ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::_secondary_data, ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::displacement_size, MeshLib::Element::getID(), NumLib::initShapeMatrices(), ProcessLib::HydroMechanics::SecondaryData< ShapeMatrixType >::N_u, and MaterialLib::Solids::selectSolidConstitutiveRelation().

Member Function Documentation

◆ assemble()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int DisplacementDim>
void ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::assemble ( double const  t,
double const  dt,
std::vector< double > const &  local_x,
std::vector< double > const &  local_x_dot,
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 772 of file TH2MFEM-impl.h.

778 {
779  auto const matrix_size = gas_pressure_size + capillary_pressure_size +
781  assert(local_x.size() == matrix_size);
782 
783  auto const gas_pressure = Eigen::Map<VectorType<gas_pressure_size> const>(
784  local_x.data() + gas_pressure_index, gas_pressure_size);
785 
786  auto const capillary_pressure =
787  Eigen::Map<VectorType<capillary_pressure_size> const>(
789 
790  auto const capillary_pressure_dot =
791  Eigen::Map<VectorType<capillary_pressure_size> const>(
792  local_x_dot.data() + capillary_pressure_index,
794 
795  // pointer to local_M_data vector
796  auto local_M =
797  MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
798  local_M_data, matrix_size, matrix_size);
799 
800  // pointer to local_K_data vector
801  auto local_K =
802  MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
803  local_K_data, matrix_size, matrix_size);
804 
805  // pointer to local_rhs_data vector
806  auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
807  local_rhs_data, matrix_size);
808 
809  // component-formulation
810  // W - liquid phase main component
811  // C - gas phase main component
812  // pointer-matrices to the mass matrix - C component equation
813  auto MCpG = local_M.template block<C_size, gas_pressure_size>(
815  auto MCpC = local_M.template block<C_size, capillary_pressure_size>(
817  auto MCT = local_M.template block<C_size, temperature_size>(
819  auto MCu = local_M.template block<C_size, displacement_size>(
821 
822  // pointer-matrices to the stiffness matrix - C component equation
823  auto LCpG = local_K.template block<C_size, gas_pressure_size>(
825  auto LCpC = local_K.template block<C_size, capillary_pressure_size>(
827  auto LCT = local_K.template block<C_size, temperature_size>(
829 
830  // pointer-matrices to the mass matrix - W component equation
831  auto MWpG = local_M.template block<W_size, gas_pressure_size>(
833  auto MWpC = local_M.template block<W_size, capillary_pressure_size>(
835  auto MWT = local_M.template block<W_size, temperature_size>(
837  auto MWu = local_M.template block<W_size, displacement_size>(
839 
840  // pointer-matrices to the stiffness matrix - W component equation
841  auto LWpG = local_K.template block<W_size, gas_pressure_size>(
843  auto LWpC = local_K.template block<W_size, capillary_pressure_size>(
845  auto LWT = local_K.template block<W_size, temperature_size>(
847 
848  // pointer-matrices to the mass matrix - temperature equation
849  auto MTu = local_M.template block<temperature_size, displacement_size>(
851 
852  // pointer-matrices to the stiffness matrix - temperature equation
853  auto KTT = local_K.template block<temperature_size, temperature_size>(
855 
856  // pointer-matrices to the stiffness matrix - displacement equation
857  auto KUpG = local_K.template block<displacement_size, gas_pressure_size>(
859  auto KUpC =
860  local_K.template block<displacement_size, capillary_pressure_size>(
862 
863  // pointer-vectors to the right hand side terms - C-component equation
864  auto fC = local_f.template segment<C_size>(C_index);
865  // pointer-vectors to the right hand side terms - W-component equation
866  auto fW = local_f.template segment<W_size>(W_index);
867  // pointer-vectors to the right hand side terms - temperature equation
868  auto fT = local_f.template segment<temperature_size>(temperature_index);
869  // pointer-vectors to the right hand side terms - displacement equation
870  auto fU = local_f.template segment<displacement_size>(displacement_index);
871 
873  pos.setElementID(_element.getID());
874 
875  unsigned const n_integration_points =
876  _integration_method.getNumberOfPoints();
877 
879  Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
880  Eigen::Map<Eigen::VectorXd const>(local_x_dot.data(),
881  local_x_dot.size()),
882  t, dt);
883 
884  for (unsigned int_point = 0; int_point < n_integration_points; int_point++)
885  {
886  pos.setIntegrationPoint(int_point);
887  auto& ip = _ip_data[int_point];
888 
889  auto const& Np = ip.N_p;
890  auto const& NT = Np;
891  auto const& Nu = ip.N_u;
892 
893  auto const& NpT = Np.transpose().eval();
894  auto const& NTT = NT.transpose().eval();
895 
896  auto const& gradNp = ip.dNdx_p;
897  auto const& gradNT = gradNp;
898  auto const& gradNu = ip.dNdx_u;
899 
900  auto const& gradNpT = gradNp.transpose().eval();
901  auto const& gradNTT = gradNT.transpose().eval();
902 
903  auto const& Nu_op = ip.N_u_op;
904  auto const& w = ip.integration_weight;
905 
906  auto const& m = Invariants::identity2;
907 
908  auto const mT = m.transpose().eval();
909 
910  auto const x_coord =
911  NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
913  _element, Nu);
914 
915  auto const Bu =
916  LinearBMatrix::computeBMatrix<DisplacementDim,
917  ShapeFunctionDisplacement::NPOINTS,
918  typename BMatricesType::BMatrixType>(
919  gradNu, Nu, x_coord, _is_axially_symmetric);
920 
921  auto const BuT = Bu.transpose().eval();
922 
923  double const pCap = Np.dot(capillary_pressure);
924 
925  GlobalDimVectorType const gradpGR = gradNp * gas_pressure;
926  GlobalDimVectorType const gradpCap = gradNp * capillary_pressure;
927 
928  double const pCap_dot = Np.dot(capillary_pressure_dot);
929  auto& beta_T_SR = ip.beta_T_SR;
930 
931  auto const I =
932  Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
933 
934  const double sD_G = ip.diffusion_coefficient_vapour;
935  const double sD_L = ip.diffusion_coefficient_solvate;
936 
937  auto const D_C_G = (sD_G * I).eval();
938  auto const D_W_G = (sD_G * I).eval();
939  auto const D_C_L = (sD_L * I).eval();
940  auto const D_W_L = (sD_L * I).eval();
941 
942  auto& k_S = ip.k_S;
943 
944  auto& s_L = ip.s_L;
945  auto const s_G = 1. - s_L;
946  auto const s_L_dot = (s_L - ip.s_L_prev) / dt;
947 
948  auto& alpha_B = ip.alpha_B;
949  auto& beta_p_SR = ip.beta_p_SR;
950 
951  auto const& b = _process_data.specific_body_force;
952 
953  // porosity
954  auto& phi = ip.phi;
955 
956  // volume fraction
957  auto const phi_G = s_G * phi;
958  auto const phi_L = s_L * phi;
959  auto const phi_S = 1. - phi;
960 
961  // solid phase density
962  auto& rho_SR = ip.rhoSR;
963  // effective density
964  auto const rho = phi_G * ip.rhoGR + phi_L * ip.rhoLR + phi_S * rho_SR;
965 
966  // abbreviations
967  const double rho_C_FR = s_G * ip.rhoCGR + s_L * ip.rhoCLR;
968  const double rho_W_FR = s_G * ip.rhoWGR + s_L * ip.rhoWLR;
969 
970  // phase specific enthalpies
971  auto& h_G = ip.h_G;
972  auto& h_L = ip.h_L;
973 
974  auto const rho_C_GR_dot = (ip.rhoCGR - ip.rhoCGR_prev) / dt;
975  auto const rho_C_LR_dot = (ip.rhoCLR - ip.rhoCLR_prev) / dt;
976  auto const rho_W_GR_dot = (ip.rhoWGR - ip.rhoWGR_prev) / dt;
977  auto const rho_W_LR_dot = (ip.rhoWLR - ip.rhoWLR_prev) / dt;
978 
979  auto const rho_h_eff = ip.rho_G_h_G + ip.rho_L_h_L + ip.rho_S_h_S;
980 
981  auto const rho_u_eff_dot = (ip.rho_u_eff - ip.rho_u_eff_prev) / dt;
982 
983  auto const k_over_mu_G = (k_S * ip.k_rel_G / ip.muGR).eval();
984  auto const k_over_mu_L = (k_S * ip.k_rel_L / ip.muLR).eval();
985 
986  GlobalDimVectorType const w_GS =
987  k_over_mu_G * ip.rhoGR * b - k_over_mu_G * gradpGR;
988 
989  GlobalDimVectorType const w_LS = k_over_mu_L * gradpCap +
990  k_over_mu_L * ip.rhoLR * b -
991  k_over_mu_L * gradpGR;
992 
993  // ---------------------------------------------------------------------
994  // C-component equation
995  // ---------------------------------------------------------------------
996 
997  MCpG.noalias() += NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * Np * w;
998  MCpC.noalias() -=
999  NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1000 
1001  if (_process_data.apply_mass_lumping)
1002  {
1003  if (pCap_dot != 0.) // avoid division by Zero
1004  {
1005  MCpC.noalias() +=
1006  NpT *
1007  (phi * (ip.rhoCLR - ip.rhoCGR) -
1008  rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1009  s_L_dot / pCap_dot * Np * w;
1010  }
1011  }
1012 
1013  MCT.noalias() -= NpT * rho_C_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1014  MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
1015 
1016  auto const advection_C_G = (ip.rhoCGR * k_over_mu_G).eval();
1017  auto const advection_C_L = (ip.rhoCLR * k_over_mu_L).eval();
1018  auto const diffusion_C_G_p =
1019  (phi_G * ip.rhoGR * D_C_G * ip.dxmCG_dpGR).eval();
1020  auto const diffusion_C_L_p =
1021  (phi_L * ip.rhoLR * D_C_L * ip.dxmCL_dpLR).eval();
1022  auto const diffusion_C_G_T =
1023  (phi_G * ip.rhoGR * D_C_G * ip.dxmCG_dT).eval();
1024  auto const diffusion_C_L_T =
1025  (phi_L * ip.rhoLR * D_C_L * ip.dxmCL_dT).eval();
1026 
1027  auto const advection_C = (advection_C_G + advection_C_L).eval();
1028  auto const diffusion_C_p = (diffusion_C_G_p + diffusion_C_L_p).eval();
1029  auto const diffusion_C_T = (diffusion_C_G_T + diffusion_C_L_T).eval();
1030 
1031  LCpG.noalias() += gradNpT * (advection_C + diffusion_C_p) * gradNp * w;
1032 
1033  LCpC.noalias() -=
1034  gradNpT * (advection_C_L + diffusion_C_L_p) * gradNp * w;
1035 
1036  LCT.noalias() += gradNpT * (diffusion_C_T)*gradNp * w;
1037 
1038  fC.noalias() += gradNpT *
1039  (advection_C_G * ip.rhoGR + advection_C_L * ip.rhoLR) *
1040  b * w;
1041 
1042  if (!_process_data.apply_mass_lumping)
1043  {
1044  fC.noalias() -= NpT *
1045  (phi * (ip.rhoCLR - ip.rhoCGR) -
1046  rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1047  s_L_dot * w;
1048  }
1049  // fC_III
1050  fC.noalias() -=
1051  NpT * phi * (s_G * rho_C_GR_dot + s_L * rho_C_LR_dot) * w;
1052 
1053  // ---------------------------------------------------------------------
1054  // W-component equation
1055  // ---------------------------------------------------------------------
1056 
1057  MWpG.noalias() += NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * Np * w;
1058  MWpC.noalias() -=
1059  NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1060 
1061  if (_process_data.apply_mass_lumping)
1062  {
1063  if (pCap_dot != 0.) // avoid division by Zero
1064  {
1065  MWpC.noalias() +=
1066  NpT *
1067  (phi * (ip.rhoWLR - ip.rhoWGR) -
1068  rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1069  s_L_dot / pCap_dot * Np * w;
1070  }
1071  }
1072 
1073  MWT.noalias() -= NpT * rho_W_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1074 
1075  MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
1076 
1077  auto const advection_W_G = (ip.rhoWGR * k_over_mu_G).eval();
1078  auto const advection_W_L = (ip.rhoWLR * k_over_mu_L).eval();
1079  auto const diffusion_W_G_p =
1080  (phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dpGR).eval();
1081  auto const diffusion_W_L_p =
1082  (phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dpLR).eval();
1083  auto const diffusion_W_G_T =
1084  (phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dT).eval();
1085  auto const diffusion_W_L_T =
1086  (phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dT).eval();
1087 
1088  auto const advection_W = (advection_W_G + advection_W_L).eval();
1089  auto const diffusion_W_p = (diffusion_W_G_p + diffusion_W_L_p).eval();
1090  auto const diffusion_W_T = (diffusion_W_G_T + diffusion_W_L_T).eval();
1091 
1092  LWpG.noalias() += gradNpT * (advection_W + diffusion_W_p) * gradNp * w;
1093 
1094  LWpC.noalias() -=
1095  gradNpT * (advection_W_L + diffusion_W_L_p) * gradNp * w;
1096 
1097  LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
1098 
1099  fW.noalias() += gradNpT *
1100  (advection_W_G * ip.rhoGR + advection_W_L * ip.rhoLR) *
1101  b * w;
1102 
1103  if (!_process_data.apply_mass_lumping)
1104  {
1105  fW.noalias() -= NpT *
1106  (phi * (ip.rhoWLR - ip.rhoWGR) -
1107  rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1108  s_L_dot * w;
1109  }
1110 
1111  fW.noalias() -=
1112  NpT * phi * (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
1113 
1114  // ---------------------------------------------------------------------
1115  // - temperature equation
1116  // ---------------------------------------------------------------------
1117 
1118  MTu.noalias() += NTT * rho_h_eff * mT * Bu * w;
1119 
1120  KTT.noalias() += gradNTT * ip.lambda * gradNT * w;
1121 
1122  fT.noalias() -= NTT * rho_u_eff_dot * w;
1123 
1124  fT.noalias() +=
1125  gradNTT * (ip.rhoGR * h_G * w_GS + ip.rhoLR * h_L * w_LS) * w;
1126 
1127  fT.noalias() +=
1128  NTT * (ip.rhoGR * w_GS.transpose() + ip.rhoLR * w_LS.transpose()) *
1129  b * w;
1130 
1131  // ---------------------------------------------------------------------
1132  // - displacement equation
1133  // ---------------------------------------------------------------------
1134 
1135  KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
1136 
1137  KUpC.noalias() += (BuT * alpha_B * s_L * m * Np) * w;
1138 
1139  fU.noalias() -= (BuT * ip.sigma_eff - Nu_op.transpose() * rho * b) * w;
1140 
1141  if (_process_data.apply_mass_lumping)
1142  {
1143  MCpG = MCpG.colwise().sum().eval().asDiagonal();
1144  MCpC = MCpC.colwise().sum().eval().asDiagonal();
1145  MWpG = MWpG.colwise().sum().eval().asDiagonal();
1146  MWpC = MWpC.colwise().sum().eval().asDiagonal();
1147  }
1148  } // int_point-loop
1149 }
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:82
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition: BMatrixPolicy.h:55
static const int temperature_index
Definition: TH2MFEM.h:409
static const int gas_pressure_index
Definition: TH2MFEM.h:405
static const int gas_pressure_size
Definition: TH2MFEM.h:406
typename ShapeMatricesTypePressure::GlobalDimVectorType GlobalDimVectorType
Definition: TH2MFEM.h:69
std::vector< ConstitutiveVariables< DisplacementDim > > updateConstitutiveVariables(Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_dot, double const t, double const dt)
Definition: TH2MFEM-impl.h:100
static const int capillary_pressure_index
Definition: TH2MFEM.h:407
static const int temperature_size
Definition: TH2MFEM.h:410
static const int displacement_index
Definition: TH2MFEM.h:411
static const int capillary_pressure_size
Definition: TH2MFEM.h:408
@ rho
density. For some models, rho substitutes p
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
Definition: LinearBMatrix.h:42
static Eigen::Matrix< double, KelvinVectorSize, 1 > const identity2
Kelvin mapping of 2nd order identity tensor.
Definition: KelvinVector.h:76

References MaterialPropertyLib::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), NumLib::interpolateXCoordinate(), ParameterLib::SpatialPosition::setElementID(), and ParameterLib::SpatialPosition::setIntegrationPoint().

◆ assembleWithJacobian()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int DisplacementDim>
void ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::assembleWithJacobian ( double const  t,
double const  dt,
std::vector< double > const &  local_x,
std::vector< double > const &  local_xdot,
const double  ,
const double  ,
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 1156 of file TH2MFEM-impl.h.

1165 {
1166  auto const matrix_size = gas_pressure_size + capillary_pressure_size +
1168  assert(local_x.size() == matrix_size);
1169 
1170  auto const temperature = Eigen::Map<VectorType<temperature_size> const>(
1171  local_x.data() + temperature_index, temperature_size);
1172 
1173  auto const gas_pressure = Eigen::Map<VectorType<gas_pressure_size> const>(
1174  local_x.data() + gas_pressure_index, gas_pressure_size);
1175 
1176  auto const capillary_pressure =
1177  Eigen::Map<VectorType<capillary_pressure_size> const>(
1179 
1180  auto const gas_pressure_dot =
1181  Eigen::Map<VectorType<gas_pressure_size> const>(
1182  local_xdot.data() + gas_pressure_index, gas_pressure_size);
1183 
1184  auto const capillary_pressure_dot =
1185  Eigen::Map<VectorType<capillary_pressure_size> const>(
1186  local_xdot.data() + capillary_pressure_index,
1188 
1189  auto const temperature_dot = Eigen::Map<VectorType<temperature_size> const>(
1190  local_xdot.data() + temperature_index, temperature_size);
1191 
1192  auto const displacement_dot =
1193  Eigen::Map<VectorType<displacement_size> const>(
1194  local_xdot.data() + displacement_index, displacement_size);
1195 
1196  auto local_Jac =
1197  MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1198  local_Jac_data, matrix_size, matrix_size);
1199 
1200  auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
1201  local_rhs_data, matrix_size);
1202 
1203  // component-formulation
1204  // W - liquid phase main component
1205  // C - gas phase main component
1206 
1207  // C component equation matrices
1208  MatrixType<C_size, gas_pressure_size> MCpG =
1209  MatrixType<C_size, gas_pressure_size>::Zero(C_size, gas_pressure_size);
1210  MatrixType<C_size, capillary_pressure_size> MCpC =
1211  MatrixType<C_size, capillary_pressure_size>::Zero(
1213  MatrixType<C_size, temperature_size> MCT =
1214  MatrixType<C_size, temperature_size>::Zero(C_size, temperature_size);
1215  MatrixType<C_size, displacement_size> MCu =
1216  MatrixType<C_size, displacement_size>::Zero(C_size, displacement_size);
1217 
1218  MatrixType<C_size, gas_pressure_size> LCpG =
1219  MatrixType<C_size, gas_pressure_size>::Zero(C_size, gas_pressure_size);
1220  MatrixType<C_size, capillary_pressure_size> LCpC =
1221  MatrixType<C_size, capillary_pressure_size>::Zero(
1223  MatrixType<C_size, temperature_size> LCT =
1224  MatrixType<C_size, temperature_size>::Zero(C_size, temperature_size);
1225 
1226  // mass matrix - W component equation
1227  MatrixType<W_size, gas_pressure_size> MWpG =
1228  MatrixType<W_size, gas_pressure_size>::Zero(W_size, gas_pressure_size);
1229  MatrixType<W_size, capillary_pressure_size> MWpC =
1230  MatrixType<W_size, capillary_pressure_size>::Zero(
1232  MatrixType<W_size, temperature_size> MWT =
1233  MatrixType<W_size, temperature_size>::Zero(W_size, temperature_size);
1234  MatrixType<W_size, displacement_size> MWu =
1235  MatrixType<W_size, displacement_size>::Zero(W_size, displacement_size);
1236 
1237  // stiffness matrix - W component equation
1238  MatrixType<W_size, gas_pressure_size> LWpG =
1239  MatrixType<W_size, gas_pressure_size>::Zero(W_size, gas_pressure_size);
1240  MatrixType<W_size, capillary_pressure_size> LWpC =
1241  MatrixType<W_size, capillary_pressure_size>::Zero(
1243  MatrixType<W_size, temperature_size> LWT =
1244  MatrixType<W_size, temperature_size>::Zero(W_size, temperature_size);
1245 
1246  // mass matrix - temperature equation
1247  MatrixType<temperature_size, displacement_size> MTu =
1248  MatrixType<temperature_size, displacement_size>::Zero(
1250 
1251  // stiffness matrix - temperature equation
1252  MatrixType<temperature_size, temperature_size> KTT =
1253  MatrixType<temperature_size, temperature_size>::Zero(temperature_size,
1255 
1256  // stiffness matrices - displacement equation coupling into pressures
1257  MatrixType<displacement_size, gas_pressure_size> KUpG =
1258  MatrixType<displacement_size, gas_pressure_size>::Zero(
1260  MatrixType<displacement_size, capillary_pressure_size> KUpC =
1261  MatrixType<displacement_size, capillary_pressure_size>::Zero(
1263 
1264  // pointer-vectors to the right hand side terms - C-component equation
1265  auto fC = local_f.template segment<C_size>(C_index);
1266  // pointer-vectors to the right hand side terms - W-component equation
1267  auto fW = local_f.template segment<W_size>(W_index);
1268  // pointer-vectors to the right hand side terms - temperature equation
1269  auto fT = local_f.template segment<temperature_size>(temperature_index);
1270  // pointer-vectors to the right hand side terms - displacement equation
1271  auto fU = local_f.template segment<displacement_size>(displacement_index);
1272 
1274  pos.setElementID(_element.getID());
1275 
1276  unsigned const n_integration_points =
1277  _integration_method.getNumberOfPoints();
1278 
1279  auto const ip_constitutive_variables = updateConstitutiveVariables(
1280  Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1281  Eigen::Map<Eigen::VectorXd const>(local_xdot.data(), local_xdot.size()),
1282  t, dt);
1283 
1284  for (unsigned int_point = 0; int_point < n_integration_points; int_point++)
1285  {
1286  pos.setIntegrationPoint(int_point);
1287  auto& ip = _ip_data[int_point];
1288  auto& ip_cv = ip_constitutive_variables[int_point];
1289 
1290  auto const& Np = ip.N_p;
1291  auto const& NT = Np;
1292  auto const& Nu = ip.N_u;
1293 
1294  auto const& NpT = Np.transpose().eval();
1295  auto const& NTT = NT.transpose().eval();
1296 
1297  auto const& gradNp = ip.dNdx_p;
1298  auto const& gradNT = gradNp;
1299  auto const& gradNu = ip.dNdx_u;
1300 
1301  auto const& gradNpT = gradNp.transpose().eval();
1302  auto const& gradNTT = gradNT.transpose().eval();
1303 
1304  auto const& Nu_op = ip.N_u_op;
1305  auto const& w = ip.integration_weight;
1306 
1307  auto const& m = Invariants::identity2;
1308  auto const mT = m.transpose().eval();
1309 
1310  auto const x_coord =
1311  NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
1313  _element, Nu);
1314 
1315  auto const Bu =
1316  LinearBMatrix::computeBMatrix<DisplacementDim,
1317  ShapeFunctionDisplacement::NPOINTS,
1318  typename BMatricesType::BMatrixType>(
1319  gradNu, Nu, x_coord, _is_axially_symmetric);
1320 
1321  auto const BuT = Bu.transpose().eval();
1322 
1323  double const div_u_dot = Invariants::trace(Bu * displacement_dot);
1324 
1325  double const pCap = Np.dot(capillary_pressure);
1326 
1327  GlobalDimVectorType const gradpGR = gradNp * gas_pressure;
1328  GlobalDimVectorType const gradpCap = gradNp * capillary_pressure;
1329  GlobalDimVectorType const gradT = gradNT * temperature;
1330 
1331  double const pGR_dot = Np.dot(gas_pressure_dot);
1332  double const pCap_dot = Np.dot(capillary_pressure_dot);
1333  double const T_dot = NT.dot(temperature_dot);
1334  auto& beta_T_SR = ip.beta_T_SR;
1335 
1336  auto const I =
1337  Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
1338 
1339  const double sD_G = ip.diffusion_coefficient_vapour;
1340  const double sD_L = ip.diffusion_coefficient_solvate;
1341 
1342  auto const D_C_G = (sD_G * I).eval();
1343  auto const D_W_G = (sD_G * I).eval();
1344  auto const D_C_L = (sD_L * I).eval();
1345  auto const D_W_L = (sD_L * I).eval();
1346 
1347  auto& k_S = ip.k_S;
1348 
1349  auto& s_L = ip.s_L;
1350  auto const s_G = 1. - s_L;
1351  auto const s_L_dot = (s_L - ip.s_L_prev) / dt;
1352 
1353  auto& alpha_B = ip.alpha_B;
1354  auto& beta_p_SR = ip.beta_p_SR;
1355 
1356  auto const& b = _process_data.specific_body_force;
1357 
1358  // porosity
1359  auto& phi = ip.phi;
1360 
1361  // volume fraction
1362  auto const phi_G = s_G * phi;
1363  auto const phi_L = s_L * phi;
1364  auto const phi_S = 1. - phi;
1365 
1366  // solid phase density
1367  auto& rho_SR = ip.rhoSR;
1368  // effective density
1369  auto const rho = phi_G * ip.rhoGR + phi_L * ip.rhoLR + phi_S * rho_SR;
1370 
1371  // abbreviations
1372  const double rho_C_FR = s_G * ip.rhoCGR + s_L * ip.rhoCLR;
1373  const double rho_W_FR = s_G * ip.rhoWGR + s_L * ip.rhoWLR;
1374 
1375  // phase specific enthalpies
1376  auto& h_G = ip.h_G;
1377  auto& h_L = ip.h_L;
1378 
1379  auto const rho_C_GR_dot = (ip.rhoCGR - ip.rhoCGR_prev) / dt;
1380  auto const rho_C_LR_dot = (ip.rhoCLR - ip.rhoCLR_prev) / dt;
1381  auto const rho_W_GR_dot = (ip.rhoWGR - ip.rhoWGR_prev) / dt;
1382  auto const rho_W_LR_dot = (ip.rhoWLR - ip.rhoWLR_prev) / dt;
1383 
1384  auto const rho_h_eff = ip.rho_G_h_G + ip.rho_L_h_L + ip.rho_S_h_S;
1385 
1386  auto const rho_u_eff_dot = (ip.rho_u_eff - ip.rho_u_eff_prev) / dt;
1387 
1388  auto const k_over_mu_G = (k_S * ip.k_rel_G / ip.muGR).eval();
1389  auto const k_over_mu_L = (k_S * ip.k_rel_L / ip.muLR).eval();
1390 
1391  GlobalDimVectorType const w_GS =
1392  k_over_mu_G * ip.rhoGR * b - k_over_mu_G * gradpGR;
1393 
1394  GlobalDimVectorType const w_LS = k_over_mu_L * gradpCap +
1395  k_over_mu_L * ip.rhoLR * b -
1396  k_over_mu_L * gradpGR;
1397 
1398  // ---------------------------------------------------------------------
1399  // C-component equation
1400  // ---------------------------------------------------------------------
1401 
1402  MCpG.noalias() += NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * Np * w;
1403  MCpC.noalias() -=
1404  NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1405 
1406  if (_process_data.apply_mass_lumping)
1407  {
1408  if (pCap_dot != 0.) // avoid division by Zero
1409  {
1410  MCpC.noalias() +=
1411  NpT *
1412  (phi * (ip.rhoCLR - ip.rhoCGR) -
1413  rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1414  s_L_dot / pCap_dot * Np * w;
1415  }
1416  }
1417 
1418  MCT.noalias() -= NpT * rho_C_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1419  // d (fC_4_MCT * T_dot)/d T
1420  local_Jac
1421  .template block<C_size, temperature_size>(C_index,
1423  .noalias() += NpT * ip_cv.dfC_4_MCT_dT * T_dot * NT * w;
1424 
1425  MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
1426  // d (fC_4_MCu * u_dot)/d T
1427  local_Jac
1428  .template block<C_size, temperature_size>(C_index,
1430  .noalias() += NpT * ip_cv.dfC_4_MCu_dT * div_u_dot * NT * w;
1431 
1432  auto const advection_C_G = (ip.rhoCGR * k_over_mu_G).eval();
1433  auto const advection_C_L = (ip.rhoCLR * k_over_mu_L).eval();
1434  auto const diffusion_C_G_p =
1435  (phi_G * ip.rhoGR * D_C_G * ip.dxmCG_dpGR).eval();
1436  auto const diffusion_C_L_p =
1437  (phi_L * ip.rhoLR * D_C_L * ip.dxmCL_dpLR).eval();
1438  auto const diffusion_C_G_T =
1439  (phi_G * ip.rhoGR * D_C_G * ip.dxmCG_dT).eval();
1440  auto const diffusion_C_L_T =
1441  (phi_L * ip.rhoLR * D_C_L * ip.dxmCL_dT).eval();
1442 
1443  auto const advection_C = (advection_C_G + advection_C_L).eval();
1444  auto const diffusion_C_p = (diffusion_C_G_p + diffusion_C_L_p).eval();
1445  auto const diffusion_C_T = (diffusion_C_G_T + diffusion_C_L_T).eval();
1446 
1447  LCpG.noalias() += gradNpT * (advection_C + diffusion_C_p) * gradNp * w;
1448 
1449  // d (fC_4_LCpG * grad p_GR)/d p_GR
1450  local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1451  gradNpT *
1452  (ip_cv.dadvection_C_dp_GR
1453  // + ip_cv.ddiffusion_C_p_dp_GR TODO (naumov)
1454  ) *
1455  gradpGR * Np * w;
1456 
1457  // d (fC_4_LCpG * grad p_GR)/d T
1458  local_Jac
1459  .template block<C_size, temperature_size>(C_index,
1461  .noalias() += gradNpT * ip_cv.dfC_4_LCpG_dT * gradpGR * NT * w;
1462 
1463  // d (fC_4_MCpG * p_GR_dot)/d p_GR
1464  local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1465  NpT * ip_cv.dfC_4_MCpG_dp_GR * pGR_dot * Np * w;
1466 
1467  // d (fC_4_MCpG * p_GR_dot)/d T
1468  local_Jac
1469  .template block<C_size, temperature_size>(C_index,
1471  .noalias() += NpT * ip_cv.dfC_4_MCpG_dT * pGR_dot * NT * w;
1472 
1473  LCpC.noalias() -=
1474  gradNpT * (advection_C_L + diffusion_C_L_p) * gradNp * w;
1475 
1476  LCT.noalias() += gradNpT * diffusion_C_T * gradNp * w;
1477 
1478  // fC_1
1479  fC.noalias() += gradNpT *
1480  (advection_C_G * ip.rhoGR + advection_C_L * ip.rhoLR) *
1481  b * w;
1482 
1483  if (!_process_data.apply_mass_lumping)
1484  {
1485  // fC_2 = \int a * s_L_dot
1486  auto const a = phi * (ip.rhoCLR - ip.rhoCGR) -
1487  rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR;
1488  fC.noalias() -= NpT * a * s_L_dot * w;
1489 
1490  local_Jac.template block<C_size, C_size>(C_index, C_index)
1491  .noalias() +=
1492  NpT *
1493  (ip_cv.dfC_2a_dp_GR * s_L_dot /*- a * (ds_L_dp_GR = 0) / dt*/) *
1494  Np * w;
1495 
1496  local_Jac.template block<C_size, W_size>(C_index, W_index)
1497  .noalias() +=
1498  NpT *
1499  (ip_cv.dfC_2a_dp_cap * s_L_dot + a * ip_cv.ds_L_dp_cap / dt) *
1500  Np * w;
1501 
1502  local_Jac
1503  .template block<C_size, temperature_size>(C_index,
1505  .noalias() += NpT * ip_cv.dfC_2a_dT * s_L_dot * NT * w;
1506  }
1507  {
1508  // fC_3 = \int phi * a
1509  double const a = s_G * rho_C_GR_dot + s_L * rho_C_LR_dot;
1510  fC.noalias() -= NpT * phi * a * w;
1511 
1512  local_Jac.template block<C_size, C_size>(C_index, C_index)
1513  .noalias() += NpT * phi * ip_cv.dfC_3a_dp_GR * Np * w;
1514 
1515  local_Jac.template block<C_size, W_size>(C_index, W_index)
1516  .noalias() += NpT * phi * ip_cv.dfC_3a_dp_cap * Np * w;
1517 
1518  local_Jac
1519  .template block<C_size, temperature_size>(C_index,
1521  .noalias() += NpT *
1522  (
1523 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
1524  ip.dphi_dT * a +
1525 #endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
1526  phi * ip_cv.dfC_3a_dT) *
1527  NT * w;
1528  }
1529  // ---------------------------------------------------------------------
1530  // W-component equation
1531  // ---------------------------------------------------------------------
1532 
1533  MWpG.noalias() += NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * Np * w;
1534  MWpC.noalias() -=
1535  NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1536 
1537  if (_process_data.apply_mass_lumping)
1538  {
1539  if (pCap_dot != 0.) // avoid division by Zero
1540  {
1541  MWpC.noalias() +=
1542  NpT *
1543  (phi * (ip.rhoWLR - ip.rhoWGR) -
1544  rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1545  s_L_dot / pCap_dot * Np * w;
1546  }
1547  }
1548 
1549  MWT.noalias() -= NpT * rho_W_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1550 
1551  MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
1552 
1553  auto const advection_W_G = (ip.rhoWGR * k_over_mu_G).eval();
1554  auto const advection_W_L = (ip.rhoWLR * k_over_mu_L).eval();
1555  auto const diffusion_W_G_p = phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dpGR;
1556  auto const diffusion_W_L_p = phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dpLR;
1557  auto const diffusion_W_G_T = phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dT;
1558  auto const diffusion_W_L_T = phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dT;
1559 
1560  auto const advection_W = advection_W_G + advection_W_L;
1561  auto const diffusion_W_p = diffusion_W_G_p + diffusion_W_L_p;
1562  auto const diffusion_W_T = diffusion_W_G_T + diffusion_W_L_T;
1563 
1564  LWpG.noalias() += gradNpT * (advection_W + diffusion_W_p) * gradNp * w;
1565 
1566  // fW_4 LWpG' parts; LWpG = \int grad (a + d) grad
1567  local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() -=
1568  gradNpT * (ip_cv.dfW_4_LWpG_a_dp_GR + ip_cv.dfW_4_LWpG_d_dp_GR) *
1569  gradpGR * Np * w;
1570 
1571  local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() -=
1572  gradNpT * (ip_cv.dfW_4_LWpG_a_dp_cap + ip_cv.dfW_4_LWpG_d_dp_cap) *
1573  gradpGR * Np * w;
1574 
1575  local_Jac
1576  .template block<W_size, temperature_size>(W_index,
1578  .noalias() -= gradNpT *
1579  (ip_cv.dfW_4_LWpG_a_dT + ip_cv.dfW_4_LWpG_d_dT) *
1580  gradpGR * NT * w;
1581 
1582  LWpC.noalias() -=
1583  gradNpT * (advection_W_L + diffusion_W_L_p) * gradNp * w;
1584 
1585  // fW_4 LWp_cap' parts; LWpC = \int grad (a + d) grad
1586  local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() -=
1587  gradNpT * (ip_cv.dfW_4_LWpC_a_dp_GR + ip_cv.dfW_4_LWpC_d_dp_GR) *
1588  gradpCap * Np * w;
1589 
1590  local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() -=
1591  gradNpT * (ip_cv.dfW_4_LWpC_a_dp_cap + ip_cv.dfW_4_LWpC_d_dp_cap) *
1592  gradpCap * Np * w;
1593 
1594  local_Jac
1595  .template block<W_size, temperature_size>(W_index,
1597  .noalias() -= gradNpT *
1598  (ip_cv.dfW_4_LWpC_a_dT + ip_cv.dfW_4_LWpC_d_dT) *
1599  gradpCap * NT * w;
1600 
1601  LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
1602 
1603  // fW_1
1604  fW.noalias() += gradNpT *
1605  (advection_W_G * ip.rhoGR + advection_W_L * ip.rhoLR) *
1606  b * w;
1607 
1608  // fW_2 = \int (f - g) * s_L_dot
1609  if (!_process_data.apply_mass_lumping)
1610  {
1611  double const f = phi * (ip.rhoWLR - ip.rhoWGR);
1612  double const g = rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR;
1613 
1614  fW.noalias() -= NpT * (f - g) * s_L_dot * w;
1615 
1616  local_Jac.template block<W_size, C_size>(W_index, C_index)
1617  .noalias() += NpT * (ip_cv.dfW_2a_dp_GR - ip_cv.dfW_2b_dp_GR) *
1618  s_L_dot * Np * w;
1619 
1620  // sign negated because of dp_cap = -dp_LR
1621  // TODO (naumov) Had to change the sign to get equal Jacobian WW
1622  // blocks in A2 Test. Where is the error?
1623  local_Jac.template block<W_size, W_size>(W_index, W_index)
1624  .noalias() +=
1625  NpT *
1626  ((ip_cv.dfW_2a_dp_cap - ip_cv.dfW_2b_dp_cap) * s_L_dot +
1627  (f - g) * ip_cv.ds_L_dp_cap / dt) *
1628  Np * w;
1629 
1630  local_Jac
1631  .template block<W_size, temperature_size>(W_index,
1633  .noalias() +=
1634  NpT * (ip_cv.dfW_2a_dT - ip_cv.dfW_2b_dT) * s_L_dot * Np * w;
1635  }
1636 
1637  // fW_3 = \int phi * a
1638  fW.noalias() -=
1639  NpT * phi * (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
1640 
1641  local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1642  NpT * phi * ip_cv.dfW_3a_dp_GR * Np * w;
1643 
1644  local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1645  NpT * phi * ip_cv.dfW_3a_dp_cap * Np * w;
1646 
1647  local_Jac
1648  .template block<W_size, temperature_size>(W_index,
1650  .noalias() +=
1651  NpT *
1652  (
1653 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
1654  ip.dphi_dT * (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) +
1655 #endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
1656  phi * ip_cv.dfW_3a_dT) *
1657  NT * w;
1658 
1659  // ---------------------------------------------------------------------
1660  // - temperature equation
1661  // ---------------------------------------------------------------------
1662 
1663  MTu.noalias() += NTT * rho_h_eff * mT * Bu * w;
1664 
1665  // dfT_4/dp_GR
1666  // d (MTu * u_dot)/dp_GR
1667  local_Jac
1668  .template block<temperature_size, C_size>(temperature_index,
1669  C_index)
1670  .noalias() += NTT * ip_cv.drho_h_eff_dp_GR * div_u_dot * NT * w;
1671 
1672  // dfT_4/dp_cap
1673  // d (MTu * u_dot)/dp_cap
1674  local_Jac
1675  .template block<temperature_size, W_size>(temperature_index,
1676  W_index)
1677  .noalias() -= NTT * ip_cv.drho_h_eff_dp_cap * div_u_dot * NT * w;
1678 
1679  // dfT_4/dT
1680  // d (MTu * u_dot)/dT
1681  local_Jac
1682  .template block<temperature_size, temperature_size>(
1684  .noalias() += NTT * ip_cv.drho_h_eff_dT * div_u_dot * NT * w;
1685 
1686  KTT.noalias() += gradNTT * ip.lambda * gradNT * w;
1687 
1688  // d KTT/dp_GR * T
1689  // TODO (naumov) always zero if lambda_xR have no derivatives wrt. p_GR.
1690  // dlambda_dp_GR =
1691  // (dphi_G_dp_GR = 0) * lambdaGR + phi_G * dlambda_GR_dp_GR +
1692  // (dphi_L_dp_GR = 0) * lambdaLR + phi_L * dlambda_LR_dp_GR +
1693  // (dphi_S_dp_GR = 0) * lambdaSR + phi_S * dlambda_SR_dp_GR +
1694  // = 0
1695  //
1696  // Since dlambda/dp_GR is 0 the derivative is omitted:
1697  // local_Jac
1698  // .template block<temperature_size, C_size>(temperature_index,
1699  // C_index)
1700  // .noalias() += gradNTT * dlambda_dp_GR * gradT * Np * w;
1701 
1702  // d KTT/dp_cap * T
1703  local_Jac
1704  .template block<temperature_size, W_size>(temperature_index,
1705  W_index)
1706  .noalias() += gradNTT * ip_cv.dlambda_dp_cap * gradT * Np * w;
1707 
1708  // d KTT/dT * T
1709  local_Jac
1710  .template block<temperature_size, temperature_size>(
1712  .noalias() += gradNTT * ip_cv.dlambda_dT * gradT * NT * w;
1713 
1714  // fT_1
1715  fT.noalias() -= NTT * rho_u_eff_dot * w;
1716 
1717  // dfT_1/dp_GR
1718  local_Jac
1719  .template block<temperature_size, C_size>(temperature_index,
1720  C_index)
1721  .noalias() += NpT / dt * ip_cv.drho_u_eff_dp_GR * Np * w;
1722 
1723  // dfT_1/dp_cap
1724  local_Jac
1725  .template block<temperature_size, W_size>(temperature_index,
1726  W_index)
1727  .noalias() += NpT / dt * ip_cv.drho_u_eff_dp_cap * Np * w;
1728 
1729  // dfT_1/dT
1730  // MTT
1731  local_Jac
1732  .template block<temperature_size, temperature_size>(
1734  .noalias() += NTT * ip_cv.drho_u_eff_dT / dt * NT * w;
1735 
1736  // fT_2
1737  fT.noalias() +=
1738  gradNTT * (ip.rhoGR * h_G * w_GS + ip.rhoLR * h_L * w_LS) * w;
1739 
1740  // dfT_2/dp_GR
1741  local_Jac
1742  .template block<temperature_size, C_size>(temperature_index,
1743  C_index)
1744  .noalias() -=
1745  // dfT_2/dp_GR first part
1746  gradNTT * ip_cv.drho_GR_h_w_eff_dp_GR_Npart * Np * w +
1747  // dfT_2/dp_GR second part
1748  gradNTT * ip_cv.drho_GR_h_w_eff_dp_GR_gradNpart * gradNp * w;
1749 
1750  // dfT_2/dp_cap
1751  local_Jac
1752  .template block<temperature_size, W_size>(temperature_index,
1753  W_index)
1754  .noalias() -=
1755  // first part of dfT_2/dp_cap
1756  gradNTT * (-ip_cv.drho_LR_h_w_eff_dp_cap_Npart) * Np * w +
1757  // second part of dfT_2/dp_cap
1758  gradNTT * (-ip_cv.drho_LR_h_w_eff_dp_cap_gradNpart) * gradNp * w;
1759 
1760  // dfT_2/dT
1761  local_Jac
1762  .template block<temperature_size, temperature_size>(
1764  .noalias() -= gradNTT * ip_cv.drho_GR_h_w_eff_dT * NT * w;
1765 
1766  // fT_3
1767  fT.noalias() +=
1768  NTT * (ip.rhoGR * w_GS.transpose() + ip.rhoLR * w_LS.transpose()) *
1769  b * w;
1770 
1771  // ---------------------------------------------------------------------
1772  // - displacement equation
1773  // ---------------------------------------------------------------------
1774 
1775  KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
1776 
1777  // dfU_2/dp_GR part i.e. part of the d(KUpG*p_GR)/dp_GR derivative is
1778  // dKUpG/dp_GR + KUpG. The former is zero, the latter is handled below.
1779 
1780  KUpC.noalias() += (BuT * alpha_B * s_L * m * Np) * w;
1781 
1782  // dfU_2/dp_LR part i.e. part of the d(KUpC*p_cap)/dp_LR derivative is
1783  // dKUpC/dp_LR + KUpC. The latter is handled below, the former here:
1784  local_Jac
1785  .template block<displacement_size, W_size>(displacement_index,
1786  W_index)
1787  .noalias() += BuT * alpha_B * ip_cv.ds_L_dp_cap * pCap * m * Np * w;
1788 
1789  local_Jac
1790  .template block<displacement_size, displacement_size>(
1792  .noalias() += BuT * ip_cv.C * Bu * w;
1793 
1794  // fU_1
1795  fU.noalias() -= (BuT * ip.sigma_eff - Nu_op.transpose() * rho * b) * w;
1796 
1797  // KuT
1798  local_Jac
1799  .template block<displacement_size, temperature_size>(
1801  .noalias() -= BuT * (ip_cv.C * ip.alpha_T_SR) * NT * w;
1802 
1803  /* TODO (naumov) Test with gravity needed to check this Jacobian part.
1804  local_Jac
1805  .template block<displacement_size, temperature_size>(
1806  displacement_index, temperature_index)
1807  .noalias() += Nu_op.transpose() * ip_cv.drho_dT * b * Nu_op * w;
1808  */
1809 
1810  if (_process_data.apply_mass_lumping)
1811  {
1812  MCpG = MCpG.colwise().sum().eval().asDiagonal();
1813  MCpC = MCpC.colwise().sum().eval().asDiagonal();
1814  MWpG = MWpG.colwise().sum().eval().asDiagonal();
1815  MWpC = MWpC.colwise().sum().eval().asDiagonal();
1816  }
1817  } // int_point-loop
1818 
1819  // --- Gas ---
1820  // fC_4
1821  fC.noalias() -= LCpG * gas_pressure + LCpC * capillary_pressure +
1822  LCT * temperature + MCpG * gas_pressure_dot +
1823  MCpC * capillary_pressure_dot + MCT * temperature_dot +
1824  MCu * displacement_dot;
1825 
1826  local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1827  LCpG + MCpG / dt;
1828  local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
1829  LCpC + MCpC / dt;
1830  local_Jac
1831  .template block<C_size, temperature_size>(C_index, temperature_index)
1832  .noalias() += LCT + MCT / dt;
1833  local_Jac
1834  .template block<C_size, displacement_size>(C_index, displacement_index)
1835  .noalias() += MCu / dt;
1836 
1837  // --- Capillary pressure ---
1838  // fW_4
1839  fW.noalias() -= LWpG * gas_pressure + LWpC * capillary_pressure +
1840  LWT * temperature + MWpG * gas_pressure_dot +
1841  MWpC * capillary_pressure_dot + MWT * temperature_dot +
1842  MWu * displacement_dot;
1843 
1844  local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1845  LWpC + MWpC / dt;
1846  local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1847  LWpG + MWpG / dt;
1848  local_Jac
1849  .template block<W_size, temperature_size>(W_index, temperature_index)
1850  .noalias() += LWT + MWT / dt;
1851  local_Jac
1852  .template block<W_size, displacement_size>(W_index, displacement_index)
1853  .noalias() += MWu / dt;
1854 
1855  // --- Temperature ---
1856  // fT_4
1857  fT.noalias() -= KTT * temperature + MTu * displacement_dot;
1858 
1859  local_Jac
1860  .template block<temperature_size, temperature_size>(temperature_index,
1862  .noalias() += KTT;
1863  local_Jac
1864  .template block<temperature_size, displacement_size>(temperature_index,
1866  .noalias() += MTu / dt;
1867 
1868  // --- Displacement ---
1869  // fU_2
1870  fU.noalias() -= KUpG * gas_pressure + KUpC * capillary_pressure;
1871 
1872  local_Jac
1873  .template block<displacement_size, C_size>(displacement_index, C_index)
1874  .noalias() += KUpG;
1875  local_Jac
1876  .template block<displacement_size, W_size>(displacement_index, W_index)
1877  .noalias() += KUpC;
1878 }
static double trace(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Trace of the corresponding tensor.

References MaterialPropertyLib::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), NumLib::interpolateXCoordinate(), ParameterLib::SpatialPosition::setElementID(), and ParameterLib::SpatialPosition::setIntegrationPoint().

◆ computeSecondaryVariableConcrete()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

2095 {
2096  auto const gas_pressure =
2097  local_x.template segment<gas_pressure_size>(gas_pressure_index);
2098  auto const capillary_pressure =
2099  local_x.template segment<capillary_pressure_size>(
2101  auto const liquid_pressure = (gas_pressure - capillary_pressure).eval();
2102 
2104  ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2105  DisplacementDim>(_element, _is_axially_symmetric, gas_pressure,
2106  *_process_data.gas_pressure_interpolated);
2107 
2109  ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2111  *_process_data.capillary_pressure_interpolated);
2112 
2114  ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2115  DisplacementDim>(_element, _is_axially_symmetric, liquid_pressure,
2116  *_process_data.liquid_pressure_interpolated);
2117 
2118  auto const temperature =
2119  local_x.template segment<temperature_size>(temperature_index);
2120 
2122  ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2123  DisplacementDim>(_element, _is_axially_symmetric, temperature,
2124  *_process_data.temperature_interpolated);
2125 
2126  unsigned const n_integration_points =
2127  _integration_method.getNumberOfPoints();
2128 
2129  double saturation_avg = 0;
2130 
2131  updateConstitutiveVariables(local_x, local_x_dot, t, dt);
2132 
2133  for (unsigned ip = 0; ip < n_integration_points; ip++)
2134  {
2135  saturation_avg += _ip_data[ip].s_L;
2136  }
2137  saturation_avg /= n_integration_points;
2138  (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
2139 }
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)

References MaterialPropertyLib::capillary_pressure, and NumLib::interpolateToHigherOrderNodes().

◆ getEpsilon()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int DisplacementDim>
std::vector<double> ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::getEpsilon ( ) const
inlineoverrideprivatevirtual

Implements ProcessLib::TH2M::LocalAssemblerInterface.

Definition at line 219 of file TH2MFEM.h.

220  {
221  constexpr int kelvin_vector_size =
223 
224  return transposeInPlace<kelvin_vector_size>(
225  [this](std::vector<double>& values)
226  { return getIntPtEpsilon(0, {}, {}, values); });
227  }
virtual std::vector< double > const & getIntPtEpsilon(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
Definition: TH2MFEM.h:229
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:23

References ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::getIntPtEpsilon(), and MathLib::KelvinVector::kelvin_vector_dimensions().

◆ getIntPtDarcyVelocityGas()

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

Implements ProcessLib::TH2M::LocalAssemblerInterface.

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

1890 {
1891  auto const num_intpts = _ip_data.size();
1892 
1893  constexpr int process_id = 0; // monolithic scheme;
1894  auto const indices =
1895  NumLib::getIndices(_element.getID(), *dof_table[process_id]);
1896  assert(!indices.empty());
1897  auto const local_x = x[process_id]->get(indices);
1898 
1899  cache.clear();
1900  auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
1901  double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
1902  cache, DisplacementDim, num_intpts);
1903 
1904  auto const pGR =
1905  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
1906  gas_pressure_size> const>(local_x.data() + gas_pressure_index,
1908  auto const pCap =
1909  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
1910  capillary_pressure_size> const>(
1912  auto const T =
1913  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
1914  temperature_size> const>(local_x.data() + temperature_index,
1916 
1917  unsigned const n_integration_points =
1918  _integration_method.getNumberOfPoints();
1919 
1921  pos.setElementID(_element.getID());
1922 
1923  auto const& medium = *_process_data.media_map->getMedium(_element.getID());
1924  auto const& gas_phase = medium.phase("Gas");
1925 
1926  MPL::VariableArray vars;
1927 
1928  for (unsigned ip = 0; ip < n_integration_points; ip++)
1929  {
1930  pos.setIntegrationPoint(ip);
1931 
1932  auto const& N_p = _ip_data[ip].N_p;
1933 
1934  vars[static_cast<int>(MPL::Variable::temperature)] =
1935  N_p.dot(T); // N_p = N_T
1936  vars[static_cast<int>(MPL::Variable::phase_pressure)] = N_p.dot(pGR);
1937  vars[static_cast<int>(MPL::Variable::capillary_pressure)] =
1938  N_p.dot(pCap);
1939 
1940  // TODO (naumov) Temporary value not used by current material
1941  // models. Need extension of secondary variables interface.
1942  double const dt = std::numeric_limits<double>::quiet_NaN();
1943 
1944  auto const mu_GR = gas_phase.property(MPL::PropertyType::viscosity)
1945  .template value<double>(vars, pos, t, dt);
1946 
1947  GlobalDimMatrixType k_S = MPL::formEigenTensor<DisplacementDim>(
1948  medium.property(MPL::PropertyType::permeability)
1949  .value(vars, pos, t, dt));
1950 
1951  auto const s_L = medium.property(MPL::PropertyType::saturation)
1952  .template value<double>(vars, pos, t, dt);
1953 
1954  vars[static_cast<int>(MPL::Variable::liquid_saturation)] = s_L;
1955 
1956  auto const k_rel =
1957  medium
1958  .property(
1960  .template value<double>(vars, pos, t, dt);
1961 
1962  auto const k_over_mu = k_S * k_rel / mu_GR;
1963 
1964  vars[static_cast<int>(MPL::Variable::molar_mass)] = 0.1;
1965  auto const rho_GR = gas_phase.property(MPL::PropertyType::density)
1966  .template value<double>(vars, pos, t, dt);
1967  auto const& b = _process_data.specific_body_force;
1968 
1969  // Compute the velocity
1970  auto const& dNdx_p = _ip_data[ip].dNdx_p;
1971  cache_matrix.col(ip).noalias() =
1972  -k_over_mu * dNdx_p * pGR + k_over_mu * rho_GR * b;
1973  }
1974 
1975  return cache;
1976 }
typename ShapeMatricesTypePressure::GlobalDimMatrixType GlobalDimMatrixType
Definition: TH2MFEM.h:66
typename ShapeMatricesTypePressure::template VectorType< N > VectorType
Definition: TH2MFEM.h:59
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
Definition: VariableType.h:108
@ relative_permeability_nonwetting_phase
Definition: PropertyType.h:77
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:32
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)

References MaterialPropertyLib::capillary_pressure, MathLib::createZeroedMatrix(), MaterialPropertyLib::density, NumLib::getIndices(), MaterialPropertyLib::molar_mass, MaterialPropertyLib::permeability, MaterialPropertyLib::relative_permeability_nonwetting_phase, MaterialPropertyLib::saturation, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), MaterialPropertyLib::temperature, and MaterialPropertyLib::viscosity.

◆ getIntPtDarcyVelocityLiquid()

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

Implements ProcessLib::TH2M::LocalAssemblerInterface.

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

1988 {
1989  auto const num_intpts = _ip_data.size();
1990 
1991  constexpr int process_id = 0; // monolithic scheme;
1992  auto const indices =
1993  NumLib::getIndices(_element.getID(), *dof_table[process_id]);
1994  assert(!indices.empty());
1995  auto const local_x = x[process_id]->get(indices);
1996 
1997  cache.clear();
1998  auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
1999  double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2000  cache, DisplacementDim, num_intpts);
2001 
2002  auto const pGR =
2003  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
2004  gas_pressure_size> const>(local_x.data() + gas_pressure_index,
2006  auto const pCap =
2007  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
2008  capillary_pressure_size> const>(
2010  auto const pLR = pGR - pCap;
2011  auto const T =
2012  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
2013  temperature_size> const>(local_x.data() + temperature_index,
2015 
2016  unsigned const n_integration_points =
2017  _integration_method.getNumberOfPoints();
2018 
2020  pos.setElementID(_element.getID());
2021 
2022  auto const& medium = *_process_data.media_map->getMedium(_element.getID());
2023  auto const& liquid_phase = medium.phase("AqueousLiquid");
2024 
2025  MPL::VariableArray vars;
2026 
2027  for (unsigned ip = 0; ip < n_integration_points; ip++)
2028  {
2029  pos.setIntegrationPoint(ip);
2030 
2031  auto const& N_p = _ip_data[ip].N_p;
2032 
2033  vars[static_cast<int>(MPL::Variable::temperature)] = N_p.dot(T);
2034  vars[static_cast<int>(MPL::Variable::phase_pressure)] = N_p.dot(pGR);
2035  vars[static_cast<int>(MPL::Variable::liquid_phase_pressure)] =
2036  N_p.dot(pLR);
2037  vars[static_cast<int>(MPL::Variable::capillary_pressure)] =
2038  N_p.dot(pCap);
2039 
2040  // TODO (naumov) Temporary value not used by current material
2041  // models. Need extension of secondary variables interface.
2042  double const dt = std::numeric_limits<double>::quiet_NaN();
2043 
2044  auto const mu_LR = liquid_phase.property(MPL::PropertyType::viscosity)
2045  .template value<double>(vars, pos, t, dt);
2046  GlobalDimMatrixType k_S = MPL::formEigenTensor<DisplacementDim>(
2047  medium.property(MPL::PropertyType::permeability)
2048  .value(vars, pos, t, dt));
2049 
2050  auto const s_L = medium.property(MPL::PropertyType::saturation)
2051  .template value<double>(vars, pos, t, dt);
2052 
2053  vars[static_cast<int>(MPL::Variable::liquid_saturation)] = s_L;
2054 
2055  auto const k_rel =
2057  .template value<double>(vars, pos, t, dt);
2058 
2059  auto const k_over_mu = k_S * k_rel / mu_LR;
2060 
2061  vars[static_cast<int>(MPL::Variable::molar_fraction)] = 1.0;
2062 
2063  auto const cCL = [&]()
2064  {
2065  if (liquid_phase.hasProperty(MPL::PropertyType::concentration))
2066  {
2067  return liquid_phase.property(MPL::PropertyType::concentration)
2068  .template value<double>(vars, pos, t, dt); // in mol*m^(-3)
2069  }
2070  return 0.;
2071  }();
2072 
2073  vars[static_cast<int>(MPL::Variable::concentration)] = cCL;
2074 
2075  auto const rho_LR = liquid_phase.property(MPL::PropertyType::density)
2076  .template value<double>(vars, pos, t, dt);
2077  auto const& b = _process_data.specific_body_force;
2078 
2079  // Compute the velocity
2080  auto const& dNdx_p = _ip_data[ip].dNdx_p;
2081  cache_matrix.col(ip).noalias() =
2082  -k_over_mu * dNdx_p * pLR + k_over_mu * rho_LR * b;
2083  }
2084 
2085  return cache;
2086 }
@ concentration
used to specify decay rate of a substance.
Definition: PropertyType.h:46

References MaterialPropertyLib::capillary_pressure, MaterialPropertyLib::concentration, MathLib::createZeroedMatrix(), MaterialPropertyLib::density, NumLib::getIndices(), MaterialPropertyLib::permeability, MaterialPropertyLib::relative_permeability, MaterialPropertyLib::saturation, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), MaterialPropertyLib::temperature, and MaterialPropertyLib::viscosity.

◆ getIntPtEnthalpyGas()

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

◆ getIntPtEnthalpyLiquid()

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

◆ getIntPtEnthalpySolid()

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

◆ getIntPtEpsilon()

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

◆ getIntPtGasDensity()

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

◆ getIntPtLiquidDensity()

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

◆ getIntPtMassFractionGas()

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

◆ getIntPtMassFractionLiquid()

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

◆ getIntPtMoleFractionGas()

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

◆ getIntPtPorosity()

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

◆ getIntPtRelativePermeabilityGas()

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

◆ getIntPtRelativePermeabilityLiquid()

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

◆ getIntPtSaturation()

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

◆ getIntPtSigma()

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

◆ getIntPtSolidDensity()

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

◆ getIntPtVapourPressure()

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

◆ getSaturation()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int DisplacementDim>
std::vector<double> ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::getSaturation ( ) const
inlineoverrideprivatevirtual

Implements ProcessLib::TH2M::LocalAssemblerInterface.

Definition at line 293 of file TH2MFEM.h.

294  {
295  std::vector<double> result;
296  getIntPtSaturation(0, {}, {}, result);
297  return result;
298  }
virtual std::vector< double > const & getIntPtSaturation(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
Definition: TH2MFEM.h:300

References ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::getIntPtSaturation().

◆ getShapeMatrix()

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

Provides the shape matrix at the given integration point.

Implements NumLib::ExtrapolatableElement.

Definition at line 162 of file TH2MFEM.h.

164  {
165  auto const& N_u = _secondary_data.N_u[integration_point];
166 
167  // assumes N is stored contiguously in memory
168  return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
169  }

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

◆ getSigma()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int DisplacementDim>
std::vector<double> ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::getSigma ( ) const
inlineoverrideprivatevirtual

Implements ProcessLib::TH2M::LocalAssemblerInterface.

Definition at line 196 of file TH2MFEM.h.

197  {
198  {
199  constexpr int kelvin_vector_size =
201  DisplacementDim);
202 
203  return transposeInPlace<kelvin_vector_size>(
204  [this](std::vector<double>& values)
205  { return getIntPtSigma(0, {}, {}, values); });
206  }
207  }
std::vector< double > const & getIntPtSigma(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
Definition: TH2MFEM.h:209

References ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::getIntPtSigma(), and MathLib::KelvinVector::kelvin_vector_dimensions().

◆ initializeConcrete()

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

Set initial stress from parameter.

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 114 of file TH2MFEM.h.

115  {
116  unsigned const n_integration_points =
117  _integration_method.getNumberOfPoints();
118 
119  for (unsigned ip = 0; ip < n_integration_points; ip++)
120  {
121  auto& ip_data = _ip_data[ip];
122 
124  if (_process_data.initial_stress != nullptr)
125  {
126  ParameterLib::SpatialPosition const x_position{
127  std::nullopt, _element.getID(), ip,
129  ShapeFunctionDisplacement,
131  _element, ip_data.N_u))};
132 
133  ip_data.sigma_eff =
135  DisplacementDim>((*_process_data.initial_stress)(
136  std::numeric_limits<
137  double>::quiet_NaN() /* time independent */,
138  x_position));
139  }
140 
141  ip_data.pushBackState();
142  }
143  }
Eigen::Matrix< double, Eigen::MatrixBase< Derived >::RowsAtCompileTime, 1 > symmetricTensorToKelvinVector(Eigen::MatrixBase< Derived > const &v)
Definition: KelvinVector.h:170
MathLib::TemplatePoint< double, 3 > Point3d
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)

References ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::_element, ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::_integration_method, ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::_ip_data, ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::_process_data, MeshLib::Element::getID(), NumLib::interpolateCoordinates(), and MathLib::KelvinVector::symmetricTensorToKelvinVector().

◆ postTimestepConcrete()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 145 of file TH2MFEM.h.

148  {
149  unsigned const n_integration_points =
150  _integration_method.getNumberOfPoints();
151 
152  for (unsigned ip = 0; ip < n_integration_points; ip++)
153  {
154  _ip_data[ip].pushBackState();
155  }
156  }

References ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::_integration_method, and ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::_ip_data.

◆ setInitialConditionsConcrete()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

748 {
749  [[maybe_unused]] auto const matrix_size =
752 
753  assert(local_x.size() == matrix_size);
754 
756  Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
757  Eigen::VectorXd::Zero(matrix_size), t, 0);
758 
759  unsigned const n_integration_points =
760  _integration_method.getNumberOfPoints();
761  for (unsigned ip = 0; ip < n_integration_points; ip++)
762  {
763  auto& ip_data = _ip_data[ip];
764  ip_data.pushBackState();
765  }
766 }

◆ setIPDataInitialConditions()

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

Implements ProcessLib::TH2M::LocalAssemblerInterface.

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

702 {
703  if (integration_order !=
704  static_cast<int>(_integration_method.getIntegrationOrder()))
705  {
706  OGS_FATAL(
707  "Setting integration point initial conditions; The integration "
708  "order of the local assembler for element {:d} is different "
709  "from the integration order in the initial condition.",
710  _element.getID());
711  }
712 
713  if (name == "sigma_ip")
714  {
715  if (_process_data.initial_stress != nullptr)
716  {
717  OGS_FATAL(
718  "Setting initial conditions for stress from integration "
719  "point data and from a parameter '{:s}' is not possible "
720  "simultaneously.",
721  _process_data.initial_stress->name);
722  }
723  return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
724  values, _ip_data, &IpData::sigma_eff);
725  }
726 
727  if (name == "saturation_ip")
728  {
730  &IpData::s_L);
731  }
732  if (name == "epsilon_ip")
733  {
734  return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
735  values, _ip_data, &IpData::eps);
736  }
737  return 0;
738 }
#define OGS_FATAL(...)
Definition: Error.h:26
std::size_t setIntegrationPointScalarData(double const *values, std::vector< IntegrationPointData, Eigen::aligned_allocator< IntegrationPointData >> &ip_data, MemberType member)

References MaterialPropertyLib::name, OGS_FATAL, and ProcessLib::setIntegrationPointScalarData().

◆ setSigma()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int DisplacementDim>
std::size_t ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::setSigma ( double const *  values)
inlineprivate

◆ updateConstitutiveVariables()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int DisplacementDim>
std::vector< ConstitutiveVariables< DisplacementDim > > ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::updateConstitutiveVariables ( Eigen::VectorXd const &  local_x,
Eigen::VectorXd const &  local_x_dot,
double const  t,
double const  dt 
)
private

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

103 {
104  [[maybe_unused]] auto const matrix_size =
107 
108  assert(local_x.size() == matrix_size);
109 
110  auto const gas_pressure =
111  local_x.template segment<gas_pressure_size>(gas_pressure_index);
112  auto const capillary_pressure =
113  local_x.template segment<capillary_pressure_size>(
115 
116  auto const temperature =
117  local_x.template segment<temperature_size>(temperature_index);
118  auto const temperature_dot =
119  local_x_dot.template segment<temperature_size>(temperature_index);
120 
121  auto const displacement =
122  local_x.template segment<displacement_size>(displacement_index);
123 
125  pos.setElementID(_element.getID());
126 
127  auto const& medium = *_process_data.media_map->getMedium(_element.getID());
128  auto const& gas_phase = medium.phase("Gas");
129  auto const& liquid_phase = medium.phase("AqueousLiquid");
130  auto const& solid_phase = medium.phase("Solid");
131 
132  unsigned const n_integration_points =
133  _integration_method.getNumberOfPoints();
134 
135  std::vector<ConstitutiveVariables<DisplacementDim>>
136  ip_constitutive_variables(n_integration_points);
137 
138  for (unsigned ip = 0; ip < n_integration_points; ip++)
139  {
140  pos.setIntegrationPoint(ip);
141  auto& ip_data = _ip_data[ip];
142  auto& ip_cv = ip_constitutive_variables[ip];
143 
144  auto const& Np = ip_data.N_p;
145  auto const& NT = Np;
146  auto const& Nu = ip_data.N_u;
147  auto const& gradNu = ip_data.dNdx_u;
148  auto const& gradNp = ip_data.dNdx_p;
149  auto const x_coord =
150  NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
152  _element, Nu);
153 
154  double const T = NT.dot(temperature);
155  double const pGR = Np.dot(gas_pressure);
156  double const pCap = Np.dot(capillary_pressure);
157  double const pLR = pGR - pCap;
158  GlobalDimVectorType const gradpGR = gradNp * gas_pressure;
159  GlobalDimVectorType const gradpCap = gradNp * capillary_pressure;
160 
161  MPL::VariableArray vars;
162  vars[static_cast<int>(MPL::Variable::temperature)] = T;
163  vars[static_cast<int>(MPL::Variable::phase_pressure)] = pGR;
164  vars[static_cast<int>(MPL::Variable::capillary_pressure)] = pCap;
165  vars[static_cast<int>(MPL::Variable::liquid_phase_pressure)] = pLR;
166 
167  // medium properties
168  auto const K_S = ip_data.solid_material.getBulkModulus(t, pos);
169 
170  ip_data.alpha_B = medium.property(MPL::PropertyType::biot_coefficient)
171  .template value<double>(vars, pos, t, dt);
172 
173  ip_data.s_L =
174  medium.property(MPL::PropertyType::saturation)
175  .template value<double>(
176  vars, pos, t, std::numeric_limits<double>::quiet_NaN());
177 
178  vars[static_cast<int>(MPL::Variable::liquid_saturation)] = ip_data.s_L;
179 
180  // intrinsic permeability
181  ip_data.k_S = MPL::formEigenTensor<DisplacementDim>(
182  medium.property(MPL::PropertyType::permeability)
183  .value(vars, pos, t, dt));
184 
185  auto const Bu =
186  LinearBMatrix::computeBMatrix<DisplacementDim,
187  ShapeFunctionDisplacement::NPOINTS,
188  typename BMatricesType::BMatrixType>(
189  gradNu, Nu, x_coord, _is_axially_symmetric);
190 
191  auto& eps = ip_data.eps;
192  eps.noalias() = Bu * displacement;
193 
194  // relative permeability
195  // Set mechanical variables for the intrinsic permeability model
196  // For stress dependent permeability.
197  {
198  // Note: if Bishop model is available, ip_data.s_L in the following
199  // computation should be replaced with the Bishop value.
200  auto const sigma_total =
201  (_ip_data[ip].sigma_eff - ip_data.alpha_B *
202  (pGR - ip_data.s_L * pCap) *
204  .eval();
205 
206  vars[static_cast<int>(MPL::Variable::total_stress)]
207  .emplace<SymmetricTensor>(
209  sigma_total));
210  }
211  // For strain dependent permeability
212  vars[static_cast<int>(MPL::Variable::volumetric_strain)] =
213  Invariants::trace(eps);
214  vars[static_cast<int>(MPL::Variable::equivalent_plastic_strain)] =
215  _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
216  vars[static_cast<int>(MPL::Variable::mechanical_strain)]
218  eps);
219 
220  ip_data.k_rel_G =
221  medium
222  .property(
224  .template value<double>(vars, pos, t, dt);
225 
226  auto const dk_rel_G_ds_L =
228  .template dValue<double>(vars, MPL::Variable::liquid_saturation,
229  pos, t, dt);
230 
231  ip_data.k_rel_L =
233  .template value<double>(vars, pos, t, dt);
234 
235  auto const dk_rel_L_ds_L =
237  .template dValue<double>(vars, MPL::Variable::liquid_saturation,
238  pos, t, dt);
239 
240  // solid phase compressibility
241  ip_data.beta_p_SR = (1. - ip_data.alpha_B) / K_S;
242 
243  // solid phase linear thermal expansion coefficient
244  ip_data.alpha_T_SR = MathLib::KelvinVector::tensorToKelvin<
246  solid_phase
247  .property(
249  .value(vars, pos, t, dt)));
250 
251  // isotropic solid phase volumetric thermal expansion coefficient
252  ip_data.beta_T_SR = Invariants::trace(ip_data.alpha_T_SR);
253 
254  double const T_dot = NT.dot(temperature_dot);
256  dthermal_strain = ip_data.alpha_T_SR * T_dot * dt;
257 
258  auto& eps_prev = ip_data.eps_prev;
259  auto& eps_m = ip_data.eps_m;
260  auto& eps_m_prev = ip_data.eps_m_prev;
261 
262  eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
265  eps_m);
266 
267  auto const rho_ref_SR =
268  solid_phase.property(MPL::PropertyType::density)
269  .template value<double>(
270  vars, pos, t, std::numeric_limits<double>::quiet_NaN());
271 
272  auto const lambdaSR = MPL::formEigenTensor<DisplacementDim>(
273  solid_phase.property(MPL::PropertyType::thermal_conductivity)
274  .value(vars, pos, t, dt));
275 
276  double const T0 = _process_data.reference_temperature(t, pos)[0];
277  double const delta_T(T - T0);
278  ip_data.thermal_volume_strain = ip_data.beta_T_SR * delta_T;
279 
280  // initial porosity
281  auto const phi_0 = medium.property(MPL::PropertyType::porosity)
282  .template value<double>(vars, pos, t, dt);
283 
284  auto const phi_S_0 = 1. - phi_0;
285 
286 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
287  auto const& m = Invariants::identity2;
288  double const div_u = m.transpose() * eps;
289 
290  const double phi_S = phi_S_0 * (1. + ip_data.thermal_volume_strain -
291  ip_data.alpha_B * div_u);
292 #else // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
293  const double phi_S = phi_S_0;
294 #endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
295 
296  // porosity
297  ip_data.phi = 1. - phi_S;
298 
299  // solid phase density
300 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
301  auto const rhoSR = rho_ref_SR * (1. - ip_data.thermal_volume_strain +
302  (ip_data.alpha_B - 1.) * div_u);
303 #else // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
304  auto const rhoSR = rho_ref_SR;
305 #endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
306 
307  auto const T_prev = T - T_dot * dt;
308  ip_cv.C = ip_data.updateConstitutiveRelation(vars, t, pos, dt, T_prev);
309 
310  // constitutive model object as specified in process creation
311  auto& ptm = *_process_data.phase_transition_model_;
312  ptm.computeConstitutiveVariables(&medium, vars, pos, t, dt);
313  auto& c = ptm.cv;
314  auto const phi_L = ip_data.s_L * ip_data.phi;
315  auto const phi_G = (1. - ip_data.s_L) * ip_data.phi;
316 
317  // TODO (Grunwald): individual volume fractions can be stored in a
318  // container of type MPL::Composition (a.k.a. std::vector<double>) which
319  // can be stored in the variable array for access in MPL properties
320  // ---
321  // MaterialPropertyLib::Composition volume_fraction{phi_G, phi_L,
322  // phi_S};
323  // vars[static_cast<int>(MPL::Variable::volume_fraction)] =
324  // volume_fraction;
325  // ---
326  // TODO (Grunwald) replace effective thermal conductivity by a more
327  // sophisticated law by allowing the law to be chosen in the project
328  // file as medium property, e.g.
329  // lambda = medium.property(MPL::PropertyType::thermal_conductivity)..
330  // where volume fraction is stored in the variable array
331 
332  auto const lambdaGR = MPL::formEigenTensor<DisplacementDim>(c.lambdaGR);
333  auto const lambdaLR = MPL::formEigenTensor<DisplacementDim>(c.lambdaLR);
334 
335  ip_data.lambda = phi_G * lambdaGR + phi_L * lambdaLR + phi_S * lambdaSR;
336 
337  auto const cpS =
338  solid_phase.property(MPL::PropertyType::specific_heat_capacity)
339  .template value<double>(vars, pos, t, dt);
340  ip_data.h_S = cpS * T;
341  auto const u_S = ip_data.h_S;
342 
343  ip_data.rho_u_eff = phi_G * c.rhoGR * c.uG + phi_L * c.rhoLR * c.uL +
344  phi_S * rhoSR * u_S;
345 
346  ip_data.rho_G_h_G = phi_G * c.rhoGR * c.hG;
347  ip_data.rho_L_h_L = phi_L * c.rhoLR * c.hL;
348  ip_data.rho_S_h_S = phi_S * rhoSR * ip_data.h_S;
349 
350  ip_data.muGR = c.muGR;
351  ip_data.muLR = c.muLR;
352 
353  ip_data.rhoGR = c.rhoGR;
354  ip_data.rhoLR = c.rhoLR;
355  ip_data.rhoSR = rhoSR;
356 
357  ip_data.rhoCGR = c.rhoCGR;
358  ip_data.rhoCLR = c.rhoCLR;
359  ip_data.rhoWGR = c.rhoWGR;
360  ip_data.rhoWLR = c.rhoWLR;
361 
362  ip_data.dxmCG_dpGR = c.dxmCG_dpGR;
363  ip_data.dxmCG_dT = c.dxmCG_dT;
364  ip_data.dxmCL_dpLR = c.dxmCL_dpLR;
365  ip_data.dxmCL_dT = c.dxmCL_dT;
366  ip_data.dxmWG_dpGR = c.dxmWG_dpGR;
367  ip_data.dxmWG_dT = c.dxmWG_dT;
368  ip_data.dxmWL_dpLR = c.dxmWL_dpLR;
369  ip_data.dxmWL_dT = c.dxmWL_dT;
370 
371  // for variable output
372  ip_data.xnCG = c.xnCG;
373  ip_data.xmCG = c.xmCG;
374  ip_data.xmWL = c.xmWL;
375 
376  ip_data.diffusion_coefficient_vapour = c.diffusion_coefficient_vapour;
377  ip_data.diffusion_coefficient_solvate = c.diffusion_coefficient_solvate;
378 
379  ip_data.h_G = c.hG;
380  ip_data.h_L = c.hL;
381  ip_data.pWGR = c.pWGR;
382 
383  // ---------------------------------------------------------------------
384  // Derivatives for Jacobian
385  // ---------------------------------------------------------------------
386  auto const drho_LR_dT =
387  liquid_phase.property(MPL::PropertyType::density)
388  .template dValue<double>(vars, MPL::Variable::temperature, pos,
389  t, dt);
390  auto const drho_SR_dT =
391  solid_phase.property(MPL::PropertyType::density)
392  .template dValue<double>(vars, MPL::Variable::temperature,
393  pos, t, dt)
394 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
395  * (1. - ip_data.thermal_volume_strain +
396  (ip_data.alpha_B - 1.) * div_u) -
397  rho_ref_SR * ip_data.beta_T_SR
398 #endif
399  ;
400 
401  // porosity
402  auto const dphi_0_dT =
403  medium[MPL::PropertyType::porosity].template dValue<double>(
404  vars, MPL::Variable::temperature, pos, t, dt);
405 
406  auto const dphi_S_0_dT = -dphi_0_dT;
407  const double dphi_S_dT = dphi_S_0_dT
408 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
409  * (1. + ip_data.thermal_volume_strain -
410  ip_data.alpha_B * div_u) +
411  phi_S_0 * ip_data.beta_T_SR
412 #endif
413  ;
414 
415  ip_cv.drho_u_eff_dT =
416  phi_G * c.drho_GR_dT * c.uG + phi_G * c.rhoGR * c.du_G_dT +
417  phi_L * drho_LR_dT * c.uL + phi_L * c.rhoLR * c.du_L_dT +
418  phi_S * drho_SR_dT * u_S + phi_S * rhoSR * cpS +
419  dphi_S_dT * rhoSR * u_S;
420 
421  ip_cv.ds_L_dp_cap =
422  medium[MPL::PropertyType::saturation].template dValue<double>(
423  vars, MPL::Variable::capillary_pressure, pos, t, dt);
424 
425  // ds_L_dp_GR = 0;
426  // ds_G_dp_GR = -ds_L_dp_GR;
427  double const ds_G_dp_cap = -ip_cv.ds_L_dp_cap;
428 
429  // dphi_G_dp_GR = -ds_L_dp_GR * ip_data.phi = 0;
430  double const dphi_G_dp_cap = -ip_cv.ds_L_dp_cap * ip_data.phi;
431  // dphi_L_dp_GR = ds_L_dp_GR * ip_data.phi = 0;
432  double const dphi_L_dp_cap = ip_cv.ds_L_dp_cap * ip_data.phi;
433 
434  auto const dlambda_GR_dT = MPL::formEigenTensor<DisplacementDim>(
435  gas_phase[MPL::PropertyType::thermal_conductivity].dValue(
436  vars, MPL::Variable::temperature, pos, t, dt));
437  auto const dlambda_LR_dT = MPL::formEigenTensor<DisplacementDim>(
438  liquid_phase[MPL::PropertyType::thermal_conductivity].dValue(
439  vars, MPL::Variable::temperature, pos, t, dt));
440  auto const dlambda_SR_dT = MPL::formEigenTensor<DisplacementDim>(
441  solid_phase[MPL::PropertyType::thermal_conductivity].dValue(
442  vars, MPL::Variable::temperature, pos, t, dt));
443 
444  ip_cv.dlambda_dp_cap =
445  dphi_G_dp_cap * lambdaGR + dphi_L_dp_cap * lambdaLR;
446 
447  ip_cv.dlambda_dT = phi_G * dlambda_GR_dT + phi_L * dlambda_LR_dT +
448  phi_S * dlambda_SR_dT + dphi_S_dT * lambdaSR;
449 
450  // From p_LR = p_GR - p_cap it follows for
451  // drho_LR/dp_GR = drho_LR/dp_LR * dp_LR/dp_GR
452  // = drho_LR/dp_LR * (dp_GR/dp_GR - dp_cap/dp_GR)
453  // = drho_LR/dp_LR * (1 - 0)
454  double const drho_LR_dp_GR = c.drho_LR_dp_LR;
455  double const drho_LR_dp_cap = -c.drho_LR_dp_LR;
456  // drho_GR_dp_cap = 0;
457 
458  ip_cv.drho_h_eff_dp_GR =
459  /*(dphi_G_dp_GR = 0) * c.rhoGR * c.hG +*/ phi_G * c.drho_GR_dp_GR *
460  c.hG +
461  /*(dphi_L_dp_GR = 0) * c.rhoLR * c.hL +*/ phi_L * drho_LR_dp_GR *
462  c.hL;
463  ip_cv.drho_h_eff_dp_cap = dphi_G_dp_cap * c.rhoGR * c.hG +
464  /*phi_G * (drho_GR_dp_cap = 0) * c.hG +*/
465  dphi_L_dp_cap * c.rhoLR * c.hL +
466  phi_L * drho_LR_dp_cap * c.hL;
467 
468  // TODO (naumov) Extend for temperature dependent porosities.
469  constexpr double dphi_G_dT = 0;
470  constexpr double dphi_L_dT = 0;
471  ip_cv.drho_h_eff_dT =
472  dphi_G_dT * c.rhoGR * c.hG + phi_G * c.drho_GR_dT * c.hG +
473  phi_G * c.rhoGR * c.dh_G_dT + dphi_L_dT * c.rhoLR * c.hL +
474  phi_L * drho_LR_dT * c.hL + phi_L * c.rhoLR * c.dh_L_dT +
475  dphi_S_dT * rhoSR * ip_data.h_S + phi_S * drho_SR_dT * ip_data.h_S +
476  phi_S * rhoSR * cpS;
477 
478  ip_cv.drho_u_eff_dp_GR =
479  /*(dphi_G_dp_GR = 0) * c.rhoGR * c.uG +*/
480  phi_G * c.drho_GR_dp_GR * c.uG + phi_G * c.rhoGR * c.du_G_dp_GR +
481  /*(dphi_L_dp_GR = 0) * c.rhoLR * c.uL +*/
482  phi_L * drho_LR_dp_GR * c.uL + phi_L * c.rhoLR * c.du_L_dp_GR;
483 
484  ip_cv.drho_u_eff_dp_cap = dphi_G_dp_cap * c.rhoGR * c.uG +
485  /*phi_G * (drho_GR_dp_cap = 0) * c.uG +*/
486  dphi_L_dp_cap * c.rhoLR * c.uL +
487  phi_L * drho_LR_dp_cap * c.uL +
488  phi_L * c.rhoLR * c.du_L_dp_cap;
489 
490  auto const& b = _process_data.specific_body_force;
491  auto const k_over_mu_G =
492  (ip_data.k_S * ip_data.k_rel_G / ip_data.muGR).eval();
493  auto const k_over_mu_L =
494  (ip_data.k_S * ip_data.k_rel_L / ip_data.muLR).eval();
495 
496  // dk_over_mu_G_dp_GR =
497  // ip_data.k_S * dk_rel_G_ds_L * (ds_L_dp_GR = 0) / ip_data.muGR =
498  // 0;
499  // dk_over_mu_L_dp_GR =
500  // ip_data.k_S * dk_rel_L_ds_L * (ds_L_dp_GR = 0) / ip_data.muLR =
501  // 0;
502  ip_cv.dk_over_mu_G_dp_cap =
503  ip_data.k_S * dk_rel_G_ds_L * ip_cv.ds_L_dp_cap / ip_data.muGR;
504  ip_cv.dk_over_mu_L_dp_cap =
505  ip_data.k_S * dk_rel_L_ds_L * ip_cv.ds_L_dp_cap / ip_data.muLR;
506 
507  GlobalDimVectorType const w_GS =
508  k_over_mu_G * c.rhoGR * b - k_over_mu_G * gradpGR;
509  GlobalDimVectorType const w_LS = k_over_mu_L * gradpCap +
510  k_over_mu_L * c.rhoLR * b -
511  k_over_mu_L * gradpGR;
512 
513  ip_cv.drho_GR_h_w_eff_dp_GR_Npart =
514  c.drho_GR_dp_GR * c.hG * w_GS +
515  c.rhoGR * c.hG * k_over_mu_G * c.drho_GR_dp_GR * b;
516  ip_cv.drho_GR_h_w_eff_dp_GR_gradNpart =
517  -c.rhoGR * c.hG * k_over_mu_G - c.rhoLR * c.hL * k_over_mu_L;
518 
519  ip_cv.drho_LR_h_w_eff_dp_cap_Npart =
520  -drho_LR_dp_cap * c.hL * w_LS -
521  c.rhoLR * c.hL * k_over_mu_L * drho_LR_dp_cap * b;
522  ip_cv.drho_LR_h_w_eff_dp_cap_gradNpart =
523  // TODO (naumov) why the minus sign??????
524  -c.rhoLR * c.hL * k_over_mu_L;
525 
526  ip_cv.drho_GR_h_w_eff_dT =
527  c.drho_GR_dT * c.hG * w_GS + c.rhoGR * c.dh_G_dT * w_GS +
528  drho_LR_dT * c.hL * w_LS + c.rhoLR * c.dh_L_dT * w_LS;
529  // TODO (naumov) + k_over_mu_G * drho_GR_dT * b + k_over_mu_L *
530  // drho_LR_dT * b
531 
532  // Derivatives of s_G * rho_C_GR_dot + s_L * rho_C_LR_dot abbreviated
533  // here with S_rho_C_eff.
534  double const s_L = ip_data.s_L;
535  double const s_G = 1. - ip_data.s_L;
536  double const rho_C_GR_dot = (ip_data.rhoCGR - ip_data.rhoCGR_prev) / dt;
537  double const rho_C_LR_dot = (ip_data.rhoCLR - ip_data.rhoCLR_prev) / dt;
538  double const rho_C_FR = s_G * ip_data.rhoCGR + s_L * ip_data.rhoCLR;
539  double const rho_W_FR = s_G * ip_data.rhoWGR + s_L * ip_data.rhoWLR;
540  // TODO (naumov) Extend for partially saturated media.
541  constexpr double drho_C_GR_dp_cap = 0;
542  ip_cv.dfC_3a_dp_GR =
543  /*(ds_G_dp_GR = 0) * rho_C_GR_dot +*/ s_G * c.drho_C_GR_dp_GR / dt +
544  /*(ds_L_dp_GR = 0) * rho_C_LR_dot +*/ s_L * c.drho_C_LR_dp_GR / dt;
545  ip_cv.dfC_3a_dp_cap =
546  ds_G_dp_cap * rho_C_GR_dot + s_G * drho_C_GR_dp_cap / dt +
547  ip_cv.ds_L_dp_cap * rho_C_LR_dot - s_L * c.drho_C_LR_dp_LR / dt;
548  ip_cv.dfC_3a_dT = s_G * c.drho_C_GR_dT / dt + s_L * c.drho_C_LR_dT / dt;
549 
550  double const drho_C_FR_dp_GR =
551  /*(ds_G_dp_GR = 0) * ip_data.rhoCGR +*/ s_G * c.drho_C_GR_dp_GR +
552  /*(ds_L_dp_GR = 0) * ip_data.rhoCLR +*/ s_L * c.drho_C_LR_dp_GR;
553  ip_cv.dfC_4_MCpG_dp_GR = drho_C_FR_dp_GR *
554  (ip_data.alpha_B - ip_data.phi) *
555  ip_data.beta_p_SR;
556 
557  double const drho_C_FR_dT = s_G * c.drho_C_GR_dT + s_L * c.drho_C_LR_dT;
558  ip_cv.dfC_4_MCpG_dT =
559  drho_C_FR_dT * (ip_data.alpha_B - ip_data.phi) * ip_data.beta_p_SR
560 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
561  - rho_C_FR * ip_data.dphi_dT * ip_data.beta_p_SR
562 #endif
563  ;
564 
565  ip_cv.dfC_4_MCT_dT =
566  drho_C_FR_dT * (ip_data.alpha_B - ip_data.phi) * ip_data.beta_T_SR
567 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
568  + rho_C_FR * (ip_data.alpha_B - ip_data.dphi_dT) * ip_data.beta_T_SR
569 #endif
570  ;
571 
572  ip_cv.dfC_4_MCu_dT = drho_C_FR_dT * ip_data.alpha_B;
573 
574  ip_cv.dfC_2a_dp_GR = -ip_data.phi * c.drho_C_GR_dp_GR -
575  drho_C_FR_dp_GR * pCap *
576  (ip_data.alpha_B - ip_data.phi) *
577  ip_data.beta_p_SR;
578 
579  double const drho_C_FR_dp_cap =
580  ds_G_dp_cap * ip_data.rhoCGR + s_G * drho_C_GR_dp_cap +
581  ip_cv.ds_L_dp_cap * ip_data.rhoCLR - s_L * c.drho_C_LR_dp_LR;
582 
583  ip_cv.dfC_2a_dp_cap =
584  ip_data.phi * (-c.drho_C_LR_dp_LR - drho_C_GR_dp_cap) -
585  drho_C_FR_dp_cap * pCap * (ip_data.alpha_B - ip_data.phi) *
586  ip_data.beta_p_SR +
587  rho_C_FR * (ip_data.alpha_B - ip_data.phi) * ip_data.beta_p_SR;
588 
589  ip_cv.dfC_2a_dT =
590 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
591  ip_data.dphi_dT * (ip_data.rhoCLR - ip_data.rhoCGR) +
592 #endif
593  ip_data.phi * (c.drho_C_LR_dT - c.drho_C_GR_dT) -
594  drho_C_FR_dT * pCap * (ip_data.alpha_B - ip_data.phi) *
595  ip_data.beta_p_SR
596 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
597  + rho_C_FR * pCap * ip_data.dphi_dT * ip_data.beta_p_SR
598 #endif
599  ;
600 
601  ip_cv.dadvection_C_dp_GR =
602  c.drho_C_GR_dp_GR * k_over_mu_G + c.drho_C_LR_dp_GR * k_over_mu_L;
603 
604  ip_cv.dfC_4_LCpG_dT =
605  c.drho_C_GR_dT * k_over_mu_G + c.drho_C_LR_dT * k_over_mu_L
606  // + ip_cv.ddiffusion_C_p_dT TODO (naumov)
607  ;
608 
609  double const drho_W_FR_dp_GR =
610  /*(ds_G_dp_GR = 0) * ip_data.rhoWGR +*/ s_G * c.drho_W_GR_dp_GR +
611  /*(ds_L_dp_GR = 0) * ip_data.rhoWLR +*/ s_L * c.drho_W_LR_dp_GR;
612  double const drho_W_FR_dp_cap =
613  ds_G_dp_cap * ip_data.rhoWGR + s_G * c.drho_W_GR_dp_cap +
614  ip_cv.ds_L_dp_cap * ip_data.rhoWLR - s_L * c.drho_W_LR_dp_LR;
615  double const drho_W_FR_dT = s_G * c.drho_W_GR_dT + s_L * c.drho_W_LR_dT;
616 
617  ip_cv.dfW_2a_dp_GR =
618  ip_data.phi * (c.drho_W_LR_dp_GR - c.drho_W_GR_dp_GR);
619  ip_cv.dfW_2b_dp_GR = drho_W_FR_dp_GR * pCap *
620  (ip_data.alpha_B - ip_data.phi) *
621  ip_data.beta_p_SR;
622  ip_cv.dfW_2a_dp_cap =
623  ip_data.phi * (-c.drho_W_LR_dp_LR - c.drho_W_GR_dp_cap);
624  ip_cv.dfW_2b_dp_cap =
625  drho_W_FR_dp_cap * pCap * (ip_data.alpha_B - ip_data.phi) *
626  ip_data.beta_p_SR +
627  rho_W_FR * (ip_data.alpha_B - ip_data.phi) * ip_data.beta_p_SR;
628 
629  ip_cv.dfW_2a_dT =
630 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
631  ip_data.dphi_dT * (ip_data.rhoWLR - ip_data.rhoWGR) +
632 #endif
633  ip_data.phi * (c.drho_W_LR_dT - c.drho_W_GR_dT);
634  ip_cv.dfW_2b_dT =
635  drho_W_FR_dT * pCap * (ip_data.alpha_B - ip_data.phi) *
636  ip_data.beta_p_SR
637 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
638  - rho_W_FR * pCap * ip_data.dphi_dT * ip_data.beta_p_SR
639 #endif
640  ;
641 
642  double const rho_W_GR_dot = (ip_data.rhoWGR - ip_data.rhoWGR_prev) / dt;
643  double const rho_W_LR_dot = (ip_data.rhoWLR - ip_data.rhoWLR_prev) / dt;
644 
645  ip_cv.dfW_3a_dp_GR =
646  /*(ds_G_dp_GR = 0) * rho_W_GR_dot +*/ s_G * c.drho_W_GR_dp_GR / dt +
647  /*(ds_L_dp_GR = 0) * rho_W_LR_dot +*/ s_L * c.drho_W_LR_dp_GR / dt;
648  ip_cv.dfW_3a_dp_cap =
649  ds_G_dp_cap * rho_W_GR_dot + s_G * c.drho_W_GR_dp_cap / dt +
650  ip_cv.ds_L_dp_cap * rho_W_LR_dot - s_L * c.drho_W_LR_dp_LR / dt;
651  ip_cv.dfW_3a_dT = s_G * c.drho_W_GR_dT / dt + s_L * c.drho_W_LR_dT / dt;
652 
653  ip_cv.dfW_4_LWpG_a_dp_GR = c.drho_W_GR_dp_GR * k_over_mu_G
654  // + rhoWGR * (dk_over_mu_G_dp_GR = 0)
655  + c.drho_W_LR_dp_GR * k_over_mu_L
656  // + rhoWLR * (dk_over_mu_L_dp_GR = 0)
657  ;
658  ip_cv.dfW_4_LWpG_a_dp_cap = -c.drho_W_LR_dp_LR * k_over_mu_L;
659  ip_cv.dfW_4_LWpG_a_dT =
660  c.drho_W_GR_dT * k_over_mu_G
661  //+ rhoWGR * (dk_over_mu_G_dT != 0 TODO for mu_G(T))
662  + c.drho_W_LR_dT * k_over_mu_L
663  //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_G(T))
664  ;
665 
666  // TODO (naumov) for dxmW*/d* != 0
667  ip_cv.dfW_4_LWpG_d_dp_GR =
668  Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
669  ip_cv.dfW_4_LWpG_d_dp_cap =
670  Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
671  ip_cv.dfW_4_LWpG_d_dT =
672  Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
673 
674  ip_cv.dfW_4_LWpC_a_dp_GR = c.drho_W_LR_dp_GR * k_over_mu_L
675  //+ rhoWLR * (dk_over_mu_L_dp_GR = 0)
676  ;
677  ip_cv.dfW_4_LWpC_a_dp_cap = -c.drho_W_LR_dp_LR * k_over_mu_L +
678  ip_data.rhoWLR * ip_cv.dk_over_mu_L_dp_cap;
679  ip_cv.dfW_4_LWpC_a_dT = c.drho_W_LR_dT * k_over_mu_L
680  //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_L(T))
681  ;
682 
683  // TODO (naumov) for dxmW*/d* != 0
684  ip_cv.dfW_4_LWpC_d_dp_GR =
685  Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
686  ip_cv.dfW_4_LWpC_d_dp_cap =
687  Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
688  ip_cv.dfW_4_LWpC_d_dT =
689  Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
690  }
691 
692  return ip_constitutive_variables;
693 }
template Eigen::Matrix< double, 3, 3 > formEigenTensor< 3 >(MaterialPropertyLib::PropertyDataType const &values)
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:48
KelvinVectorType< DisplacementDim > tensorToKelvin(Eigen::Matrix< double, 3, 3 > const &m)

References MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::c, MaterialPropertyLib::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), MaterialPropertyLib::density, MaterialPropertyLib::formEigenTensor< 3 >(), NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::mechanical_strain, MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, MaterialPropertyLib::relative_permeability, MaterialPropertyLib::relative_permeability_nonwetting_phase, MaterialPropertyLib::saturation, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), MaterialPropertyLib::specific_heat_capacity, MaterialPropertyLib::temperature, MathLib::KelvinVector::tensorToKelvin(), MaterialPropertyLib::thermal_conductivity, and MaterialPropertyLib::thermal_expansivity.

Member Data Documentation

◆ _element

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int DisplacementDim>
MeshLib::Element const& ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::_element
private

◆ _integration_method

◆ _ip_data

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

Definition at line 394 of file TH2MFEM.h.

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

◆ _is_axially_symmetric

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int DisplacementDim>
bool const ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::_is_axially_symmetric
private

Definition at line 398 of file TH2MFEM.h.

◆ _process_data

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int DisplacementDim>
TH2MProcessData<DisplacementDim>& ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::_process_data
private

◆ _secondary_data

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

◆ C_index

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

Definition at line 415 of file TH2MFEM.h.

◆ C_size

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

Definition at line 416 of file TH2MFEM.h.

◆ capillary_pressure_index

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

Definition at line 407 of file TH2MFEM.h.

◆ capillary_pressure_size

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

Definition at line 408 of file TH2MFEM.h.

◆ displacement_index

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

Definition at line 411 of file TH2MFEM.h.

◆ displacement_size

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

Definition at line 412 of file TH2MFEM.h.

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

◆ gas_pressure_index

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

Definition at line 405 of file TH2MFEM.h.

◆ gas_pressure_size

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

Definition at line 406 of file TH2MFEM.h.

◆ KelvinVectorSize

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int DisplacementDim>
int const ProcessLib::TH2M::TH2MLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::KelvinVectorSize
static
Initial value:

Definition at line 71 of file TH2MFEM.h.

◆ temperature_index

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

Definition at line 409 of file TH2MFEM.h.

◆ temperature_size

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

Definition at line 410 of file TH2MFEM.h.

◆ W_index

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

Definition at line 417 of file TH2MFEM.h.

◆ W_size

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

Definition at line 418 of file TH2MFEM.h.


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