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 52 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
using ShapeMatricesTypePressure
using GlobalDimMatrixType
using BMatricesType
using KelvinVectorType = typename BMatricesType::KelvinVectorType
using IpData
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)
void setInitialConditionsConcrete (Eigen::VectorXd const local_x, double const t, 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 > &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_b_data, std::vector< double > &local_Jac_data) override
void initializeConcrete () 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.
Public Member Functions inherited from ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >
 LocalAssemblerInterface (MeshLib::Element const &e, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, RichardsMechanicsProcessData< DisplacementDim > &process_data)
std::size_t setIPDataInitialConditions (std::string_view name, double const *values, int const integration_order)
std::vector< double > getMaterialStateVariableInternalState (std::function< std::span< double >(typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables &)> const &get_values_span, int const &n_components) const
unsigned getNumberOfIntegrationPoints () const
int getMaterialID () const
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt (unsigned integration_point) const
void postTimestepConcrete (Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, int const) override final
Public Member Functions inherited from ProcessLib::LocalAssemblerInterface
virtual ~LocalAssemblerInterface ()=default
virtual void setInitialConditions (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, int const process_id)
virtual void initialize (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
virtual void preAssemble (double const, double const, std::vector< double > const &)
virtual void assembleForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
virtual void computeSecondaryVariable (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, int const process_id)
virtual void preTimestep (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
virtual void postTimestep (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
void postNonLinearSolver (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
virtual Eigen::Vector3d getFlux (MathLib::Point3d const &, double const, std::vector< double > const &) const
virtual Eigen::Vector3d getFlux (MathLib::Point3d const &, double const, std::vector< std::vector< double > > const &) const
 Fits to staggered scheme.
virtual std::optional< VectorSegment > getVectorDeformationSegment () const
Public Member Functions inherited from NumLib::ExtrapolatableElement
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_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_b_data, std::vector< double > &local_Jac_data)

Static Private Member Functions

static void assembleWithJacobianEvalConstitutiveSetting (double const t, double const dt, ParameterLib::SpatialPosition const &x_position, IpData &ip_data, MPL::VariableArray &variables, MPL::VariableArray &variables_prev, MPL::Medium const *const medium, TemperatureData const T_data, CapillaryPressureData< DisplacementDim > const &p_cap_data, ConstitutiveData< DisplacementDim > &CD, StatefulData< DisplacementDim > &SD, StatefulDataPrev< DisplacementDim > const &SD_prev, std::optional< MicroPorosityParameters > const &micro_porosity_parameters, MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material, ProcessLib::ThermoRichardsMechanics::MaterialStateData< DisplacementDim > &material_state_data)
static constexpr auto localDOF (auto const &x)

Private Attributes

std::vector< IpData, Eigen::aligned_allocator< IpData > > ip_data_
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

Additional Inherited Members

Static Public Member Functions inherited from ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >
static auto getReflectionDataForOutput ()
Protected Attributes inherited from ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >
RichardsMechanicsProcessData< DisplacementDim > & process_data_
NumLib::GenericIntegrationMethod const & integration_method_
MeshLib::Element const & element_
bool const is_axially_symmetric_
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material_
std::vector< StatefulData< DisplacementDim > > current_states_
std::vector< StatefulDataPrev< DisplacementDim > > prev_states_
std::vector< ProcessLib::ThermoRichardsMechanics::MaterialStateData< DisplacementDim > > material_states_
std::vector< OutputData< DisplacementDim > > output_data_

Member Typedef Documentation

◆ BMatricesType

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

Definition at line 64 of file RichardsMechanicsFEM.h.

◆ GlobalDimMatrixType

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

Definition at line 61 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 75 of file RichardsMechanicsFEM.h.

◆ IpData

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::IpData
Initial value:
IntegrationPointData<BMatricesType, ShapeMatricesTypeDisplacement,
ShapeMatricesTypePressure, DisplacementDim,
ShapeFunctionDisplacement::NPOINTS>
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::ShapeMatricesTypeDisplacement
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
Definition HydroMechanicsFEM.h:170
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::ShapeMatricesTypePressure
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
Definition HydroMechanicsFEM.h:174
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::BMatricesType
BMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > BMatricesType
Definition HydroMechanicsFEM.h:386

Definition at line 68 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 66 of file RichardsMechanicsFEM.h.

◆ ShapeMatricesTypeDisplacement

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ShapeMatricesTypeDisplacement
Initial value:
ShapeMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>
ShapeMatrixPolicyType
EigenFixedShapeMatrixPolicy< ShapeFunction, GlobalDim > ShapeMatrixPolicyType
Definition ShapeMatrixPolicy.h:141

Definition at line 56 of file RichardsMechanicsFEM.h.

◆ ShapeMatricesTypePressure

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

Definition at line 58 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 77 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

References RichardsMechanicsLocalAssembler().

Referenced by RichardsMechanicsLocalAssembler(), RichardsMechanicsLocalAssembler(), and assembleWithJacobianEvalConstitutiveSetting().

◆ RichardsMechanicsLocalAssembler() [2/3]

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

References RichardsMechanicsLocalAssembler().

◆ 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 129 of file RichardsMechanicsFEM-impl.h.

136 : LocalAssemblerInterface<DisplacementDim>{
137 e, integration_method, is_axially_symmetric, process_data}
138{
139 unsigned const n_integration_points =
140 this->integration_method_.getNumberOfPoints();
141
142 ip_data_.resize(n_integration_points);
143 secondary_data_.N_u.resize(n_integration_points);
144
145 auto const shape_matrices_u =
146 NumLib::initShapeMatrices<ShapeFunctionDisplacement,
147 ShapeMatricesTypeDisplacement,
148 DisplacementDim>(e, is_axially_symmetric,
149 this->integration_method_);
150
151 auto const shape_matrices_p =
152 NumLib::initShapeMatrices<ShapeFunctionPressure,
153 ShapeMatricesTypePressure, DisplacementDim>(
154 e, is_axially_symmetric, this->integration_method_);
155
156 auto const& medium =
157 this->process_data_.media_map.getMedium(this->element_.getID());
158
159 for (unsigned ip = 0; ip < n_integration_points; ip++)
160 {
161 auto& ip_data = ip_data_[ip];
162 auto const& sm_u = shape_matrices_u[ip];
163 ip_data_[ip].integration_weight =
164 this->integration_method_.getWeightedPoint(ip).getWeight() *
165 sm_u.integralMeasure * sm_u.detJ;
166
167 ip_data.N_u = sm_u.N;
168 ip_data.dNdx_u = sm_u.dNdx;
169
170 ParameterLib::SpatialPosition x_position = {
171 std::nullopt, this->element_.getID(),
172 MathLib::Point3d(
173 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
174 ShapeMatricesTypeDisplacement>(
175 this->element_, ip_data.N_u))};
176
177 ip_data.N_p = shape_matrices_p[ip].N;
178 ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
179
180 // Initial porosity. Could be read from integration point data or mesh.
181 auto& porosity =
182 std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(
183 this->current_states_[ip])
184 .phi;
185 porosity = medium->property(MPL::porosity)
186 .template initialValue<double>(
187 x_position,
188 std::numeric_limits<
189 double>::quiet_NaN() /* t independent */);
190
191 auto& transport_porosity =
192 std::get<
193 ProcessLib::ThermoRichardsMechanics::TransportPorosityData>(
194 this->current_states_[ip])
195 .phi;
196 transport_porosity = porosity;
197 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
198 {
199 transport_porosity =
200 medium->property(MPL::transport_porosity)
201 .template initialValue<double>(
202 x_position,
203 std::numeric_limits<
204 double>::quiet_NaN() /* t independent */);
205 }
206
207 secondary_data_.N_u[ip] = shape_matrices_u[ip].N;
208 }
209}
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler::ShapeMatricesTypeDisplacement
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
Definition RichardsMechanicsFEM.h:56
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler::secondary_data_
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > secondary_data_
Definition RichardsMechanicsFEM.h:251
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler::ShapeMatricesTypePressure
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
Definition RichardsMechanicsFEM.h:58
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler::ip_data_
std::vector< IpData, Eigen::aligned_allocator< IpData > > ip_data_
Definition RichardsMechanicsFEM.h:247
ProcessLib::RichardsMechanics::LocalAssemblerInterface::current_states_
std::vector< StatefulData< DisplacementDim > > current_states_
Definition LocalAssemblerInterface.h:185
ProcessLib::RichardsMechanics::LocalAssemblerInterface::LocalAssemblerInterface
LocalAssemblerInterface(MeshLib::Element const &e, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, RichardsMechanicsProcessData< DisplacementDim > &process_data)
Definition LocalAssemblerInterface.h:32
ProcessLib::RichardsMechanics::LocalAssemblerInterface::process_data_
RichardsMechanicsProcessData< DisplacementDim > & process_data_
Definition LocalAssemblerInterface.h:178
ProcessLib::RichardsMechanics::LocalAssemblerInterface::integration_method_
NumLib::GenericIntegrationMethod const & integration_method_
Definition LocalAssemblerInterface.h:179

References ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::LocalAssemblerInterface().

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 376 of file RichardsMechanicsFEM-impl.h.

382{
383 assert(local_x.size() == pressure_size + displacement_size);
384
385 auto const [p_L, u] = localDOF(local_x);
386 auto const [p_L_prev, u_prev] = localDOF(local_x_prev);
387
388 auto K = MathLib::createZeroedMatrix<
389 typename ShapeMatricesTypeDisplacement::template MatrixType<
390 displacement_size + pressure_size,
391 displacement_size + pressure_size>>(
392 local_K_data, displacement_size + pressure_size,
393 displacement_size + pressure_size);
394
395 auto M = MathLib::createZeroedMatrix<
396 typename ShapeMatricesTypeDisplacement::template MatrixType<
397 displacement_size + pressure_size,
398 displacement_size + pressure_size>>(
399 local_M_data, displacement_size + pressure_size,
400 displacement_size + pressure_size);
401
402 auto rhs = MathLib::createZeroedVector<
403 typename ShapeMatricesTypeDisplacement::template VectorType<
404 displacement_size + pressure_size>>(
405 local_rhs_data, displacement_size + pressure_size);
406
407 auto const& identity2 = MathLib::KelvinVector::Invariants<
408 MathLib::KelvinVector::kelvin_vector_dimensions(
409 DisplacementDim)>::identity2;
410
411 auto const& medium =
412 this->process_data_.media_map.getMedium(this->element_.getID());
413 auto const& liquid_phase = medium->phase("AqueousLiquid");
414 auto const& solid_phase = medium->phase("Solid");
415 MPL::VariableArray variables;
416 MPL::VariableArray variables_prev;
417
418 ParameterLib::SpatialPosition x_position;
419 x_position.setElementID(this->element_.getID());
420
421 unsigned const n_integration_points =
422 this->integration_method_.getNumberOfPoints();
423 for (unsigned ip = 0; ip < n_integration_points; ip++)
424 {
425 auto const& w = ip_data_[ip].integration_weight;
426
427 auto const& N_u = ip_data_[ip].N_u;
428 auto const& dNdx_u = ip_data_[ip].dNdx_u;
429
430 auto const& N_p = ip_data_[ip].N_p;
431 auto const& dNdx_p = ip_data_[ip].dNdx_p;
432
433 x_position = {
434 std::nullopt, this->element_.getID(),
435 MathLib::Point3d(
436 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
437 ShapeMatricesTypeDisplacement>(
438 this->element_, N_u))};
439 auto const x_coord = x_position.getCoordinates().value()[0];
440
441 auto const B =
442 LinearBMatrix::computeBMatrix<DisplacementDim,
443 ShapeFunctionDisplacement::NPOINTS,
444 typename BMatricesType::BMatrixType>(
445 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
446
447 auto& eps =
448 std::get<StrainData<DisplacementDim>>(this->current_states_[ip]);
449 eps.eps.noalias() = B * u;
450
451 auto& S_L =
452 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(
453 this->current_states_[ip])
454 .S_L;
455 auto const S_L_prev =
456 std::get<
457 PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
458 this->prev_states_[ip])
459 ->S_L;
460
461 double p_cap_ip;
462 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
463
464 double p_cap_prev_ip;
465 NumLib::shapeFunctionInterpolate(-p_L_prev, N_p, p_cap_prev_ip);
466
467 variables.capillary_pressure = p_cap_ip;
468 variables.liquid_phase_pressure = -p_cap_ip;
469 // setting pG to 1 atm
470 // TODO : rewrite equations s.t. p_L = pG-p_cap
471 variables.gas_phase_pressure = 1.0e5;
472
473 auto const temperature =
474 medium->property(MPL::PropertyType::reference_temperature)
475 .template value<double>(variables, x_position, t, dt);
476 variables.temperature = temperature;
477
478 auto const alpha =
479 medium->property(MPL::PropertyType::biot_coefficient)
480 .template value<double>(variables, x_position, t, dt);
481 auto& SD = this->current_states_[ip];
482 variables.stress =
483 std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
484 DisplacementDim>>(SD)
485 .sigma_eff;
486 // Set mechanical strain temporary to compute tangent stiffness.
487 variables.mechanical_strain
488 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
489 eps.eps);
490 auto const C_el = ip_data_[ip].computeElasticTangentStiffness(
491 variables, t, x_position, dt, this->solid_material_,
492 *this->material_states_[ip].material_state_variables);
493
494 auto const beta_SR = (1 - alpha) / this->solid_material_.getBulkModulus(
495 t, x_position, &C_el);
496 variables.grain_compressibility = beta_SR;
497
498 auto const rho_LR =
499 liquid_phase.property(MPL::PropertyType::density)
500 .template value<double>(variables, x_position, t, dt);
501 variables.density = rho_LR;
502
503 auto const& b = this->process_data_.specific_body_force;
504
505 S_L = medium->property(MPL::PropertyType::saturation)
506 .template value<double>(variables, x_position, t, dt);
507 variables.liquid_saturation = S_L;
508 variables_prev.liquid_saturation = S_L_prev;
509
510 // tangent derivative for Jacobian
511 double const dS_L_dp_cap =
512 medium->property(MPL::PropertyType::saturation)
513 .template dValue<double>(variables,
514 MPL::Variable::capillary_pressure,
515 x_position, t, dt);
516 // secant derivative from time discretization for storage
517 // use tangent, if secant is not available
518 double const DeltaS_L_Deltap_cap =
519 (p_cap_ip == p_cap_prev_ip)
520 ? dS_L_dp_cap
521 : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
522
523 auto const chi = [medium, x_position, t, dt](double const S_L)
524 {
525 MPL::VariableArray vs;
526 vs.liquid_saturation = S_L;
527 return medium->property(MPL::PropertyType::bishops_effective_stress)
528 .template value<double>(vs, x_position, t, dt);
529 };
530 double const chi_S_L = chi(S_L);
531 double const chi_S_L_prev = chi(S_L_prev);
532
533 double const p_FR = -chi_S_L * p_cap_ip;
534 variables.effective_pore_pressure = p_FR;
535 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
536
537 // Set volumetric strain rate for the general case without swelling.
538 variables.volumetric_strain = Invariants::trace(eps.eps);
539 variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
540
541 auto& phi = std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(
542 this->current_states_[ip])
543 .phi;
544 { // Porosity update
545 auto const phi_prev = std::get<PrevState<
546 ProcessLib::ThermoRichardsMechanics::PorosityData>>(
547 this->prev_states_[ip])
548 ->phi;
549 variables_prev.porosity = phi_prev;
550 phi = medium->property(MPL::PropertyType::porosity)
551 .template value<double>(variables, variables_prev,
552 x_position, t, dt);
553 variables.porosity = phi;
554 }
555
556 if (alpha < phi)
557 {
558 OGS_FATAL(
559 "RichardsMechanics: Biot-coefficient {} is smaller than "
560 "porosity {} in element/integration point {}/{}.",
561 alpha, phi, this->element_.getID(), ip);
562 }
563
564 // Swelling and possibly volumetric strain rate update.
565 {
566 auto& sigma_sw =
567 std::get<ProcessLib::ThermoRichardsMechanics::
568 ConstitutiveStress_StrainTemperature::
569 SwellingDataStateful<DisplacementDim>>(
570 this->current_states_[ip])
571 .sigma_sw;
572 auto const& sigma_sw_prev = std::get<PrevState<
573 ProcessLib::ThermoRichardsMechanics::
574 ConstitutiveStress_StrainTemperature::SwellingDataStateful<
575 DisplacementDim>>>(this->prev_states_[ip])
576 ->sigma_sw;
577
578 // If there is swelling, compute it. Update volumetric strain rate,
579 // s.t. it corresponds to the mechanical part only.
580 sigma_sw = sigma_sw_prev;
581 if (solid_phase.hasProperty(
582 MPL::PropertyType::swelling_stress_rate))
583 {
584 auto const sigma_sw_dot =
585 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
586 MPL::formEigenTensor<3>(
587 solid_phase[MPL::PropertyType::swelling_stress_rate]
588 .value(variables, variables_prev, x_position, t,
589 dt)));
590 sigma_sw += sigma_sw_dot * dt;
591
592 variables.volumetric_mechanical_strain =
593 variables.volumetric_strain +
594 identity2.transpose() * C_el.inverse() * sigma_sw;
595 variables_prev.volumetric_mechanical_strain =
596 variables_prev.volumetric_strain +
597 identity2.transpose() * C_el.inverse() * sigma_sw_prev;
598 }
599 else
600 {
601 variables.volumetric_mechanical_strain =
602 variables.volumetric_strain;
603 variables_prev.volumetric_mechanical_strain =
604 variables_prev.volumetric_strain;
605 }
606
607 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
608 {
609 auto& transport_porosity =
610 std::get<ProcessLib::ThermoRichardsMechanics::
611 TransportPorosityData>(
612 this->current_states_[ip])
613 .phi;
614 auto const transport_porosity_prev =
615 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
616 TransportPorosityData>>(
617 this->prev_states_[ip])
618 ->phi;
619 variables_prev.transport_porosity = transport_porosity_prev;
620
621 transport_porosity =
622 medium->property(MPL::PropertyType::transport_porosity)
623 .template value<double>(variables, variables_prev,
624 x_position, t, dt);
625 variables.transport_porosity = transport_porosity;
626 }
627 else
628 {
629 variables.transport_porosity = phi;
630 }
631 }
632
633 double const k_rel =
634 medium->property(MPL::PropertyType::relative_permeability)
635 .template value<double>(variables, x_position, t, dt);
636 auto const mu =
637 liquid_phase.property(MPL::PropertyType::viscosity)
638 .template value<double>(variables, x_position, t, dt);
639
640 auto const& sigma_sw =
641 std::get<ProcessLib::ThermoRichardsMechanics::
642 ConstitutiveStress_StrainTemperature::
643 SwellingDataStateful<DisplacementDim>>(
644 this->current_states_[ip])
645 .sigma_sw;
646 auto const& sigma_eff =
647 std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
648 DisplacementDim>>(this->current_states_[ip])
649 .sigma_eff;
650
651 // Set mechanical variables for the intrinsic permeability model
652 // For stress dependent permeability.
653 {
654 auto const sigma_total =
655 (sigma_eff - alpha * p_FR * identity2).eval();
656
657 // For stress dependent permeability.
658 variables.total_stress.emplace<SymmetricTensor>(
659 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
660 sigma_total));
661 }
662
663 variables.equivalent_plastic_strain =
664 this->material_states_[ip]
665 .material_state_variables->getEquivalentPlasticStrain();
666
667 auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
668 medium->property(MPL::PropertyType::permeability)
669 .value(variables, x_position, t, dt));
670
671 GlobalDimMatrixType const rho_K_over_mu =
672 K_intrinsic * rho_LR * k_rel / mu;
673
674 //
675 // displacement equation, displacement part
676 //
677 {
678 auto& eps_m = std::get<ProcessLib::ConstitutiveRelations::
679 MechanicalStrainData<DisplacementDim>>(
680 this->current_states_[ip])
681 .eps_m;
682 eps_m.noalias() =
683 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
684 ? eps.eps + C_el.inverse() * sigma_sw
685 : eps.eps;
686 variables.mechanical_strain.emplace<
687 MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
688 eps_m);
689 }
690
691 {
692 auto& SD = this->current_states_[ip];
693 auto const& SD_prev = this->prev_states_[ip];
694 auto& sigma_eff =
695 std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
696 DisplacementDim>>(SD);
697 auto const& sigma_eff_prev =
698 std::get<PrevState<ProcessLib::ConstitutiveRelations::
699 EffectiveStressData<DisplacementDim>>>(
700 SD_prev);
701 auto const& eps_m =
702 std::get<ProcessLib::ConstitutiveRelations::
703 MechanicalStrainData<DisplacementDim>>(SD);
704 auto& eps_m_prev =
705 std::get<PrevState<ProcessLib::ConstitutiveRelations::
706 MechanicalStrainData<DisplacementDim>>>(
707 SD_prev);
708
709 ip_data_[ip].updateConstitutiveRelation(
710 variables, t, x_position, dt, temperature, sigma_eff,
711 sigma_eff_prev, eps_m, eps_m_prev, this->solid_material_,
712 this->material_states_[ip].material_state_variables);
713 }
714
715 // p_SR
716 variables.solid_grain_pressure =
717 p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
718 auto const rho_SR =
719 solid_phase.property(MPL::PropertyType::density)
720 .template value<double>(variables, x_position, t, dt);
721
722 //
723 // displacement equation, displacement part
724 //
725 double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
726 rhs.template segment<displacement_size>(displacement_index).noalias() -=
727 (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
728
729 //
730 // pressure equation, pressure part.
731 //
732 auto const beta_LR =
733 1 / rho_LR *
734 liquid_phase.property(MPL::PropertyType::density)
735 .template dValue<double>(variables,
736 MPL::Variable::liquid_phase_pressure,
737 x_position, t, dt);
738
739 double const a0 = S_L * (alpha - phi) * beta_SR;
740 // Volumetric average specific storage of the solid and fluid phases.
741 double const specific_storage =
742 DeltaS_L_Deltap_cap * (p_cap_ip * a0 - phi) +
743 S_L * (phi * beta_LR + a0);
744 M.template block<pressure_size, pressure_size>(pressure_index,
745 pressure_index)
746 .noalias() += N_p.transpose() * rho_LR * specific_storage * N_p * w;
747
748 K.template block<pressure_size, pressure_size>(pressure_index,
749 pressure_index)
750 .noalias() += dNdx_p.transpose() * rho_K_over_mu * dNdx_p * w;
751
752 rhs.template segment<pressure_size>(pressure_index).noalias() +=
753 dNdx_p.transpose() * rho_LR * rho_K_over_mu * b * w;
754
755 //
756 // displacement equation, pressure part
757 //
758 K.template block<displacement_size, pressure_size>(displacement_index,
759 pressure_index)
760 .noalias() -= B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
761
762 //
763 // pressure equation, displacement part.
764 //
765 M.template block<pressure_size, displacement_size>(pressure_index,
766 displacement_index)
767 .noalias() += N_p.transpose() * S_L * rho_LR * alpha *
768 identity2.transpose() * B * w;
769 }
770
771 if (this->process_data_.apply_mass_lumping)
772 {
773 auto Mpp = M.template block<pressure_size, pressure_size>(
774 pressure_index, pressure_index);
775 Mpp = Mpp.colwise().sum().eval().asDiagonal();
776 }
777}
OGS_FATAL
#define OGS_FATAL(...)
Definition Error.h:26
ProcessLib::BMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim >::BMatrixType
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition BMatrixPolicy.h:55
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler::SymmetricTensor
Eigen::Matrix< double, KelvinVectorSize, 1 > SymmetricTensor
Definition RichardsMechanicsFEM.h:77
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler::GlobalDimMatrixType
typename ShapeMatricesTypePressure::GlobalDimMatrixType GlobalDimMatrixType
Definition RichardsMechanicsFEM.h:61
MaterialPropertyLib::formEigenTensor
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Definition FormEigenTensor.cpp:126
MathLib::KelvinVector::kelvinVectorToSymmetricTensor
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Definition KelvinVector.cpp:142
MathLib::KelvinVector::kelvin_vector_dimensions
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition KelvinVector.h:24
MathLib::KelvinVector::tensorToKelvin
KelvinVectorType< DisplacementDim > tensorToKelvin(Eigen::Matrix< double, 3, 3 > const &m)
NumLib::detail::shapeFunctionInterpolate
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
Definition Interpolation.h:27
MathLib::KelvinVector::Invariants::trace
static double trace(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Trace of the corresponding tensor.
Definition KelvinVector-impl.h:62
ProcessLib::RichardsMechanics::LocalAssemblerInterface::prev_states_
std::vector< StatefulDataPrev< DisplacementDim > > prev_states_
Definition LocalAssemblerInterface.h:186
ProcessLib::RichardsMechanics::LocalAssemblerInterface::solid_material_
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material_
Definition LocalAssemblerInterface.h:183
ProcessLib::RichardsMechanics::LocalAssemblerInterface::material_states_
std::vector< ProcessLib::ThermoRichardsMechanics::MaterialStateData< DisplacementDim > > material_states_
Definition LocalAssemblerInterface.h:189

References assemble(), MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::bishops_effective_stress, MaterialPropertyLib::capillary_pressure, MaterialPropertyLib::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::current_states_, MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::density, displacement_index, displacement_size, MaterialPropertyLib::VariableArray::effective_pore_pressure, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::element_, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::gas_phase_pressure, ParameterLib::SpatialPosition::getCoordinates(), MaterialPropertyLib::VariableArray::grain_compressibility, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::integration_method_, NumLib::interpolateCoordinates(), ip_data_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::is_axially_symmetric_, MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, localDOF(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::material_states_, MaterialPropertyLib::VariableArray::mechanical_strain, N_u_op, OGS_FATAL, MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, pressure_index, pressure_size, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::prev_states_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::process_data_, MaterialPropertyLib::reference_temperature, MaterialPropertyLib::relative_permeability, MaterialPropertyLib::saturation, ParameterLib::SpatialPosition::setElementID(), NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::VariableArray::solid_grain_pressure, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::solid_material_, MaterialPropertyLib::VariableArray::stress, MaterialPropertyLib::swelling_stress_rate, MaterialPropertyLib::VariableArray::temperature, MathLib::KelvinVector::tensorToKelvin(), MaterialPropertyLib::VariableArray::total_stress, MathLib::KelvinVector::Invariants< KelvinVectorSize >::trace(), MaterialPropertyLib::transport_porosity, MaterialPropertyLib::VariableArray::transport_porosity, MaterialPropertyLib::viscosity, MaterialPropertyLib::VariableArray::volumetric_mechanical_strain, and MaterialPropertyLib::VariableArray::volumetric_strain.

Referenced by assemble().

◆ 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 > & local_rhs_data,
std::vector< double > & local_Jac_data )
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

1101{
1102 assert(local_x.size() == pressure_size + displacement_size);
1103
1104 auto const [p_L, u] = localDOF(local_x);
1105 auto const [p_L_prev, u_prev] = localDOF(local_x_prev);
1106
1107 auto local_Jac = MathLib::createZeroedMatrix<
1108 typename ShapeMatricesTypeDisplacement::template MatrixType<
1109 displacement_size + pressure_size,
1110 displacement_size + pressure_size>>(
1111 local_Jac_data, displacement_size + pressure_size,
1112 displacement_size + pressure_size);
1113
1114 auto local_rhs = MathLib::createZeroedVector<
1115 typename ShapeMatricesTypeDisplacement::template VectorType<
1116 displacement_size + pressure_size>>(
1117 local_rhs_data, displacement_size + pressure_size);
1118
1119 auto const& identity2 = MathLib::KelvinVector::Invariants<
1120 MathLib::KelvinVector::kelvin_vector_dimensions(
1121 DisplacementDim)>::identity2;
1122
1123 typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
1124 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
1125 pressure_size);
1126
1127 typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_p =
1128 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
1129 pressure_size);
1130
1131 typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_S_Jpp =
1132 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
1133 pressure_size);
1134
1135 typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_S =
1136 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
1137 pressure_size);
1138
1139 typename ShapeMatricesTypeDisplacement::template MatrixType<
1140 displacement_size, pressure_size>
1141 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
1142 displacement_size, pressure_size>::Zero(displacement_size,
1143 pressure_size);
1144
1145 typename ShapeMatricesTypeDisplacement::template MatrixType<
1146 pressure_size, displacement_size>
1147 Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
1148 pressure_size, displacement_size>::Zero(pressure_size,
1149 displacement_size);
1150
1151 auto const& medium =
1152 this->process_data_.media_map.getMedium(this->element_.getID());
1153 auto const& liquid_phase = medium->phase("AqueousLiquid");
1154 auto const& solid_phase = medium->phase("Solid");
1155 MPL::VariableArray variables;
1156 MPL::VariableArray variables_prev;
1157
1158 unsigned const n_integration_points =
1159 this->integration_method_.getNumberOfPoints();
1160 for (unsigned ip = 0; ip < n_integration_points; ip++)
1161 {
1162 ConstitutiveData<DisplacementDim> CD;
1163 auto& SD = this->current_states_[ip];
1164 auto const& SD_prev = this->prev_states_[ip];
1165 [[maybe_unused]] auto models = createConstitutiveModels(
1166 this->process_data_, this->solid_material_);
1167
1168 auto const& w = ip_data_[ip].integration_weight;
1169
1170 auto const& N_u = ip_data_[ip].N_u;
1171 auto const& dNdx_u = ip_data_[ip].dNdx_u;
1172
1173 auto const& N_p = ip_data_[ip].N_p;
1174 auto const& dNdx_p = ip_data_[ip].dNdx_p;
1175
1176 ParameterLib::SpatialPosition x_position = {
1177 std::nullopt, this->element_.getID(),
1178 MathLib::Point3d(
1179 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
1180 ShapeMatricesTypeDisplacement>(
1181 this->element_, N_u))};
1182 auto const x_coord = x_position.getCoordinates().value()[0];
1183
1184 auto const B =
1185 LinearBMatrix::computeBMatrix<DisplacementDim,
1186 ShapeFunctionDisplacement::NPOINTS,
1187 typename BMatricesType::BMatrixType>(
1188 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
1189
1190 double p_cap_ip;
1191 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
1192
1193 double p_cap_prev_ip;
1194 NumLib::shapeFunctionInterpolate(-p_L_prev, N_p, p_cap_prev_ip);
1195
1196 variables.capillary_pressure = p_cap_ip;
1197 variables.liquid_phase_pressure = -p_cap_ip;
1198 // setting pG to 1 atm
1199 // TODO : rewrite equations s.t. p_L = pG-p_cap
1200 variables.gas_phase_pressure = 1.0e5;
1201
1202 auto const temperature =
1203 medium->property(MPL::PropertyType::reference_temperature)
1204 .template value<double>(variables, x_position, t, dt);
1205 variables.temperature = temperature;
1206
1207 std::get<StrainData<DisplacementDim>>(SD).eps.noalias() = B * u;
1208
1209 assembleWithJacobianEvalConstitutiveSetting(
1210 t, dt, x_position, ip_data_[ip], variables, variables_prev, medium,
1211 TemperatureData{temperature},
1212 CapillaryPressureData<DisplacementDim>{
1213 p_cap_ip, p_cap_prev_ip,
1214 Eigen::Vector<double, DisplacementDim>::Zero()},
1215 CD, SD, SD_prev, this->process_data_.micro_porosity_parameters,
1216 this->solid_material_, this->material_states_[ip]);
1217
1218 {
1219 auto const& C = *std::get<StiffnessTensor<DisplacementDim>>(CD);
1220 local_Jac
1221 .template block<displacement_size, displacement_size>(
1222 displacement_index, displacement_index)
1223 .noalias() += B.transpose() * C * B * w;
1224 }
1225
1226 auto const& b = this->process_data_.specific_body_force;
1227
1228 {
1229 auto const& sigma_eff =
1230 std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
1231 DisplacementDim>>(this->current_states_[ip])
1232 .sigma_eff;
1233 double const rho = *std::get<Density>(CD);
1234 local_rhs.template segment<displacement_size>(displacement_index)
1235 .noalias() -= (B.transpose() * sigma_eff -
1236 N_u_op(N_u).transpose() * rho * b) *
1237 w;
1238 }
1239
1240 //
1241 // displacement equation, pressure part
1242 //
1243
1244 double const alpha =
1245 *std::get<ProcessLib::ThermoRichardsMechanics::BiotData>(CD);
1246 double const dS_L_dp_cap =
1247 std::get<ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv>(
1248 CD)
1249 .dS_L_dp_cap;
1250
1251 {
1252 double const chi_S_L =
1253 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
1254 .chi_S_L;
1255 Kup.noalias() +=
1256 B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
1257 double const dchi_dS_L =
1258 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
1259 .dchi_dS_L;
1260
1261 local_Jac
1262 .template block<displacement_size, pressure_size>(
1263 displacement_index, pressure_index)
1264 .noalias() -= B.transpose() * alpha *
1265 (chi_S_L + dchi_dS_L * p_cap_ip * dS_L_dp_cap) *
1266 identity2 * N_p * w;
1267 }
1268
1269 double const phi =
1270 std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(CD).phi;
1271 double const rho_LR = *std::get<LiquidDensity>(CD);
1272 local_Jac
1273 .template block<displacement_size, pressure_size>(
1274 displacement_index, pressure_index)
1275 .noalias() +=
1276 N_u_op(N_u).transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w;
1277
1278 // For the swelling stress with double structure model the corresponding
1279 // Jacobian u-p entry would be required, but it does not improve
1280 // convergence and sometimes worsens it:
1281 // if (medium->hasProperty(MPL::PropertyType::saturation_micro))
1282 // {
1283 // -B.transpose() *
1284 // dsigma_sw_dS_L_m* dS_L_m_dp_cap_m*(p_L_m - p_L_m_prev) /
1285 // (p_cap_ip - p_cap_prev_ip) * N_p* w;
1286 // }
1287 if (!medium->hasProperty(MPL::PropertyType::saturation_micro) &&
1288 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
1289 {
1290 using DimMatrix = Eigen::Matrix<double, 3, 3>;
1291 auto const dsigma_sw_dS_L =
1292 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
1293 solid_phase
1294 .property(MPL::PropertyType::swelling_stress_rate)
1295 .template dValue<DimMatrix>(
1296 variables, variables_prev,
1297 MPL::Variable::liquid_saturation, x_position, t,
1298 dt));
1299 local_Jac
1300 .template block<displacement_size, pressure_size>(
1301 displacement_index, pressure_index)
1302 .noalias() +=
1303 B.transpose() * dsigma_sw_dS_L * dS_L_dp_cap * N_p * w;
1304 }
1305 //
1306 // pressure equation, displacement part.
1307 //
1308 double const S_L =
1309 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(
1310 this->current_states_[ip])
1311 .S_L;
1312 if (this->process_data_.explicit_hm_coupling_in_unsaturated_zone)
1313 {
1314 double const chi_S_L_prev = std::get<PrevState<
1315 ProcessLib::ThermoRichardsMechanics::BishopsData>>(CD)
1316 ->chi_S_L;
1317 Kpu.noalias() += N_p.transpose() * chi_S_L_prev * rho_LR * alpha *
1318 identity2.transpose() * B * w;
1319 }
1320 else
1321 {
1322 Kpu.noalias() += N_p.transpose() * S_L * rho_LR * alpha *
1323 identity2.transpose() * B * w;
1324 }
1325
1326 //
1327 // pressure equation, pressure part.
1328 //
1329
1330 double const k_rel =
1331 std::get<ProcessLib::ThermoRichardsMechanics::PermeabilityData<
1332 DisplacementDim>>(CD)
1333 .k_rel;
1334 auto const& K_intrinsic =
1335 std::get<ProcessLib::ThermoRichardsMechanics::PermeabilityData<
1336 DisplacementDim>>(CD)
1337 .Ki;
1338 double const mu =
1339 *std::get<ProcessLib::ThermoRichardsMechanics::LiquidViscosityData>(
1340 CD);
1341
1342 GlobalDimMatrixType const rho_Ki_over_mu = K_intrinsic * rho_LR / mu;
1343
1344 laplace_p.noalias() +=
1345 dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
1346
1347 auto const beta_LR =
1348 1 / rho_LR *
1349 liquid_phase.property(MPL::PropertyType::density)
1350 .template dValue<double>(variables,
1351 MPL::Variable::liquid_phase_pressure,
1352 x_position, t, dt);
1353
1354 double const beta_SR =
1355 std::get<
1356 ProcessLib::ThermoRichardsMechanics::SolidCompressibilityData>(
1357 CD)
1358 .beta_SR;
1359 double const a0 = (alpha - phi) * beta_SR;
1360 double const specific_storage_a_p = S_L * (phi * beta_LR + S_L * a0);
1361 double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
1362
1363 double const dspecific_storage_a_p_dp_cap =
1364 dS_L_dp_cap * (phi * beta_LR + 2 * S_L * a0);
1365 double const dspecific_storage_a_S_dp_cap =
1366 -a0 * (S_L + p_cap_ip * dS_L_dp_cap);
1367
1368 storage_p_a_p.noalias() +=
1369 N_p.transpose() * rho_LR * specific_storage_a_p * N_p * w;
1370
1371 double const DeltaS_L_Deltap_cap =
1372 std::get<SaturationSecantDerivative>(CD).DeltaS_L_Deltap_cap;
1373 storage_p_a_S.noalias() -= N_p.transpose() * rho_LR *
1374 specific_storage_a_S * DeltaS_L_Deltap_cap *
1375 N_p * w;
1376
1377 local_Jac
1378 .template block<pressure_size, pressure_size>(pressure_index,
1379 pressure_index)
1380 .noalias() += N_p.transpose() * (p_cap_ip - p_cap_prev_ip) / dt *
1381 rho_LR * dspecific_storage_a_p_dp_cap * N_p * w;
1382
1383 double const S_L_prev =
1384 std::get<
1385 PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
1386 this->prev_states_[ip])
1387 ->S_L;
1388 storage_p_a_S_Jpp.noalias() -=
1389 N_p.transpose() * rho_LR *
1390 ((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap +
1391 specific_storage_a_S * dS_L_dp_cap) /
1392 dt * N_p * w;
1393
1394 if (!this->process_data_.explicit_hm_coupling_in_unsaturated_zone)
1395 {
1396 local_Jac
1397 .template block<pressure_size, pressure_size>(pressure_index,
1398 pressure_index)
1399 .noalias() -= N_p.transpose() * rho_LR * dS_L_dp_cap * alpha *
1400 identity2.transpose() * B * (u - u_prev) / dt *
1401 N_p * w;
1402 }
1403
1404 double const dk_rel_dS_l =
1405 medium->property(MPL::PropertyType::relative_permeability)
1406 .template dValue<double>(variables,
1407 MPL::Variable::liquid_saturation,
1408 x_position, t, dt);
1409 typename ShapeMatricesTypeDisplacement::GlobalDimVectorType const
1410 grad_p_cap = -dNdx_p * p_L;
1411 local_Jac
1412 .template block<pressure_size, pressure_size>(pressure_index,
1413 pressure_index)
1414 .noalias() += dNdx_p.transpose() * rho_Ki_over_mu * grad_p_cap *
1415 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
1416
1417 local_Jac
1418 .template block<pressure_size, pressure_size>(pressure_index,
1419 pressure_index)
1420 .noalias() += dNdx_p.transpose() * rho_LR * rho_Ki_over_mu * b *
1421 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
1422
1423 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
1424 dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
1425
1426 if (medium->hasProperty(MPL::PropertyType::saturation_micro))
1427 {
1428 double const alpha_bar =
1429 this->process_data_.micro_porosity_parameters
1430 ->mass_exchange_coefficient;
1431 auto const p_L_m =
1432 *std::get<MicroPressure>(this->current_states_[ip]);
1433 local_rhs.template segment<pressure_size>(pressure_index)
1434 .noalias() -=
1435 N_p.transpose() * alpha_bar / mu * (-p_cap_ip - p_L_m) * w;
1436
1437 local_Jac
1438 .template block<pressure_size, pressure_size>(pressure_index,
1439 pressure_index)
1440 .noalias() += N_p.transpose() * alpha_bar / mu * N_p * w;
1441 if (p_cap_ip != p_cap_prev_ip)
1442 {
1443 auto const p_L_m_prev = **std::get<PrevState<MicroPressure>>(
1444 this->prev_states_[ip]);
1445 local_Jac
1446 .template block<pressure_size, pressure_size>(
1447 pressure_index, pressure_index)
1448 .noalias() += N_p.transpose() * alpha_bar / mu *
1449 (p_L_m - p_L_m_prev) /
1450 (p_cap_ip - p_cap_prev_ip) * N_p * w;
1451 }
1452 }
1453 }
1454
1455 if (this->process_data_.apply_mass_lumping)
1456 {
1457 storage_p_a_p = storage_p_a_p.colwise().sum().eval().asDiagonal();
1458 storage_p_a_S = storage_p_a_S.colwise().sum().eval().asDiagonal();
1459 storage_p_a_S_Jpp =
1460 storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
1461 }
1462
1463 // pressure equation, pressure part.
1464 local_Jac
1465 .template block<pressure_size, pressure_size>(pressure_index,
1466 pressure_index)
1467 .noalias() += laplace_p + storage_p_a_p / dt + storage_p_a_S_Jpp;
1468
1469 // pressure equation, displacement part.
1470 local_Jac
1471 .template block<pressure_size, displacement_size>(pressure_index,
1472 displacement_index)
1473 .noalias() = Kpu / dt;
1474
1475 // pressure equation
1476 local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
1477 laplace_p * p_L +
1478 (storage_p_a_p + storage_p_a_S) * (p_L - p_L_prev) / dt +
1479 Kpu * (u - u_prev) / dt;
1480
1481 // displacement equation
1482 local_rhs.template segment<displacement_size>(displacement_index)
1483 .noalias() += Kup * p_L;
1484}
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler::assembleWithJacobianEvalConstitutiveSetting
static void assembleWithJacobianEvalConstitutiveSetting(double const t, double const dt, ParameterLib::SpatialPosition const &x_position, IpData &ip_data, MPL::VariableArray &variables, MPL::VariableArray &variables_prev, MPL::Medium const *const medium, TemperatureData const T_data, CapillaryPressureData< DisplacementDim > const &p_cap_data, ConstitutiveData< DisplacementDim > &CD, StatefulData< DisplacementDim > &SD, StatefulDataPrev< DisplacementDim > const &SD_prev, std::optional< MicroPorosityParameters > const &micro_porosity_parameters, MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material, ProcessLib::ThermoRichardsMechanics::MaterialStateData< DisplacementDim > &material_state_data)
Definition RichardsMechanicsFEM-impl.h:783
ProcessLib::RichardsMechanics::createConstitutiveModels
ConstitutiveModels< DisplacementDim > createConstitutiveModels(TRMProcessData const &process_data, MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material)
Definition ConstitutiveModels.h:24
EigenFixedShapeMatrixPolicy::NodalMatrixType
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
Definition ShapeMatrixPolicy.h:74
EigenFixedShapeMatrixPolicy::GlobalDimVectorType
VectorType< GlobalDim > GlobalDimVectorType
Definition ShapeMatrixPolicy.h:85

References assembleWithJacobianEvalConstitutiveSetting(), MaterialPropertyLib::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), ProcessLib::RichardsMechanics::createConstitutiveModels(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::current_states_, MaterialPropertyLib::density, displacement_index, displacement_size, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::element_, MaterialPropertyLib::VariableArray::gas_phase_pressure, ParameterLib::SpatialPosition::getCoordinates(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::integration_method_, NumLib::interpolateCoordinates(), ip_data_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::is_axially_symmetric_, MathLib::KelvinVector::kelvin_vector_dimensions(), MaterialPropertyLib::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::liquid_saturation, localDOF(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::material_states_, N_u_op, pressure_index, pressure_size, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::prev_states_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::process_data_, MaterialPropertyLib::reference_temperature, MaterialPropertyLib::relative_permeability, MaterialPropertyLib::saturation_micro, NumLib::detail::shapeFunctionInterpolate(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::solid_material_, MaterialPropertyLib::swelling_stress_rate, MaterialPropertyLib::VariableArray::temperature, and MathLib::KelvinVector::tensorToKelvin().

◆ assembleWithJacobianEvalConstitutiveSetting()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianEvalConstitutiveSetting ( double const t,
double const dt,
ParameterLib::SpatialPosition const & x_position,
IpData & ip_data,
MPL::VariableArray & variables,
MPL::VariableArray & variables_prev,
MPL::Medium const *const medium,
TemperatureData const T_data,
CapillaryPressureData< DisplacementDim > const & p_cap_data,
ConstitutiveData< DisplacementDim > & CD,
StatefulData< DisplacementDim > & SD,
StatefulDataPrev< DisplacementDim > const & SD_prev,
std::optional< MicroPorosityParameters > const & micro_porosity_parameters,
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material,
ProcessLib::ThermoRichardsMechanics::MaterialStateData< DisplacementDim > & material_state_data )
staticprivate

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

800{
801 auto const& liquid_phase = medium->phase("AqueousLiquid");
802 auto const& solid_phase = medium->phase("Solid");
803
804 auto const& identity2 = MathLib::KelvinVector::Invariants<
805 MathLib::KelvinVector::kelvin_vector_dimensions(
806 DisplacementDim)>::identity2;
807
808 double const temperature = T_data();
809 double const p_cap_ip = p_cap_data.p_cap;
810 double const p_cap_prev_ip = p_cap_data.p_cap_prev;
811
812 auto const& eps = std::get<StrainData<DisplacementDim>>(SD);
813 auto& S_L =
814 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(SD).S_L;
815 auto const S_L_prev =
816 std::get<
817 PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
818 SD_prev)
819 ->S_L;
820 auto const alpha =
821 medium->property(MPL::PropertyType::biot_coefficient)
822 .template value<double>(variables, x_position, t, dt);
823 *std::get<ProcessLib::ThermoRichardsMechanics::BiotData>(CD) = alpha;
824
825 variables.stress =
826 std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
827 DisplacementDim>>(SD)
828 .sigma_eff;
829 // Set mechanical strain temporary to compute tangent stiffness.
830 variables.mechanical_strain
831 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
832 eps.eps);
833 auto const C_el = ip_data.computeElasticTangentStiffness(
834 variables, t, x_position, dt, solid_material,
835 *material_state_data.material_state_variables);
836
837 auto const beta_SR =
838 (1 - alpha) / solid_material.getBulkModulus(t, x_position, &C_el);
839 variables.grain_compressibility = beta_SR;
840 std::get<ProcessLib::ThermoRichardsMechanics::SolidCompressibilityData>(CD)
841 .beta_SR = beta_SR;
842
843 auto const rho_LR =
844 liquid_phase.property(MPL::PropertyType::density)
845 .template value<double>(variables, x_position, t, dt);
846 variables.density = rho_LR;
847 *std::get<LiquidDensity>(CD) = rho_LR;
848
849 S_L = medium->property(MPL::PropertyType::saturation)
850 .template value<double>(variables, x_position, t, dt);
851 variables.liquid_saturation = S_L;
852 variables_prev.liquid_saturation = S_L_prev;
853
854 // tangent derivative for Jacobian
855 double const dS_L_dp_cap =
856 medium->property(MPL::PropertyType::saturation)
857 .template dValue<double>(variables,
858 MPL::Variable::capillary_pressure,
859 x_position, t, dt);
860 std::get<ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv>(CD)
861 .dS_L_dp_cap = dS_L_dp_cap;
862 // secant derivative from time discretization for storage
863 // use tangent, if secant is not available
864 double const DeltaS_L_Deltap_cap =
865 (p_cap_ip == p_cap_prev_ip)
866 ? dS_L_dp_cap
867 : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
868 std::get<SaturationSecantDerivative>(CD).DeltaS_L_Deltap_cap =
869 DeltaS_L_Deltap_cap;
870
871 auto const chi = [medium, x_position, t, dt](double const S_L)
872 {
873 MPL::VariableArray vs;
874 vs.liquid_saturation = S_L;
875 return medium->property(MPL::PropertyType::bishops_effective_stress)
876 .template value<double>(vs, x_position, t, dt);
877 };
878 double const chi_S_L = chi(S_L);
879 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD).chi_S_L =
880 chi_S_L;
881 double const chi_S_L_prev = chi(S_L_prev);
882 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::BishopsData>>(CD)
883 ->chi_S_L = chi_S_L_prev;
884
885 auto const dchi_dS_L =
886 medium->property(MPL::PropertyType::bishops_effective_stress)
887 .template dValue<double>(
888 variables, MPL::Variable::liquid_saturation, x_position, t, dt);
889 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD).dchi_dS_L =
890 dchi_dS_L;
891
892 double const p_FR = -chi_S_L * p_cap_ip;
893 variables.effective_pore_pressure = p_FR;
894 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
895
896 // Set volumetric strain rate for the general case without swelling.
897 variables.volumetric_strain = Invariants::trace(eps.eps);
898 // TODO (CL) changed that, using eps_prev for the moment, not B * u_prev
899 // variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
900 variables_prev.volumetric_strain = Invariants::trace(
901 std::get<PrevState<StrainData<DisplacementDim>>>(SD_prev)->eps);
902
903 auto& phi =
904 std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(SD).phi;
905 { // Porosity update
906 auto const phi_prev =
907 std::get<
908 PrevState<ProcessLib::ThermoRichardsMechanics::PorosityData>>(
909 SD_prev)
910 ->phi;
911 variables_prev.porosity = phi_prev;
912 phi = medium->property(MPL::PropertyType::porosity)
913 .template value<double>(variables, variables_prev, x_position,
914 t, dt);
915 variables.porosity = phi;
916 }
917 std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(CD).phi = phi;
918
919 if (alpha < phi)
920 {
921 auto const eid =
922 x_position.getElementID()
923 ? static_cast<std::ptrdiff_t>(*x_position.getElementID())
924 : static_cast<std::ptrdiff_t>(-1);
925 OGS_FATAL(
926 "RichardsMechanics: Biot-coefficient {} is smaller than porosity "
927 "{} in element {}.",
928 alpha, phi, eid);
929 }
930
931 auto const mu = liquid_phase.property(MPL::PropertyType::viscosity)
932 .template value<double>(variables, x_position, t, dt);
933 *std::get<ProcessLib::ThermoRichardsMechanics::LiquidViscosityData>(CD) =
934 mu;
935
936 {
937 // Swelling and possibly volumetric strain rate update.
938 auto& sigma_sw =
939 std::get<ProcessLib::ThermoRichardsMechanics::
940 ConstitutiveStress_StrainTemperature::
941 SwellingDataStateful<DisplacementDim>>(SD);
942 auto const& sigma_sw_prev =
943 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
944 ConstitutiveStress_StrainTemperature::
945 SwellingDataStateful<DisplacementDim>>>(
946 SD_prev);
947 auto const transport_porosity_prev = std::get<PrevState<
948 ProcessLib::ThermoRichardsMechanics::TransportPorosityData>>(
949 SD_prev);
950 auto const phi_prev = std::get<
951 PrevState<ProcessLib::ThermoRichardsMechanics::PorosityData>>(
952 SD_prev);
953 auto& transport_porosity = std::get<
954 ProcessLib::ThermoRichardsMechanics::TransportPorosityData>(SD);
955 auto& p_L_m = std::get<MicroPressure>(SD);
956 auto const p_L_m_prev = std::get<PrevState<MicroPressure>>(SD_prev);
957 auto& S_L_m = std::get<MicroSaturation>(SD);
958 auto const S_L_m_prev = std::get<PrevState<MicroSaturation>>(SD_prev);
959
960 updateSwellingStressAndVolumetricStrain<DisplacementDim>(
961 *medium, solid_phase, C_el, rho_LR, mu, micro_porosity_parameters,
962 alpha, phi, p_cap_ip, variables, variables_prev, x_position, t, dt,
963 sigma_sw, sigma_sw_prev, transport_porosity_prev, phi_prev,
964 transport_porosity, p_L_m_prev, S_L_m_prev, p_L_m, S_L_m);
965 }
966
967 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
968 {
969 if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
970 {
971 auto& transport_porosity =
972 std::get<
973 ProcessLib::ThermoRichardsMechanics::TransportPorosityData>(
974 SD)
975 .phi;
976 auto const transport_porosity_prev = std::get<PrevState<
977 ProcessLib::ThermoRichardsMechanics::TransportPorosityData>>(
978 SD_prev)
979 ->phi;
980 variables_prev.transport_porosity = transport_porosity_prev;
981
982 transport_porosity =
983 medium->property(MPL::PropertyType::transport_porosity)
984 .template value<double>(variables, variables_prev,
985 x_position, t, dt);
986 variables.transport_porosity = transport_porosity;
987 }
988 }
989 else
990 {
991 variables.transport_porosity = phi;
992 }
993
994 // Set mechanical variables for the intrinsic permeability model
995 // For stress dependent permeability.
996 {
997 // TODO mechanical constitutive relation will be evaluated afterwards
998 auto const sigma_total =
999 (std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
1000 DisplacementDim>>(SD)
1001 .sigma_eff +
1002 alpha * p_FR * identity2)
1003 .eval();
1004 // For stress dependent permeability.
1005 variables.total_stress.emplace<SymmetricTensor>(
1006 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_total));
1007 }
1008
1009 variables.equivalent_plastic_strain =
1010 material_state_data.material_state_variables
1011 ->getEquivalentPlasticStrain();
1012
1013 double const k_rel =
1014 medium->property(MPL::PropertyType::relative_permeability)
1015 .template value<double>(variables, x_position, t, dt);
1016
1017 auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
1018 medium->property(MPL::PropertyType::permeability)
1019 .value(variables, x_position, t, dt));
1020
1021 std::get<
1022 ProcessLib::ThermoRichardsMechanics::PermeabilityData<DisplacementDim>>(
1023 CD)
1024 .k_rel = k_rel;
1025 std::get<
1026 ProcessLib::ThermoRichardsMechanics::PermeabilityData<DisplacementDim>>(
1027 CD)
1028 .Ki = K_intrinsic;
1029
1030 //
1031 // displacement equation, displacement part
1032 //
1033
1034 {
1035 auto& sigma_sw =
1036 std::get<ProcessLib::ThermoRichardsMechanics::
1037 ConstitutiveStress_StrainTemperature::
1038 SwellingDataStateful<DisplacementDim>>(SD)
1039 .sigma_sw;
1040
1041 auto& eps_m =
1042 std::get<ProcessLib::ConstitutiveRelations::MechanicalStrainData<
1043 DisplacementDim>>(SD)
1044 .eps_m;
1045 eps_m.noalias() =
1046 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
1047 ? eps.eps + C_el.inverse() * sigma_sw
1048 : eps.eps;
1049 variables.mechanical_strain
1050 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
1051 eps_m);
1052 }
1053
1054 {
1055 auto& sigma_eff =
1056 std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
1057 DisplacementDim>>(SD);
1058 auto const& sigma_eff_prev =
1059 std::get<PrevState<ProcessLib::ConstitutiveRelations::
1060 EffectiveStressData<DisplacementDim>>>(
1061 SD_prev);
1062 auto const& eps_m =
1063 std::get<ProcessLib::ConstitutiveRelations::MechanicalStrainData<
1064 DisplacementDim>>(SD);
1065 auto& eps_m_prev =
1066 std::get<PrevState<ProcessLib::ConstitutiveRelations::
1067 MechanicalStrainData<DisplacementDim>>>(
1068 SD_prev);
1069
1070 auto C = ip_data.updateConstitutiveRelation(
1071 variables, t, x_position, dt, temperature, sigma_eff,
1072 sigma_eff_prev, eps_m, eps_m_prev, solid_material,
1073 material_state_data.material_state_variables);
1074
1075 *std::get<StiffnessTensor<DisplacementDim>>(CD) = std::move(C);
1076 }
1077
1078 // p_SR
1079 variables.solid_grain_pressure =
1080 p_FR - std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
1081 DisplacementDim>>(SD)
1082 .sigma_eff.dot(identity2) /
1083 (3 * (1 - phi));
1084 auto const rho_SR =
1085 solid_phase.property(MPL::PropertyType::density)
1086 .template value<double>(variables, x_position, t, dt);
1087
1088 double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
1089 *std::get<Density>(CD) = rho;
1090}
ProcessLib::RichardsMechanics::updateSwellingStressAndVolumetricStrain
void updateSwellingStressAndVolumetricStrain(MaterialPropertyLib::Medium const &medium, MaterialPropertyLib::Phase const &solid_phase, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > const &C_el, double const rho_LR, double const mu, std::optional< MicroPorosityParameters > micro_porosity_parameters, double const alpha, double const phi, double const p_cap_ip, MPL::VariableArray &variables, MPL::VariableArray &variables_prev, ParameterLib::SpatialPosition const &x_position, double const t, double const dt, ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature::SwellingDataStateful< DisplacementDim > &sigma_sw, PrevState< ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature::SwellingDataStateful< DisplacementDim > > const &sigma_sw_prev, PrevState< ProcessLib::ThermoRichardsMechanics::TransportPorosityData > const phi_M_prev, PrevState< ProcessLib::ThermoRichardsMechanics::PorosityData > const phi_prev, ProcessLib::ThermoRichardsMechanics::TransportPorosityData &phi_M, PrevState< MicroPressure > const p_L_m_prev, PrevState< MicroSaturation > const S_L_m_prev, MicroPressure &p_L_m, MicroSaturation &S_L_m)
Definition RichardsMechanicsFEM-impl.h:35

References RichardsMechanicsLocalAssembler(), MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::bishops_effective_stress, MaterialPropertyLib::capillary_pressure, MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::effective_pore_pressure, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialLib::Solids::MechanicsBase< DisplacementDim >::getBulkModulus(), ParameterLib::SpatialPosition::getElementID(), MaterialPropertyLib::VariableArray::grain_compressibility, MaterialPropertyLib::Medium::hasProperty(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::liquid_saturation, MaterialPropertyLib::VariableArray::liquid_saturation, ProcessLib::ThermoRichardsMechanics::MaterialStateData< DisplacementDim >::material_state_variables, MaterialPropertyLib::VariableArray::mechanical_strain, OGS_FATAL, ProcessLib::RichardsMechanics::CapillaryPressureData< DisplacementDim >::p_cap, ProcessLib::RichardsMechanics::CapillaryPressureData< DisplacementDim >::p_cap_prev, MaterialPropertyLib::permeability, MaterialPropertyLib::Medium::phase(), MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, MaterialPropertyLib::Medium::property(), MaterialPropertyLib::relative_permeability, MaterialPropertyLib::saturation, MaterialPropertyLib::saturation_micro, ProcessLib::ConstitutiveRelations::EffectiveStressData< DisplacementDim >::sigma_eff, MaterialPropertyLib::VariableArray::solid_grain_pressure, MaterialPropertyLib::VariableArray::stress, MaterialPropertyLib::swelling_stress_rate, MaterialPropertyLib::VariableArray::total_stress, MathLib::KelvinVector::Invariants< KelvinVectorSize >::trace(), MaterialPropertyLib::transport_porosity, MaterialPropertyLib::VariableArray::transport_porosity, ProcessLib::RichardsMechanics::updateSwellingStressAndVolumetricStrain(), MaterialPropertyLib::Property::value(), MaterialPropertyLib::viscosity, and MaterialPropertyLib::VariableArray::volumetric_strain.

Referenced by assembleWithJacobian().

◆ 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_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_b_dataRight hand side vector of an element.
local_Jac_dataElement Jacobian matrix for the Newton-Raphson method.

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

1510{
1511 OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
1512}

References OGS_FATAL.

Referenced by assembleWithJacobianForStaggeredScheme().

◆ 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_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_b_dataRight hand side vector of an element.
local_Jac_dataElement Jacobian matrix for the Newton-Raphson method.

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

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

References OGS_FATAL.

Referenced by assembleWithJacobianForStaggeredScheme().

◆ 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_b_data,
std::vector< double > & local_Jac_data )
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

1524{
1525 // For the equations with pressure
1526 if (process_id == 0)
1527 {
1528 assembleWithJacobianForPressureEquations(t, dt, local_x, local_x_prev,
1529 local_b_data, local_Jac_data);
1530 return;
1531 }
1532
1533 // For the equations with deformation
1534 assembleWithJacobianForDeformationEquations(t, dt, local_x, local_x_prev,
1535 local_b_data, local_Jac_data);
1536}
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler::assembleWithJacobianForPressureEquations
void assembleWithJacobianForPressureEquations(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
Definition RichardsMechanicsFEM-impl.h:1490
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler::assembleWithJacobianForDeformationEquations
void assembleWithJacobianForDeformationEquations(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
Definition RichardsMechanicsFEM-impl.h:1504

References assembleWithJacobianForDeformationEquations(), and assembleWithJacobianForPressureEquations().

◆ 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 1541 of file RichardsMechanicsFEM-impl.h.

1545{
1546 auto const [p_L, u] = localDOF(local_x);
1547 auto const [p_L_prev, u_prev] = localDOF(local_x_prev);
1548
1549 auto const& identity2 = MathLib::KelvinVector::Invariants<
1550 MathLib::KelvinVector::kelvin_vector_dimensions(
1551 DisplacementDim)>::identity2;
1552
1553 auto const& medium =
1554 this->process_data_.media_map.getMedium(this->element_.getID());
1555 auto const& liquid_phase = medium->phase("AqueousLiquid");
1556 auto const& solid_phase = medium->phase("Solid");
1557 MPL::VariableArray variables;
1558 MPL::VariableArray variables_prev;
1559
1560 unsigned const n_integration_points =
1561 this->integration_method_.getNumberOfPoints();
1562
1563 double saturation_avg = 0;
1564 double porosity_avg = 0;
1565
1566 using KV = MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
1567 KV sigma_avg = KV::Zero();
1568
1569 for (unsigned ip = 0; ip < n_integration_points; ip++)
1570 {
1571 auto const& N_p = ip_data_[ip].N_p;
1572 auto const& N_u = ip_data_[ip].N_u;
1573 auto const& dNdx_u = ip_data_[ip].dNdx_u;
1574
1575 ParameterLib::SpatialPosition x_position = {
1576 std::nullopt, this->element_.getID(),
1577 MathLib::Point3d(
1578 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
1579 ShapeMatricesTypeDisplacement>(
1580 this->element_, N_u))};
1581 auto const x_coord = x_position.getCoordinates().value()[0];
1582
1583 auto const B =
1584 LinearBMatrix::computeBMatrix<DisplacementDim,
1585 ShapeFunctionDisplacement::NPOINTS,
1586 typename BMatricesType::BMatrixType>(
1587 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
1588
1589 double p_cap_ip;
1590 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
1591
1592 double p_cap_prev_ip;
1593 NumLib::shapeFunctionInterpolate(-p_L_prev, N_p, p_cap_prev_ip);
1594
1595 variables.capillary_pressure = p_cap_ip;
1596 variables.liquid_phase_pressure = -p_cap_ip;
1597 // setting pG to 1 atm
1598 // TODO : rewrite equations s.t. p_L = pG-p_cap
1599 variables.gas_phase_pressure = 1.0e5;
1600
1601 auto const temperature =
1602 medium->property(MPL::PropertyType::reference_temperature)
1603 .template value<double>(variables, x_position, t, dt);
1604 variables.temperature = temperature;
1605
1606 auto& eps =
1607 std::get<StrainData<DisplacementDim>>(this->current_states_[ip])
1608 .eps;
1609 eps.noalias() = B * u;
1610 auto& S_L =
1611 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(
1612 this->current_states_[ip])
1613 .S_L;
1614 auto const S_L_prev =
1615 std::get<
1616 PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
1617 this->prev_states_[ip])
1618 ->S_L;
1619 S_L = medium->property(MPL::PropertyType::saturation)
1620 .template value<double>(variables, x_position, t, dt);
1621 variables.liquid_saturation = S_L;
1622 variables_prev.liquid_saturation = S_L_prev;
1623
1624 auto const chi = [medium, x_position, t, dt](double const S_L)
1625 {
1626 MPL::VariableArray vs;
1627 vs.liquid_saturation = S_L;
1628 return medium->property(MPL::PropertyType::bishops_effective_stress)
1629 .template value<double>(vs, x_position, t, dt);
1630 };
1631 double const chi_S_L = chi(S_L);
1632 double const chi_S_L_prev = chi(S_L_prev);
1633
1634 auto const alpha =
1635 medium->property(MPL::PropertyType::biot_coefficient)
1636 .template value<double>(variables, x_position, t, dt);
1637 auto& SD = this->current_states_[ip];
1638 variables.stress =
1639 std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
1640 DisplacementDim>>(SD)
1641 .sigma_eff;
1642 // Set mechanical strain temporary to compute tangent stiffness.
1643 variables.mechanical_strain
1644 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
1645 eps);
1646 auto const C_el = ip_data_[ip].computeElasticTangentStiffness(
1647 variables, t, x_position, dt, this->solid_material_,
1648 *this->material_states_[ip].material_state_variables);
1649
1650 auto const beta_SR = (1 - alpha) / this->solid_material_.getBulkModulus(
1651 t, x_position, &C_el);
1652 variables.grain_compressibility = beta_SR;
1653
1654 variables.effective_pore_pressure = -chi_S_L * p_cap_ip;
1655 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
1656
1657 // Set volumetric strain rate for the general case without swelling.
1658 variables.volumetric_strain = Invariants::trace(eps);
1659 variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
1660
1661 auto& phi = std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(
1662 this->current_states_[ip])
1663 .phi;
1664 { // Porosity update
1665 auto const phi_prev = std::get<PrevState<
1666 ProcessLib::ThermoRichardsMechanics::PorosityData>>(
1667 this->prev_states_[ip])
1668 ->phi;
1669 variables_prev.porosity = phi_prev;
1670 phi = medium->property(MPL::PropertyType::porosity)
1671 .template value<double>(variables, variables_prev,
1672 x_position, t, dt);
1673 variables.porosity = phi;
1674 }
1675
1676 auto const rho_LR =
1677 liquid_phase.property(MPL::PropertyType::density)
1678 .template value<double>(variables, x_position, t, dt);
1679 variables.density = rho_LR;
1680 auto const mu =
1681 liquid_phase.property(MPL::PropertyType::viscosity)
1682 .template value<double>(variables, x_position, t, dt);
1683
1684 {
1685 // Swelling and possibly volumetric strain rate update.
1686 auto& sigma_sw =
1687 std::get<ProcessLib::ThermoRichardsMechanics::
1688 ConstitutiveStress_StrainTemperature::
1689 SwellingDataStateful<DisplacementDim>>(
1690 this->current_states_[ip]);
1691 auto const& sigma_sw_prev = std::get<
1692 PrevState<ProcessLib::ThermoRichardsMechanics::
1693 ConstitutiveStress_StrainTemperature::
1694 SwellingDataStateful<DisplacementDim>>>(
1695 this->prev_states_[ip]);
1696 auto const transport_porosity_prev = std::get<PrevState<
1697 ProcessLib::ThermoRichardsMechanics::TransportPorosityData>>(
1698 this->prev_states_[ip]);
1699 auto const phi_prev = std::get<
1700 PrevState<ProcessLib::ThermoRichardsMechanics::PorosityData>>(
1701 this->prev_states_[ip]);
1702 auto& transport_porosity = std::get<
1703 ProcessLib::ThermoRichardsMechanics::TransportPorosityData>(
1704 this->current_states_[ip]);
1705 auto& p_L_m = std::get<MicroPressure>(this->current_states_[ip]);
1706 auto const p_L_m_prev =
1707 std::get<PrevState<MicroPressure>>(this->prev_states_[ip]);
1708 auto& S_L_m = std::get<MicroSaturation>(this->current_states_[ip]);
1709 auto const S_L_m_prev =
1710 std::get<PrevState<MicroSaturation>>(this->prev_states_[ip]);
1711
1712 updateSwellingStressAndVolumetricStrain<DisplacementDim>(
1713 *medium, solid_phase, C_el, rho_LR, mu,
1714 this->process_data_.micro_porosity_parameters, alpha, phi,
1715 p_cap_ip, variables, variables_prev, x_position, t, dt,
1716 sigma_sw, sigma_sw_prev, transport_porosity_prev, phi_prev,
1717 transport_porosity, p_L_m_prev, S_L_m_prev, p_L_m, S_L_m);
1718 }
1719
1720 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
1721 {
1722 if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
1723 {
1724 auto& transport_porosity =
1725 std::get<ProcessLib::ThermoRichardsMechanics::
1726 TransportPorosityData>(
1727 this->current_states_[ip])
1728 .phi;
1729 auto const transport_porosity_prev =
1730 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
1731 TransportPorosityData>>(
1732 this->prev_states_[ip])
1733 ->phi;
1734
1735 variables_prev.transport_porosity = transport_porosity_prev;
1736
1737 transport_porosity =
1738 medium->property(MPL::PropertyType::transport_porosity)
1739 .template value<double>(variables, variables_prev,
1740 x_position, t, dt);
1741 variables.transport_porosity = transport_porosity;
1742 }
1743 }
1744 else
1745 {
1746 variables.transport_porosity = phi;
1747 }
1748
1749 auto const& sigma_eff =
1750 std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
1751 DisplacementDim>>(this->current_states_[ip])
1752 .sigma_eff;
1753
1754 // Set mechanical variables for the intrinsic permeability model
1755 // For stress dependent permeability.
1756 {
1757 auto const sigma_total =
1758 (sigma_eff + alpha * chi_S_L * identity2 * p_cap_ip).eval();
1759 // For stress dependent permeability.
1760 variables.total_stress.emplace<SymmetricTensor>(
1761 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
1762 sigma_total));
1763 }
1764
1765 variables.equivalent_plastic_strain =
1766 this->material_states_[ip]
1767 .material_state_variables->getEquivalentPlasticStrain();
1768
1769 auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
1770 medium->property(MPL::PropertyType::permeability)
1771 .value(variables, x_position, t, dt));
1772
1773 double const k_rel =
1774 medium->property(MPL::PropertyType::relative_permeability)
1775 .template value<double>(variables, x_position, t, dt);
1776
1777 GlobalDimMatrixType const K_over_mu = k_rel * K_intrinsic / mu;
1778
1779 double const p_FR = -chi_S_L * p_cap_ip;
1780 // p_SR
1781 variables.solid_grain_pressure =
1782 p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
1783 auto const rho_SR =
1784 solid_phase.property(MPL::PropertyType::density)
1785 .template value<double>(variables, x_position, t, dt);
1786 *std::get<DrySolidDensity>(this->output_data_[ip]) = (1 - phi) * rho_SR;
1787
1788 {
1789 auto& SD = this->current_states_[ip];
1790 auto const& sigma_sw =
1791 std::get<ProcessLib::ThermoRichardsMechanics::
1792 ConstitutiveStress_StrainTemperature::
1793 SwellingDataStateful<DisplacementDim>>(SD)
1794 .sigma_sw;
1795 auto& eps_m =
1796 std::get<ProcessLib::ConstitutiveRelations::
1797 MechanicalStrainData<DisplacementDim>>(SD)
1798 .eps_m;
1799 eps_m.noalias() =
1800 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
1801 ? eps + C_el.inverse() * sigma_sw
1802 : eps;
1803 variables.mechanical_strain.emplace<
1804 MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
1805 eps_m);
1806 }
1807
1808 {
1809 auto& SD = this->current_states_[ip];
1810 auto const& SD_prev = this->prev_states_[ip];
1811 auto& sigma_eff =
1812 std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
1813 DisplacementDim>>(SD);
1814 auto const& sigma_eff_prev =
1815 std::get<PrevState<ProcessLib::ConstitutiveRelations::
1816 EffectiveStressData<DisplacementDim>>>(
1817 SD_prev);
1818 auto const& eps_m =
1819 std::get<ProcessLib::ConstitutiveRelations::
1820 MechanicalStrainData<DisplacementDim>>(SD);
1821 auto const& eps_m_prev =
1822 std::get<PrevState<ProcessLib::ConstitutiveRelations::
1823 MechanicalStrainData<DisplacementDim>>>(
1824 SD_prev);
1825
1826 ip_data_[ip].updateConstitutiveRelation(
1827 variables, t, x_position, dt, temperature, sigma_eff,
1828 sigma_eff_prev, eps_m, eps_m_prev, this->solid_material_,
1829 this->material_states_[ip].material_state_variables);
1830 }
1831
1832 auto const& b = this->process_data_.specific_body_force;
1833
1834 // Compute the velocity
1835 auto const& dNdx_p = ip_data_[ip].dNdx_p;
1836 std::get<
1837 ProcessLib::ThermoRichardsMechanics::DarcyLawData<DisplacementDim>>(
1838 this->output_data_[ip])
1839 ->noalias() = -K_over_mu * dNdx_p * p_L + rho_LR * K_over_mu * b;
1840
1841 saturation_avg += S_L;
1842 porosity_avg += phi;
1843 sigma_avg += sigma_eff;
1844 }
1845 saturation_avg /= n_integration_points;
1846 porosity_avg /= n_integration_points;
1847 sigma_avg /= n_integration_points;
1848
1849 (*this->process_data_.element_saturation)[this->element_.getID()] =
1850 saturation_avg;
1851 (*this->process_data_.element_porosity)[this->element_.getID()] =
1852 porosity_avg;
1853
1854 Eigen::Map<KV>(
1855 &(*this->process_data_.element_stresses)[this->element_.getID() *
1856 KV::RowsAtCompileTime]) =
1857 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_avg);
1858
1859 NumLib::interpolateToHigherOrderNodes<
1860 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1861 DisplacementDim>(this->element_, this->is_axially_symmetric_, p_L,
1862 *this->process_data_.pressure_interpolated);
1863}
MeshLib::Element::getID
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:91
ProcessLib::RichardsMechanics::LocalAssemblerInterface::output_data_
std::vector< OutputData< DisplacementDim > > output_data_
Definition LocalAssemblerInterface.h:190

References MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::bishops_effective_stress, MaterialPropertyLib::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::current_states_, MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::effective_pore_pressure, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::element_, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::gas_phase_pressure, ParameterLib::SpatialPosition::getCoordinates(), MaterialPropertyLib::VariableArray::grain_compressibility, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::integration_method_, NumLib::interpolateCoordinates(), NumLib::interpolateToHigherOrderNodes(), ip_data_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::is_axially_symmetric_, MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, localDOF(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::material_states_, MaterialPropertyLib::VariableArray::mechanical_strain, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::output_data_, MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::prev_states_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::process_data_, MaterialPropertyLib::reference_temperature, MaterialPropertyLib::relative_permeability, MaterialPropertyLib::saturation, MaterialPropertyLib::saturation_micro, NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::VariableArray::solid_grain_pressure, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::solid_material_, MaterialPropertyLib::VariableArray::stress, MaterialPropertyLib::swelling_stress_rate, MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, MathLib::KelvinVector::Invariants< KelvinVectorSize >::trace(), MaterialPropertyLib::transport_porosity, MaterialPropertyLib::VariableArray::transport_porosity, ProcessLib::RichardsMechanics::updateSwellingStressAndVolumetricStrain(), MaterialPropertyLib::viscosity, and MaterialPropertyLib::VariableArray::volumetric_strain.

◆ 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 165 of file RichardsMechanicsFEM.h.

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

References secondary_data_.

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

118 {
119 unsigned const n_integration_points =
120 this->integration_method_.getNumberOfPoints();
121
122 for (unsigned ip = 0; ip < n_integration_points; ip++)
123 {
124 auto& SD = this->current_states_[ip];
125 auto& ip_data = ip_data_[ip];
126
127 ParameterLib::SpatialPosition const x_position{
128 std::nullopt, this->element_.getID(),
129 MathLib::Point3d(NumLib::interpolateCoordinates<
130 ShapeFunctionDisplacement,
131 ShapeMatricesTypeDisplacement>(this->element_,
132 ip_data.N_u))};
133
135 if (this->process_data_.initial_stress.value)
136 {
137 std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
138 DisplacementDim>>(SD)
139 .sigma_eff =
140 MathLib::KelvinVector::symmetricTensorToKelvinVector<
141 DisplacementDim>(
142 // The data in process_data_.initial_stress.value can
143 // be total stress or effective stress.
144 (*this->process_data_.initial_stress.value)(
145 std::numeric_limits<
146 double>::quiet_NaN() /* time independent */,
147 x_position));
148 }
149
150 double const t = 0; // TODO (naumov) pass t from top
151 this->solid_material_.initializeInternalStateVariables(
152 t, x_position,
153 *this->material_states_[ip].material_state_variables);
154
155 this->material_states_[ip].pushBackState();
156
157 this->prev_states_[ip] = SD;
158 }
159 }

References ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::current_states_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::element_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::integration_method_, NumLib::interpolateCoordinates(), ip_data_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::material_states_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::prev_states_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::process_data_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::solid_material_, and MathLib::KelvinVector::symmetricTensorToKelvinVector().

◆ localDOF()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
constexpr auto ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::localDOF ( auto const & x)
inlinestaticconstexprprivate

Definition at line 240 of file RichardsMechanicsFEM.h.

241 {
242 return NumLib::localDOF<
243 ShapeFunctionPressure,
244 NumLib::Vectorial<ShapeFunctionDisplacement, DisplacementDim>>(x);
245 }

References NumLib::localDOF().

Referenced by assemble(), assembleWithJacobian(), computeSecondaryVariableConcrete(), and setInitialConditionsConcrete().

◆ setInitialConditionsConcrete()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::setInitialConditionsConcrete ( Eigen::VectorXd const local_x,
double const t,
int const process_id )
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

218{
219 assert(local_x.size() == pressure_size + displacement_size);
220
221 auto const [p_L, u] = localDOF(local_x);
222
223 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
224 auto const& medium =
225 this->process_data_.media_map.getMedium(this->element_.getID());
226 MPL::VariableArray variables;
227
228 auto const& solid_phase = medium->phase("Solid");
229
230 auto const& identity2 = MathLib::KelvinVector::Invariants<
231 MathLib::KelvinVector::kelvin_vector_dimensions(
232 DisplacementDim)>::identity2;
233
234 unsigned const n_integration_points =
235 this->integration_method_.getNumberOfPoints();
236 for (unsigned ip = 0; ip < n_integration_points; ip++)
237 {
238 auto const& N_p = ip_data_[ip].N_p;
239
240 ParameterLib::SpatialPosition x_position = {
241 std::nullopt, this->element_.getID(),
242 MathLib::Point3d(
243 NumLib::interpolateCoordinates<ShapeFunctionPressure,
244 ShapeMatricesTypePressure>(
245 this->element_, N_p))};
246
247 double p_cap_ip;
248 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
249
250 variables.capillary_pressure = p_cap_ip;
251 variables.liquid_phase_pressure = -p_cap_ip;
252 // setting pG to 1 atm
253 // TODO : rewrite equations s.t. p_L = pG-p_cap
254 variables.gas_phase_pressure = 1.0e5;
255
256 {
257 auto& p_L_m = std::get<MicroPressure>(this->current_states_[ip]);
258 auto& p_L_m_prev =
259 std::get<PrevState<MicroPressure>>(this->prev_states_[ip]);
260 **p_L_m_prev = -p_cap_ip;
261 *p_L_m = -p_cap_ip;
262 }
263
264 auto const temperature =
265 medium->property(MPL::PropertyType::reference_temperature)
266 .template value<double>(variables, x_position, t, dt);
267 variables.temperature = temperature;
268
269 auto& S_L_prev =
270 std::get<
271 PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
272 this->prev_states_[ip])
273 ->S_L;
274 S_L_prev = medium->property(MPL::PropertyType::saturation)
275 .template value<double>(variables, x_position, t, dt);
276
277 if (this->process_data_.initial_stress.isTotalStress())
278 {
279 auto const alpha_b =
280 medium->property(MPL::PropertyType::biot_coefficient)
281 .template value<double>(variables, x_position, t, dt);
282
283 variables.liquid_saturation = S_L_prev;
284 double const chi_S_L =
285 medium->property(MPL::PropertyType::bishops_effective_stress)
286 .template value<double>(variables, x_position, t, dt);
287
288 // Initial stresses are total stress, which were assigned to
289 // sigma_eff in
290 // RichardsMechanicsLocalAssembler::initializeConcrete().
291 auto& sigma_eff =
292 std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
293 DisplacementDim>>(this->current_states_[ip]);
294
295 auto& sigma_eff_prev =
296 std::get<PrevState<ProcessLib::ConstitutiveRelations::
297 EffectiveStressData<DisplacementDim>>>(
298 this->prev_states_[ip]);
299
300 // Reset sigma_eff to effective stress
301 sigma_eff.sigma_eff.noalias() +=
302 chi_S_L * alpha_b * (-p_cap_ip) * identity2;
303 sigma_eff_prev->sigma_eff = sigma_eff.sigma_eff;
304 }
305
306 if (medium->hasProperty(MPL::PropertyType::saturation_micro))
307 {
308 MPL::VariableArray vars;
309 vars.capillary_pressure = p_cap_ip;
310
311 auto& S_L_m = std::get<MicroSaturation>(this->current_states_[ip]);
312 auto& S_L_m_prev =
313 std::get<PrevState<MicroSaturation>>(this->prev_states_[ip]);
314
315 *S_L_m = medium->property(MPL::PropertyType::saturation_micro)
316 .template value<double>(vars, x_position, t, dt);
317 *S_L_m_prev = S_L_m;
318 }
319
320 // Set eps_m_prev from potentially non-zero eps and sigma_sw from
321 // restart.
322 auto& SD = this->current_states_[ip];
323 variables.stress =
324 std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
325 DisplacementDim>>(SD)
326 .sigma_eff;
327
328 auto const& N_u = ip_data_[ip].N_u;
329 auto const& dNdx_u = ip_data_[ip].dNdx_u;
330 auto const x_coord =
331 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
332 ShapeMatricesTypeDisplacement>(
333 this->element_, N_u);
334 auto const B =
335 LinearBMatrix::computeBMatrix<DisplacementDim,
336 ShapeFunctionDisplacement::NPOINTS,
337 typename BMatricesType::BMatrixType>(
338 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
339 auto& eps =
340 std::get<StrainData<DisplacementDim>>(this->current_states_[ip])
341 .eps;
342 eps.noalias() = B * u;
343
344 // Set mechanical strain temporary to compute tangent stiffness.
345 variables.mechanical_strain
346 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
347 eps);
348
349 auto const C_el = ip_data_[ip].computeElasticTangentStiffness(
350 variables, t, x_position, dt, this->solid_material_,
351 *this->material_states_[ip].material_state_variables);
352
353 auto const& sigma_sw =
354 std::get<ProcessLib::ThermoRichardsMechanics::
355 ConstitutiveStress_StrainTemperature::
356 SwellingDataStateful<DisplacementDim>>(
357 this->current_states_[ip])
358 .sigma_sw;
359 auto& eps_m_prev =
360 std::get<PrevState<ProcessLib::ConstitutiveRelations::
361 MechanicalStrainData<DisplacementDim>>>(
362 this->prev_states_[ip])
363 ->eps_m;
364
365 eps_m_prev.noalias() =
366 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
367 ? eps + C_el.inverse() * sigma_sw
368 : eps;
369 }
370}

References MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::bishops_effective_stress, MaterialPropertyLib::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::current_states_, displacement_size, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::element_, MaterialPropertyLib::VariableArray::gas_phase_pressure, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::integration_method_, NumLib::interpolateCoordinates(), NumLib::interpolateXCoordinate(), ip_data_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::is_axially_symmetric_, MathLib::KelvinVector::kelvin_vector_dimensions(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, localDOF(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::material_states_, MaterialPropertyLib::VariableArray::mechanical_strain, pressure_size, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::prev_states_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::process_data_, MaterialPropertyLib::reference_temperature, MaterialPropertyLib::saturation, MaterialPropertyLib::saturation_micro, NumLib::detail::shapeFunctionInterpolate(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::solid_material_, MaterialPropertyLib::VariableArray::stress, MaterialPropertyLib::swelling_stress_rate, and MaterialPropertyLib::VariableArray::temperature.

Member Data Documentation

◆ 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 255 of file RichardsMechanicsFEM.h.

Referenced by assemble(), and assembleWithJacobian().

◆ 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 256 of file RichardsMechanicsFEM.h.

Referenced by assemble(), assembleWithJacobian(), and setInitialConditionsConcrete().

◆ ip_data_

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
std::vector<IpData, Eigen::aligned_allocator<IpData> > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ip_data_
private

Definition at line 247 of file RichardsMechanicsFEM.h.

Referenced by assemble(), assembleWithJacobian(), computeSecondaryVariableConcrete(), initializeConcrete(), and setInitialConditionsConcrete().

◆ KelvinVectorSize

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
int const ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::KelvinVectorSize
static
Initial value:
=
MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim)

Definition at line 73 of file RichardsMechanicsFEM.h.

◆ N_u_op

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
auto& ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::N_u_op
staticconstexpr
Initial value:
= MathLib::eigenBlockMatrixView<
DisplacementDim,
typename ShapeMatricesTypeDisplacement::NodalRowVectorType>
MathLib::eigenBlockMatrixView
constexpr Eigen::CwiseNullaryOp< EigenBlockMatrixViewFunctor< D, M >, typename EigenBlockMatrixViewFunctor< D, M >::Matrix > eigenBlockMatrixView(const Eigen::MatrixBase< M > &matrix)
Definition EigenBlockMatrixView.h:49
EigenFixedShapeMatrixPolicy::NodalRowVectorType
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
Definition ShapeMatrixPolicy.h:78

Definition at line 79 of file RichardsMechanicsFEM.h.

Referenced by assemble(), and assembleWithJacobian().

◆ pressure_index

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

Definition at line 253 of file RichardsMechanicsFEM.h.

Referenced by assemble(), and assembleWithJacobian().

◆ 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 254 of file RichardsMechanicsFEM.h.

Referenced by assemble(), assembleWithJacobian(), and setInitialConditionsConcrete().

◆ secondary_data_

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

Definition at line 251 of file RichardsMechanicsFEM.h.

Referenced by getShapeMatrix().

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