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 45 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.
int getNumberOfVectorElementsForDeformation () const override
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.
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 57 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:77

Definition at line 54 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 68 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:163
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::ShapeMatricesTypePressure
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
Definition HydroMechanicsFEM.h:167
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::BMatricesType
BMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > BMatricesType
Definition HydroMechanicsFEM.h:384

Definition at line 61 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 59 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:134

Definition at line 49 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 51 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 70 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 121 of file RichardsMechanicsFEM-impl.h.

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

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

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

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

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

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

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

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

1507{
1508 OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
1509}

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

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

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

1521{
1522 // For the equations with pressure
1523 if (process_id == 0)
1524 {
1525 assembleWithJacobianForPressureEquations(t, dt, local_x, local_x_prev,
1526 local_b_data, local_Jac_data);
1527 return;
1528 }
1529
1530 // For the equations with deformation
1531 assembleWithJacobianForDeformationEquations(t, dt, local_x, local_x_prev,
1532 local_b_data, local_Jac_data);
1533}
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:1487
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:1501

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

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

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.

◆ getNumberOfVectorElementsForDeformation()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 167 of file RichardsMechanicsFEM.h.

168 {
169 return displacement_size;
170 }

References displacement_size.

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

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

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

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

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

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

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

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

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, 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::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 253 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 254 of file RichardsMechanicsFEM.h.

Referenced by assemble(), assembleWithJacobian(), getNumberOfVectorElementsForDeformation(), 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 245 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 66 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:43
EigenFixedShapeMatrixPolicy::NodalRowVectorType
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
Definition ShapeMatrixPolicy.h:71

Definition at line 72 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 251 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 252 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 249 of file RichardsMechanicsFEM.h.

Referenced by getShapeMatrix().

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