OGS
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > Class Template Reference

Detailed Description

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
class ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >

Definition at line 49 of file RichardsMechanicsFEM.h.

#include <RichardsMechanicsFEM.h>

Inheritance diagram for ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >:
[legend]
Collaboration diagram for ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >:
[legend]

Public Types

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

Public Member Functions

 RichardsMechanicsLocalAssembler (RichardsMechanicsLocalAssembler const &)=delete
 
 RichardsMechanicsLocalAssembler (RichardsMechanicsLocalAssembler &&)=delete
 
 RichardsMechanicsLocalAssembler (MeshLib::Element const &e, std::size_t const, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, RichardsMechanicsProcessData< DisplacementDim > &process_data)
 
std::size_t setIPDataInitialConditions (std::string_view const name, double const *values, int const integration_order) override
 
void setInitialConditionsConcrete (std::vector< double > const &local_x, double const t, bool const use_monolithic_scheme, int const process_id) override
 
void assemble (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_rhs_data) override
 
void assembleWithJacobian (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &, std::vector< double > &, std::vector< double > &local_rhs_data, std::vector< double > &local_Jac_data) override
 
void assembleWithJacobianForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
 
void initializeConcrete () override
 
void postTimestepConcrete (Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, bool const, int const) override
 
void computeSecondaryVariableConcrete (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev) override
 
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix (const unsigned integration_point) const override
 Provides the shape matrix at the given integration point.
 
std::vector< double > 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 > getMicroSaturation () const override
 
std::vector< double > const & getIntPtMicroSaturation (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 > getMicroPressure () const override
 
std::vector< double > const & getIntPtMicroPressure (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
 
int getMaterialID () const override
 
std::vector< double > getMaterialStateVariableInternalState (std::function< std::span< double >(typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables &)> const &get_values_span, int const &n_components) 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
 
virtual std::size_t setIPDataInitialConditions (std::string_view const name, double const *values, int const integration_order)=0
 
virtual std::vector< double > getSigma () const =0
 
virtual 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 =0
 
virtual std::vector< double > getSwellingStress () const =0
 
virtual 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 =0
 
virtual std::vector< double > getEpsilon () const =0
 
virtual 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 =0
 
virtual 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 =0
 
virtual std::vector< double > getSaturation () const =0
 
virtual 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 =0
 
virtual std::vector< double > getMicroSaturation () const =0
 
virtual std::vector< double > const & getIntPtMicroSaturation (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > getMicroPressure () const =0
 
virtual std::vector< double > const & getIntPtMicroPressure (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > getPorosity () const =0
 
virtual std::vector< double > const & getIntPtPorosity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > getTransportPorosity () const =0
 
virtual 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 =0
 
virtual 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 =0
 
virtual std::vector< double > getMaterialStateVariableInternalState (std::function< std::span< double >(typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables &)> const &get_values_span, int const &n_components) const =0
 
virtual unsigned getNumberOfIntegrationPoints () const =0
 
virtual int getMaterialID () const =0
 
virtual MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt (unsigned) const =0
 
- Public Member Functions inherited from ProcessLib::LocalAssemblerInterface
virtual ~LocalAssemblerInterface ()=default
 
void setInitialConditions (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, bool const use_monolithic_scheme, int const process_id)
 
virtual void initialize (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
 
virtual void preAssemble (double const, double const, std::vector< double > const &)
 
virtual void assemble (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
 
virtual void assembleForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
 
virtual void assembleWithJacobian (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
 
virtual void assembleWithJacobianForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
 
virtual void computeSecondaryVariable (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, int const process_id)
 
virtual void preTimestep (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
 
virtual void postTimestep (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, bool const use_monolithic_scheme, int const process_id)
 
void postNonLinearSolver (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, bool const use_monolithic_scheme, int const process_id)
 
virtual Eigen::Vector3d getFlux (MathLib::Point3d const &, double const, std::vector< double > const &) const
 
virtual Eigen::Vector3d getFlux (MathLib::Point3d const &, double const, std::vector< std::vector< double > > const &) const
 Fits to staggered scheme.
 
- Public Member Functions inherited from NumLib::ExtrapolatableElement
virtual Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix (const unsigned integration_point) const =0
 Provides the shape matrix at the given integration point.
 
virtual ~ExtrapolatableElement ()=default
 

Static Public Attributes

static int const KelvinVectorSize
 
static constexpr auto & N_u_op
 

Private Member Functions

void assembleWithJacobianForDeformationEquations (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
 
void assembleWithJacobianForPressureEquations (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
 
unsigned getNumberOfIntegrationPoints () const override
 
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt (unsigned integration_point) const override
 

Private Attributes

RichardsMechanicsProcessData< DisplacementDim > & _process_data
 
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
 
NumLib::GenericIntegrationMethod const & _integration_method
 
MeshLib::Element const & _element
 
bool const _is_axially_symmetric
 
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
 

Static Private Attributes

static const int pressure_index = 0
 
static const int pressure_size = ShapeFunctionPressure::NPOINTS
 
static const int displacement_index = ShapeFunctionPressure::NPOINTS
 
static const int displacement_size
 

Member Typedef Documentation

◆ BMatricesType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::BMatricesType = BMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>

Definition at line 61 of file RichardsMechanicsFEM.h.

◆ GlobalDimMatrixType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::GlobalDimMatrixType = typename ShapeMatricesTypePressure::GlobalDimMatrixType

Definition at line 58 of file RichardsMechanicsFEM.h.

◆ Invariants

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

Definition at line 72 of file RichardsMechanicsFEM.h.

◆ IpData

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::IpData = IntegrationPointData<BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, ShapeFunctionDisplacement::NPOINTS>

Definition at line 65 of file RichardsMechanicsFEM.h.

◆ KelvinVectorType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::KelvinVectorType = typename BMatricesType::KelvinVectorType

Definition at line 63 of file RichardsMechanicsFEM.h.

◆ ShapeMatricesTypeDisplacement

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ShapeMatricesTypeDisplacement = ShapeMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>

Definition at line 53 of file RichardsMechanicsFEM.h.

◆ ShapeMatricesTypePressure

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ShapeMatricesTypePressure = ShapeMatrixPolicyType<ShapeFunctionPressure, DisplacementDim>

Definition at line 55 of file RichardsMechanicsFEM.h.

◆ SymmetricTensor

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

Definition at line 74 of file RichardsMechanicsFEM.h.

Constructor & Destructor Documentation

◆ RichardsMechanicsLocalAssembler() [1/3]

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

◆ RichardsMechanicsLocalAssembler() [2/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::RichardsMechanicsLocalAssembler ( RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > &&  )
delete

◆ RichardsMechanicsLocalAssembler() [3/3]

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

Definition at line 117 of file RichardsMechanicsFEM-impl.h.

124 : _process_data(process_data),
125 _integration_method(integration_method),
126 _element(e),
127 _is_axially_symmetric(is_axially_symmetric)
128{
129 unsigned const n_integration_points =
131
132 _ip_data.reserve(n_integration_points);
133 _secondary_data.N_u.resize(n_integration_points);
134
135 auto const shape_matrices_u =
136 NumLib::initShapeMatrices<ShapeFunctionDisplacement,
138 DisplacementDim>(e, is_axially_symmetric,
140
141 auto const shape_matrices_p =
142 NumLib::initShapeMatrices<ShapeFunctionPressure,
143 ShapeMatricesTypePressure, DisplacementDim>(
144 e, is_axially_symmetric, _integration_method);
145
146 auto const& solid_material =
148 _process_data.solid_materials, _process_data.material_ids,
149 e.getID());
150
151 auto const& medium = _process_data.media_map.getMedium(_element.getID());
152
154 x_position.setElementID(_element.getID());
155 for (unsigned ip = 0; ip < n_integration_points; ip++)
156 {
157 x_position.setIntegrationPoint(ip);
158 _ip_data.emplace_back(solid_material);
159 auto& ip_data = _ip_data[ip];
160 auto const& sm_u = shape_matrices_u[ip];
161 _ip_data[ip].integration_weight =
163 sm_u.integralMeasure * sm_u.detJ;
164
165 ip_data.N_u = sm_u.N;
166 ip_data.dNdx_u = sm_u.dNdx;
167
168 ip_data.N_p = shape_matrices_p[ip].N;
169 ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
170
171 // Initial porosity. Could be read from integration point data or mesh.
172 ip_data.porosity =
173 medium->property(MPL::porosity)
174 .template initialValue<double>(
175 x_position,
176 std::numeric_limits<
177 double>::quiet_NaN() /* t independent */);
178
179 ip_data.transport_porosity = ip_data.porosity;
180 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
181 {
182 ip_data.transport_porosity =
183 medium->property(MPL::transport_porosity)
184 .template initialValue<double>(
185 x_position,
186 std::numeric_limits<
187 double>::quiet_NaN() /* t independent */);
188 }
189
190 _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
191 }
192}
double getWeight() const
Definition: WeightedPoint.h:80
std::size_t getID() const
Returns the ID of the element.
Definition: Element.h:89
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
RichardsMechanicsProcessData< DisplacementDim > & _process_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
auto & selectSolidConstitutiveRelation(SolidMaterialsMap 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::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_element, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_integration_method, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_process_data, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_secondary_data, MeshLib::Element::getID(), NumLib::GenericIntegrationMethod::getNumberOfPoints(), MathLib::WeightedPoint::getWeight(), NumLib::GenericIntegrationMethod::getWeightedPoint(), NumLib::initShapeMatrices(), ProcessLib::HydroMechanics::SecondaryData< ShapeMatrixType >::N_u, MaterialPropertyLib::porosity, MaterialLib::Solids::selectSolidConstitutiveRelation(), ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), and MaterialPropertyLib::transport_porosity.

Member Function Documentation

◆ assemble()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 358 of file RichardsMechanicsFEM-impl.h.

364{
365 assert(local_x.size() == pressure_size + displacement_size);
366
367 auto p_L =
368 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
369 pressure_size> const>(local_x.data() + pressure_index,
371
372 auto u =
373 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
374 displacement_size> const>(local_x.data() + displacement_index,
376
377 auto p_L_prev =
378 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
379 pressure_size> const>(local_x_prev.data() + pressure_index,
381
382 auto u_prev =
383 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
384 displacement_size> const>(local_x_prev.data() + displacement_index,
386
388 typename ShapeMatricesTypeDisplacement::template MatrixType<
391 local_K_data, displacement_size + pressure_size,
393
395 typename ShapeMatricesTypeDisplacement::template MatrixType<
398 local_M_data, displacement_size + pressure_size,
400
402 typename ShapeMatricesTypeDisplacement::template VectorType<
404 local_rhs_data, displacement_size + pressure_size);
405
406 auto const& identity2 = MathLib::KelvinVector::Invariants<
408 DisplacementDim)>::identity2;
409
410 auto const& medium = _process_data.media_map.getMedium(_element.getID());
411 auto const& liquid_phase = medium->phase("AqueousLiquid");
412 auto const& solid_phase = medium->phase("Solid");
413 MPL::VariableArray variables;
414 MPL::VariableArray variables_prev;
415
417 x_position.setElementID(_element.getID());
418
419 unsigned const n_integration_points =
421 for (unsigned ip = 0; ip < n_integration_points; ip++)
422 {
423 x_position.setIntegrationPoint(ip);
424 auto const& w = _ip_data[ip].integration_weight;
425
426 auto const& N_u = _ip_data[ip].N_u;
427 auto const& dNdx_u = _ip_data[ip].dNdx_u;
428
429 auto const& N_p = _ip_data[ip].N_p;
430 auto const& dNdx_p = _ip_data[ip].dNdx_p;
431
432 auto const x_coord =
433 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
435 _element, N_u);
436 auto const B =
437 LinearBMatrix::computeBMatrix<DisplacementDim,
438 ShapeFunctionDisplacement::NPOINTS,
440 dNdx_u, N_u, x_coord, _is_axially_symmetric);
441
442 auto& eps = _ip_data[ip].eps;
443 auto& eps_m = _ip_data[ip].eps_m;
444 eps.noalias() = B * u;
445
446 auto& S_L = _ip_data[ip].saturation;
447 auto const S_L_prev = _ip_data[ip].saturation_prev;
448
449 double p_cap_ip;
450 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
451
452 double p_cap_prev_ip;
453 NumLib::shapeFunctionInterpolate(-p_L_prev, N_p, p_cap_prev_ip);
454
455 variables.capillary_pressure = p_cap_ip;
456 variables.liquid_phase_pressure = -p_cap_ip;
457 // setting pG to 1 atm
458 // TODO : rewrite equations s.t. p_L = pG-p_cap
459 variables.gas_phase_pressure = 1.0e5;
460
461 auto const temperature =
462 medium->property(MPL::PropertyType::reference_temperature)
463 .template value<double>(variables, x_position, t, dt);
464 variables.temperature = temperature;
465
466 auto const alpha =
467 medium->property(MPL::PropertyType::biot_coefficient)
468 .template value<double>(variables, x_position, t, dt);
469 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
470 t, x_position, dt, temperature);
471
472 auto const beta_SR =
473 (1 - alpha) /
474 _ip_data[ip].solid_material.getBulkModulus(t, x_position, &C_el);
475 variables.grain_compressibility = beta_SR;
476
477 auto const rho_LR =
478 liquid_phase.property(MPL::PropertyType::density)
479 .template value<double>(variables, x_position, t, dt);
480 variables.density = rho_LR;
481
482 auto const& b = _process_data.specific_body_force;
483
484 S_L = medium->property(MPL::PropertyType::saturation)
485 .template value<double>(variables, x_position, t, dt);
486 variables.liquid_saturation = S_L;
487 variables_prev.liquid_saturation = S_L_prev;
488
489 // tangent derivative for Jacobian
490 double const dS_L_dp_cap =
491 medium->property(MPL::PropertyType::saturation)
492 .template dValue<double>(variables,
493 MPL::Variable::capillary_pressure,
494 x_position, t, dt);
495 // secant derivative from time discretization for storage
496 // use tangent, if secant is not available
497 double const DeltaS_L_Deltap_cap =
498 (p_cap_ip == p_cap_prev_ip)
499 ? dS_L_dp_cap
500 : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
501
502 auto const chi = [medium, x_position, t, dt](double const S_L)
503 {
505 vs.liquid_saturation = S_L;
506 return medium->property(MPL::PropertyType::bishops_effective_stress)
507 .template value<double>(vs, x_position, t, dt);
508 };
509 double const chi_S_L = chi(S_L);
510 double const chi_S_L_prev = chi(S_L_prev);
511
512 double const p_FR = -chi_S_L * p_cap_ip;
513 variables.effective_pore_pressure = p_FR;
514 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
515
516 // Set volumetric strain rate for the general case without swelling.
517 variables.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
518 variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
519
520 auto& phi = _ip_data[ip].porosity;
521 { // Porosity update
522 variables_prev.porosity = _ip_data[ip].porosity_prev;
523 phi = medium->property(MPL::PropertyType::porosity)
524 .template value<double>(variables, variables_prev,
525 x_position, t, dt);
526 variables.porosity = phi;
527 }
528
529 if (alpha < phi)
530 {
531 OGS_FATAL(
532 "RichardsMechanics: Biot-coefficient {} is smaller than "
533 "porosity {} in element/integration point {}/{}.",
534 alpha, phi, _element.getID(), ip);
535 }
536
537 // Swelling and possibly volumetric strain rate update.
538 {
539 auto& sigma_sw = _ip_data[ip].sigma_sw;
540 auto const& sigma_sw_prev = _ip_data[ip].sigma_sw_prev;
541
542 // If there is swelling, compute it. Update volumetric strain rate,
543 // s.t. it corresponds to the mechanical part only.
544 sigma_sw = sigma_sw_prev;
545 if (solid_phase.hasProperty(
546 MPL::PropertyType::swelling_stress_rate))
547 {
548 auto const sigma_sw_dot =
549 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
551 solid_phase[MPL::PropertyType::swelling_stress_rate]
552 .value(variables, variables_prev, x_position, t,
553 dt)));
554 sigma_sw += sigma_sw_dot * dt;
555
556 // !!! Misusing volumetric strain for mechanical volumetric
557 // strain just to update the transport porosity !!!
558 variables.volumetric_strain +=
559 identity2.transpose() * C_el.inverse() * sigma_sw;
560 variables_prev.volumetric_strain +=
561 identity2.transpose() * C_el.inverse() * sigma_sw_prev;
562 }
563
564 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
565 {
566 variables_prev.transport_porosity =
567 _ip_data[ip].transport_porosity_prev;
568
569 _ip_data[ip].transport_porosity =
570 medium->property(MPL::PropertyType::transport_porosity)
571 .template value<double>(variables, variables_prev,
572 x_position, t, dt);
573 variables.transport_porosity = _ip_data[ip].transport_porosity;
574 }
575 else
576 {
577 variables.transport_porosity = phi;
578 }
579 }
580
581 double const k_rel =
582 medium->property(MPL::PropertyType::relative_permeability)
583 .template value<double>(variables, x_position, t, dt);
584 auto const mu =
585 liquid_phase.property(MPL::PropertyType::viscosity)
586 .template value<double>(variables, x_position, t, dt);
587
588 auto const& sigma_sw = _ip_data[ip].sigma_sw;
589 auto const& sigma_eff = _ip_data[ip].sigma_eff;
590
591 // Set mechanical variables for the intrinsic permeability model
592 // For stress dependent permeability.
593 {
594 auto const sigma_total =
595 (sigma_eff - alpha * p_FR * identity2).eval();
596
597 // For stress dependent permeability.
598 variables.total_stress.emplace<SymmetricTensor>(
600 sigma_total));
601 }
602
603 variables.equivalent_plastic_strain =
604 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
605
606 auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
607 medium->property(MPL::PropertyType::permeability)
608 .value(variables, x_position, t, dt));
609
610 GlobalDimMatrixType const rho_K_over_mu =
611 K_intrinsic * rho_LR * k_rel / mu;
612
613 //
614 // displacement equation, displacement part
615 //
616 eps_m.noalias() =
617 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
618 ? eps + C_el.inverse() * sigma_sw
619 : eps;
620 variables.mechanical_strain
622 eps_m);
623
624 _ip_data[ip].updateConstitutiveRelation(variables, t, x_position, dt,
625 temperature);
626
627 // p_SR
628 variables.solid_grain_pressure =
629 p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
630 auto const rho_SR =
631 solid_phase.property(MPL::PropertyType::density)
632 .template value<double>(variables, x_position, t, dt);
633
634 //
635 // displacement equation, displacement part
636 //
637 double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
638 rhs.template segment<displacement_size>(displacement_index).noalias() -=
639 (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
640
641 //
642 // pressure equation, pressure part.
643 //
644 auto const beta_LR =
645 1 / rho_LR *
646 liquid_phase.property(MPL::PropertyType::density)
647 .template dValue<double>(variables,
648 MPL::Variable::liquid_phase_pressure,
649 x_position, t, dt);
650
651 double const a0 = S_L * (alpha - phi) * beta_SR;
652 // Volumetric average specific storage of the solid and fluid phases.
653 double const specific_storage =
654 DeltaS_L_Deltap_cap * (p_cap_ip * a0 - phi) +
655 S_L * (phi * beta_LR + a0);
656 M.template block<pressure_size, pressure_size>(pressure_index,
658 .noalias() += N_p.transpose() * rho_LR * specific_storage * N_p * w;
659
660 K.template block<pressure_size, pressure_size>(pressure_index,
662 .noalias() += dNdx_p.transpose() * rho_K_over_mu * dNdx_p * w;
663
664 rhs.template segment<pressure_size>(pressure_index).noalias() +=
665 dNdx_p.transpose() * rho_LR * rho_K_over_mu * b * w;
666
667 //
668 // displacement equation, pressure part
669 //
670 K.template block<displacement_size, pressure_size>(displacement_index,
672 .noalias() -= B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
673
674 //
675 // pressure equation, displacement part.
676 //
677 M.template block<pressure_size, displacement_size>(pressure_index,
679 .noalias() += N_p.transpose() * S_L * rho_LR * alpha *
680 identity2.transpose() * B * w;
681 }
682
683 if (_process_data.apply_mass_lumping)
684 {
685 auto Mpp = M.template block<pressure_size, pressure_size>(
687 Mpp = Mpp.colwise().sum().eval().asDiagonal();
688 }
689}
#define OGS_FATAL(...)
Definition: Error.h:26
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > mechanical_strain
Definition: VariableType.h:182
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > total_stress
Definition: VariableType.h:201
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition: BMatrixPolicy.h:55
typename ShapeMatricesTypePressure::GlobalDimMatrixType GlobalDimMatrixType
Eigen::Matrix< double, KelvinVectorSize, 1 > SymmetricTensor
@ rho
density. For some models, rho substitutes p
template Eigen::Matrix< double, 3, 3 > formEigenTensor< 3 >(MaterialPropertyLib::PropertyDataType const &values)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:49
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:24
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:32
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
Definition: Interpolation.h:27
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
auto eval(Function &f, Tuples &... ts) -> typename detail::GetFunctionReturnType< decltype(&Function::eval)>::type
Definition: Apply.h:267
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 double trace(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Trace of the corresponding tensor.

References MaterialPropertyLib::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::effective_pore_pressure, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor< 3 >(), MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::VariableArray::grain_compressibility, NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, MaterialPropertyLib::VariableArray::mechanical_strain, OGS_FATAL, MaterialPropertyLib::VariableArray::porosity, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::VariableArray::solid_grain_pressure, MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, MaterialPropertyLib::VariableArray::transport_porosity, and MaterialPropertyLib::VariableArray::volumetric_strain.

◆ assembleWithJacobian()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 694 of file RichardsMechanicsFEM-impl.h.

702{
703 assert(local_x.size() == pressure_size + displacement_size);
704
705 auto p_L =
706 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
707 pressure_size> const>(local_x.data() + pressure_index,
709
710 auto u =
711 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
712 displacement_size> const>(local_x.data() + displacement_index,
714
715 auto p_L_prev =
716 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
717 pressure_size> const>(local_x_prev.data() + pressure_index,
719 auto u_prev =
720 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
721 displacement_size> const>(local_x_prev.data() + displacement_index,
723
724 auto local_Jac = MathLib::createZeroedMatrix<
725 typename ShapeMatricesTypeDisplacement::template MatrixType<
728 local_Jac_data, displacement_size + pressure_size,
730
731 auto local_rhs = MathLib::createZeroedVector<
732 typename ShapeMatricesTypeDisplacement::template VectorType<
734 local_rhs_data, displacement_size + pressure_size);
735
736 auto const& identity2 = MathLib::KelvinVector::Invariants<
738 DisplacementDim)>::identity2;
739
741 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
743
744 typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_p =
745 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
747
748 typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_S_Jpp =
749 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
751
752 typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_S =
753 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
755
756 typename ShapeMatricesTypeDisplacement::template MatrixType<
758 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
761
762 typename ShapeMatricesTypeDisplacement::template MatrixType<
764 Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
767
768 auto const& medium = _process_data.media_map.getMedium(_element.getID());
769 auto const& liquid_phase = medium->phase("AqueousLiquid");
770 auto const& solid_phase = medium->phase("Solid");
771 MPL::VariableArray variables;
772 MPL::VariableArray variables_prev;
773
775 x_position.setElementID(_element.getID());
776
777 unsigned const n_integration_points =
779 for (unsigned ip = 0; ip < n_integration_points; ip++)
780 {
781 x_position.setIntegrationPoint(ip);
782 auto const& w = _ip_data[ip].integration_weight;
783
784 auto const& N_u = _ip_data[ip].N_u;
785 auto const& dNdx_u = _ip_data[ip].dNdx_u;
786
787 auto const& N_p = _ip_data[ip].N_p;
788 auto const& dNdx_p = _ip_data[ip].dNdx_p;
789
790 auto const x_coord =
791 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
793 _element, N_u);
794 auto const B =
795 LinearBMatrix::computeBMatrix<DisplacementDim,
796 ShapeFunctionDisplacement::NPOINTS,
798 dNdx_u, N_u, x_coord, _is_axially_symmetric);
799
800 double p_cap_ip;
801 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
802
803 double p_cap_prev_ip;
804 NumLib::shapeFunctionInterpolate(-p_L_prev, N_p, p_cap_prev_ip);
805
806 variables.capillary_pressure = p_cap_ip;
807 variables.liquid_phase_pressure = -p_cap_ip;
808 // setting pG to 1 atm
809 // TODO : rewrite equations s.t. p_L = pG-p_cap
810 variables.gas_phase_pressure = 1.0e5;
811
812 auto const temperature =
813 medium->property(MPL::PropertyType::reference_temperature)
814 .template value<double>(variables, x_position, t, dt);
815 variables.temperature = temperature;
816
817 auto& eps = _ip_data[ip].eps;
818 auto& eps_m = _ip_data[ip].eps_m;
819 eps.noalias() = B * u;
820 auto const& sigma_eff = _ip_data[ip].sigma_eff;
821 auto& S_L = _ip_data[ip].saturation;
822 auto const S_L_prev = _ip_data[ip].saturation_prev;
823 auto const alpha =
824 medium->property(MPL::PropertyType::biot_coefficient)
825 .template value<double>(variables, x_position, t, dt);
826
827 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
828 t, x_position, dt, temperature);
829
830 auto const beta_SR =
831 (1 - alpha) /
832 _ip_data[ip].solid_material.getBulkModulus(t, x_position, &C_el);
833 variables.grain_compressibility = beta_SR;
834
835 auto const rho_LR =
836 liquid_phase.property(MPL::PropertyType::density)
837 .template value<double>(variables, x_position, t, dt);
838 variables.density = rho_LR;
839
840 auto const& b = _process_data.specific_body_force;
841
842 S_L = medium->property(MPL::PropertyType::saturation)
843 .template value<double>(variables, x_position, t, dt);
844 variables.liquid_saturation = S_L;
845 variables_prev.liquid_saturation = S_L_prev;
846
847 // tangent derivative for Jacobian
848 double const dS_L_dp_cap =
849 medium->property(MPL::PropertyType::saturation)
850 .template dValue<double>(variables,
851 MPL::Variable::capillary_pressure,
852 x_position, t, dt);
853 // secant derivative from time discretization for storage
854 // use tangent, if secant is not available
855 double const DeltaS_L_Deltap_cap =
856 (p_cap_ip == p_cap_prev_ip)
857 ? dS_L_dp_cap
858 : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
859
860 auto const chi = [medium, x_position, t, dt](double const S_L)
861 {
863 vs.liquid_saturation = S_L;
864 return medium->property(MPL::PropertyType::bishops_effective_stress)
865 .template value<double>(vs, x_position, t, dt);
866 };
867 double const chi_S_L = chi(S_L);
868 double const chi_S_L_prev = chi(S_L_prev);
869
870 double const p_FR = -chi_S_L * p_cap_ip;
871 variables.effective_pore_pressure = p_FR;
872 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
873
874 // Set volumetric strain rate for the general case without swelling.
875 variables.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
876 variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
877
878 auto& phi = _ip_data[ip].porosity;
879 { // Porosity update
880
881 variables_prev.porosity = _ip_data[ip].porosity_prev;
882 phi = medium->property(MPL::PropertyType::porosity)
883 .template value<double>(variables, variables_prev,
884 x_position, t, dt);
885 variables.porosity = phi;
886 }
887
888 if (alpha < phi)
889 {
890 OGS_FATAL(
891 "RichardsMechanics: Biot-coefficient {} is smaller than "
892 "porosity {} in element/integration point {}/{}.",
893 alpha, phi, _element.getID(), ip);
894 }
895
896 auto const mu =
897 liquid_phase.property(MPL::PropertyType::viscosity)
898 .template value<double>(variables, x_position, t, dt);
899
900 // Swelling and possibly volumetric strain rate update.
901 updateSwellingStressAndVolumetricStrain<DisplacementDim>(
902 _ip_data[ip], *medium, solid_phase, C_el, rho_LR, mu,
903 _process_data.micro_porosity_parameters, alpha, phi, p_cap_ip,
904 variables, variables_prev, x_position, t, dt);
905
906 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
907 {
908 if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
909 {
910 variables_prev.transport_porosity =
911 _ip_data[ip].transport_porosity_prev;
912
913 _ip_data[ip].transport_porosity =
914 medium->property(MPL::PropertyType::transport_porosity)
915 .template value<double>(variables, variables_prev,
916 x_position, t, dt);
917 variables.transport_porosity = _ip_data[ip].transport_porosity;
918 }
919 }
920 else
921 {
922 variables.transport_porosity = phi;
923 }
924
925 double const k_rel =
926 medium->property(MPL::PropertyType::relative_permeability)
927 .template value<double>(variables, x_position, t, dt);
928
929 // Set mechanical variables for the intrinsic permeability model
930 // For stress dependent permeability.
931 {
932 auto const sigma_total =
933 (_ip_data[ip].sigma_eff + alpha * p_FR * identity2).eval();
934 // For stress dependent permeability.
935 variables.total_stress.emplace<SymmetricTensor>(
937 sigma_total));
938 }
939
940 variables.equivalent_plastic_strain =
941 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
942
943 auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
944 medium->property(MPL::PropertyType::permeability)
945 .value(variables, x_position, t, dt));
946
947 GlobalDimMatrixType const rho_Ki_over_mu = K_intrinsic * rho_LR / mu;
948
949 //
950 // displacement equation, displacement part
951 //
952 auto& sigma_sw = _ip_data[ip].sigma_sw;
953
954 eps_m.noalias() =
955 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
956 ? eps + C_el.inverse() * sigma_sw
957 : eps;
958 variables.mechanical_strain
960 eps_m);
961
962 auto C = _ip_data[ip].updateConstitutiveRelation(
963 variables, t, x_position, dt, temperature);
964
965 local_Jac
966 .template block<displacement_size, displacement_size>(
968 .noalias() += B.transpose() * C * B * w;
969
970 // p_SR
971 variables.solid_grain_pressure =
972 p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
973 auto const rho_SR =
974 solid_phase.property(MPL::PropertyType::density)
975 .template value<double>(variables, x_position, t, dt);
976
977 double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
978 local_rhs.template segment<displacement_size>(displacement_index)
979 .noalias() -=
980 (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
981
982 //
983 // displacement equation, pressure part
984 //
985 Kup.noalias() += B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
986
987 auto const dchi_dS_L =
988 medium->property(MPL::PropertyType::bishops_effective_stress)
989 .template dValue<double>(variables,
990 MPL::Variable::liquid_saturation,
991 x_position, t, dt);
992 local_Jac
993 .template block<displacement_size, pressure_size>(
995 .noalias() -= B.transpose() * alpha *
996 (chi_S_L + dchi_dS_L * p_cap_ip * dS_L_dp_cap) *
997 identity2 * N_p * w;
998
999 local_Jac
1000 .template block<displacement_size, pressure_size>(
1002 .noalias() +=
1003 N_u_op(N_u).transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w;
1004
1005 // For the swelling stress with double structure model the corresponding
1006 // Jacobian u-p entry would be required, but it does not improve
1007 // convergence and sometimes worsens it:
1008 // if (medium->hasProperty(MPL::PropertyType::saturation_micro))
1009 // {
1010 // -B.transpose() *
1011 // dsigma_sw_dS_L_m* dS_L_m_dp_cap_m*(p_L_m - p_L_m_prev) /
1012 // (p_cap_ip - p_cap_prev_ip) * N_p* w;
1013 // }
1014 if (!medium->hasProperty(MPL::PropertyType::saturation_micro) &&
1015 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
1016 {
1017 using DimMatrix = Eigen::Matrix<double, 3, 3>;
1018 auto const dsigma_sw_dS_L =
1019 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
1020 solid_phase
1021 .property(MPL::PropertyType::swelling_stress_rate)
1022 .template dValue<DimMatrix>(
1023 variables, variables_prev,
1024 MPL::Variable::liquid_saturation, x_position, t,
1025 dt));
1026 local_Jac
1027 .template block<displacement_size, pressure_size>(
1029 .noalias() +=
1030 B.transpose() * dsigma_sw_dS_L * dS_L_dp_cap * N_p * w;
1031 }
1032 //
1033 // pressure equation, displacement part.
1034 //
1035 if (_process_data.explicit_hm_coupling_in_unsaturated_zone)
1036 {
1037 Kpu.noalias() += N_p.transpose() * chi_S_L_prev * rho_LR * alpha *
1038 identity2.transpose() * B * w;
1039 }
1040 else
1041 {
1042 Kpu.noalias() += N_p.transpose() * S_L * rho_LR * alpha *
1043 identity2.transpose() * B * w;
1044 }
1045
1046 //
1047 // pressure equation, pressure part.
1048 //
1049 laplace_p.noalias() +=
1050 dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
1051
1052 auto const beta_LR =
1053 1 / rho_LR *
1054 liquid_phase.property(MPL::PropertyType::density)
1055 .template dValue<double>(variables,
1056 MPL::Variable::liquid_phase_pressure,
1057 x_position, t, dt);
1058
1059 double const a0 = (alpha - phi) * beta_SR;
1060 double const specific_storage_a_p = S_L * (phi * beta_LR + S_L * a0);
1061 double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
1062
1063 double const dspecific_storage_a_p_dp_cap =
1064 dS_L_dp_cap * (phi * beta_LR + 2 * S_L * a0);
1065 double const dspecific_storage_a_S_dp_cap =
1066 -a0 * (S_L + p_cap_ip * dS_L_dp_cap);
1067
1068 storage_p_a_p.noalias() +=
1069 N_p.transpose() * rho_LR * specific_storage_a_p * N_p * w;
1070
1071 storage_p_a_S.noalias() -= N_p.transpose() * rho_LR *
1072 specific_storage_a_S * DeltaS_L_Deltap_cap *
1073 N_p * w;
1074
1075 local_Jac
1076 .template block<pressure_size, pressure_size>(pressure_index,
1078 .noalias() += N_p.transpose() * (p_cap_ip - p_cap_prev_ip) / dt *
1079 rho_LR * dspecific_storage_a_p_dp_cap * N_p * w;
1080
1081 storage_p_a_S_Jpp.noalias() -=
1082 N_p.transpose() * rho_LR *
1083 ((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap +
1084 specific_storage_a_S * dS_L_dp_cap) /
1085 dt * N_p * w;
1086
1087 if (!_process_data.explicit_hm_coupling_in_unsaturated_zone)
1088 {
1089 local_Jac
1090 .template block<pressure_size, pressure_size>(pressure_index,
1092 .noalias() -= N_p.transpose() * rho_LR * dS_L_dp_cap * alpha *
1093 identity2.transpose() * B * (u - u_prev) / dt *
1094 N_p * w;
1095 }
1096
1097 double const dk_rel_dS_l =
1098 medium->property(MPL::PropertyType::relative_permeability)
1099 .template dValue<double>(variables,
1100 MPL::Variable::liquid_saturation,
1101 x_position, t, dt);
1103 grad_p_cap = -dNdx_p * p_L;
1104 local_Jac
1105 .template block<pressure_size, pressure_size>(pressure_index,
1107 .noalias() += dNdx_p.transpose() * rho_Ki_over_mu * grad_p_cap *
1108 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
1109
1110 local_Jac
1111 .template block<pressure_size, pressure_size>(pressure_index,
1113 .noalias() += dNdx_p.transpose() * rho_LR * rho_Ki_over_mu * b *
1114 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
1115
1116 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
1117 dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
1118
1119 if (medium->hasProperty(MPL::PropertyType::saturation_micro))
1120 {
1121 double const alpha_bar = _process_data.micro_porosity_parameters
1122 ->mass_exchange_coefficient;
1123 auto const p_L_m = _ip_data[ip].liquid_pressure_m;
1124 local_rhs.template segment<pressure_size>(pressure_index)
1125 .noalias() -=
1126 N_p.transpose() * alpha_bar / mu * (-p_cap_ip - p_L_m) * w;
1127
1128 local_Jac
1129 .template block<pressure_size, pressure_size>(pressure_index,
1131 .noalias() += N_p.transpose() * alpha_bar / mu * N_p * w;
1132 if (p_cap_ip != p_cap_prev_ip)
1133 {
1134 double const p_L_m_prev = _ip_data[ip].liquid_pressure_m_prev;
1135 local_Jac
1136 .template block<pressure_size, pressure_size>(
1138 .noalias() += N_p.transpose() * alpha_bar / mu *
1139 (p_L_m - p_L_m_prev) /
1140 (p_cap_ip - p_cap_prev_ip) * N_p * w;
1141 }
1142 }
1143 }
1144
1145 if (_process_data.apply_mass_lumping)
1146 {
1147 storage_p_a_p = storage_p_a_p.colwise().sum().eval().asDiagonal();
1148 storage_p_a_S = storage_p_a_S.colwise().sum().eval().asDiagonal();
1149 storage_p_a_S_Jpp =
1150 storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
1151 }
1152
1153 // pressure equation, pressure part.
1154 local_Jac
1155 .template block<pressure_size, pressure_size>(pressure_index,
1157 .noalias() += laplace_p + storage_p_a_p / dt + storage_p_a_S_Jpp;
1158
1159 // pressure equation, displacement part.
1160 local_Jac
1161 .template block<pressure_size, displacement_size>(pressure_index,
1163 .noalias() = Kpu / dt;
1164
1165 // pressure equation
1166 local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
1167 laplace_p * p_L +
1168 (storage_p_a_p + storage_p_a_S) * (p_L - p_L_prev) / dt +
1169 Kpu * (u - u_prev) / dt;
1170
1171 // displacement equation
1172 local_rhs.template segment<displacement_size>(displacement_index)
1173 .noalias() += Kup * p_L;
1174}
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
VectorType< GlobalDim > GlobalDimVectorType

References MaterialPropertyLib::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::effective_pore_pressure, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::VariableArray::grain_compressibility, NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, MaterialPropertyLib::VariableArray::mechanical_strain, OGS_FATAL, MaterialPropertyLib::VariableArray::porosity, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::VariableArray::solid_grain_pressure, MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, MaterialPropertyLib::VariableArray::transport_porosity, and MaterialPropertyLib::VariableArray::volumetric_strain.

◆ assembleWithJacobianForDeformationEquations()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianForDeformationEquations ( double const  t,
double const  dt,
Eigen::VectorXd const &  local_x,
Eigen::VectorXd const &  local_x_prev,
std::vector< double > &  local_M_data,
std::vector< double > &  local_K_data,
std::vector< double > &  local_b_data,
std::vector< double > &  local_Jac_data 
)
private

Assemble local matrices and vectors arise from the linearized discretized weak form of the residual of the momentum balance equation,

\[ \nabla (\sigma - \alpha_b p \mathrm{I}) = f \]

where \( \sigma\) is the effective stress tensor, \(p\) is the pore pressure, \(\alpha_b\) is the Biot constant, \(\mathrm{I}\) is the identity tensor, and \(f\) is the body force.

Parameters
tTime
dtTime increment
local_xNodal values of \(x\) of an element.
local_x_prevNodal values of \(x_{prev}\) of an element.
local_M_dataMass matrix of an element, which takes the form of \( \int N^T N\mathrm{d}\Omega\). Not used.
local_K_dataLaplacian matrix of an element, which takes the form of \( \int (\nabla N)^T K \nabla N\mathrm{d}\Omega\). Not used.
local_b_dataRight hand side vector of an element.
local_Jac_dataElement Jacobian matrix for the Newton-Raphson method.

Definition at line 1485 of file RichardsMechanicsFEM-impl.h.

1494{
1495 OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
1496}

References OGS_FATAL.

◆ assembleWithJacobianForPressureEquations()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianForPressureEquations ( double const  t,
double const  dt,
Eigen::VectorXd const &  local_x,
Eigen::VectorXd const &  local_x_prev,
std::vector< double > &  local_M_data,
std::vector< double > &  local_K_data,
std::vector< double > &  local_b_data,
std::vector< double > &  local_Jac_data 
)
private

Assemble local matrices and vectors arise from the linearized discretized weak form of the residual of the mass balance equation of single phase flow,

\[ \alpha \cdot{p} - \nabla (K (\nabla p + \rho g \nabla z) + \alpha_b \nabla \cdot \dot{u} = Q \]

where \( alpha\) is a coefficient may depend on storage or the fluid density change, \( \rho\) is the fluid density, \(g\) is the gravitational acceleration, \(z\) is the vertical coordinate, \(u\) is the displacement, and \(Q\) is the source/sink term.

Parameters
tTime
dtTime increment
local_xNodal values of \(x\) of an element.
local_x_prevNodal values of \(x_{prev}\) of an element.
local_M_dataMass matrix of an element, which takes the form of \( \int N^T N\mathrm{d}\Omega\). Not used.
local_K_dataLaplacian matrix of an element, which takes the form of \( \int (\nabla N)^T K \nabla N\mathrm{d}\Omega\). Not used.
local_b_dataRight hand side vector of an element.
local_Jac_dataElement Jacobian matrix for the Newton-Raphson method.

Definition at line 1469 of file RichardsMechanicsFEM-impl.h.

1478{
1479 OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
1480}

References OGS_FATAL.

◆ assembleWithJacobianForStaggeredScheme()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianForStaggeredScheme ( double const  t,
double const  dt,
Eigen::VectorXd const &  local_x,
Eigen::VectorXd const &  local_x_prev,
int const  process_id,
std::vector< double > &  local_M_data,
std::vector< double > &  local_K_data,
std::vector< double > &  local_b_data,
std::vector< double > &  local_Jac_data 
)
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 1501 of file RichardsMechanicsFEM-impl.h.

1507{
1508 // For the equations with pressure
1509 if (process_id == 0)
1510 {
1511 assembleWithJacobianForPressureEquations(t, dt, local_x, local_x_prev,
1512 local_M_data, local_K_data,
1513 local_b_data, local_Jac_data);
1514 return;
1515 }
1516
1517 // For the equations with deformation
1518 assembleWithJacobianForDeformationEquations(t, dt, local_x, local_x_prev,
1519 local_M_data, local_K_data,
1520 local_b_data, local_Jac_data);
1521}
void assembleWithJacobianForDeformationEquations(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
void assembleWithJacobianForPressureEquations(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)

◆ computeSecondaryVariableConcrete()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 1526 of file RichardsMechanicsFEM-impl.h.

1530{
1531 auto p_L = local_x.template segment<pressure_size>(pressure_index);
1532 auto u = local_x.template segment<displacement_size>(displacement_index);
1533
1534 auto p_L_prev =
1535 local_x_prev.template segment<pressure_size>(pressure_index);
1536 auto u_prev =
1537 local_x_prev.template segment<displacement_size>(displacement_index);
1538
1539 auto const& identity2 = MathLib::KelvinVector::Invariants<
1541 DisplacementDim)>::identity2;
1542
1543 auto const& medium = _process_data.media_map.getMedium(_element.getID());
1544 auto const& liquid_phase = medium->phase("AqueousLiquid");
1545 auto const& solid_phase = medium->phase("Solid");
1546 MPL::VariableArray variables;
1547 MPL::VariableArray variables_prev;
1548
1550 x_position.setElementID(_element.getID());
1551
1552 unsigned const n_integration_points =
1554
1555 double saturation_avg = 0;
1556 double porosity_avg = 0;
1557
1559 KV sigma_avg = KV::Zero();
1560
1561 for (unsigned ip = 0; ip < n_integration_points; ip++)
1562 {
1563 x_position.setIntegrationPoint(ip);
1564 auto const& N_p = _ip_data[ip].N_p;
1565 auto const& N_u = _ip_data[ip].N_u;
1566 auto const& dNdx_u = _ip_data[ip].dNdx_u;
1567
1568 auto const x_coord =
1569 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
1571 _element, N_u);
1572 auto const B =
1573 LinearBMatrix::computeBMatrix<DisplacementDim,
1574 ShapeFunctionDisplacement::NPOINTS,
1576 dNdx_u, N_u, x_coord, _is_axially_symmetric);
1577
1578 double p_cap_ip;
1579 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
1580
1581 double p_cap_prev_ip;
1582 NumLib::shapeFunctionInterpolate(-p_L_prev, N_p, p_cap_prev_ip);
1583
1584 variables.capillary_pressure = p_cap_ip;
1585 variables.liquid_phase_pressure = -p_cap_ip;
1586 // setting pG to 1 atm
1587 // TODO : rewrite equations s.t. p_L = pG-p_cap
1588 variables.gas_phase_pressure = 1.0e5;
1589
1590 auto const temperature =
1591 medium->property(MPL::PropertyType::reference_temperature)
1592 .template value<double>(variables, x_position, t, dt);
1593 variables.temperature = temperature;
1594
1595 auto& eps = _ip_data[ip].eps;
1596 eps.noalias() = B * u;
1597 auto& eps_m = _ip_data[ip].eps_m;
1598 auto& S_L = _ip_data[ip].saturation;
1599 auto const S_L_prev = _ip_data[ip].saturation_prev;
1600 S_L = medium->property(MPL::PropertyType::saturation)
1601 .template value<double>(variables, x_position, t, dt);
1602 variables.liquid_saturation = S_L;
1603 variables_prev.liquid_saturation = S_L_prev;
1604
1605 auto const chi = [medium, x_position, t, dt](double const S_L)
1606 {
1608 vs.liquid_saturation = S_L;
1609 return medium->property(MPL::PropertyType::bishops_effective_stress)
1610 .template value<double>(vs, x_position, t, dt);
1611 };
1612 double const chi_S_L = chi(S_L);
1613 double const chi_S_L_prev = chi(S_L_prev);
1614
1615 auto const alpha =
1616 medium->property(MPL::PropertyType::biot_coefficient)
1617 .template value<double>(variables, x_position, t, dt);
1618
1619 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
1620 t, x_position, dt, temperature);
1621
1622 auto const beta_SR =
1623 (1 - alpha) /
1624 _ip_data[ip].solid_material.getBulkModulus(t, x_position, &C_el);
1625 variables.grain_compressibility = beta_SR;
1626
1627 variables.effective_pore_pressure = -chi_S_L * p_cap_ip;
1628 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
1629
1630 // Set volumetric strain rate for the general case without swelling.
1631 variables.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
1632 variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
1633
1634 auto& phi = _ip_data[ip].porosity;
1635 { // Porosity update
1636 variables_prev.porosity = _ip_data[ip].porosity_prev;
1637 phi = medium->property(MPL::PropertyType::porosity)
1638 .template value<double>(variables, variables_prev,
1639 x_position, t, dt);
1640 variables.porosity = phi;
1641 }
1642
1643 auto const rho_LR =
1644 liquid_phase.property(MPL::PropertyType::density)
1645 .template value<double>(variables, x_position, t, dt);
1646 variables.density = rho_LR;
1647 auto const mu =
1648 liquid_phase.property(MPL::PropertyType::viscosity)
1649 .template value<double>(variables, x_position, t, dt);
1650
1651 // Swelling and possibly volumetric strain rate update.
1652 updateSwellingStressAndVolumetricStrain<DisplacementDim>(
1653 _ip_data[ip], *medium, solid_phase, C_el, rho_LR, mu,
1654 _process_data.micro_porosity_parameters, alpha, phi, p_cap_ip,
1655 variables, variables_prev, x_position, t, dt);
1656
1657 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
1658 {
1659 if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
1660 {
1661 variables_prev.transport_porosity =
1662 _ip_data[ip].transport_porosity_prev;
1663
1664 _ip_data[ip].transport_porosity =
1665 medium->property(MPL::PropertyType::transport_porosity)
1666 .template value<double>(variables, variables_prev,
1667 x_position, t, dt);
1668 variables.transport_porosity = _ip_data[ip].transport_porosity;
1669 }
1670 }
1671 else
1672 {
1673 variables.transport_porosity = phi;
1674 }
1675
1676 // Set mechanical variables for the intrinsic permeability model
1677 // For stress dependent permeability.
1678 {
1679 auto const sigma_total = (_ip_data[ip].sigma_eff +
1680 alpha * chi_S_L * identity2 * p_cap_ip)
1681 .eval();
1682 // For stress dependent permeability.
1683 variables.total_stress.emplace<SymmetricTensor>(
1685 sigma_total));
1686 }
1687
1688 variables.equivalent_plastic_strain =
1689 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
1690
1691 auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
1692 medium->property(MPL::PropertyType::permeability)
1693 .value(variables, x_position, t, dt));
1694
1695 double const k_rel =
1696 medium->property(MPL::PropertyType::relative_permeability)
1697 .template value<double>(variables, x_position, t, dt);
1698
1699 GlobalDimMatrixType const K_over_mu = k_rel * K_intrinsic / mu;
1700
1701 auto const& sigma_eff = _ip_data[ip].sigma_eff;
1702 double const p_FR = -chi_S_L * p_cap_ip;
1703 // p_SR
1704 variables.solid_grain_pressure =
1705 p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
1706 auto const rho_SR =
1707 solid_phase.property(MPL::PropertyType::density)
1708 .template value<double>(variables, x_position, t, dt);
1709 _ip_data[ip].dry_density_solid = (1 - phi) * rho_SR;
1710
1711 auto& sigma_sw = _ip_data[ip].sigma_sw;
1712
1713 eps_m.noalias() =
1714 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
1715 ? eps + C_el.inverse() * sigma_sw
1716 : eps;
1717 variables.mechanical_strain
1719 eps_m);
1720
1721 _ip_data[ip].updateConstitutiveRelation(variables, t, x_position, dt,
1722 temperature);
1723
1724 auto const& b = _process_data.specific_body_force;
1725
1726 // Compute the velocity
1727 auto const& dNdx_p = _ip_data[ip].dNdx_p;
1728 _ip_data[ip].v_darcy.noalias() =
1729 -K_over_mu * dNdx_p * p_L + rho_LR * K_over_mu * b;
1730
1731 saturation_avg += S_L;
1732 porosity_avg += phi;
1733 sigma_avg += sigma_eff;
1734 }
1735 saturation_avg /= n_integration_points;
1736 porosity_avg /= n_integration_points;
1737 sigma_avg /= n_integration_points;
1738
1739 (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
1740 (*_process_data.element_porosity)[_element.getID()] = porosity_avg;
1741
1742 Eigen::Map<KV>(&(*_process_data.element_stresses)[_element.getID() *
1743 KV::RowsAtCompileTime]) =
1745
1747 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1748 DisplacementDim>(_element, _is_axially_symmetric, p_L,
1749 *_process_data.pressure_interpolated);
1750}
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::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::effective_pore_pressure, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::VariableArray::grain_compressibility, NumLib::interpolateToHigherOrderNodes(), NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, MaterialPropertyLib::VariableArray::mechanical_strain, MaterialPropertyLib::VariableArray::porosity, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::VariableArray::solid_grain_pressure, MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, MaterialPropertyLib::VariableArray::transport_porosity, and MaterialPropertyLib::VariableArray::volumetric_strain.

◆ getEpsilon()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getEpsilon
overridevirtual

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

Definition at line 1251 of file RichardsMechanicsFEM-impl.h.

1252{
1253 constexpr int kelvin_vector_size =
1255
1256 return transposeInPlace<kelvin_vector_size>(
1257 [this](std::vector<double>& values)
1258 { return getIntPtEpsilon(0, {}, {}, values); });
1259}
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 ShapeFunctionPressure , int DisplacementDim>
std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, 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::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1304 of file RichardsMechanicsFEM-impl.h.

1310{
1311 unsigned const n_integration_points =
1313
1314 cache.clear();
1315 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
1316 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
1317 cache, DisplacementDim, n_integration_points);
1318
1319 for (unsigned ip = 0; ip < n_integration_points; ip++)
1320 {
1321 cache_matrix.col(ip).noalias() = _ip_data[ip].v_darcy;
1322 }
1323
1324 return cache;
1325}

References MathLib::createZeroedMatrix().

◆ getIntPtDryDensitySolid()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, 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::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1455 of file RichardsMechanicsFEM-impl.h.

1461{
1464}
std::vector< double > const & getIntegrationPointScalarData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)

References ProcessLib::getIntegrationPointScalarData().

◆ getIntPtEpsilon()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, 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::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1264 of file RichardsMechanicsFEM-impl.h.

1270{
1271 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
1272 _ip_data, &IpData::eps, cache);
1273}

◆ getIntPtMicroPressure()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtMicroPressure ( const double  t,
std::vector< GlobalVector * > const &  x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  dof_table,
std::vector< double > &  cache 
) const
overridevirtual

◆ getIntPtMicroSaturation()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtMicroSaturation ( const double  t,
std::vector< GlobalVector * > const &  x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  dof_table,
std::vector< double > &  cache 
) const
overridevirtual

◆ getIntPtPorosity()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, 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 ShapeFunctionPressure , int DisplacementDim>
std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, 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 ShapeFunctionPressure , int DisplacementDim>
std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, 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::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1193 of file RichardsMechanicsFEM-impl.h.

1199{
1200 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
1201 _ip_data, &IpData::sigma_eff, cache);
1202}

◆ getIntPtSwellingStress()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, 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::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1221 of file RichardsMechanicsFEM-impl.h.

1227{
1228 constexpr int kelvin_vector_size =
1230 auto const n_integration_points = _ip_data.size();
1231
1232 cache.clear();
1233 auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
1234 double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
1235 cache, kelvin_vector_size, n_integration_points);
1236
1237 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1238 {
1239 auto const& sigma_sw = _ip_data[ip].sigma_sw;
1240 cache_mat.col(ip) =
1242 }
1243
1244 return cache;
1245}

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

◆ getIntPtTransportPorosity()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtTransportPorosity ( const double  t,
std::vector< GlobalVector * > const &  x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const &  dof_table,
std::vector< double > &  cache 
) const
overridevirtual

◆ getMaterialID()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
int ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getMaterialID
overridevirtual

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

Definition at line 1279 of file RichardsMechanicsFEM-impl.h.

1280{
1281 return _process_data.material_ids == nullptr
1282 ? 0
1283 : (*_process_data.material_ids)[_element.getID()];
1284}

◆ getMaterialStateVariableInternalState()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getMaterialStateVariableInternalState ( std::function< std::span< double >(typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables &)> const &  get_values_span,
int const &  n_components 
) const
overridevirtual

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

Definition at line 1289 of file RichardsMechanicsFEM-impl.h.

1295{
1297 _ip_data, &IpData::material_state_variables, get_values_span,
1298 n_components);
1299}
std::vector< double > getIntegrationPointDataMaterialStateVariables(IntegrationPointDataVector const &ip_data_vector, MemberType member, std::function< std::span< double >(MaterialStateVariables &)> get_values_span, int const n_components)
std::unique_ptr< typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables > material_state_variables

References ProcessLib::getIntegrationPointDataMaterialStateVariables().

◆ getMaterialStateVariablesAt()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getMaterialStateVariablesAt ( unsigned  integration_point) const
overrideprivatevirtual

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

Definition at line 1766 of file RichardsMechanicsFEM-impl.h.

1768{
1769 return *_ip_data[integration_point].material_state_variables;
1770}

◆ getMicroPressure()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getMicroPressure
overridevirtual

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

Definition at line 1381 of file RichardsMechanicsFEM-impl.h.

1382{
1383 std::vector<double> result;
1384 getIntPtMicroPressure(0, {}, {}, result);
1385 return result;
1386}
std::vector< double > const & getIntPtMicroPressure(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

◆ getMicroSaturation()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getMicroSaturation
overridevirtual

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

Definition at line 1356 of file RichardsMechanicsFEM-impl.h.

1357{
1358 std::vector<double> result;
1359 getIntPtMicroSaturation(0, {}, {}, result);
1360 return result;
1361}
std::vector< double > const & getIntPtMicroSaturation(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

◆ getNumberOfIntegrationPoints()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
unsigned ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getNumberOfIntegrationPoints
overrideprivatevirtual

◆ getPorosity()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getPorosity
overridevirtual

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

Definition at line 1406 of file RichardsMechanicsFEM-impl.h.

1407{
1408 std::vector<double> result;
1409 getIntPtPorosity(0, {}, {}, result);
1410 return result;
1411}
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 ShapeFunctionPressure , int DisplacementDim>
std::vector< double > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getSaturation
overridevirtual

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

Definition at line 1331 of file RichardsMechanicsFEM-impl.h.

1332{
1333 std::vector<double> result;
1334 getIntPtSaturation(0, {}, {}, result);
1335 return result;
1336}
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 ShapeFunctionPressure , int DisplacementDim>
Eigen::Map< const Eigen::RowVectorXd > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getShapeMatrix ( const unsigned  integration_point) const
inlineoverridevirtual

Provides the shape matrix at the given integration point.

Implements NumLib::ExtrapolatableElement.

Definition at line 178 of file RichardsMechanicsFEM.h.

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

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

◆ getSigma()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getSigma
overridevirtual

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

Definition at line 1180 of file RichardsMechanicsFEM-impl.h.

1181{
1182 constexpr int kelvin_vector_size =
1184
1185 return transposeInPlace<kelvin_vector_size>(
1186 [this](std::vector<double>& values)
1187 { return getIntPtSigma(0, {}, {}, values); });
1188}
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 ShapeFunctionPressure , int DisplacementDim>
std::vector< double > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getSwellingStress
overridevirtual

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

Definition at line 1208 of file RichardsMechanicsFEM-impl.h.

1209{
1210 constexpr int kelvin_vector_size =
1212
1213 return transposeInPlace<kelvin_vector_size>(
1214 [this](std::vector<double>& values)
1215 { return getIntPtSwellingStress(0, {}, {}, values); });
1216}
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 ShapeFunctionPressure , int DisplacementDim>
std::vector< double > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getTransportPorosity
overridevirtual

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

Definition at line 1431 of file RichardsMechanicsFEM-impl.h.

1432{
1433 std::vector<double> result;
1434 getIntPtTransportPorosity(0, {}, {}, result);
1435 return result;
1436}
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 ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::initializeConcrete ( )
inlineoverridevirtual

Set initial stress from parameter.

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 124 of file RichardsMechanicsFEM.h.

125 {
126 unsigned const n_integration_points =
128
129 for (unsigned ip = 0; ip < n_integration_points; ip++)
130 {
131 auto& ip_data = _ip_data[ip];
132
133 ParameterLib::SpatialPosition const x_position{
134 std::nullopt, _element.getID(), ip,
137 ShapeFunctionDisplacement,
139
141 if (_process_data.initial_stress != nullptr)
142 {
143 ip_data.sigma_eff =
145 DisplacementDim>((*_process_data.initial_stress)(
146 std::numeric_limits<
147 double>::quiet_NaN() /* time independent */,
148 x_position));
149 }
150
151 double const t = 0; // TODO (naumov) pass t from top
152 ip_data.solid_material.initializeInternalStateVariables(
153 t, x_position, *ip_data.material_state_variables);
154
155 ip_data.pushBackState();
156 }
157 }
Eigen::Matrix< double, Eigen::MatrixBase< Derived >::RowsAtCompileTime, 1 > symmetricTensorToKelvinVector(Eigen::MatrixBase< Derived > const &v)
Definition: KelvinVector.h:201
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)

References ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_element, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_integration_method, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_process_data, MeshLib::Element::getID(), NumLib::GenericIntegrationMethod::getNumberOfPoints(), NumLib::interpolateCoordinates(), and MathLib::KelvinVector::symmetricTensorToKelvinVector().

◆ postTimestepConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::postTimestepConcrete ( Eigen::VectorXd const &  ,
Eigen::VectorXd const &  ,
double const  ,
double const  ,
bool const  ,
int const   
)
inlineoverridevirtual

◆ setInitialConditionsConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, 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 280 of file RichardsMechanicsFEM-impl.h.

285{
286 assert(local_x.size() == pressure_size + displacement_size);
287
288 auto p_L =
289 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
290 pressure_size> const>(local_x.data() + pressure_index,
292
293 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
294 auto const& medium = _process_data.media_map.getMedium(_element.getID());
295 MPL::VariableArray variables;
296
298 x_position.setElementID(_element.getID());
299
300 auto const& solid_phase = medium->phase("Solid");
301
302 unsigned const n_integration_points =
304 for (unsigned ip = 0; ip < n_integration_points; ip++)
305 {
306 x_position.setIntegrationPoint(ip);
307
308 auto const& N_p = _ip_data[ip].N_p;
309
310 double p_cap_ip;
311 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
312
313 variables.capillary_pressure = p_cap_ip;
314 variables.liquid_phase_pressure = -p_cap_ip;
315 // setting pG to 1 atm
316 // TODO : rewrite equations s.t. p_L = pG-p_cap
317 variables.gas_phase_pressure = 1.0e5;
318 _ip_data[ip].liquid_pressure_m_prev = -p_cap_ip;
319 _ip_data[ip].liquid_pressure_m = -p_cap_ip;
320
321 auto const temperature =
322 medium->property(MPL::PropertyType::reference_temperature)
323 .template value<double>(variables, x_position, t, dt);
324 variables.temperature = temperature;
325
326 _ip_data[ip].saturation_prev =
327 medium->property(MPL::PropertyType::saturation)
328 .template value<double>(variables, x_position, t, dt);
329
330 if (medium->hasProperty(MPL::PropertyType::saturation_micro))
331 {
333 vars.capillary_pressure = p_cap_ip;
334 double const S_L_m =
335 medium->property(MPL::PropertyType::saturation_micro)
336 .template value<double>(vars, x_position, t, dt);
337 _ip_data[ip].saturation_m_prev = S_L_m;
338 }
339
340 // Set eps_m_prev from potentially non-zero eps and sigma_sw from
341 // restart.
342 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
343 t, x_position, dt, temperature);
344 auto& eps = _ip_data[ip].eps;
345 auto& sigma_sw = _ip_data[ip].sigma_sw;
346
347 _ip_data[ip].eps_m_prev.noalias() =
348 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
349 ? eps + C_el.inverse() * sigma_sw
350 : eps;
351 }
352}

References MaterialPropertyLib::VariableArray::capillary_pressure, MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::VariableArray::liquid_phase_pressure, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), NumLib::detail::shapeFunctionInterpolate(), and MaterialPropertyLib::VariableArray::temperature.

◆ setIPDataInitialConditions()

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

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

Definition at line 198 of file RichardsMechanicsFEM-impl.h.

201{
202 if (integration_order !=
203 static_cast<int>(_integration_method.getIntegrationOrder()))
204 {
205 OGS_FATAL(
206 "Setting integration point initial conditions; The integration "
207 "order of the local assembler for element {:d} is different "
208 "from the integration order in the initial condition.",
209 _element.getID());
210 }
211
212 if (name == "sigma")
213 {
214 if (_process_data.initial_stress != nullptr)
215 {
216 OGS_FATAL(
217 "Setting initial conditions for stress from integration "
218 "point data and from a parameter '{:s}' is not possible "
219 "simultaneously.",
220 _process_data.initial_stress->name);
221 }
222 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
223 values, _ip_data, &IpData::sigma_eff);
224 }
225
226 if (name == "saturation")
227 {
230 }
231 if (name == "porosity")
232 {
235 }
236 if (name == "transport_porosity")
237 {
240 }
241 if (name == "swelling_stress")
242 {
243 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
244 values, _ip_data, &IpData::sigma_sw);
245 }
246 if (name == "epsilon")
247 {
248 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
249 values, _ip_data, &IpData::eps);
250 }
251 if (name.starts_with("material_state_variable_"))
252 {
253 name.remove_prefix(24);
254
255 // Using first ip data for solid material. TODO (naumov) move solid
256 // material into element, store only material state in IPs.
257 auto const& internal_variables =
258 _ip_data[0].solid_material.getInternalVariables();
259 if (auto const iv = std::find_if(
260 begin(internal_variables), end(internal_variables),
261 [&name](auto const& iv) { return iv.name == name; });
262 iv != end(internal_variables))
263 {
264 DBUG("Setting material state variable '{:s}'", name);
267 iv->reference);
268 }
269
270 ERR("Could not find variable {:s} in solid material model's internal "
271 "variables.",
272 name);
273 }
274 return 0;
275}
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:30
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:45
std::size_t setIntegrationPointDataMaterialStateVariables(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType member, std::function< std::span< double >(MaterialStateVariables &)> get_values_span)
std::size_t setIntegrationPointScalarData(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType IpData::*const member)

References DBUG(), ERR(), OGS_FATAL, ProcessLib::setIntegrationPointDataMaterialStateVariables(), and ProcessLib::setIntegrationPointScalarData().

Member Data Documentation

◆ _element

◆ _integration_method

◆ _ip_data

◆ _is_axially_symmetric

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
bool const ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_is_axially_symmetric
private

Definition at line 339 of file RichardsMechanicsFEM.h.

◆ _process_data

◆ _secondary_data

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

◆ displacement_index

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
const int ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::displacement_index = ShapeFunctionPressure::NPOINTS
staticprivate

Definition at line 346 of file RichardsMechanicsFEM.h.

◆ displacement_size

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

Definition at line 347 of file RichardsMechanicsFEM.h.

◆ KelvinVectorSize

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
int const ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::KelvinVectorSize
static
Initial value:

Definition at line 70 of file RichardsMechanicsFEM.h.

◆ N_u_op

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

Definition at line 76 of file RichardsMechanicsFEM.h.

◆ pressure_index

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
const int ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::pressure_index = 0
staticprivate

Definition at line 344 of file RichardsMechanicsFEM.h.

◆ pressure_size

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
const int ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::pressure_size = ShapeFunctionPressure::NPOINTS
staticprivate

Definition at line 345 of file RichardsMechanicsFEM.h.


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