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

Detailed Description

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

Definition at line 47 of file RichardsMechanicsFEM.h.

#include <RichardsMechanicsFEM.h>

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

Public Types

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

Public Member Functions

 RichardsMechanicsLocalAssembler (RichardsMechanicsLocalAssembler const &)=delete
 
 RichardsMechanicsLocalAssembler (RichardsMechanicsLocalAssembler &&)=delete
 
 RichardsMechanicsLocalAssembler (MeshLib::Element const &e, std::size_t const, bool const is_axially_symmetric, unsigned const integration_order, RichardsMechanicsProcessData< DisplacementDim > &process_data)
 
std::size_t setIPDataInitialConditions (std::string const &name, double const *values, int const integration_order) override
 
void setInitialConditionsConcrete (std::vector< double > const &local_x, double const t, bool const use_monolithic_scheme, int const process_id) override
 
void assemble (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_xdot, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_rhs_data) override
 
void assembleWithJacobian (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_xdot, const double, const double, std::vector< double > &, std::vector< double > &, std::vector< double > &local_rhs_data, std::vector< double > &local_Jac_data) override
 
void assembleWithJacobianForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, const double dxdot_dx, const double dx_dx, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
 
void initializeConcrete () override
 
void postTimestepConcrete (Eigen::VectorXd const &, double const, double const) override
 
void computeSecondaryVariableConcrete (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_dot) override
 
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix (const unsigned integration_point) const override
 Provides the shape matrix at the given integration point. More...
 
std::vector< double > getSigma () const override
 
std::vector< double > const & getIntPtDarcyVelocity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::vector< double > getSaturation () const override
 
std::vector< double > const & getIntPtSaturation (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::vector< double > getMicroSaturation () const override
 
std::vector< double > const & getIntPtMicroSaturation (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::vector< double > getMicroPressure () const override
 
std::vector< double > const & getIntPtMicroPressure (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::vector< double > getPorosity () const override
 
std::vector< double > const & getIntPtPorosity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::vector< double > getTransportPorosity () const override
 
std::vector< double > const & getIntPtTransportPorosity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::vector< double > const & getIntPtSigma (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::vector< double > getSwellingStress () const override
 
std::vector< double > const & getIntPtSwellingStress (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::vector< double > getEpsilon () const override
 
std::vector< double > const & getIntPtEpsilon (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::vector< double > getMaterialStateVariableInternalState (std::function< BaseLib::DynamicSpan< double >(typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables &)> const &get_values_span, int const &n_components) const override
 
std::vector< double > const & getIntPtDryDensitySolid (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
- Public Member Functions inherited from ProcessLib::LocalAssemblerInterface
virtual ~LocalAssemblerInterface ()=default
 
void setInitialConditions (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, bool const use_monolithic_scheme, int const process_id)
 
virtual void initialize (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
 
virtual void preAssemble (double const, double const, std::vector< double > const &)
 
virtual void assembleForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
 
virtual void computeSecondaryVariable (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_dot, int const process_id)
 
virtual void preTimestep (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
 
virtual void postTimestep (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double const dt)
 
void postNonLinearSolver (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, GlobalVector const &xdot, double const t, double const dt, bool const use_monolithic_scheme, int const process_id)
 
virtual std::vector< double > interpolateNodalValuesToIntegrationPoints (std::vector< double > const &)
 
virtual Eigen::Vector3d getFlux (MathLib::Point3d const &, double const, std::vector< double > const &) const
 
virtual Eigen::Vector3d getFlux (MathLib::Point3d const &, double const, std::vector< std::vector< double >> const &) const
 Fits to staggered scheme. More...
 
- Public Member Functions inherited from NumLib::ExtrapolatableElement
virtual ~ExtrapolatableElement ()=default
 

Static Public Attributes

static int const KelvinVectorSize
 

Private Member Functions

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

Private Attributes

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

Static Private Attributes

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

Member Typedef Documentation

◆ BMatricesType

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

Definition at line 59 of file RichardsMechanicsFEM.h.

◆ GlobalDimMatrixType

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

Definition at line 56 of file RichardsMechanicsFEM.h.

◆ Invariants

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

Definition at line 70 of file RichardsMechanicsFEM.h.

◆ IpData

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

Definition at line 63 of file RichardsMechanicsFEM.h.

◆ KelvinVectorType

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

Definition at line 61 of file RichardsMechanicsFEM.h.

◆ ShapeMatricesTypeDisplacement

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

Definition at line 51 of file RichardsMechanicsFEM.h.

◆ ShapeMatricesTypePressure

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

Definition at line 53 of file RichardsMechanicsFEM.h.

◆ SymmetricTensor

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

Definition at line 72 of file RichardsMechanicsFEM.h.

Constructor & Destructor Documentation

◆ RichardsMechanicsLocalAssembler() [1/3]

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

◆ RichardsMechanicsLocalAssembler() [2/3]

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

◆ RichardsMechanicsLocalAssembler() [3/3]

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

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

130  : _process_data(process_data),
131  _integration_method(integration_order),
132  _element(e),
133  _is_axially_symmetric(is_axially_symmetric)
134 {
135  unsigned const n_integration_points =
136  _integration_method.getNumberOfPoints();
137 
138  _ip_data.reserve(n_integration_points);
139  _secondary_data.N_u.resize(n_integration_points);
140 
141  auto const shape_matrices_u =
142  NumLib::initShapeMatrices<ShapeFunctionDisplacement,
144  DisplacementDim>(e, is_axially_symmetric,
146 
147  auto const shape_matrices_p =
148  NumLib::initShapeMatrices<ShapeFunctionPressure,
149  ShapeMatricesTypePressure, DisplacementDim>(
150  e, is_axially_symmetric, _integration_method);
151 
152  auto const& solid_material =
154  _process_data.solid_materials, _process_data.material_ids,
155  e.getID());
156 
157  auto const& medium = _process_data.media_map->getMedium(_element.getID());
158 
160  x_position.setElementID(_element.getID());
161  for (unsigned ip = 0; ip < n_integration_points; ip++)
162  {
163  x_position.setIntegrationPoint(ip);
164  _ip_data.emplace_back(solid_material);
165  auto& ip_data = _ip_data[ip];
166  auto const& sm_u = shape_matrices_u[ip];
167  _ip_data[ip].integration_weight =
168  _integration_method.getWeightedPoint(ip).getWeight() *
169  sm_u.integralMeasure * sm_u.detJ;
170 
171  ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
172  DisplacementDim, displacement_size>::Zero(DisplacementDim,
174  for (int i = 0; i < DisplacementDim; ++i)
175  {
176  ip_data.N_u_op
177  .template block<1, displacement_size / DisplacementDim>(
178  i, i * displacement_size / DisplacementDim)
179  .noalias() = sm_u.N;
180  }
181 
182  ip_data.N_u = sm_u.N;
183  ip_data.dNdx_u = sm_u.dNdx;
184 
185  ip_data.N_p = shape_matrices_p[ip].N;
186  ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
187 
188  // Initial porosity. Could be read from integration point data or mesh.
189  ip_data.porosity =
190  medium->property(MPL::porosity)
191  .template initialValue<double>(
192  x_position,
193  std::numeric_limits<
194  double>::quiet_NaN() /* t independent */);
195 
196  ip_data.transport_porosity = ip_data.porosity;
197  if (medium->hasProperty(MPL::PropertyType::transport_porosity))
198  {
199  ip_data.transport_porosity =
200  medium->property(MPL::transport_porosity)
201  .template initialValue<double>(
202  x_position,
203  std::numeric_limits<
204  double>::quiet_NaN() /* t independent */);
205  }
206 
207  _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
208  }
209 }
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:82
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
RichardsMechanicsProcessData< DisplacementDim > & _process_data
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
MechanicsBase< DisplacementDim > & selectSolidConstitutiveRelation(std::map< int, std::unique_ptr< MechanicsBase< DisplacementDim >>> const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u

References ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::_element, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::_integration_method, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::_ip_data, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::_process_data, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::_secondary_data, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::displacement_size, MeshLib::Element::getID(), NumLib::initShapeMatrices(), ProcessLib::HydroMechanics::SecondaryData< ShapeMatrixType >::N_u, MaterialPropertyLib::porosity, MaterialLib::Solids::selectSolidConstitutiveRelation(), ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), and MaterialPropertyLib::transport_porosity.

Member Function Documentation

◆ assemble()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

382 {
383  assert(local_x.size() == pressure_size + displacement_size);
384 
385  auto p_L =
386  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
387  pressure_size> const>(local_x.data() + pressure_index,
388  pressure_size);
389 
390  auto u =
391  Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
392  displacement_size> const>(local_x.data() + displacement_index,
394 
395  auto p_L_dot =
396  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
397  pressure_size> const>(local_xdot.data() + pressure_index,
398  pressure_size);
399 
400  auto u_dot =
401  Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
402  displacement_size> const>(local_xdot.data() + displacement_index,
404 
406  typename ShapeMatricesTypeDisplacement::template MatrixType<
409  local_K_data, displacement_size + pressure_size,
411 
413  typename ShapeMatricesTypeDisplacement::template MatrixType<
416  local_M_data, displacement_size + pressure_size,
418 
419  auto rhs = MathLib::createZeroedVector<
420  typename ShapeMatricesTypeDisplacement::template VectorType<
422  local_rhs_data, displacement_size + pressure_size);
423 
424  auto const& identity2 = MathLib::KelvinVector::Invariants<
426  DisplacementDim)>::identity2;
427 
428  auto const& medium = _process_data.media_map->getMedium(_element.getID());
429  auto const& liquid_phase = medium->phase("AqueousLiquid");
430  auto const& solid_phase = medium->phase("Solid");
431  MPL::VariableArray variables;
432  MPL::VariableArray variables_prev;
433 
435  x_position.setElementID(_element.getID());
436 
437  unsigned const n_integration_points =
438  _integration_method.getNumberOfPoints();
439  for (unsigned ip = 0; ip < n_integration_points; ip++)
440  {
441  x_position.setIntegrationPoint(ip);
442  auto const& w = _ip_data[ip].integration_weight;
443 
444  auto const& N_u_op = _ip_data[ip].N_u_op;
445 
446  auto const& N_u = _ip_data[ip].N_u;
447  auto const& dNdx_u = _ip_data[ip].dNdx_u;
448 
449  auto const& N_p = _ip_data[ip].N_p;
450  auto const& dNdx_p = _ip_data[ip].dNdx_p;
451 
452  auto const x_coord =
453  NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
455  _element, N_u);
456  auto const B =
457  LinearBMatrix::computeBMatrix<DisplacementDim,
458  ShapeFunctionDisplacement::NPOINTS,
459  typename BMatricesType::BMatrixType>(
460  dNdx_u, N_u, x_coord, _is_axially_symmetric);
461 
462  auto& eps = _ip_data[ip].eps;
463  auto& eps_m = _ip_data[ip].eps_m;
464  eps.noalias() = B * u;
465 
466  auto& S_L = _ip_data[ip].saturation;
467  auto const S_L_prev = _ip_data[ip].saturation_prev;
468 
469  double p_cap_ip;
470  NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
471 
472  double p_cap_dot_ip;
473  NumLib::shapeFunctionInterpolate(-p_L_dot, N_p, p_cap_dot_ip);
474 
475  variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
476  p_cap_ip;
477  variables[static_cast<int>(MPL::Variable::phase_pressure)] = -p_cap_ip;
478 
479  auto const temperature =
481  .template value<double>(variables, x_position, t, dt);
482  variables[static_cast<int>(MPL::Variable::temperature)] = temperature;
483 
484  auto const alpha =
485  medium->property(MPL::PropertyType::biot_coefficient)
486  .template value<double>(variables, x_position, t, dt);
487  auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
488  t, x_position, dt, temperature);
489 
490  auto const beta_SR =
491  (1 - alpha) /
492  _ip_data[ip].solid_material.getBulkModulus(t, x_position, &C_el);
493  variables[static_cast<int>(MPL::Variable::grain_compressibility)] =
494  beta_SR;
495 
496  auto const rho_LR =
497  liquid_phase.property(MPL::PropertyType::density)
498  .template value<double>(variables, x_position, t, dt);
499 
500  auto const& b = _process_data.specific_body_force;
501 
502  S_L = medium->property(MPL::PropertyType::saturation)
503  .template value<double>(variables, x_position, t, dt);
504  variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
505  variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
506  S_L_prev;
507 
508  // tangent derivative for Jacobian
509  double const dS_L_dp_cap =
510  medium->property(MPL::PropertyType::saturation)
511  .template dValue<double>(variables,
513  x_position, t, dt);
514  // secant derivative from time discretization for storage
515  // use tangent, if secant is not available
516  double const DeltaS_L_Deltap_cap =
517  (p_cap_dot_ip == 0) ? dS_L_dp_cap
518  : (S_L - S_L_prev) / (dt * p_cap_dot_ip);
519 
520  auto const chi = [medium, x_position, t, dt](double const S_L)
521  {
523  vs[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
524  return medium->property(MPL::PropertyType::bishops_effective_stress)
525  .template value<double>(vs, x_position, t, dt);
526  };
527  double const chi_S_L = chi(S_L);
528  double const chi_S_L_prev = chi(S_L_prev);
529 
530  double const p_FR = -chi_S_L * p_cap_ip;
531  variables[static_cast<int>(MPL::Variable::effective_pore_pressure)] =
532  p_FR;
533  variables_prev[static_cast<int>(
534  MPL::Variable::effective_pore_pressure)] =
535  -chi_S_L_prev * (p_cap_ip - p_cap_dot_ip * dt);
536 
537  // Set volumetric strain rate for the general case without swelling.
538  variables[static_cast<int>(MPL::Variable::volumetric_strain)]
539  .emplace<double>(Invariants::trace(_ip_data[ip].eps));
540  variables_prev[static_cast<int>(MPL::Variable::volumetric_strain)]
541  .emplace<double>(Invariants::trace(B * (u - u_dot * dt)));
542 
543  auto& phi = _ip_data[ip].porosity;
544  { // Porosity update
545  variables_prev[static_cast<int>(MPL::Variable::porosity)] =
546  _ip_data[ip].porosity_prev;
547  phi = medium->property(MPL::PropertyType::porosity)
548  .template value<double>(variables, variables_prev,
549  x_position, t, dt);
550  variables[static_cast<int>(MPL::Variable::porosity)] = phi;
551  }
552 
553  if (alpha < phi)
554  {
555  OGS_FATAL(
556  "RichardsMechanics: Biot-coefficient {} is smaller than "
557  "porosity {} in element/integration point {}/{}.",
558  alpha, phi, _element.getID(), ip);
559  }
560 
561  // Swelling and possibly volumetric strain rate update.
562  {
563  auto& sigma_sw = _ip_data[ip].sigma_sw;
564  auto const& sigma_sw_prev = _ip_data[ip].sigma_sw_prev;
565 
566  // If there is swelling, compute it. Update volumetric strain rate,
567  // s.t. it corresponds to the mechanical part only.
568  sigma_sw = sigma_sw_prev;
569  if (solid_phase.hasProperty(
571  {
572  using DimMatrix = Eigen::Matrix<double, 3, 3>;
573  auto const sigma_sw_dot =
574  MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
575  solid_phase
577  .template value<DimMatrix>(
578  variables, variables_prev, x_position, t, dt));
579  sigma_sw += sigma_sw_dot * dt;
580 
581  // !!! Misusing volumetric strain for mechanical volumetric
582  // strain just to update the transport porosity !!!
583  std::get<double>(variables[static_cast<int>(
584  MPL::Variable::volumetric_strain)]) +=
585  identity2.transpose() * C_el.inverse() * sigma_sw;
586  std::get<double>(variables_prev[static_cast<int>(
587  MPL::Variable::volumetric_strain)]) +=
588  identity2.transpose() * C_el.inverse() * sigma_sw_prev;
589  }
590 
591  if (medium->hasProperty(MPL::PropertyType::transport_porosity))
592  {
593  variables_prev[static_cast<int>(
595  _ip_data[ip].transport_porosity_prev;
596 
597  _ip_data[ip].transport_porosity =
598  medium->property(MPL::PropertyType::transport_porosity)
599  .template value<double>(variables, variables_prev,
600  x_position, t, dt);
601  variables[static_cast<int>(MPL::Variable::transport_porosity)] =
602  _ip_data[ip].transport_porosity;
603  }
604  else
605  {
606  variables[static_cast<int>(MPL::Variable::transport_porosity)] =
607  phi;
608  }
609  }
610 
611  double const k_rel =
613  .template value<double>(variables, x_position, t, dt);
614  auto const mu =
615  liquid_phase.property(MPL::PropertyType::viscosity)
616  .template value<double>(variables, x_position, t, dt);
617 
618  auto const& sigma_sw = _ip_data[ip].sigma_sw;
619  auto const& sigma_eff = _ip_data[ip].sigma_eff;
620 
621  // Set mechanical variables for the intrinsic permeability model
622  // For stress dependent permeability.
623  {
624  auto const sigma_total =
625  (sigma_eff - alpha * p_FR * identity2).eval();
626 
627  // For stress dependent permeability.
628  variables[static_cast<int>(MPL::Variable::total_stress)]
629  .emplace<SymmetricTensor>(
631  sigma_total));
632  }
633 
634  variables[static_cast<int>(
636  _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
637 
638  auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
639  medium->property(MPL::PropertyType::permeability)
640  .value(variables, x_position, t, dt));
641 
642  GlobalDimMatrixType const rho_K_over_mu =
643  K_intrinsic * rho_LR * k_rel / mu;
644 
645  //
646  // displacement equation, displacement part
647  //
648  eps_m.noalias() =
649  solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
650  ? eps + C_el.inverse() * sigma_sw
651  : eps;
652  variables[static_cast<int>(MPL::Variable::mechanical_strain)]
654  eps_m);
655 
656  _ip_data[ip].updateConstitutiveRelation(variables, t, x_position, dt,
657  temperature);
658 
659  // p_SR
660  variables[static_cast<int>(MPL::Variable::solid_grain_pressure)] =
661  p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
662  auto const rho_SR =
663  solid_phase.property(MPL::PropertyType::density)
664  .template value<double>(variables, x_position, t, dt);
665 
666  //
667  // displacement equation, displacement part
668  //
669  double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
670  rhs.template segment<displacement_size>(displacement_index).noalias() -=
671  (B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w;
672 
673  //
674  // pressure equation, pressure part.
675  //
676  auto const beta_LR = 1 / rho_LR *
677  liquid_phase.property(MPL::PropertyType::density)
678  .template dValue<double>(
679  variables, MPL::Variable::phase_pressure,
680  x_position, t, dt);
681 
682  double const a0 = S_L * (alpha - phi) * beta_SR;
683  // Volumetric average specific storage of the solid and fluid phases.
684  double const specific_storage =
685  DeltaS_L_Deltap_cap * (p_cap_ip * a0 - phi) +
686  S_L * (phi * beta_LR + a0);
687  M.template block<pressure_size, pressure_size>(pressure_index,
689  .noalias() += N_p.transpose() * rho_LR * specific_storage * N_p * w;
690 
691  K.template block<pressure_size, pressure_size>(pressure_index,
693  .noalias() += dNdx_p.transpose() * rho_K_over_mu * dNdx_p * w;
694 
695  rhs.template segment<pressure_size>(pressure_index).noalias() +=
696  dNdx_p.transpose() * rho_LR * rho_K_over_mu * b * w;
697 
698  //
699  // displacement equation, pressure part
700  //
701  K.template block<displacement_size, pressure_size>(displacement_index,
703  .noalias() -= B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
704 
705  //
706  // pressure equation, displacement part.
707  //
708  M.template block<pressure_size, displacement_size>(pressure_index,
710  .noalias() += N_p.transpose() * S_L * rho_LR * alpha *
711  identity2.transpose() * B * w;
712  }
713 
714  if (_process_data.apply_mass_lumping)
715  {
716  auto Mpp = M.template block<pressure_size, pressure_size>(
718  Mpp = Mpp.colwise().sum().eval().asDiagonal();
719  }
720 }
#define OGS_FATAL(...)
Definition: Error.h:26
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition: BMatrixPolicy.h:55
typename ShapeMatricesTypePressure::GlobalDimMatrixType GlobalDimMatrixType
@ rho
density. For some models, rho substitutes p
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
Definition: VariableType.h:108
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:48
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:23
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:32
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
void shapeFunctionInterpolate(const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
Definition: Interpolation.h:79
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
Definition: LinearBMatrix.h:42
static double trace(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Trace of the corresponding tensor.

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

◆ assembleWithJacobian()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

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

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

◆ assembleWithJacobianForDeformationEquations()

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

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

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

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

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

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

1535 {
1536  OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
1537 }

References OGS_FATAL.

◆ assembleWithJacobianForPressureEquations()

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

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

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

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

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

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

1518 {
1519  OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
1520 }

References OGS_FATAL.

◆ assembleWithJacobianForStaggeredScheme()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

1550 {
1551  // For the equations with pressure
1552  if (process_id == 0)
1553  {
1555  t, dt, local_x, local_xdot, dxdot_dx, dx_dx, local_M_data,
1556  local_K_data, local_b_data, local_Jac_data);
1557  return;
1558  }
1559 
1560  // For the equations with deformation
1562  t, dt, local_x, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
1563  local_b_data, local_Jac_data);
1564 }
void assembleWithJacobianForDeformationEquations(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, const double dxdot_dx, const double dx_dx, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
void assembleWithJacobianForPressureEquations(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, const double dxdot_dx, const double dx_dx, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)

◆ computeSecondaryVariableConcrete()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

1574 {
1575  auto p_L = local_x.template segment<pressure_size>(pressure_index);
1576  auto u = local_x.template segment<displacement_size>(displacement_index);
1577 
1578  auto p_L_dot = local_x_dot.template segment<pressure_size>(pressure_index);
1579  auto u_dot =
1580  local_x_dot.template segment<displacement_size>(displacement_index);
1581 
1582  auto const& identity2 = MathLib::KelvinVector::Invariants<
1584  DisplacementDim)>::identity2;
1585 
1586  auto const& medium = _process_data.media_map->getMedium(_element.getID());
1587  auto const& liquid_phase = medium->phase("AqueousLiquid");
1588  auto const& solid_phase = medium->phase("Solid");
1589  MPL::VariableArray variables;
1590  MPL::VariableArray variables_prev;
1591 
1592  ParameterLib::SpatialPosition x_position;
1593  x_position.setElementID(_element.getID());
1594 
1595  unsigned const n_integration_points =
1596  _integration_method.getNumberOfPoints();
1597 
1598  double saturation_avg = 0;
1599  double porosity_avg = 0;
1600 
1602  KV sigma_avg = KV::Zero();
1603 
1604  for (unsigned ip = 0; ip < n_integration_points; ip++)
1605  {
1606  x_position.setIntegrationPoint(ip);
1607  auto const& N_p = _ip_data[ip].N_p;
1608  auto const& N_u = _ip_data[ip].N_u;
1609  auto const& dNdx_u = _ip_data[ip].dNdx_u;
1610 
1611  auto const x_coord =
1612  NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
1614  _element, N_u);
1615  auto const B =
1616  LinearBMatrix::computeBMatrix<DisplacementDim,
1617  ShapeFunctionDisplacement::NPOINTS,
1618  typename BMatricesType::BMatrixType>(
1619  dNdx_u, N_u, x_coord, _is_axially_symmetric);
1620 
1621  double p_cap_ip;
1622  NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
1623 
1624  double p_cap_dot_ip;
1625  NumLib::shapeFunctionInterpolate(-p_L_dot, N_p, p_cap_dot_ip);
1626 
1627  variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
1628  p_cap_ip;
1629  variables[static_cast<int>(MPL::Variable::phase_pressure)] = -p_cap_ip;
1630 
1631  auto const temperature =
1633  .template value<double>(variables, x_position, t, dt);
1634  variables[static_cast<int>(MPL::Variable::temperature)] = temperature;
1635 
1636  auto& eps = _ip_data[ip].eps;
1637  eps.noalias() = B * u;
1638  auto& eps_m = _ip_data[ip].eps_m;
1639  auto& S_L = _ip_data[ip].saturation;
1640  auto const S_L_prev = _ip_data[ip].saturation_prev;
1641  S_L = medium->property(MPL::PropertyType::saturation)
1642  .template value<double>(variables, x_position, t, dt);
1643  variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
1644  variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
1645  S_L_prev;
1646 
1647  auto const chi = [medium, x_position, t, dt](double const S_L)
1648  {
1649  MPL::VariableArray vs;
1650  vs.fill(std::numeric_limits<double>::quiet_NaN());
1651  vs[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
1652  return medium->property(MPL::PropertyType::bishops_effective_stress)
1653  .template value<double>(vs, x_position, t, dt);
1654  };
1655  double const chi_S_L = chi(S_L);
1656  double const chi_S_L_prev = chi(S_L_prev);
1657 
1658  auto const alpha =
1659  medium->property(MPL::PropertyType::biot_coefficient)
1660  .template value<double>(variables, x_position, t, dt);
1661 
1662  auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
1663  t, x_position, dt, temperature);
1664 
1665  auto const beta_SR =
1666  (1 - alpha) /
1667  _ip_data[ip].solid_material.getBulkModulus(t, x_position, &C_el);
1668  variables[static_cast<int>(MPL::Variable::grain_compressibility)] =
1669  beta_SR;
1670 
1671  variables[static_cast<int>(MPL::Variable::effective_pore_pressure)] =
1672  -chi_S_L * p_cap_ip;
1673  variables_prev[static_cast<int>(
1674  MPL::Variable::effective_pore_pressure)] =
1675  -chi_S_L_prev * (p_cap_ip - p_cap_dot_ip * dt);
1676 
1677  // Set volumetric strain rate for the general case without swelling.
1678  variables[static_cast<int>(MPL::Variable::volumetric_strain)]
1679  .emplace<double>(Invariants::trace(_ip_data[ip].eps));
1680  variables_prev[static_cast<int>(MPL::Variable::volumetric_strain)]
1681  .emplace<double>(Invariants::trace(B * (u - u_dot * dt)));
1682 
1683  auto& phi = _ip_data[ip].porosity;
1684  { // Porosity update
1685  variables_prev[static_cast<int>(MPL::Variable::porosity)] =
1686  _ip_data[ip].porosity_prev;
1687  phi = medium->property(MPL::PropertyType::porosity)
1688  .template value<double>(variables, variables_prev,
1689  x_position, t, dt);
1690  variables[static_cast<int>(MPL::Variable::porosity)] = phi;
1691  }
1692 
1693  auto const mu =
1694  liquid_phase.property(MPL::PropertyType::viscosity)
1695  .template value<double>(variables, x_position, t, dt);
1696  auto const rho_LR =
1697  liquid_phase.property(MPL::PropertyType::density)
1698  .template value<double>(variables, x_position, t, dt);
1699 
1700  // Swelling and possibly volumetric strain rate update.
1701  updateSwellingStressAndVolumetricStrain<DisplacementDim>(
1702  _ip_data[ip], *medium, solid_phase, C_el, rho_LR, mu,
1703  _process_data.micro_porosity_parameters, alpha, phi, p_cap_ip,
1704  variables, variables_prev, x_position, t, dt);
1705 
1706  if (medium->hasProperty(MPL::PropertyType::transport_porosity))
1707  {
1708  if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
1709  {
1710  variables_prev[static_cast<int>(
1712  _ip_data[ip].transport_porosity_prev;
1713 
1714  _ip_data[ip].transport_porosity =
1715  medium->property(MPL::PropertyType::transport_porosity)
1716  .template value<double>(variables, variables_prev,
1717  x_position, t, dt);
1718  variables[static_cast<int>(MPL::Variable::transport_porosity)] =
1719  _ip_data[ip].transport_porosity;
1720  }
1721  }
1722  else
1723  {
1724  variables[static_cast<int>(MPL::Variable::transport_porosity)] =
1725  phi;
1726  }
1727 
1728  // Set mechanical variables for the intrinsic permeability model
1729  // For stress dependent permeability.
1730  {
1731  auto const sigma_total = (_ip_data[ip].sigma_eff +
1732  alpha * chi_S_L * identity2 * p_cap_ip)
1733  .eval();
1734  // For stress dependent permeability.
1735  variables[static_cast<int>(MPL::Variable::total_stress)]
1736  .emplace<SymmetricTensor>(
1738  sigma_total));
1739  }
1740 
1741  variables[static_cast<int>(
1743  _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
1744 
1745  auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
1746  medium->property(MPL::PropertyType::permeability)
1747  .value(variables, x_position, t, dt));
1748 
1749  double const k_rel =
1751  .template value<double>(variables, x_position, t, dt);
1752 
1753  GlobalDimMatrixType const K_over_mu = k_rel * K_intrinsic / mu;
1754 
1755  auto const& sigma_eff = _ip_data[ip].sigma_eff;
1756  double const p_FR = -chi_S_L * p_cap_ip;
1757  // p_SR
1758  variables[static_cast<int>(MPL::Variable::solid_grain_pressure)] =
1759  p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
1760  auto const rho_SR =
1761  solid_phase.property(MPL::PropertyType::density)
1762  .template value<double>(variables, x_position, t, dt);
1763  _ip_data[ip].dry_density_solid = (1 - phi) * rho_SR;
1764 
1765  auto& sigma_sw = _ip_data[ip].sigma_sw;
1766 
1767  eps_m.noalias() =
1768  solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
1769  ? eps + C_el.inverse() * sigma_sw
1770  : eps;
1771  variables[static_cast<int>(MPL::Variable::mechanical_strain)]
1773  eps_m);
1774 
1775  _ip_data[ip].updateConstitutiveRelation(variables, t, x_position, dt,
1776  temperature);
1777 
1778  auto const& b = _process_data.specific_body_force;
1779 
1780  // Compute the velocity
1781  auto const& dNdx_p = _ip_data[ip].dNdx_p;
1782  _ip_data[ip].v_darcy.noalias() =
1783  -K_over_mu * dNdx_p * p_L + rho_LR * K_over_mu * b;
1784 
1785  saturation_avg += S_L;
1786  porosity_avg += phi;
1787  sigma_avg += sigma_eff;
1788  }
1789  saturation_avg /= n_integration_points;
1790  porosity_avg /= n_integration_points;
1791  sigma_avg /= n_integration_points;
1792 
1793  (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
1794  (*_process_data.element_porosity)[_element.getID()] = porosity_avg;
1795 
1796  Eigen::Map<KV>(&(*_process_data.element_stresses)[_element.getID() *
1797  KV::RowsAtCompileTime]) =
1799 
1801  ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1802  DisplacementDim>(_element, _is_axially_symmetric, p_L,
1803  *_process_data.pressure_interpolated);
1804 }
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)

References MaterialPropertyLib::alpha, MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::bishops_effective_stress, MaterialPropertyLib::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), MaterialPropertyLib::density, MaterialPropertyLib::equivalent_plastic_strain, NumLib::interpolateToHigherOrderNodes(), NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, MaterialPropertyLib::reference_temperature, MaterialPropertyLib::relative_permeability, MaterialPropertyLib::saturation, MaterialPropertyLib::saturation_micro, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), NumLib::shapeFunctionInterpolate(), MaterialPropertyLib::swelling_stress_rate, MaterialPropertyLib::temperature, MaterialPropertyLib::transport_porosity, and MaterialPropertyLib::viscosity.

◆ getEpsilon()

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

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

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

1293 {
1294  constexpr int kelvin_vector_size =
1296 
1297  return transposeInPlace<kelvin_vector_size>(
1298  [this](std::vector<double>& values)
1299  { return getIntPtEpsilon(0, {}, {}, values); });
1300 }
std::vector< double > const & getIntPtEpsilon(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

References MathLib::KelvinVector::kelvin_vector_dimensions().

◆ getIntPtDarcyVelocity()

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

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

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

1343 {
1344  unsigned const n_integration_points =
1345  _integration_method.getNumberOfPoints();
1346 
1347  cache.clear();
1348  auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
1349  double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
1350  cache, DisplacementDim, n_integration_points);
1351 
1352  for (unsigned ip = 0; ip < n_integration_points; ip++)
1353  {
1354  cache_matrix.col(ip).noalias() = _ip_data[ip].v_darcy;
1355  }
1356 
1357  return cache;
1358 }

References MathLib::createZeroedMatrix().

◆ getIntPtDryDensitySolid()

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

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

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

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

References ProcessLib::getIntegrationPointScalarData().

◆ getIntPtEpsilon()

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

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

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

1312 {
1313  return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
1314  _ip_data, &IpData::eps, cache);
1315 }

◆ getIntPtMicroPressure()

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

◆ getIntPtMicroSaturation()

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

◆ getIntPtPorosity()

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

◆ getIntPtSaturation()

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

◆ getIntPtSigma()

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

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

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

1239 {
1240  return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
1241  _ip_data, &IpData::sigma_eff, cache);
1242 }

◆ getIntPtSwellingStress()

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

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

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

1268 {
1269  constexpr int kelvin_vector_size =
1271  auto const n_integration_points = _ip_data.size();
1272 
1273  cache.clear();
1274  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
1275  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
1276  cache, kelvin_vector_size, n_integration_points);
1277 
1278  for (unsigned ip = 0; ip < n_integration_points; ++ip)
1279  {
1280  auto const& sigma_sw = _ip_data[ip].sigma_sw;
1281  cache_mat.col(ip) =
1283  }
1284 
1285  return cache;
1286 }

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

◆ getIntPtTransportPorosity()

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

◆ getMaterialStateVariableInternalState()

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

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

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

1327 {
1329  _ip_data, &IpData::material_state_variables, get_values_span,
1330  n_components);
1331 }
std::vector< double > getIntegrationPointDataMaterialStateVariables(std::vector< IntegrationPointData, Eigen::aligned_allocator< IntegrationPointData >> const &ip_data_vector, MemberType member, std::function< BaseLib::DynamicSpan< double >(MaterialStateVariables &)> get_values_span, int const n_components)
std::unique_ptr< typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables > material_state_variables

References ProcessLib::getIntegrationPointDataMaterialStateVariables().

◆ getMaterialStateVariablesAt()

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

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

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

1823 {
1824  return *_ip_data[integration_point].material_state_variables;
1825 }

◆ getMicroPressure()

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

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

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

1417 {
1418  std::vector<double> result;
1419  getIntPtMicroPressure(0, {}, {}, result);
1420  return result;
1421 }
std::vector< double > const & getIntPtMicroPressure(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

◆ getMicroSaturation()

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

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

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

1391 {
1392  std::vector<double> result;
1393  getIntPtMicroSaturation(0, {}, {}, result);
1394  return result;
1395 }
std::vector< double > const & getIntPtMicroSaturation(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

◆ getNumberOfIntegrationPoints()

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

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

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

1811 {
1812  return _integration_method.getNumberOfPoints();
1813 }

◆ getPorosity()

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

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

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

1443 {
1444  std::vector<double> result;
1445  getIntPtPorosity(0, {}, {}, result);
1446  return result;
1447 }
std::vector< double > const & getIntPtPorosity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

◆ getSaturation()

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

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

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

1365 {
1366  std::vector<double> result;
1367  getIntPtSaturation(0, {}, {}, result);
1368  return result;
1369 }
std::vector< double > const & getIntPtSaturation(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

◆ getShapeMatrix()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int DisplacementDim>
Eigen::Map<const Eigen::RowVectorXd> ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::getShapeMatrix ( const unsigned  integration_point) const
inlineoverridevirtual

Provides the shape matrix at the given integration point.

Implements NumLib::ExtrapolatableElement.

Definition at line 168 of file RichardsMechanicsFEM.h.

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

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

◆ getSigma()

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

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

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

1220 {
1221  constexpr int kelvin_vector_size =
1223 
1224  return transposeInPlace<kelvin_vector_size>(
1225  [this](std::vector<double>& values)
1226  { return getIntPtSigma(0, {}, {}, values); });
1227 }
std::vector< double > const & getIntPtSigma(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

References MathLib::KelvinVector::kelvin_vector_dimensions().

◆ getSwellingStress()

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

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

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

1249 {
1250  constexpr int kelvin_vector_size =
1252 
1253  return transposeInPlace<kelvin_vector_size>(
1254  [this](std::vector<double>& values)
1255  { return getIntPtSwellingStress(0, {}, {}, values); });
1256 }
std::vector< double > const & getIntPtSwellingStress(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

References MathLib::KelvinVector::kelvin_vector_dimensions().

◆ getTransportPorosity()

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

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

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

1469 {
1470  std::vector<double> result;
1471  getIntPtTransportPorosity(0, {}, {}, result);
1472  return result;
1473 }
std::vector< double > const & getIntPtTransportPorosity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

◆ initializeConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int DisplacementDim>
void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::initializeConcrete ( )
inlineoverridevirtual

Set initial stress from parameter.

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 120 of file RichardsMechanicsFEM.h.

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

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

◆ postTimestepConcrete()

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

◆ setInitialConditionsConcrete()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

304 {
305  assert(local_x.size() == pressure_size + displacement_size);
306 
307  auto p_L =
308  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
309  pressure_size> const>(local_x.data() + pressure_index,
310  pressure_size);
311 
312  constexpr double dt = std::numeric_limits<double>::quiet_NaN();
313  auto const& medium = _process_data.media_map->getMedium(_element.getID());
314  MPL::VariableArray variables;
315 
317  x_position.setElementID(_element.getID());
318 
319  auto const& solid_phase = medium->phase("Solid");
320 
321  unsigned const n_integration_points =
322  _integration_method.getNumberOfPoints();
323  for (unsigned ip = 0; ip < n_integration_points; ip++)
324  {
325  x_position.setIntegrationPoint(ip);
326 
327  auto const& N_p = _ip_data[ip].N_p;
328 
329  double p_cap_ip;
330  NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
331 
332  variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
333  p_cap_ip;
334  variables[static_cast<int>(MPL::Variable::phase_pressure)] = -p_cap_ip;
335  _ip_data[ip].liquid_pressure_m_prev = -p_cap_ip;
336  _ip_data[ip].liquid_pressure_m = -p_cap_ip;
337 
338  auto const temperature =
340  .template value<double>(variables, x_position, t, dt);
341  variables[static_cast<int>(MPL::Variable::temperature)] = temperature;
342 
343  _ip_data[ip].saturation_prev =
344  medium->property(MPL::PropertyType::saturation)
345  .template value<double>(variables, x_position, t, dt);
346 
347  if (medium->hasProperty(MPL::PropertyType::saturation_micro))
348  {
349  MPL::VariableArray vars;
350  vars[static_cast<int>(MPL::Variable::capillary_pressure)] =
351  p_cap_ip;
352  double const S_L_m =
353  medium->property(MPL::PropertyType::saturation_micro)
354  .template value<double>(vars, x_position, t, dt);
355  _ip_data[ip].saturation_m_prev = S_L_m;
356  }
357 
358  // Set eps_m_prev from potentially non-zero eps and sigma_sw from
359  // restart.
360  auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
361  t, x_position, dt, temperature);
362  auto& eps = _ip_data[ip].eps;
363  auto& sigma_sw = _ip_data[ip].sigma_sw;
364 
365  _ip_data[ip].eps_m_prev.noalias() =
366  solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
367  ? eps + C_el.inverse() * sigma_sw
368  : eps;
369  }
370 }

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

◆ setIPDataInitialConditions()

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

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

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

218 {
219  if (integration_order !=
220  static_cast<int>(_integration_method.getIntegrationOrder()))
221  {
222  OGS_FATAL(
223  "Setting integration point initial conditions; The integration "
224  "order of the local assembler for element {:d} is different "
225  "from the integration order in the initial condition.",
226  _element.getID());
227  }
228 
229  if (name == "sigma_ip")
230  {
231  if (_process_data.initial_stress != nullptr)
232  {
233  OGS_FATAL(
234  "Setting initial conditions for stress from integration "
235  "point data and from a parameter '{:s}' is not possible "
236  "simultaneously.",
237  _process_data.initial_stress->name);
238  }
239  return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
240  values, _ip_data, &IpData::sigma_eff);
241  }
242 
243  if (name == "saturation_ip")
244  {
247  }
248  if (name == "porosity_ip")
249  {
252  }
253  if (name == "transport_porosity_ip")
254  {
257  }
258  if (name == "swelling_stress_ip")
259  {
260  return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
261  values, _ip_data, &IpData::sigma_sw);
262  }
263  if (name == "epsilon_ip")
264  {
265  return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
266  values, _ip_data, &IpData::eps);
267  }
268  if (name.starts_with("material_state_variable_") && name.ends_with("_ip"))
269  {
270  std::string const variable_name = name.substr(24, name.size() - 24 - 3);
271  DBUG("Setting material state variable '{:s}'", variable_name);
272 
273  // Using first ip data for solid material. TODO (naumov) move solid
274  // material into element, store only material state in IPs.
275  auto const& internal_variables =
276  _ip_data[0].solid_material.getInternalVariables();
277  if (auto const iv =
278  std::find_if(begin(internal_variables), end(internal_variables),
279  [&variable_name](auto const& iv)
280  { return iv.name == variable_name; });
281  iv != end(internal_variables))
282  {
285  iv->reference);
286  }
287 
288  ERR("Could not find variable {:s} in solid material model's internal "
289  "variables.",
290  variable_name);
291  }
292  return 0;
293 }
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:42
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:27
std::size_t setIntegrationPointScalarData(double const *values, std::vector< IntegrationPointData, Eigen::aligned_allocator< IntegrationPointData >> &ip_data, MemberType member)
std::size_t setIntegrationPointDataMaterialStateVariables(double const *values, std::vector< IntegrationPointData, Eigen::aligned_allocator< IntegrationPointData >> &ip_data_vector, MemberType member, std::function< BaseLib::DynamicSpan< double >(MaterialStateVariables &)> get_values_span)

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

Member Data Documentation

◆ _element

◆ _integration_method

◆ _ip_data

◆ _is_axially_symmetric

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

Definition at line 333 of file RichardsMechanicsFEM.h.

◆ _process_data

◆ _secondary_data

◆ displacement_index

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

Definition at line 340 of file RichardsMechanicsFEM.h.

◆ displacement_size

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int DisplacementDim>
const int ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, DisplacementDim >::displacement_size
staticprivate

◆ KelvinVectorSize

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

Definition at line 68 of file RichardsMechanicsFEM.h.

◆ pressure_index

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

Definition at line 338 of file RichardsMechanicsFEM.h.

◆ pressure_size

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

Definition at line 339 of file RichardsMechanicsFEM.h.


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