OGS
ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim > Class Template Reference

Detailed Description

template<typename ShapeFunctionDisplacement, typename ShapeFunction, typename IntegrationMethod, int DisplacementDim>
class ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >

Definition at line 46 of file ThermoRichardsMechanicsFEM.h.

#include <ThermoRichardsMechanicsFEM.h>

Inheritance diagram for ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >:
[legend]
Collaboration diagram for ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >:
[legend]

Public Types

using ShapeMatricesTypeDisplacement = ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim >
 
using ShapeMatricesType = ShapeMatrixPolicyType< ShapeFunction, DisplacementDim >
 
using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType
 
using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType
 
using BMatricesType = BMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim >
 
using KelvinVectorType = typename BMatricesType::KelvinVectorType
 
using IpData = IntegrationPointData< BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesType, DisplacementDim, ShapeFunctionDisplacement::NPOINTS >
 
using Invariants = MathLib::KelvinVector::Invariants< KelvinVectorSize >
 
using SymmetricTensor = Eigen::Matrix< double, KelvinVectorSize, 1 >
 

Public Member Functions

 ThermoRichardsMechanicsLocalAssembler (ThermoRichardsMechanicsLocalAssembler const &)=delete
 
 ThermoRichardsMechanicsLocalAssembler (ThermoRichardsMechanicsLocalAssembler &&)=delete
 
 ThermoRichardsMechanicsLocalAssembler (MeshLib::Element const &e, std::size_t const, bool const is_axially_symmetric, unsigned const integration_order, ThermoRichardsMechanicsProcessData< DisplacementDim > &process_data)
 
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 use_monolithic_scheme, int const process_id) 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 > getSigma () const override
 
std::vector< double > const & getIntPtDarcyVelocity (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 > getSaturation () const override
 
std::vector< double > const & getIntPtSaturation (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 > getPorosity () const override
 
std::vector< double > const & getIntPtPorosity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::vector< double > getTransportPorosity () const override
 
std::vector< double > const & getIntPtTransportPorosity (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 & getIntPtSigma (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 > getSwellingStress () const override
 
std::vector< double > const & getIntPtSwellingStress (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 > getEpsilon () const override
 
std::vector< double > const & getIntPtEpsilon (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 & getIntPtDryDensitySolid (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
- Public Member Functions inherited from ProcessLib::LocalAssemblerInterface
virtual ~LocalAssemblerInterface ()=default
 
void setInitialConditions (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, bool const use_monolithic_scheme, int const process_id)
 
virtual void initialize (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
 
virtual void preAssemble (double const, double const, std::vector< double > const &)
 
virtual void assemble (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_xdot, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
 
virtual void assembleForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_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 Member Functions

unsigned getNumberOfIntegrationPoints () const override
 
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt (unsigned integration_point) const override
 

Private Attributes

ThermoRichardsMechanicsProcessData< 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::ShapeTypesecondary_data_
 

Static Private Attributes

static const int temperature_index = 0
 
static const int temperature_size = ShapeFunction::NPOINTS
 
static const int pressure_index = temperature_size
 
static const int pressure_size = ShapeFunction::NPOINTS
 
static const int displacement_index = 2 * ShapeFunction::NPOINTS
 
static const int displacement_size
 

Member Typedef Documentation

◆ BMatricesType

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
using ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::BMatricesType = BMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>

Definition at line 60 of file ThermoRichardsMechanicsFEM.h.

◆ GlobalDimMatrixType

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
using ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType

Definition at line 57 of file ThermoRichardsMechanicsFEM.h.

◆ GlobalDimVectorType

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
using ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType

Definition at line 58 of file ThermoRichardsMechanicsFEM.h.

◆ Invariants

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
using ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>

Definition at line 71 of file ThermoRichardsMechanicsFEM.h.

◆ IpData

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
using ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::IpData = IntegrationPointData<BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesType, DisplacementDim, ShapeFunctionDisplacement::NPOINTS>

Definition at line 64 of file ThermoRichardsMechanicsFEM.h.

◆ KelvinVectorType

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
using ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::KelvinVectorType = typename BMatricesType::KelvinVectorType

Definition at line 62 of file ThermoRichardsMechanicsFEM.h.

◆ ShapeMatricesType

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
using ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, DisplacementDim>

Definition at line 54 of file ThermoRichardsMechanicsFEM.h.

◆ ShapeMatricesTypeDisplacement

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
using ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::ShapeMatricesTypeDisplacement = ShapeMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>

Definition at line 50 of file ThermoRichardsMechanicsFEM.h.

◆ SymmetricTensor

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
using ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>

Definition at line 73 of file ThermoRichardsMechanicsFEM.h.

Constructor & Destructor Documentation

◆ ThermoRichardsMechanicsLocalAssembler() [1/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::ThermoRichardsMechanicsLocalAssembler ( ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim > const &  )
delete

◆ ThermoRichardsMechanicsLocalAssembler() [2/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::ThermoRichardsMechanicsLocalAssembler ( ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim > &&  )
delete

◆ ThermoRichardsMechanicsLocalAssembler() [3/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::ThermoRichardsMechanicsLocalAssembler ( MeshLib::Element const &  e,
std::size_t const  ,
bool const  is_axially_symmetric,
unsigned const  integration_order,
ThermoRichardsMechanicsProcessData< DisplacementDim > &  process_data 
)

Definition at line 33 of file ThermoRichardsMechanicsFEM-impl.h.

40  : process_data_(process_data),
41  integration_method_(integration_order),
42  element_(e),
43  is_axially_symmetric_(is_axially_symmetric)
44 {
45  unsigned const n_integration_points =
46  integration_method_.getNumberOfPoints();
47 
48  ip_data_.reserve(n_integration_points);
49  secondary_data_.N_u.resize(n_integration_points);
50 
51  auto const shape_matrices_u =
52  NumLib::initShapeMatrices<ShapeFunctionDisplacement,
54  DisplacementDim>(e, is_axially_symmetric,
56 
57  auto const shape_matrices =
59  DisplacementDim>(e, is_axially_symmetric,
61 
62  auto const& solid_material =
64  process_data_.solid_materials, process_data_.material_ids,
65  e.getID());
66 
67  auto const& medium = process_data_.media_map->getMedium(element_.getID());
68 
70  x_position.setElementID(element_.getID());
71  for (unsigned ip = 0; ip < n_integration_points; ip++)
72  {
73  x_position.setIntegrationPoint(ip);
74  ip_data_.emplace_back(solid_material);
75  auto& ip_data = ip_data_[ip];
76  auto const& sm_u = shape_matrices_u[ip];
77  ip_data_[ip].integration_weight =
78  integration_method_.getWeightedPoint(ip).getWeight() *
79  sm_u.integralMeasure * sm_u.detJ;
80 
81  ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
82  DisplacementDim, displacement_size>::Zero(DisplacementDim,
84  for (int i = 0; i < DisplacementDim; ++i)
85  {
86  ip_data.N_u_op
87  .template block<1, displacement_size / DisplacementDim>(
88  i, i * displacement_size / DisplacementDim)
89  .noalias() = sm_u.N;
90  }
91 
92  ip_data.N_u = sm_u.N;
93  ip_data.dNdx_u = sm_u.dNdx;
94 
95  // ip_data.N_p and ip_data.dNdx_p are used for both p and T variables
96  ip_data.N_p = shape_matrices[ip].N;
97  ip_data.dNdx_p = shape_matrices[ip].dNdx;
98 
99  // Initial porosity. Could be read from integration point data or mesh.
100  ip_data.porosity =
101  medium->property(MPL::porosity)
102  .template initialValue<double>(
103  x_position,
104  std::numeric_limits<
105  double>::quiet_NaN() /* t independent */);
106 
107  ip_data.transport_porosity = ip_data.porosity;
108  if (medium->hasProperty(MPL::PropertyType::transport_porosity))
109  {
110  ip_data.transport_porosity =
111  medium->property(MPL::transport_porosity)
112  .template initialValue<double>(
113  x_position,
114  std::numeric_limits<
115  double>::quiet_NaN() /* t independent */);
116  }
117 
118  secondary_data_.N_u[ip] = shape_matrices_u[ip].N;
119  }
120 }
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)
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
ThermoRichardsMechanicsProcessData< DisplacementDim > & process_data_
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > secondary_data_
std::vector< IpData, Eigen::aligned_allocator< IpData > > ip_data_
ShapeMatrixPolicyType< ShapeFunction, DisplacementDim > ShapeMatricesType
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::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::displacement_size, ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::element_, MeshLib::Element::getID(), NumLib::initShapeMatrices(), ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::integration_method_, ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::ip_data_, ProcessLib::HydroMechanics::SecondaryData< ShapeMatrixType >::N_u, MaterialPropertyLib::porosity, ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::process_data_, ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::secondary_data_, MaterialLib::Solids::selectSolidConstitutiveRelation(), ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), and MaterialPropertyLib::transport_porosity.

Member Function Documentation

◆ assembleWithJacobian()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
void ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, 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 
)
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 253 of file ThermoRichardsMechanicsFEM-impl.h.

262 {
263  auto const local_matrix_dim =
265  assert(local_x.size() == local_matrix_dim);
266 
267  auto const T = Eigen::Map<typename ShapeMatricesType::template VectorType<
268  temperature_size> const>(local_x.data() + temperature_index,
270  auto const p_L = Eigen::Map<
271  typename ShapeMatricesType::template VectorType<pressure_size> const>(
272  local_x.data() + pressure_index, pressure_size);
273 
274  auto const u =
275  Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
276  displacement_size> const>(local_x.data() + displacement_index,
278 
279  auto const T_dot =
280  Eigen::Map<typename ShapeMatricesType::template VectorType<
281  temperature_size> const>(local_xdot.data() + temperature_index,
283  auto const p_L_dot = Eigen::Map<
284  typename ShapeMatricesType::template VectorType<pressure_size> const>(
285  local_xdot.data() + pressure_index, pressure_size);
286  auto const u_dot =
287  Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
288  displacement_size> const>(local_xdot.data() + displacement_index,
290 
291  auto local_Jac = MathLib::createZeroedMatrix<
292  typename ShapeMatricesTypeDisplacement::template MatrixType<
293  local_matrix_dim, local_matrix_dim>>(
294  local_Jac_data, local_matrix_dim, local_matrix_dim);
295 
296  auto local_rhs =
297  MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
298  template VectorType<local_matrix_dim>>(
299  local_rhs_data, local_matrix_dim);
300 
301  auto const& identity2 = MathLib::KelvinVector::Invariants<
303  DisplacementDim)>::identity2;
304 
305  typename ShapeMatricesType::NodalMatrixType M_TT =
306  ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
308  typename ShapeMatricesType::NodalMatrixType K_TT =
309  ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
311  typename ShapeMatricesType::NodalMatrixType K_Tp =
312  ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
313  pressure_size);
314  typename ShapeMatricesType::NodalMatrixType M_pT =
315  ShapeMatricesType::NodalMatrixType::Zero(pressure_size,
317  typename ShapeMatricesType::NodalMatrixType laplace_p =
318  ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
319 
320  typename ShapeMatricesType::NodalMatrixType storage_p_a_p =
321  ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
322 
323  typename ShapeMatricesType::NodalMatrixType storage_p_a_S_Jpp =
324  ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
325 
326  typename ShapeMatricesType::NodalMatrixType storage_p_a_S =
327  ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
328 
329  typename ShapeMatricesTypeDisplacement::template MatrixType<
331  Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
334 
335  auto const& medium = process_data_.media_map->getMedium(element_.getID());
336  auto const& liquid_phase = medium->phase("AqueousLiquid");
337  auto const& solid_phase = medium->phase("Solid");
338  MPL::VariableArray variables;
339  MPL::VariableArray variables_prev;
340 
342  x_position.setElementID(element_.getID());
343 
344  unsigned const n_integration_points =
345  integration_method_.getNumberOfPoints();
346  for (unsigned ip = 0; ip < n_integration_points; ip++)
347  {
348  x_position.setIntegrationPoint(ip);
349  auto& ip_data = ip_data_[ip];
350  auto const& w = ip_data.integration_weight;
351 
352  auto const& N_u_op = ip_data.N_u_op;
353 
354  auto const& N_u = ip_data.N_u;
355  auto const& dNdx_u = ip_data.dNdx_u;
356 
357  // N and dNdx are used for both p and T variables
358  auto const& N = ip_data.N_p;
359  auto const& dNdx = ip_data.dNdx_p;
360 
361  auto const x_coord =
362  NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
364  element_, N_u);
365  auto const B =
366  LinearBMatrix::computeBMatrix<DisplacementDim,
367  ShapeFunctionDisplacement::NPOINTS,
368  typename BMatricesType::BMatrixType>(
369  dNdx_u, N_u, x_coord, is_axially_symmetric_);
370 
371  double T_ip;
373  double T_dot_ip;
374  NumLib::shapeFunctionInterpolate(T_dot, N, T_dot_ip);
375  double const dT = T_dot_ip * dt;
376 
377  double p_cap_ip;
378  NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
379 
380  double p_cap_dot_ip;
381  NumLib::shapeFunctionInterpolate(-p_L_dot, N, p_cap_dot_ip);
382 
383  variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
384  p_cap_ip;
385  variables[static_cast<int>(MPL::Variable::phase_pressure)] = -p_cap_ip;
386  variables[static_cast<int>(MPL::Variable::temperature)] = T_ip;
387 
388  auto& eps = ip_data.eps;
389  auto& eps_m = ip_data.eps_m;
390  eps.noalias() = B * u;
391  auto const& sigma_eff = ip_data.sigma_eff;
392  auto& S_L = ip_data.saturation;
393  auto const S_L_prev = ip_data.saturation_prev;
394  auto const alpha =
395  medium->property(MPL::PropertyType::biot_coefficient)
396  .template value<double>(variables, x_position, t, dt);
397 
398  double const T_ip_prev = T_ip - dT;
399  auto const C_el = ip_data.computeElasticTangentStiffness(
400  t, x_position, dt, T_ip_prev, T_ip);
401 
402  auto const beta_SR =
403  (1 - alpha) /
404  ip_data.solid_material.getBulkModulus(t, x_position, &C_el);
405  variables[static_cast<int>(MPL::Variable::grain_compressibility)] =
406  beta_SR;
407 
408  auto const rho_LR =
409  liquid_phase.property(MPL::PropertyType::density)
410  .template value<double>(variables, x_position, t, dt);
411  auto const& b = process_data_.specific_body_force;
412 
413  S_L = medium->property(MPL::PropertyType::saturation)
414  .template value<double>(variables, x_position, t, dt);
415  variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
416  variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
417  S_L_prev;
418 
419  // tangent derivative for Jacobian
420  double const dS_L_dp_cap =
421  medium->property(MPL::PropertyType::saturation)
422  .template dValue<double>(variables,
424  x_position, t, dt);
425  // secant derivative from time discretization for storage
426  // use tangent, if secant is not available
427  double const DeltaS_L_Deltap_cap =
428  (p_cap_dot_ip == 0) ? dS_L_dp_cap
429  : (S_L - S_L_prev) / (dt * p_cap_dot_ip);
430 
431  auto const chi = [medium, x_position, t, dt](double const S_L)
432  {
434  vs[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
435  return medium->property(MPL::PropertyType::bishops_effective_stress)
436  .template value<double>(vs, x_position, t, dt);
437  };
438  double const chi_S_L = chi(S_L);
439  double const chi_S_L_prev = chi(S_L_prev);
440 
441  variables[static_cast<int>(MPL::Variable::effective_pore_pressure)] =
442  -chi_S_L * p_cap_ip;
443  variables_prev[static_cast<int>(
444  MPL::Variable::effective_pore_pressure)] =
445  -chi_S_L_prev * (p_cap_ip - p_cap_dot_ip * dt);
446 
447  // Set volumetric strain rate for the general case without swelling.
448  variables[static_cast<int>(MPL::Variable::volumetric_strain)]
449  .emplace<double>(Invariants::trace(eps));
450  variables_prev[static_cast<int>(MPL::Variable::volumetric_strain)]
451  .emplace<double>(Invariants::trace(B * (u - u_dot * dt)));
452 
453  auto& phi = ip_data.porosity;
454  { // Porosity update
455 
456  variables_prev[static_cast<int>(MPL::Variable::porosity)] =
457  ip_data.porosity_prev;
458  phi = medium->property(MPL::PropertyType::porosity)
459  .template value<double>(variables, variables_prev,
460  x_position, t, dt);
461  variables[static_cast<int>(MPL::Variable::porosity)] = phi;
462  }
463 
464  if (alpha < phi)
465  {
466  OGS_FATAL(
467  "ThermoRichardsMechanics: Biot-coefficient {} is smaller than "
468  "porosity {} in element/integration point {}/{}.",
469  alpha, phi, element_.getID(), ip);
470  }
471 
472  // Swelling and possibly volumetric strain rate update.
473  if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
474  {
475  auto& sigma_sw = ip_data.sigma_sw;
476  auto const& sigma_sw_prev = ip_data.sigma_sw_prev;
477 
478  // If there is swelling, compute it. Update volumetric strain rate,
479  // s.t. it corresponds to the mechanical part only.
480  sigma_sw = sigma_sw_prev;
481 
482  using DimMatrix = Eigen::Matrix<double, 3, 3>;
483  auto const sigma_sw_dot =
484  MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
485  solid_phase
487  .template value<DimMatrix>(variables, variables_prev,
488  x_position, t, dt));
489  sigma_sw += sigma_sw_dot * dt;
490 
491  // !!! Misusing volumetric strain for mechanical volumetric
492  // strain just to update the transport porosity !!!
493  std::get<double>(variables[static_cast<int>(
494  MPL::Variable::volumetric_strain)]) +=
495  identity2.transpose() * C_el.inverse() * sigma_sw;
496  std::get<double>(variables_prev[static_cast<int>(
497  MPL::Variable::volumetric_strain)]) +=
498  identity2.transpose() * C_el.inverse() * sigma_sw_prev;
499  }
500 
501  if (solid_phase.hasProperty(MPL::PropertyType::transport_porosity))
502  {
503  variables_prev[static_cast<int>(
505  ip_data.transport_porosity_prev;
506 
507  ip_data.transport_porosity =
508  solid_phase.property(MPL::PropertyType::transport_porosity)
509  .template value<double>(variables, variables_prev,
510  x_position, t, dt);
511  variables[static_cast<int>(MPL::Variable::transport_porosity)] =
512  ip_data.transport_porosity;
513  }
514  else
515  {
516  variables[static_cast<int>(MPL::Variable::transport_porosity)] =
517  phi;
518  }
519 
520  //
521  // displacement equation, displacement part
522  //
523 
524  // Consider also anisotropic thermal expansion.
526  DisplacementDim> const solid_linear_thermal_expansivity_vector =
527  MPL::formKelvinVectorFromThermalExpansivity<DisplacementDim>(
528  solid_phase
529  .property(
531  .value(variables, x_position, t, dt));
532 
534  dthermal_strain = solid_linear_thermal_expansivity_vector * dT;
535 
536  auto& eps_prev = ip_data.eps_prev;
537  auto& eps_m_prev = ip_data.eps_m_prev;
538  eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
539 
540  if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
541  {
542  eps_m.noalias() +=
543  C_el.inverse() * (ip_data.sigma_sw - ip_data.sigma_sw_prev);
544  }
545 
546  variables[static_cast<int>(
549  eps_m);
550 
551  auto C = ip_data.updateConstitutiveRelation(variables, t, x_position,
552  dt, T_ip_prev);
553 
554  local_Jac
555  .template block<displacement_size, displacement_size>(
557  .noalias() += B.transpose() * C * B * w;
558 
559  double const p_FR = -chi_S_L * p_cap_ip;
560  // p_SR
561  variables[static_cast<int>(MPL::Variable::solid_grain_pressure)] =
562  p_FR - Invariants::trace(sigma_eff) / (3 * (1 - phi));
563  auto const rho_SR =
564  solid_phase.property(MPL::PropertyType::density)
565  .template value<double>(variables, x_position, t, dt);
566 
567  double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
568 
569  auto const sigma_total =
570  (sigma_eff + alpha * chi_S_L * identity2 * p_cap_ip).eval();
571 
572  local_rhs.template segment<displacement_size>(displacement_index)
573  .noalias() -=
574  (B.transpose() * sigma_total - N_u_op.transpose() * rho * b) * w;
575 
576  //
577  // displacement equation, pressure part
578  //
579  auto const dchi_dS_L =
581  .template dValue<double>(variables,
582  MPL::Variable::liquid_saturation,
583  x_position, t, dt);
584  local_Jac
585  .template block<displacement_size, pressure_size>(
587  .noalias() -= B.transpose() * alpha *
588  (chi_S_L + dchi_dS_L * p_cap_ip * dS_L_dp_cap) *
589  identity2 * N * w;
590 
591  local_Jac
592  .template block<displacement_size, pressure_size>(
594  .noalias() +=
595  N_u_op.transpose() * phi * rho_LR * dS_L_dp_cap * b * N * w;
596 
597  if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
598  {
599  using DimMatrix = Eigen::Matrix<double, 3, 3>;
600  auto const dsigma_sw_dS_L =
601  MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
602  solid_phase
604  .template dValue<DimMatrix>(
605  variables, variables_prev,
606  MPL::Variable::liquid_saturation, x_position, t,
607  dt));
608  local_Jac
609  .template block<displacement_size, pressure_size>(
611  .noalias() +=
612  B.transpose() * dsigma_sw_dS_L * dS_L_dp_cap * N * w;
613  }
614 
615  //
616  // pressure equation, displacement part.
617  //
618  Kpu.noalias() += N.transpose() * S_L * rho_LR * alpha *
619  identity2.transpose() * B * w;
620 
621  //
622  // pressure equation, pressure part.
623  //
624  double const k_rel =
626  .template value<double>(variables, x_position, t, dt);
627  auto const mu =
628  liquid_phase.property(MPL::PropertyType::viscosity)
629  .template value<double>(variables, x_position, t, dt);
630 
631  // Set mechanical variables for the intrinsic permeability model
632  // For stress dependent permeability.
633 
634  // For stress dependent permeability.
635  variables[static_cast<int>(MPL::Variable::total_stress)]
636  .emplace<SymmetricTensor>(
638  sigma_total));
639 
640  variables[static_cast<int>(
642  ip_data.material_state_variables->getEquivalentPlasticStrain();
643 
644  auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
645  medium->property(MPL::PropertyType::permeability)
646  .value(variables, x_position, t, dt));
647 
648  GlobalDimMatrixType const Ki_over_mu = K_intrinsic / mu;
649  GlobalDimMatrixType const rho_Ki_over_mu = rho_LR * Ki_over_mu;
650 
651  laplace_p.noalias() +=
652  dNdx.transpose() * k_rel * rho_Ki_over_mu * dNdx * w;
653 
654  auto const beta_LR = 1 / rho_LR *
655  liquid_phase.property(MPL::PropertyType::density)
656  .template dValue<double>(
657  variables, MPL::Variable::phase_pressure,
658  x_position, t, dt);
659 
660  const double alphaB_minus_phi = alpha - phi;
661  double const a0 = alphaB_minus_phi * beta_SR;
662  double const specific_storage_a_p = S_L * (phi * beta_LR + S_L * a0);
663  double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
664 
665  // Note: d beta_LR/d p is omitted because it is a small value.
666  double const dspecific_storage_a_p_dp_cap =
667  dS_L_dp_cap * (phi * beta_LR + 2 * S_L * a0);
668  double const dspecific_storage_a_S_dp_cap =
669  -a0 * (S_L + p_cap_ip * dS_L_dp_cap);
670 
671  storage_p_a_p.noalias() +=
672  N.transpose() * rho_LR * specific_storage_a_p * N * w;
673 
674  storage_p_a_S.noalias() -= N.transpose() * rho_LR *
675  specific_storage_a_S * DeltaS_L_Deltap_cap *
676  N * w;
677 
678  local_Jac
679  .template block<pressure_size, pressure_size>(pressure_index,
681  .noalias() += N.transpose() * p_cap_dot_ip * rho_LR *
682  dspecific_storage_a_p_dp_cap * N * w;
683 
684  storage_p_a_S_Jpp.noalias() -=
685  N.transpose() * rho_LR *
686  ((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap +
687  specific_storage_a_S * dS_L_dp_cap) /
688  dt * N * w;
689 
690  local_Jac
691  .template block<pressure_size, pressure_size>(pressure_index,
693  .noalias() -= N.transpose() * rho_LR * dS_L_dp_cap * alpha *
694  identity2.transpose() * B * u_dot * N * w;
695 
696  double const dk_rel_dS_L =
698  .template dValue<double>(variables,
699  MPL::Variable::liquid_saturation,
700  x_position, t, dt);
701  GlobalDimVectorType const grad_p_cap = -dNdx * p_L;
702  local_Jac
703  .template block<pressure_size, pressure_size>(pressure_index,
705  .noalias() += dNdx.transpose() * rho_Ki_over_mu * grad_p_cap *
706  dk_rel_dS_L * dS_L_dp_cap * N * w;
707 
708  local_Jac
709  .template block<pressure_size, pressure_size>(pressure_index,
711  .noalias() += dNdx.transpose() * rho_LR * rho_Ki_over_mu * b *
712  dk_rel_dS_L * dS_L_dp_cap * N * w;
713 
714  local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
715  dNdx.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
716 
717  //
718  // pressure equation, temperature part, thermal expansion.
719  //
720  {
721  double const fluid_volumetric_thermal_expansion =
723  liquid_phase, variables, rho_LR, x_position, t, dt);
724 
725  const double eff_thermal_expansion =
726  alphaB_minus_phi *
727  Invariants::trace(solid_linear_thermal_expansivity_vector) +
728  fluid_volumetric_thermal_expansion;
729 
730  M_pT.noalias() -=
731  N.transpose() * S_L * rho_LR * eff_thermal_expansion * N * w;
732  }
733 
734  //
735  // temperature equation.
736  //
737  {
738  auto const specific_heat_capacity_fluid =
739  liquid_phase
741  .template value<double>(variables, x_position, t, dt);
742 
743  auto const specific_heat_capacity_solid =
744  solid_phase
747  .template value<double>(variables, x_position, t, dt);
748 
749  M_TT.noalias() +=
750  w *
751  (rho_SR * specific_heat_capacity_solid * (1 - phi) +
752  (S_L * rho_LR * specific_heat_capacity_fluid) * phi) *
753  N.transpose() * N;
754 
755  auto const thermal_conductivity =
756  MaterialPropertyLib::formEigenTensor<DisplacementDim>(
757  medium
760  .value(variables, x_position, t, dt));
761 
762  GlobalDimVectorType const velocity_L = GlobalDimVectorType(
763  -Ki_over_mu * k_rel * (dNdx * p_L - rho_LR * b));
764 
765  K_TT.noalias() += (dNdx.transpose() * thermal_conductivity * dNdx +
766  N.transpose() * velocity_L.transpose() * dNdx *
767  rho_LR * specific_heat_capacity_fluid) *
768  w;
769 
770  //
771  // temperature equation, pressure part
772  //
773  K_Tp.noalias() -= rho_LR * specific_heat_capacity_fluid *
774  N.transpose() * (dNdx * T).transpose() * k_rel *
775  Ki_over_mu * dNdx * w;
776  K_Tp.noalias() -= rho_LR * specific_heat_capacity_fluid *
777  N.transpose() * velocity_L.dot(dNdx * T) / k_rel *
778  dk_rel_dS_L * dS_L_dp_cap * N * w;
779  }
780  }
781 
782  if (process_data_.apply_mass_lumping)
783  {
784  storage_p_a_p = storage_p_a_p.colwise().sum().eval().asDiagonal();
785  storage_p_a_S = storage_p_a_S.colwise().sum().eval().asDiagonal();
786  storage_p_a_S_Jpp =
787  storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
788  }
789 
790  //
791  // -- Jacobian
792  //
793  // temperature equation.
794  local_Jac
795  .template block<temperature_size, temperature_size>(temperature_index,
797  .noalias() += M_TT / dt + K_TT;
798  // temperature equation, pressure part
799  local_Jac
800  .template block<temperature_size, pressure_size>(temperature_index,
802  .noalias() += K_Tp;
803 
804  // pressure equation, pressure part.
805  local_Jac
806  .template block<pressure_size, pressure_size>(pressure_index,
808  .noalias() += laplace_p + storage_p_a_p / dt + storage_p_a_S_Jpp;
809 
810  // pressure equation, temperature part (contributed by thermal expansion).
811  local_Jac
812  .template block<pressure_size, temperature_size>(pressure_index,
814  .noalias() += M_pT / dt;
815 
816  // pressure equation, displacement part.
817  local_Jac
818  .template block<pressure_size, displacement_size>(pressure_index,
820  .noalias() = Kpu / dt;
821 
822  //
823  // -- Residual
824  //
825  // temperature equation
826  local_rhs.template segment<temperature_size>(temperature_index).noalias() -=
827  M_TT * T_dot + K_TT * T;
828 
829  // pressure equation
830  local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
831  laplace_p * p_L + (storage_p_a_p + storage_p_a_S) * p_L_dot +
832  Kpu * u_dot + M_pT * T_dot;
833 }
#define OGS_FATAL(...)
Definition: Error.h:26
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition: BMatrixPolicy.h:55
@ rho
density. For some models, rho substitutes p
double getLiquidThermalExpansivity(Phase const &phase, VariableArray const &vars, const double density, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
Definition: VariableType.h:108
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
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:23
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:32
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
void shapeFunctionInterpolate(const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
Definition: Interpolation.h:79
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
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
static double trace(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Trace of the corresponding tensor.

References MaterialPropertyLib::alpha, MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::bishops_effective_stress, MaterialPropertyLib::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, MaterialPropertyLib::equivalent_plastic_strain, MaterialPropertyLib::getLiquidThermalExpansivity(), NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::mechanical_strain, OGS_FATAL, MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, MaterialPropertyLib::relative_permeability, MaterialPropertyLib::saturation, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), NumLib::shapeFunctionInterpolate(), MaterialPropertyLib::specific_heat_capacity, MaterialPropertyLib::swelling_stress_rate, MaterialPropertyLib::temperature, MaterialPropertyLib::thermal_conductivity, MaterialPropertyLib::thermal_expansivity, MaterialPropertyLib::transport_porosity, and MaterialPropertyLib::viscosity.

◆ computeSecondaryVariableConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
void ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::computeSecondaryVariableConcrete ( double const  t,
double const  dt,
Eigen::VectorXd const &  local_x,
Eigen::VectorXd const &  local_x_dot 
)
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 1061 of file ThermoRichardsMechanicsFEM-impl.h.

1065 {
1066  auto const T =
1067  local_x.template segment<temperature_size>(temperature_index);
1068  auto const p_L = local_x.template segment<pressure_size>(pressure_index);
1069  auto const u =
1070  local_x.template segment<displacement_size>(displacement_index);
1071 
1072  auto const T_dot =
1073  local_x_dot.template segment<temperature_size>(temperature_index);
1074  auto const p_L_dot =
1075  local_x_dot.template segment<pressure_size>(pressure_index);
1076  auto const u_dot =
1077  local_x_dot.template segment<displacement_size>(displacement_index);
1078 
1079  auto const& identity2 = MathLib::KelvinVector::Invariants<
1081  DisplacementDim)>::identity2;
1082 
1083  auto const& medium = process_data_.media_map->getMedium(element_.getID());
1084  auto const& liquid_phase = medium->phase("AqueousLiquid");
1085  auto const& solid_phase = medium->phase("Solid");
1086  MPL::VariableArray variables;
1087  MPL::VariableArray variables_prev;
1088 
1089  ParameterLib::SpatialPosition x_position;
1090  x_position.setElementID(element_.getID());
1091 
1092  unsigned const n_integration_points =
1093  integration_method_.getNumberOfPoints();
1094 
1095  double saturation_avg = 0;
1096  double porosity_avg = 0;
1097 
1099  KV sigma_avg = KV::Zero();
1100 
1101  for (unsigned ip = 0; ip < n_integration_points; ip++)
1102  {
1103  x_position.setIntegrationPoint(ip);
1104 
1105  auto& ip_data = ip_data_[ip];
1106  // N is used for both p and T variables
1107  auto const& N = ip_data.N_p;
1108  auto const& N_u = ip_data.N_u;
1109  auto const& dNdx_u = ip_data.dNdx_u;
1110 
1111  auto const x_coord =
1112  NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
1114  element_, N_u);
1115  auto const B =
1116  LinearBMatrix::computeBMatrix<DisplacementDim,
1117  ShapeFunctionDisplacement::NPOINTS,
1118  typename BMatricesType::BMatrixType>(
1119  dNdx_u, N_u, x_coord, is_axially_symmetric_);
1120 
1121  double T_ip;
1123  double T_dot_ip;
1124  NumLib::shapeFunctionInterpolate(T_dot, N, T_dot_ip);
1125  double const dT = T_dot_ip * dt;
1126 
1127  double p_cap_ip;
1128  NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
1129 
1130  double p_cap_dot_ip;
1131  NumLib::shapeFunctionInterpolate(-p_L_dot, N, p_cap_dot_ip);
1132 
1133  variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
1134  p_cap_ip;
1135  variables[static_cast<int>(MPL::Variable::phase_pressure)] = -p_cap_ip;
1136 
1137  variables[static_cast<int>(MPL::Variable::temperature)] = T_ip;
1138 
1139  auto& eps = ip_data.eps;
1140  eps.noalias() = B * u;
1141  auto& eps_m = ip_data.eps_m;
1142  auto& S_L = ip_data.saturation;
1143  auto const S_L_prev = ip_data.saturation_prev;
1144  S_L = medium->property(MPL::PropertyType::saturation)
1145  .template value<double>(variables, x_position, t, dt);
1146  variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
1147  variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
1148  S_L_prev;
1149 
1150  auto const chi = [medium, x_position, t, dt](double const S_L)
1151  {
1152  MPL::VariableArray vs;
1153  vs.fill(std::numeric_limits<double>::quiet_NaN());
1154  vs[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
1155  return medium->property(MPL::PropertyType::bishops_effective_stress)
1156  .template value<double>(vs, x_position, t, dt);
1157  };
1158  double const chi_S_L = chi(S_L);
1159  double const chi_S_L_prev = chi(S_L_prev);
1160 
1161  auto const alpha =
1162  medium->property(MPL::PropertyType::biot_coefficient)
1163  .template value<double>(variables, x_position, t, dt);
1164 
1165  double const T_ip_prev = T_ip - dT;
1166  auto const C_el = ip_data.computeElasticTangentStiffness(
1167  t, x_position, dt, T_ip_prev, T_ip);
1168 
1169  auto const beta_SR =
1170  (1 - alpha) /
1171  ip_data.solid_material.getBulkModulus(t, x_position, &C_el);
1172  variables[static_cast<int>(MPL::Variable::grain_compressibility)] =
1173  beta_SR;
1174 
1175  variables[static_cast<int>(MPL::Variable::effective_pore_pressure)] =
1176  -chi_S_L * p_cap_ip;
1177  variables_prev[static_cast<int>(
1178  MPL::Variable::effective_pore_pressure)] =
1179  -chi_S_L_prev * (p_cap_ip - p_cap_dot_ip * dt);
1180 
1181  // Set volumetric strain rate for the general case without swelling.
1182  variables[static_cast<int>(MPL::Variable::volumetric_strain)]
1183  .emplace<double>(Invariants::trace(eps));
1184  variables_prev[static_cast<int>(MPL::Variable::volumetric_strain)]
1185  .emplace<double>(Invariants::trace(B * (u - u_dot * dt)));
1186 
1187  auto& phi = ip_data.porosity;
1188  { // Porosity update
1189  variables_prev[static_cast<int>(MPL::Variable::porosity)] =
1190  ip_data.porosity_prev;
1191  phi = medium->property(MPL::PropertyType::porosity)
1192  .template value<double>(variables, variables_prev,
1193  x_position, t, dt);
1194  variables[static_cast<int>(MPL::Variable::porosity)] = phi;
1195  }
1196 
1197  // Swelling and possibly volumetric strain rate update.
1198  if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
1199  {
1200  auto& sigma_sw = ip_data.sigma_sw;
1201  auto const& sigma_sw_prev = ip_data.sigma_sw_prev;
1202 
1203  // If there is swelling, compute it. Update volumetric strain rate,
1204  // s.t. it corresponds to the mechanical part only.
1205  sigma_sw = sigma_sw_prev;
1206 
1207  using DimMatrix = Eigen::Matrix<double, 3, 3>;
1208  auto const sigma_sw_dot =
1209  MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
1210  solid_phase
1212  .template value<DimMatrix>(variables, variables_prev,
1213  x_position, t, dt));
1214  sigma_sw += sigma_sw_dot * dt;
1215 
1216  // !!! Misusing volumetric strain for mechanical volumetric
1217  // strain just to update the transport porosity !!!
1218  std::get<double>(variables[static_cast<int>(
1219  MPL::Variable::volumetric_strain)]) +=
1220  identity2.transpose() * C_el.inverse() * sigma_sw;
1221  std::get<double>(variables_prev[static_cast<int>(
1222  MPL::Variable::volumetric_strain)]) +=
1223  identity2.transpose() * C_el.inverse() * sigma_sw_prev;
1224  }
1225 
1226  if (solid_phase.hasProperty(MPL::PropertyType::transport_porosity))
1227  {
1228  variables_prev[static_cast<int>(
1230  ip_data.transport_porosity_prev;
1231 
1232  ip_data.transport_porosity =
1233  solid_phase.property(MPL::PropertyType::transport_porosity)
1234  .template value<double>(variables, variables_prev,
1235  x_position, t, dt);
1236  variables[static_cast<int>(MPL::Variable::transport_porosity)] =
1237  ip_data.transport_porosity;
1238  }
1239  else
1240  {
1241  variables[static_cast<int>(MPL::Variable::transport_porosity)] =
1242  phi;
1243  }
1244 
1245  auto const mu =
1246  liquid_phase.property(MPL::PropertyType::viscosity)
1247  .template value<double>(variables, x_position, t, dt);
1248  auto const rho_LR =
1249  liquid_phase.property(MPL::PropertyType::density)
1250  .template value<double>(variables, x_position, t, dt);
1251 
1252  // Set mechanical variables for the intrinsic permeability model
1253  // For stress dependent permeability.
1254  {
1255  auto const sigma_total =
1256  (ip_data.sigma_eff + alpha * chi_S_L * identity2 * p_cap_ip)
1257  .eval();
1258  // For stress dependent permeability.
1259  variables[static_cast<int>(MPL::Variable::total_stress)]
1260  .emplace<SymmetricTensor>(
1262  sigma_total));
1263  }
1264 
1265  variables[static_cast<int>(
1267  ip_data.material_state_variables->getEquivalentPlasticStrain();
1268  auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
1269  medium->property(MPL::PropertyType::permeability)
1270  .value(variables, x_position, t, dt));
1271 
1272  double const k_rel =
1274  .template value<double>(variables, x_position, t, dt);
1275 
1276  GlobalDimMatrixType const K_over_mu = k_rel * K_intrinsic / mu;
1277 
1278  auto const& sigma_eff = ip_data.sigma_eff;
1279  double const p_FR = -chi_S_L * p_cap_ip;
1280  // p_SR
1281  variables[static_cast<int>(MPL::Variable::solid_grain_pressure)] =
1282  p_FR - Invariants::trace(sigma_eff) / (3 * (1 - phi));
1283  auto const rho_SR =
1284  solid_phase.property(MPL::PropertyType::density)
1285  .template value<double>(variables, x_position, t, dt);
1286  ip_data.dry_density_solid = (1 - phi) * rho_SR;
1287 
1289  DisplacementDim> const solid_linear_thermal_expansivity_vector =
1290  MPL::formKelvinVectorFromThermalExpansivity<DisplacementDim>(
1291  solid_phase
1292  .property(
1294  .value(variables, x_position, t, dt));
1295 
1297  dthermal_strain = solid_linear_thermal_expansivity_vector * dT;
1298 
1299  auto& eps_prev = ip_data.eps_prev;
1300  auto& eps_m_prev = ip_data.eps_m_prev;
1301  eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
1302 
1303  if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
1304  {
1305  eps_m.noalias() -=
1306  -C_el.inverse() * (ip_data.sigma_sw - ip_data.sigma_sw_prev);
1307  }
1308 
1309  variables[static_cast<int>(
1312  eps_m);
1313 
1314  ip_data.updateConstitutiveRelation(variables, t, x_position, dt,
1315  T_ip_prev);
1316 
1317  auto const& b = process_data_.specific_body_force;
1318 
1319  // Compute the velocity
1320  auto const& dNdx = ip_data.dNdx_p;
1321  ip_data.v_darcy.noalias() =
1322  -K_over_mu * dNdx * p_L + rho_LR * K_over_mu * b;
1323 
1324  saturation_avg += S_L;
1325  porosity_avg += phi;
1326  sigma_avg += sigma_eff;
1327  }
1328  saturation_avg /= n_integration_points;
1329  porosity_avg /= n_integration_points;
1330  sigma_avg /= n_integration_points;
1331 
1332  (*process_data_.element_saturation)[element_.getID()] = saturation_avg;
1333  (*process_data_.element_porosity)[element_.getID()] = porosity_avg;
1334 
1335  Eigen::Map<KV>(&(*process_data_.element_stresses)[element_.getID() *
1336  KV::RowsAtCompileTime]) =
1338 
1340  ShapeFunction, typename ShapeFunctionDisplacement::MeshElement,
1341  DisplacementDim>(element_, is_axially_symmetric_, p_L,
1342  *process_data_.pressure_interpolated);
1344  ShapeFunction, typename ShapeFunctionDisplacement::MeshElement,
1345  DisplacementDim>(element_, is_axially_symmetric_, T,
1346  *process_data_.temperature_interpolated);
1347 }
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::alpha, MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::bishops_effective_stress, MaterialPropertyLib::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), MaterialPropertyLib::density, MaterialPropertyLib::equivalent_plastic_strain, NumLib::interpolateToHigherOrderNodes(), NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::mechanical_strain, MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, MaterialPropertyLib::relative_permeability, MaterialPropertyLib::saturation, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), NumLib::shapeFunctionInterpolate(), MaterialPropertyLib::swelling_stress_rate, MaterialPropertyLib::temperature, MaterialPropertyLib::thermal_expansivity, MaterialPropertyLib::transport_porosity, and MaterialPropertyLib::viscosity.

◆ getEpsilon()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
std::vector< double > ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::getEpsilon
overridevirtual

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

Definition at line 912 of file ThermoRichardsMechanicsFEM-impl.h.

913 {
914  constexpr int kelvin_vector_size =
916 
917  return transposeInPlace<kelvin_vector_size>(
918  [this](std::vector<double>& values)
919  { return getIntPtEpsilon(0, {}, {}, values); });
920 }
std::vector< double > const & getIntPtEpsilon(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

References MathLib::KelvinVector::kelvin_vector_dimensions().

◆ getIntPtDarcyVelocity()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
std::vector< double > const & ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::getIntPtDarcyVelocity ( const double  t,
std::vector< GlobalVector * > const &  x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  dof_table,
std::vector< double > &  cache 
) const
overridevirtual

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

Definition at line 941 of file ThermoRichardsMechanicsFEM-impl.h.

947 {
948  unsigned const n_integration_points =
949  integration_method_.getNumberOfPoints();
950 
951  cache.clear();
952  auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
953  double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
954  cache, DisplacementDim, n_integration_points);
955 
956  for (unsigned ip = 0; ip < n_integration_points; ip++)
957  {
958  cache_matrix.col(ip).noalias() = ip_data_[ip].v_darcy;
959  }
960 
961  return cache;
962 }

References MathLib::createZeroedMatrix().

◆ getIntPtDryDensitySolid()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
std::vector< double > const & ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::getIntPtDryDensitySolid ( const double  t,
std::vector< GlobalVector * > const &  x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  dof_table,
std::vector< double > &  cache 
) const
overridevirtual

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

Definition at line 1046 of file ThermoRichardsMechanicsFEM-impl.h.

1052 {
1055 }
std::vector< double > const & getIntegrationPointScalarData(std::vector< IntegrationPointData, Eigen::aligned_allocator< IntegrationPointData >> const &ip_data, MemberType member, std::vector< double > &cache)

References ProcessLib::getIntegrationPointScalarData().

◆ getIntPtEpsilon()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
std::vector< double > const & ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::getIntPtEpsilon ( const double  t,
std::vector< GlobalVector * > const &  x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  dof_table,
std::vector< double > &  cache 
) const
overridevirtual

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

Definition at line 926 of file ThermoRichardsMechanicsFEM-impl.h.

932 {
933  return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
934  ip_data_, &IpData::eps, cache);
935 }

◆ getIntPtPorosity()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
std::vector< double > const & ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::getIntPtPorosity ( const double  t,
std::vector< GlobalVector * > const &  x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  dof_table,
std::vector< double > &  cache 
) const
overridevirtual

◆ getIntPtSaturation()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
std::vector< double > const & ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::getIntPtSaturation ( const double  t,
std::vector< GlobalVector * > const &  x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  dof_table,
std::vector< double > &  cache 
) const
overridevirtual

◆ getIntPtSigma()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
std::vector< double > const & ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::getIntPtSigma ( const double  t,
std::vector< GlobalVector * > const &  x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  dof_table,
std::vector< double > &  cache 
) const
overridevirtual

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

Definition at line 853 of file ThermoRichardsMechanicsFEM-impl.h.

859 {
860  return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
861  ip_data_, &IpData::sigma_eff, cache);
862 }

◆ getIntPtSwellingStress()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
std::vector< double > const & ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::getIntPtSwellingStress ( const double  t,
std::vector< GlobalVector * > const &  x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  dof_table,
std::vector< double > &  cache 
) const
overridevirtual

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

Definition at line 882 of file ThermoRichardsMechanicsFEM-impl.h.

888 {
889  constexpr int kelvin_vector_size =
891  auto const n_integration_points = ip_data_.size();
892 
893  cache.clear();
894  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
895  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
896  cache, kelvin_vector_size, n_integration_points);
897 
898  for (unsigned ip = 0; ip < n_integration_points; ++ip)
899  {
900  auto const& sigma_sw = ip_data_[ip].sigma_sw;
901  cache_mat.col(ip) =
903  }
904 
905  return cache;
906 }

References MathLib::createZeroedMatrix(), MathLib::KelvinVector::kelvin_vector_dimensions(), and MathLib::KelvinVector::kelvinVectorToSymmetricTensor().

◆ getIntPtTransportPorosity()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
std::vector< double > const & ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::getIntPtTransportPorosity ( const double  t,
std::vector< GlobalVector * > const &  x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  dof_table,
std::vector< double > &  cache 
) const
overridevirtual

◆ getMaterialStateVariablesAt()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::getMaterialStateVariablesAt ( unsigned  integration_point) const
overrideprivatevirtual

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

Definition at line 1363 of file ThermoRichardsMechanicsFEM-impl.h.

1365 {
1366  return *ip_data_[integration_point].material_state_variables;
1367 }

◆ getNumberOfIntegrationPoints()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
unsigned ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::getNumberOfIntegrationPoints
overrideprivatevirtual

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

Definition at line 1353 of file ThermoRichardsMechanicsFEM-impl.h.

1354 {
1355  return integration_method_.getNumberOfPoints();
1356 }

◆ getPorosity()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
std::vector< double > ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::getPorosity
overridevirtual

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

Definition at line 994 of file ThermoRichardsMechanicsFEM-impl.h.

995 {
996  std::vector<double> result;
997  getIntPtPorosity(0, {}, {}, result);
998  return result;
999 }
std::vector< double > const & getIntPtPorosity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

◆ getSaturation()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
std::vector< double > ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::getSaturation
overridevirtual

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

Definition at line 968 of file ThermoRichardsMechanicsFEM-impl.h.

969 {
970  std::vector<double> result;
971  getIntPtSaturation(0, {}, {}, result);
972  return result;
973 }
std::vector< double > const & getIntPtSaturation(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

◆ getShapeMatrix()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
Eigen::Map<const Eigen::RowVectorXd> ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::getShapeMatrix ( const unsigned  integration_point) const
inlineoverridevirtual

Provides the shape matrix at the given integration point.

Implements NumLib::ExtrapolatableElement.

Definition at line 155 of file ThermoRichardsMechanicsFEM.h.

157  {
158  auto const& N_u = secondary_data_.N_u[integration_point];
159 
160  // assumes N is stored contiguously in memory
161  return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
162  }

References ProcessLib::HydroMechanics::SecondaryData< ShapeMatrixType >::N_u, and ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::secondary_data_.

◆ getSigma()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
std::vector< double > ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::getSigma
overridevirtual

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

Definition at line 839 of file ThermoRichardsMechanicsFEM-impl.h.

840 {
841  constexpr int kelvin_vector_size =
843 
844  return transposeInPlace<kelvin_vector_size>(
845  [this](std::vector<double>& values)
846  { return getIntPtSigma(0, {}, {}, values); });
847 }
std::vector< double > const & getIntPtSigma(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

References MathLib::KelvinVector::kelvin_vector_dimensions().

◆ getSwellingStress()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
std::vector< double > ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::getSwellingStress
overridevirtual

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

Definition at line 868 of file ThermoRichardsMechanicsFEM-impl.h.

869 {
870  constexpr int kelvin_vector_size =
872 
873  return transposeInPlace<kelvin_vector_size>(
874  [this](std::vector<double>& values)
875  { return getIntPtSwellingStress(0, {}, {}, values); });
876 }
std::vector< double > const & getIntPtSwellingStress(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

References MathLib::KelvinVector::kelvin_vector_dimensions().

◆ getTransportPorosity()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
std::vector< double > ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::getTransportPorosity
overridevirtual

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

Definition at line 1020 of file ThermoRichardsMechanicsFEM-impl.h.

1021 {
1022  std::vector<double> result;
1023  getIntPtTransportPorosity(0, {}, {}, result);
1024  return result;
1025 }
std::vector< double > const & getIntPtTransportPorosity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

◆ initializeConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
void ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::initializeConcrete ( )
inlineoverridevirtual

Set initial stress from parameter.

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 107 of file ThermoRichardsMechanicsFEM.h.

108  {
109  unsigned const n_integration_points =
110  integration_method_.getNumberOfPoints();
111 
112  for (unsigned ip = 0; ip < n_integration_points; ip++)
113  {
114  auto& ip_data = ip_data_[ip];
115 
117  if (process_data_.initial_stress != nullptr)
118  {
119  ParameterLib::SpatialPosition const x_position{
120  std::nullopt, element_.getID(), ip,
122  ShapeFunctionDisplacement,
124  element_, ip_data.N_u))};
125 
126  ip_data.sigma_eff =
128  DisplacementDim>((*process_data_.initial_stress)(
129  std::numeric_limits<
130  double>::quiet_NaN() /* time independent */,
131  x_position));
132  }
133 
134  ip_data.pushBackState();
135  }
136  }
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::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::element_, MeshLib::Element::getID(), ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::integration_method_, NumLib::interpolateCoordinates(), ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::ip_data_, ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::process_data_, and MathLib::KelvinVector::symmetricTensorToKelvinVector().

◆ postTimestepConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
void ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::postTimestepConcrete ( Eigen::VectorXd const &  ,
double const  ,
double const   
)
inlineoverridevirtual

◆ setInitialConditionsConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
void ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::setInitialConditionsConcrete ( std::vector< double > const &  local_x,
double const  t,
bool const  use_monolithic_scheme,
int const  process_id 
)
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 186 of file ThermoRichardsMechanicsFEM-impl.h.

191 {
192  assert(local_x.size() ==
194 
195  auto const p_L = Eigen::Map<
196  typename ShapeMatricesType::template VectorType<pressure_size> const>(
197  local_x.data() + pressure_index, pressure_size);
198 
199  auto const T = Eigen::Map<typename ShapeMatricesType::template VectorType<
200  temperature_size> const>(local_x.data() + temperature_index,
202 
203  constexpr double dt = std::numeric_limits<double>::quiet_NaN();
204  auto const& medium = process_data_.media_map->getMedium(element_.getID());
205  MPL::VariableArray variables;
206 
208  x_position.setElementID(element_.getID());
209 
210  auto const& solid_phase = medium->phase("Solid");
211 
212  unsigned const n_integration_points =
213  integration_method_.getNumberOfPoints();
214  for (unsigned ip = 0; ip < n_integration_points; ip++)
215  {
216  x_position.setIntegrationPoint(ip);
217 
218  // N is used for both T and p variables.
219  auto const& N = ip_data_[ip].N_p;
220 
221  double p_cap_ip;
222  NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
223 
224  variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
225  p_cap_ip;
226  variables[static_cast<int>(MPL::Variable::phase_pressure)] = -p_cap_ip;
227 
228  double T_ip;
230  variables[static_cast<int>(MPL::Variable::temperature)] = T_ip;
231 
232  ip_data_[ip].saturation_prev =
233  medium->property(MPL::PropertyType::saturation)
234  .template value<double>(variables, x_position, t, dt);
235 
236  // Set eps_m_prev from potentially non-zero eps and sigma_sw from
237  // restart.
238  auto const C_el = ip_data_[ip].computeElasticTangentStiffness(
239  t, x_position, dt, T_ip, T_ip);
240  auto& eps = ip_data_[ip].eps;
241  auto& sigma_sw = ip_data_[ip].sigma_sw;
242  ip_data_[ip].eps_m_prev.noalias() =
243  solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
244  ? eps + C_el.inverse() * sigma_sw
245  : eps;
246  }
247 }

References MaterialPropertyLib::capillary_pressure, MaterialPropertyLib::saturation, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), NumLib::shapeFunctionInterpolate(), MaterialPropertyLib::swelling_stress_rate, and MaterialPropertyLib::temperature.

◆ setIPDataInitialConditions()

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

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

Definition at line 126 of file ThermoRichardsMechanicsFEM-impl.h.

129 {
130  if (integration_order !=
131  static_cast<int>(integration_method_.getIntegrationOrder()))
132  {
133  OGS_FATAL(
134  "Setting integration point initial conditions; The integration "
135  "order of the local assembler for element {:d} is different "
136  "from the integration order in the initial condition.",
137  element_.getID());
138  }
139 
140  if (name == "sigma_ip")
141  {
142  if (process_data_.initial_stress != nullptr)
143  {
144  OGS_FATAL(
145  "Setting initial conditions for stress from integration "
146  "point data and from a parameter '{:s}' is not possible "
147  "simultaneously.",
148  process_data_.initial_stress->name);
149  }
150  return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
151  values, ip_data_, &IpData::sigma_eff);
152  }
153 
154  if (name == "saturation_ip")
155  {
158  }
159  if (name == "porosity_ip")
160  {
163  }
164  if (name == "transport_porosity_ip")
165  {
168  }
169  if (name == "swelling_stress_ip")
170  {
171  return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
172  values, ip_data_, &IpData::sigma_sw);
173  }
174  if (name == "epsilon_ip")
175  {
176  return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
177  values, ip_data_, &IpData::eps);
178  }
179  return 0;
180 }
std::size_t setIntegrationPointScalarData(double const *values, std::vector< IntegrationPointData, Eigen::aligned_allocator< IntegrationPointData >> &ip_data, MemberType member)

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

Member Data Documentation

◆ displacement_index

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
const int ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::displacement_index = 2 * ShapeFunction::NPOINTS
staticprivate

Definition at line 242 of file ThermoRichardsMechanicsFEM.h.

◆ displacement_size

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
const int ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::displacement_size
staticprivate

◆ element_

◆ integration_method_

◆ ip_data_

◆ is_axially_symmetric_

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
bool const ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::is_axially_symmetric_
private

Definition at line 233 of file ThermoRichardsMechanicsFEM.h.

◆ KelvinVectorSize

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
int const ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::KelvinVectorSize
static
Initial value:

Definition at line 69 of file ThermoRichardsMechanicsFEM.h.

◆ pressure_index

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
const int ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::pressure_index = temperature_size
staticprivate

Definition at line 240 of file ThermoRichardsMechanicsFEM.h.

◆ pressure_size

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
const int ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::pressure_size = ShapeFunction::NPOINTS
staticprivate

Definition at line 241 of file ThermoRichardsMechanicsFEM.h.

◆ process_data_

◆ secondary_data_

◆ temperature_index

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
const int ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::temperature_index = 0
staticprivate

Definition at line 238 of file ThermoRichardsMechanicsFEM.h.

◆ temperature_size

template<typename ShapeFunctionDisplacement , typename ShapeFunction , typename IntegrationMethod , int DisplacementDim>
const int ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, IntegrationMethod, DisplacementDim >::temperature_size = ShapeFunction::NPOINTS
staticprivate

Definition at line 239 of file ThermoRichardsMechanicsFEM.h.


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