Loading [MathJax]/extensions/tex2jax.js
OGS
List of all members | Public Types | Public Member Functions | Static Public Attributes | Private Types | Private Member Functions | Private Attributes | Static Private Attributes
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > Class Template Reference

Detailed Description

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
class ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >

Definition at line 166 of file HydroMechanicsFEM.h.

#include <HydroMechanicsFEM.h>

Inheritance diagram for ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >:
[legend]
Collaboration diagram for ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >:
[legend]

Public Types

using ShapeMatricesTypeDisplacement
 
using ShapeMatricesTypePressure
 
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>
 
using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>
 

Public Member Functions

 HydroMechanicsLocalAssembler (HydroMechanicsLocalAssembler const &)=delete
 
 HydroMechanicsLocalAssembler (HydroMechanicsLocalAssembler &&)=delete
 
 HydroMechanicsLocalAssembler (MeshLib::Element const &e, std::size_t const, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, HydroMechanicsProcessData< DisplacementDim > &process_data)
 
std::size_t setIPDataInitialConditions (std::string_view const name, double const *values, int const integration_order) override
 Returns number of read integration points.
 
void assemble (double const, double const, std::vector< double > const &, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
 
void assembleWithJacobian (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_rhs_data, std::vector< double > &local_Jac_data) override
 
void assembleWithJacobianForStaggeredScheme (const double t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
 
void initializeConcrete () override
 
void postTimestepConcrete (Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, int const process_id) override
 
void computeSecondaryVariableConcrete (double const t, double const dt, Eigen::VectorXd const &local_xs, Eigen::VectorXd const &local_x_prev) override
 
void postNonLinearSolverConcrete (Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, int const process_id) override
 
void setInitialConditionsConcrete (Eigen::VectorXd const local_x, double const t, int const process_id) override
 
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix (const unsigned integration_point) const override
 Provides the shape matrix at the given integration point.
 
std::vector< double > getSigma () const override
 
std::vector< double > getEpsilon () const override
 
std::vector< double > getStrainRateVariable () 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 > const & getIntPtSigma (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
std::vector< double > const & getIntPtEpsilon (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
- Public Member Functions inherited from ProcessLib::HydroMechanics::LocalAssemblerInterface< DisplacementDim >
- Public Member Functions inherited from ProcessLib::LocalAssemblerInterface
virtual ~LocalAssemblerInterface ()=default
 
virtual void setInitialConditions (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, int const process_id)
 
virtual void initialize (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
 
virtual void preAssemble (double const, double const, std::vector< double > const &)
 
virtual void assembleForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
 
virtual void computeSecondaryVariable (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, int const process_id)
 
virtual void preTimestep (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
 
virtual void postTimestep (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
 
void postNonLinearSolver (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
 
virtual Eigen::Vector3d getFlux (MathLib::Point3d const &, double const, std::vector< double > const &) const
 
virtual Eigen::Vector3d getFlux (MathLib::Point3d const &, double const, std::vector< std::vector< double > > const &) const
 Fits to staggered scheme.
 
virtual std::optional< VectorSegment > getVectorDeformationSegment () const
 
- Public Member Functions inherited from NumLib::ExtrapolatableElement
virtual ~ExtrapolatableElement ()=default
 

Static Public Attributes

static int const KelvinVectorSize
 
static constexpr auto & N_u_op
 

Private Types

using BMatricesType
 
using IpData
 

Private Member Functions

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

Private Attributes

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

Static Private Attributes

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

Member Typedef Documentation

◆ BMatricesType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::BMatricesType
private
Initial value:
BMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>

Definition at line 386 of file HydroMechanicsFEM.h.

◆ Invariants

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

Definition at line 179 of file HydroMechanicsFEM.h.

◆ IpData

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

Definition at line 388 of file HydroMechanicsFEM.h.

◆ ShapeMatricesTypeDisplacement

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ShapeMatricesTypeDisplacement
Initial value:
ShapeMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>

Definition at line 170 of file HydroMechanicsFEM.h.

◆ ShapeMatricesTypePressure

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

Definition at line 174 of file HydroMechanicsFEM.h.

◆ SymmetricTensor

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

Definition at line 181 of file HydroMechanicsFEM.h.

Constructor & Destructor Documentation

◆ HydroMechanicsLocalAssembler() [1/3]

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

◆ HydroMechanicsLocalAssembler() [2/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::HydroMechanicsLocalAssembler ( HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > && )
delete

◆ HydroMechanicsLocalAssembler() [3/3]

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

Definition at line 37 of file HydroMechanicsFEM-impl.h.

44 : _process_data(process_data),
45 _integration_method(integration_method),
46 _element(e),
47 _is_axially_symmetric(is_axially_symmetric)
48{
49 unsigned const n_integration_points =
50 _integration_method.getNumberOfPoints();
51
52 _ip_data.reserve(n_integration_points);
53 _secondary_data.N_u.resize(n_integration_points);
54
55 auto const shape_matrices_u =
56 NumLib::initShapeMatrices<ShapeFunctionDisplacement,
57 ShapeMatricesTypeDisplacement,
58 DisplacementDim>(e, is_axially_symmetric,
59 _integration_method);
60
61 auto const shape_matrices_p =
62 NumLib::initShapeMatrices<ShapeFunctionPressure,
63 ShapeMatricesTypePressure, DisplacementDim>(
64 e, is_axially_symmetric, _integration_method);
65
66 auto const& solid_material =
67 MaterialLib::Solids::selectSolidConstitutiveRelation(
68 _process_data.solid_materials, _process_data.material_ids,
69 e.getID());
70
71 for (unsigned ip = 0; ip < n_integration_points; ip++)
72 {
73 _ip_data.emplace_back(solid_material);
74 auto& ip_data = _ip_data[ip];
75 auto const& sm_u = shape_matrices_u[ip];
76 ip_data.integration_weight =
77 _integration_method.getWeightedPoint(ip).getWeight() *
78 sm_u.integralMeasure * sm_u.detJ;
79
80 // Initialize current time step values
81 static const int kelvin_vector_size =
82 MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
83 ip_data.sigma_eff.setZero(kelvin_vector_size);
84 ip_data.eps.setZero(kelvin_vector_size);
85
86 // Previous time step values are not initialized and are set later.
87 ip_data.eps_prev.resize(kelvin_vector_size);
88 ip_data.sigma_eff_prev.resize(kelvin_vector_size);
89
90 ip_data.N_u = sm_u.N;
91 ip_data.dNdx_u = sm_u.dNdx;
92
93 ip_data.N_p = shape_matrices_p[ip].N;
94 ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
95
96 _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
97 }
98}
MathLib::WeightedPoint::getWeight
constexpr double getWeight() const
Definition WeightedPoint.h:82
NumLib::GenericIntegrationMethod::getWeightedPoint
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
Definition GenericIntegrationMethod.h:36
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_integration_method
NumLib::GenericIntegrationMethod const & _integration_method
Definition HydroMechanicsFEM.h:394
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_process_data
HydroMechanicsProcessData< DisplacementDim > & _process_data
Definition HydroMechanicsFEM.h:384
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_secondary_data
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
Definition HydroMechanicsFEM.h:399
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_ip_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
Definition HydroMechanicsFEM.h:392
MaterialLib::Solids::selectSolidConstitutiveRelation
auto & selectSolidConstitutiveRelation(SolidMaterialsMap const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
Definition SelectSolidConstitutiveRelation.h:36
MathLib::KelvinVector::kelvin_vector_dimensions
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition KelvinVector.h:24
NumLib::initShapeMatrices
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)
Definition InitShapeMatrices.h:54
ProcessLib::HydroMechanics::SecondaryData::N_u
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u
Definition HydroMechanicsFEM.h:161

References ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_integration_method, ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data, ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_process_data, ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_secondary_data, MeshLib::Element::getID(), NumLib::GenericIntegrationMethod::getNumberOfPoints(), MathLib::WeightedPoint::getWeight(), NumLib::GenericIntegrationMethod::getWeightedPoint(), NumLib::initShapeMatrices(), MathLib::KelvinVector::kelvin_vector_dimensions(), ProcessLib::HydroMechanics::SecondaryData< ShapeMatrixType >::N_u, and MaterialLib::Solids::selectSolidConstitutiveRelation().

Member Function Documentation

◆ assemble()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assemble ( double const ,
double const ,
std::vector< double > const & ,
std::vector< double > const & ,
std::vector< double > & ,
std::vector< double > & ,
std::vector< double > &  )
inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 203 of file HydroMechanicsFEM.h.

209 {
210 OGS_FATAL(
211 "HydroMechanicsLocalAssembler: assembly without jacobian is not "
212 "implemented.");
213 }
OGS_FATAL
#define OGS_FATAL(...)
Definition Error.h:26

References OGS_FATAL.

◆ assembleWithJacobian()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 103 of file HydroMechanicsFEM-impl.h.

109{
110 assert(local_x.size() == pressure_size + displacement_size);
111
112 auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
113 pressure_size> const>(local_x.data() + pressure_index, pressure_size);
114
115 auto u =
116 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
117 displacement_size> const>(local_x.data() + displacement_index,
118 displacement_size);
119
120 auto p_prev =
121 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
122 pressure_size> const>(local_x_prev.data() + pressure_index,
123 pressure_size);
124 auto u_prev =
125 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
126 displacement_size> const>(local_x_prev.data() + displacement_index,
127 displacement_size);
128
129 auto local_Jac = MathLib::createZeroedMatrix<
130 typename ShapeMatricesTypeDisplacement::template MatrixType<
131 displacement_size + pressure_size,
132 displacement_size + pressure_size>>(
133 local_Jac_data, displacement_size + pressure_size,
134 displacement_size + pressure_size);
135
136 auto local_rhs = MathLib::createZeroedVector<
137 typename ShapeMatricesTypeDisplacement::template VectorType<
138 displacement_size + pressure_size>>(
139 local_rhs_data, displacement_size + pressure_size);
140
141 typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
142 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
143 pressure_size);
144
145 typename ShapeMatricesTypePressure::NodalMatrixType storage_p =
146 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
147 pressure_size);
148
149 typename ShapeMatricesTypePressure::NodalMatrixType add_p_derivative =
150 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
151 pressure_size);
152
153 typename ShapeMatricesTypeDisplacement::template MatrixType<
154 displacement_size, pressure_size>
155 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
156 displacement_size, pressure_size>::Zero(displacement_size,
157 pressure_size);
158
159 typename ShapeMatricesTypeDisplacement::template MatrixType<
160 pressure_size, displacement_size>
161 Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
162 pressure_size, displacement_size>::Zero(pressure_size,
163 displacement_size);
164
165 typename ShapeMatricesTypeDisplacement::template MatrixType<
166 pressure_size, displacement_size>
167 Kpu_k = ShapeMatricesTypeDisplacement::template MatrixType<
168 pressure_size, displacement_size>::Zero(pressure_size,
169 displacement_size);
170
171 auto const& solid_material =
172 MaterialLib::Solids::selectSolidConstitutiveRelation(
173 _process_data.solid_materials, _process_data.material_ids,
174 _element.getID());
175
176 ParameterLib::SpatialPosition x_position;
177 x_position.setElementID(_element.getID());
178
179 unsigned const n_integration_points =
180 _integration_method.getNumberOfPoints();
181
182 auto const& b = _process_data.specific_body_force;
183 auto const& medium = _process_data.media_map.getMedium(_element.getID());
184 auto const& solid = medium->phase("Solid");
185 auto const& fluid = fluidPhase(*medium);
186 MPL::VariableArray vars;
187 auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
188 : vars.liquid_phase_pressure;
189
190 auto const T_ref =
191 medium->property(MPL::PropertyType::reference_temperature)
192 .template value<double>(vars, x_position, t, dt);
193 vars.temperature = T_ref;
194
195 auto const& identity2 = Invariants::identity2;
196
197 for (unsigned ip = 0; ip < n_integration_points; ip++)
198 {
199 auto const& w = _ip_data[ip].integration_weight;
200
201 auto const& N_u = _ip_data[ip].N_u;
202 auto const& dNdx_u = _ip_data[ip].dNdx_u;
203
204 auto const& N_p = _ip_data[ip].N_p;
205 auto const& dNdx_p = _ip_data[ip].dNdx_p;
206
207 auto const x_coord =
208 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
209 ShapeMatricesTypeDisplacement>(
210 _element, N_u);
211 auto const B =
212 LinearBMatrix::computeBMatrix<DisplacementDim,
213 ShapeFunctionDisplacement::NPOINTS,
214 typename BMatricesType::BMatrixType>(
215 dNdx_u, N_u, x_coord, _is_axially_symmetric);
216
217 auto& eps = _ip_data[ip].eps;
218 eps.noalias() = B * u;
219 auto const& sigma_eff = _ip_data[ip].sigma_eff;
220
221 double const p_int_pt = N_p.dot(p);
222
223 phase_pressure = p_int_pt;
224
225 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
226 t, x_position, dt, T_ref);
227 auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
228
229 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
230 .template value<double>(vars, x_position, t, dt);
231
232 auto const rho_sr =
233 solid.property(MPL::PropertyType::density)
234 .template value<double>(vars, x_position, t, dt);
235 auto const porosity =
236 medium->property(MPL::PropertyType::porosity)
237 .template value<double>(vars, x_position, t, dt);
238
239 auto const [rho_fr, mu] =
240 MaterialPropertyLib::getFluidDensityAndViscosity(t, dt, x_position,
241 fluid, vars);
242
243 auto const beta_p =
244 fluid.property(MPL::PropertyType::density)
245 .template dValue<double>(vars, _process_data.phase_variable,
246 x_position, t, dt) /
247 rho_fr;
248
249 // Set mechanical variables for the intrinsic permeability model
250 // For stress dependent permeability.
251 {
252 auto const sigma_total =
253 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
254
255 vars.total_stress.emplace<SymmetricTensor>(
256 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
257 sigma_total));
258 }
259 // For strain dependent permeability
260 vars.volumetric_strain = Invariants::trace(eps);
261 vars.equivalent_plastic_strain =
262 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
263 vars.mechanical_strain
264 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
265 eps);
266
267 auto const K = MPL::formEigenTensor<DisplacementDim>(
268 medium->property(MPL::PropertyType::permeability)
269 .value(vars, x_position, t, dt));
270 auto const dkde = MPL::formEigenTensor<
271 MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim)>(
272 (*medium)[MPL::PropertyType::permeability].dValue(
273 vars, MPL::Variable::mechanical_strain, x_position, t, dt));
274
275 auto const K_over_mu = K / mu;
276
277 auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
278 dt, u, T_ref);
279
280 //
281 // displacement equation, displacement part
282 //
283 local_Jac
284 .template block<displacement_size, displacement_size>(
285 displacement_index, displacement_index)
286 .noalias() += B.transpose() * C * B * w;
287
288 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
289 local_rhs.template segment<displacement_size>(displacement_index)
290 .noalias() -=
291 (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
292
293 //
294 // displacement equation, pressure part
295 //
296 Kup.noalias() += B.transpose() * alpha * identity2 * N_p * w;
297
298 //
299 // pressure equation, pressure part.
300 //
301 laplace_p.noalias() +=
302 rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
303
304 storage_p.noalias() +=
305 rho_fr * N_p.transpose() * N_p * w *
306 ((alpha - porosity) * (1.0 - alpha) / K_S + porosity * beta_p);
307
308 // density dependence on pressure evaluated for Darcy-term,
309 // for laplace and storage terms this dependence is neglected
310 add_p_derivative.noalias() += rho_fr * beta_p * dNdx_p.transpose() *
311 K_over_mu *
312 (dNdx_p * p - 2.0 * rho_fr * b) * N_p * w;
313
314 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
315 dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
316
317 //
318 // pressure equation, displacement part.
319 //
320 Kpu.noalias() +=
321 rho_fr * alpha * N_p.transpose() * identity2.transpose() * B * w;
322
323 Kpu_k.noalias() +=
324 dNdx_p.transpose() *
325 MathLib::KelvinVector::liftVectorToKelvin<DisplacementDim>(
326 dNdx_p * p - rho_fr * b) *
327 dkde * B * rho_fr / mu * w;
328 }
329 // displacement equation, pressure part
330 local_Jac
331 .template block<displacement_size, pressure_size>(displacement_index,
332 pressure_index)
333 .noalias() = -Kup;
334
335 if (_process_data.mass_lumping)
336 {
337 storage_p = storage_p.colwise().sum().eval().asDiagonal();
338
339 if constexpr (pressure_size == displacement_size)
340 {
341 Kpu = Kpu.colwise().sum().eval().asDiagonal();
342 Kpu_k = Kpu_k.colwise().sum().eval().asDiagonal();
343 }
344 }
345
346 // pressure equation, pressure part.
347 local_Jac
348 .template block<pressure_size, pressure_size>(pressure_index,
349 pressure_index)
350 .noalias() += laplace_p + storage_p / dt + add_p_derivative;
351
352 // pressure equation, displacement part.
353 local_Jac
354 .template block<pressure_size, displacement_size>(pressure_index,
355 displacement_index)
356 .noalias() += Kpu / dt + Kpu_k;
357
358 // pressure equation
359 local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
360 laplace_p * p + storage_p * (p - p_prev) / dt + Kpu * (u - u_prev) / dt;
361
362 // displacement equation
363 local_rhs.template segment<displacement_size>(displacement_index)
364 .noalias() += Kup * p;
365}
MeshLib::Element::getID
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:89
ParameterLib::SpatialPosition::setElementID
void setElementID(std::size_t element_id)
Definition SpatialPosition.h:76
ProcessLib::BMatrixPolicyType::BMatrixType
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition BMatrixPolicy.h:55
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::SymmetricTensor
Eigen::Matrix< double, KelvinVectorSize, 1 > SymmetricTensor
Definition HydroMechanicsFEM.h:181
MaterialLib::Fluid::PropertyVariableType::rho
@ rho
density. For some models, rho substitutes p
MaterialPropertyLib::getFluidDensityAndViscosity
std::tuple< double, double > getFluidDensityAndViscosity(double const t, double const dt, ParameterLib::SpatialPosition const &pos, MaterialPropertyLib::Phase const &fluid_phase, MaterialPropertyLib::VariableArray &vars)
It computes fluid density and viscosity for single phase flow model.
Definition GetFluidDensityAndViscosity.cpp:42
MaterialPropertyLib::formEigenTensor
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Definition FormEigenTensor.cpp:115
MaterialPropertyLib::fluidPhase
Phase const & fluidPhase(Medium const &medium)
Returns a gas or aqueous liquid phase of the given medium.
Definition Medium.cpp:100
MathLib::KelvinVector::kelvinVectorToSymmetricTensor
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Definition KelvinVector.cpp:142
MathLib::KelvinVector::KelvinVectorType
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition KelvinVector.h:47
MathLib::KelvinVector::liftVectorToKelvin
Eigen::Matrix< double, DisplacementDim, kelvin_vector_dimensions(DisplacementDim)> liftVectorToKelvin(Eigen::Matrix< double, DisplacementDim, 1 > const &v)
MathLib::createZeroedVector
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Definition EigenMapTools.h:153
MathLib::createZeroedMatrix
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition EigenMapTools.h:32
NumLib::interpolateXCoordinate
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
Definition InitShapeMatrices.h:91
ProcessLib::Graph::eval
auto eval(Function &f, Tuples &... ts) -> typename detail::GetFunctionReturnType< decltype(&Function::eval)>::type
Definition Apply.h:274
ProcessLib::LinearBMatrix::computeBMatrix
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:44
EigenFixedShapeMatrixPolicy::NodalMatrixType
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
Definition ShapeMatrixPolicy.h:74
MathLib::KelvinVector::Invariants::identity2
static Eigen::Matrix< double, KelvinVectorSize, 1 > const identity2
Kelvin mapping of 2nd order identity tensor.
Definition KelvinVector.h:107
MathLib::KelvinVector::Invariants::trace
static double trace(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Trace of the corresponding tensor.
Definition KelvinVector-impl.h:62

References ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::getFluidDensityAndViscosity(), NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MathLib::KelvinVector::liftVectorToKelvin(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::mechanical_strain, MaterialLib::Solids::selectSolidConstitutiveRelation(), ParameterLib::SpatialPosition::setElementID(), MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, and MaterialPropertyLib::VariableArray::volumetric_strain.

◆ assembleWithJacobianForDeformationEquations()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianForDeformationEquations ( const double t,
double const dt,
Eigen::VectorXd const & local_x,
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 of all coupled processes.
local_b_dataRight hand side vector of an element.
local_Jac_dataElement Jacobian matrix for the Newton-Raphson method.

Definition at line 630 of file HydroMechanicsFEM-impl.h.

634{
635 auto const p = local_x.template segment<pressure_size>(pressure_index);
636 auto const u =
637 local_x.template segment<displacement_size>(displacement_index);
638
639 auto local_Jac = MathLib::createZeroedMatrix<
640 typename ShapeMatricesTypeDisplacement::template MatrixType<
641 displacement_size, displacement_size>>(
642 local_Jac_data, displacement_size, displacement_size);
643
644 auto local_rhs =
645 MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
646 template VectorType<displacement_size>>(
647 local_b_data, displacement_size);
648
649 ParameterLib::SpatialPosition x_position;
650 x_position.setElementID(_element.getID());
651
652 auto const& medium = _process_data.media_map.getMedium(_element.getID());
653 auto const& solid = medium->phase("Solid");
654 auto const& fluid = fluidPhase(*medium);
655 MPL::VariableArray vars;
656 auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
657 : vars.liquid_phase_pressure;
658
659 auto const T_ref =
660 medium->property(MPL::PropertyType::reference_temperature)
661 .template value<double>(vars, x_position, t, dt);
662 vars.temperature = T_ref;
663
664 int const n_integration_points = _integration_method.getNumberOfPoints();
665 for (int ip = 0; ip < n_integration_points; ip++)
666 {
667 auto const& w = _ip_data[ip].integration_weight;
668
669 auto const& N_u = _ip_data[ip].N_u;
670 auto const& dNdx_u = _ip_data[ip].dNdx_u;
671
672 auto const& N_p = _ip_data[ip].N_p;
673
674 auto const x_coord =
675 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
676 ShapeMatricesTypeDisplacement>(
677 _element, N_u);
678 auto const B =
679 LinearBMatrix::computeBMatrix<DisplacementDim,
680 ShapeFunctionDisplacement::NPOINTS,
681 typename BMatricesType::BMatrixType>(
682 dNdx_u, N_u, x_coord, _is_axially_symmetric);
683
684 auto& eps = _ip_data[ip].eps;
685 auto const& sigma_eff = _ip_data[ip].sigma_eff;
686
687 phase_pressure = N_p.dot(p);
688
689 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
690 .template value<double>(vars, x_position, t, dt);
691 auto const rho_sr =
692 solid.property(MPL::PropertyType::density)
693 .template value<double>(vars, x_position, t, dt);
694 auto const porosity =
695 medium->property(MPL::PropertyType::porosity)
696 .template value<double>(vars, x_position, t, dt);
697
698 auto const rho_fr = MaterialPropertyLib::getFluidDensity(
699 t, dt, x_position, fluid, vars);
700
701 auto const& b = _process_data.specific_body_force;
702 auto const& identity2 = MathLib::KelvinVector::Invariants<
703 MathLib::KelvinVector::kelvin_vector_dimensions(
704 DisplacementDim)>::identity2;
705
706 eps.noalias() = B * u;
707 vars.mechanical_strain
708 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
709 eps);
710
711 auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
712 dt, u, T_ref);
713
714 local_Jac.noalias() += B.transpose() * C * B * w;
715
716 double p_at_xi = 0.;
717 NumLib::shapeFunctionInterpolate(p, N_p, p_at_xi);
718
719 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
720 local_rhs.noalias() -=
721 (B.transpose() * (sigma_eff - alpha * identity2 * p_at_xi) -
722 N_u_op(N_u).transpose() * rho * b) *
723 w;
724 }
725}
MaterialPropertyLib::getFluidDensity
double getFluidDensity(double const t, double const dt, ParameterLib::SpatialPosition const &pos, Phase const &fluid_phase, VariableArray &vars)
It computes fluid density for single phase flow model.
Definition GetFluidDensityAndViscosity.cpp:21
NumLib::detail::shapeFunctionInterpolate
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
Definition Interpolation.h:27

References ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::getFluidDensity(), NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::mechanical_strain, ParameterLib::SpatialPosition::setElementID(), NumLib::detail::shapeFunctionInterpolate(), and MaterialPropertyLib::VariableArray::temperature.

◆ assembleWithJacobianForPressureEquations()

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

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

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

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

Parameters
tTime
dtTime increment
local_xNodal values of \(x\) of an element of all coupled processes.
local_x_prevNodal values of \(x^{t-1}\) of an element of all coupled processes.
local_b_dataRight hand side vector of an element.
local_Jac_dataElement Jacobian matrix for the Newton-Raphson method.

Definition at line 459 of file HydroMechanicsFEM-impl.h.

464{
465 auto local_rhs =
466 MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
467 template VectorType<pressure_size>>(
468 local_b_data, pressure_size);
469
470 ParameterLib::SpatialPosition pos;
471 pos.setElementID(this->_element.getID());
472
473 auto const p = local_x.template segment<pressure_size>(pressure_index);
474
475 auto const p_prev =
476 local_x_prev.template segment<pressure_size>(pressure_index);
477
478 auto local_Jac = MathLib::createZeroedMatrix<
479 typename ShapeMatricesTypeDisplacement::template MatrixType<
480 pressure_size, pressure_size>>(local_Jac_data, pressure_size,
481 pressure_size);
482
483 typename ShapeMatricesTypePressure::NodalMatrixType storage =
484 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
485 pressure_size);
486
487 typename ShapeMatricesTypePressure::NodalMatrixType laplace =
488 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
489 pressure_size);
490
491 typename ShapeMatricesTypePressure::NodalMatrixType add_p_derivative =
492 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
493 pressure_size);
494
495 auto const& solid_material =
496 MaterialLib::Solids::selectSolidConstitutiveRelation(
497 _process_data.solid_materials, _process_data.material_ids,
498 _element.getID());
499
500 ParameterLib::SpatialPosition x_position;
501 x_position.setElementID(_element.getID());
502
503 auto const& medium = _process_data.media_map.getMedium(_element.getID());
504 auto const& fluid = fluidPhase(*medium);
505 MPL::VariableArray vars;
506 auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
507 : vars.liquid_phase_pressure;
508
509 auto const T_ref =
510 medium->property(MPL::PropertyType::reference_temperature)
511 .template value<double>(vars, x_position, t, dt);
512 vars.temperature = T_ref;
513
514 auto const& identity2 = Invariants::identity2;
515
516 auto const staggered_scheme =
517 std::get<Staggered>(_process_data.coupling_scheme);
518 auto const fixed_stress_stabilization_parameter =
519 staggered_scheme.fixed_stress_stabilization_parameter;
520 auto const fixed_stress_over_time_step =
521 staggered_scheme.fixed_stress_over_time_step;
522
523 int const n_integration_points = _integration_method.getNumberOfPoints();
524 for (int ip = 0; ip < n_integration_points; ip++)
525 {
526 auto const& w = _ip_data[ip].integration_weight;
527
528 auto const& N_p = _ip_data[ip].N_p;
529 auto const& dNdx_p = _ip_data[ip].dNdx_p;
530
531 double const p_int_pt = N_p.dot(p);
532
533 phase_pressure = p_int_pt;
534
535 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
536 t, x_position, dt, T_ref);
537 auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
538
539 auto const alpha_b =
540 medium->property(MPL::PropertyType::biot_coefficient)
541 .template value<double>(vars, x_position, t, dt);
542
543 // Set mechanical variables for the intrinsic permeability model
544 // For stress dependent permeability.
545 auto const sigma_total =
546 (_ip_data[ip].sigma_eff - alpha_b * p_int_pt * identity2).eval();
547 vars.total_stress.emplace<SymmetricTensor>(
548 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_total));
549
550 // For strain dependent permeability
551 vars.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
552 vars.equivalent_plastic_strain =
553 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
554 vars.mechanical_strain
555 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
556 _ip_data[ip].eps);
557
558 auto const K = MPL::formEigenTensor<DisplacementDim>(
559 medium->property(MPL::PropertyType::permeability)
560 .value(vars, x_position, t, dt));
561 auto const porosity =
562 medium->property(MPL::PropertyType::porosity)
563 .template value<double>(vars, x_position, t, dt);
564
565 auto const [rho_fr, mu] =
566 MaterialPropertyLib::getFluidDensityAndViscosity(t, dt, x_position,
567 fluid, vars);
568
569 auto const beta_p =
570 fluid.property(MPL::PropertyType::density)
571 .template dValue<double>(vars, _process_data.phase_variable,
572 x_position, t, dt) /
573 rho_fr;
574
575 auto const K_over_mu = K / mu;
576
577 laplace.noalias() +=
578 rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
579
580 // Artificial compressibility from the fixed stress splitting:
581 auto const beta_FS =
582 fixed_stress_stabilization_parameter * alpha_b * alpha_b / K_S;
583
584 storage.noalias() += rho_fr * N_p.transpose() * N_p * w *
585 ((alpha_b - porosity) * (1.0 - alpha_b) / K_S +
586 porosity * beta_p + beta_FS);
587
588 auto const& b = _process_data.specific_body_force;
589
590 // bodyforce-driven Darcy flow
591 local_rhs.noalias() +=
592 dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
593
594 // density dependence on pressure evaluated for Darcy-term,
595 // for laplace and storage terms this dependence is neglected (as is
596 // done for monolithic too)
597 add_p_derivative.noalias() += rho_fr * beta_p * dNdx_p.transpose() *
598 K_over_mu *
599 (dNdx_p * p - 2.0 * rho_fr * b) * N_p * w;
600
601 if (!fixed_stress_over_time_step)
602 {
603 auto const& eps = _ip_data[ip].eps;
604 auto const& eps_prev = _ip_data[ip].eps_prev;
605 const double eps_v_dot =
606 (Invariants::trace(eps) - Invariants::trace(eps_prev)) / dt;
607
608 // Constant portion of strain rate term:
609 double const strain_rate_b =
610 alpha_b * eps_v_dot -
611 beta_FS * _ip_data[ip].strain_rate_variable;
612
613 local_rhs.noalias() -= strain_rate_b * rho_fr * N_p * w;
614 }
615 else
616 {
617 // Constant portion of strain rate term:
618 local_rhs.noalias() -=
619 alpha_b * _ip_data[ip].strain_rate_variable * rho_fr * N_p * w;
620 }
621 }
622 local_Jac.noalias() = laplace + storage / dt + add_p_derivative;
623
624 local_rhs.noalias() -= laplace * p + storage * (p - p_prev) / dt;
625}

References MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::getFluidDensityAndViscosity(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::mechanical_strain, MaterialLib::Solids::selectSolidConstitutiveRelation(), ParameterLib::SpatialPosition::setElementID(), MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, and MaterialPropertyLib::VariableArray::volumetric_strain.

◆ assembleWithJacobianForStaggeredScheme()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 730 of file HydroMechanicsFEM-impl.h.

737{
738 // For the equations with pressure
739 if (process_id == _process_data.hydraulic_process_id)
740 {
741 assembleWithJacobianForPressureEquations(t, dt, local_x, local_x_prev,
742 local_b_data, local_Jac_data);
743 return;
744 }
745
746 // For the equations with deformation
747 assembleWithJacobianForDeformationEquations(t, dt, local_x, local_b_data,
748 local_Jac_data);
749}
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::assembleWithJacobianForPressureEquations
void assembleWithJacobianForPressureEquations(const double t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
Definition HydroMechanicsFEM-impl.h:460
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::assembleWithJacobianForDeformationEquations
void assembleWithJacobianForDeformationEquations(const double t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
Definition HydroMechanicsFEM-impl.h:631

◆ computeSecondaryVariableConcrete()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 1079 of file HydroMechanicsFEM-impl.h.

1083{
1084 auto const p = local_x.template segment<pressure_size>(pressure_index);
1085
1086 NumLib::interpolateToHigherOrderNodes<
1087 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1088 DisplacementDim>(_element, _is_axially_symmetric, p,
1089 *_process_data.pressure_interpolated);
1090
1091 int const elem_id = _element.getID();
1092 ParameterLib::SpatialPosition x_position;
1093 x_position.setElementID(elem_id);
1094 unsigned const n_integration_points =
1095 _integration_method.getNumberOfPoints();
1096
1097 auto const& medium = _process_data.media_map.getMedium(elem_id);
1098 MPL::VariableArray vars;
1099 auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
1100 : vars.liquid_phase_pressure;
1101
1102 SymmetricTensor k_sum = SymmetricTensor::Zero(KelvinVectorSize);
1103 auto sigma_eff_sum = MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
1104 Eigen::Matrix<double, 3, 3>::Zero());
1105
1106 auto const& identity2 = Invariants::identity2;
1107
1108 for (unsigned ip = 0; ip < n_integration_points; ip++)
1109 {
1110 auto const& eps = _ip_data[ip].eps;
1111 sigma_eff_sum += _ip_data[ip].sigma_eff;
1112
1113 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
1114 .template value<double>(vars, x_position, t, dt);
1115 double const p_int_pt = _ip_data[ip].N_p.dot(p);
1116
1117 phase_pressure = p_int_pt;
1118
1119 // Set mechanical variables for the intrinsic permeability model
1120 // For stress dependent permeability.
1121 {
1122 auto const sigma_total =
1123 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
1124 vars.total_stress.emplace<SymmetricTensor>(
1125 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
1126 sigma_total));
1127 }
1128 // For strain dependent permeability
1129 vars.volumetric_strain = Invariants::trace(eps);
1130 vars.equivalent_plastic_strain =
1131 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
1132 vars.mechanical_strain
1133 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
1134 eps);
1135
1136 k_sum += MPL::getSymmetricTensor<DisplacementDim>(
1137 medium->property(MPL::PropertyType::permeability)
1138 .value(vars, x_position, t, dt));
1139 }
1140
1141 Eigen::Map<Eigen::VectorXd>(
1142 &(*_process_data.permeability)[elem_id * KelvinVectorSize],
1143 KelvinVectorSize) = k_sum / n_integration_points;
1144
1145 Eigen::Matrix<double, 3, 3, 0, 3, 3> const sigma_avg =
1146 MathLib::KelvinVector::kelvinVectorToTensor(sigma_eff_sum) /
1147 n_integration_points;
1148
1149 Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> e_s(sigma_avg);
1150
1151 Eigen::Map<Eigen::Vector3d>(
1152 &(*_process_data.principal_stress_values)[elem_id * 3], 3) =
1153 e_s.eigenvalues();
1154
1155 auto eigen_vectors = e_s.eigenvectors();
1156
1157 for (auto i = 0; i < 3; i++)
1158 {
1159 Eigen::Map<Eigen::Vector3d>(
1160 &(*_process_data.principal_stress_vector[i])[elem_id * 3], 3) =
1161 eigen_vectors.col(i);
1162 }
1163}
MaterialPropertyLib::getSymmetricTensor
SymmetricTensor< GlobalDim > getSymmetricTensor(MaterialPropertyLib::PropertyDataType const &values)
Definition GetSymmetricTensor.cpp:95
MathLib::KelvinVector::kelvinVectorToTensor
Eigen::Matrix< double, 3, 3 > kelvinVectorToTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Definition KelvinVector.cpp:64
MathLib::KelvinVector::tensorToKelvin
KelvinVectorType< DisplacementDim > tensorToKelvin(Eigen::Matrix< double, 3, 3 > const &m)
NumLib::interpolateToHigherOrderNodes
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
Definition Interpolation.h:107

References MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::getSymmetricTensor(), NumLib::interpolateToHigherOrderNodes(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MathLib::KelvinVector::kelvinVectorToTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::mechanical_strain, ParameterLib::SpatialPosition::setElementID(), MathLib::KelvinVector::tensorToKelvin(), MaterialPropertyLib::VariableArray::total_stress, and MaterialPropertyLib::VariableArray::volumetric_strain.

◆ getEpsilon()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getEpsilon ( ) const
overridevirtual

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

Definition at line 1035 of file HydroMechanicsFEM-impl.h.

1036{
1037 auto const kelvin_vector_size =
1038 MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
1039 unsigned const n_integration_points =
1040 _integration_method.getNumberOfPoints();
1041
1042 std::vector<double> ip_epsilon_values;
1043 auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
1044 double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
1045 ip_epsilon_values, n_integration_points, kelvin_vector_size);
1046
1047 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1048 {
1049 auto const& eps = _ip_data[ip].eps;
1050 cache_mat.row(ip) =
1051 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(eps);
1052 }
1053
1054 return ip_epsilon_values;
1055}

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

◆ getIntPtDarcyVelocity()

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

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

Definition at line 370 of file HydroMechanicsFEM-impl.h.

376{
377 int const hydraulic_process_id = _process_data.hydraulic_process_id;
378 auto const indices =
379 NumLib::getIndices(_element.getID(), *dof_table[hydraulic_process_id]);
380 assert(!indices.empty());
381 auto const local_x = x[hydraulic_process_id]->get(indices);
382
383 unsigned const n_integration_points =
384 _integration_method.getNumberOfPoints();
385 cache.clear();
386 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
387 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
388 cache, DisplacementDim, n_integration_points);
389
390 auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
391 pressure_size> const>(local_x.data() + pressure_index, pressure_size);
392
393 ParameterLib::SpatialPosition x_position;
394 x_position.setElementID(_element.getID());
395
396 auto const& medium = _process_data.media_map.getMedium(_element.getID());
397 auto const& fluid = fluidPhase(*medium);
398 MPL::VariableArray vars;
399 auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
400 : vars.liquid_phase_pressure;
401
402 // TODO (naumov) Temporary value not used by current material models. Need
403 // extension of secondary variables interface.
404 double const dt = std::numeric_limits<double>::quiet_NaN();
405 vars.temperature =
406 medium->property(MPL::PropertyType::reference_temperature)
407 .template value<double>(vars, x_position, t, dt);
408
409 auto const& identity2 = Invariants::identity2;
410
411 for (unsigned ip = 0; ip < n_integration_points; ip++)
412 {
413 double const p_int_pt = _ip_data[ip].N_p.dot(p);
414
415 phase_pressure = p_int_pt;
416
417 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
418 .template value<double>(vars, x_position, t, dt);
419
420 // Set mechanical variables for the intrinsic permeability model
421 // For stress dependent permeability.
422 auto const sigma_total =
423 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
424 vars.total_stress.emplace<SymmetricTensor>(
425 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_total));
426
427 // For strain dependent permeability
428 vars.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
429 vars.equivalent_plastic_strain =
430 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
431 vars.mechanical_strain
432 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
433 _ip_data[ip].eps);
434
435 auto const K = MPL::formEigenTensor<DisplacementDim>(
436 medium->property(MPL::PropertyType::permeability)
437 .value(vars, x_position, t, dt));
438
439 auto const [rho_fr, mu] =
440 MaterialPropertyLib::getFluidDensityAndViscosity(t, dt, x_position,
441 fluid, vars);
442
443 auto const K_over_mu = K / mu;
444
445 auto const& b = _process_data.specific_body_force;
446
447 // Compute the velocity
448 auto const& dNdx_p = _ip_data[ip].dNdx_p;
449 cache_matrix.col(ip).noalias() =
450 -K_over_mu * dNdx_p * p + K_over_mu * rho_fr * b;
451 }
452
453 return cache;
454}
NumLib::getIndices
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
Definition DOFTableUtil.cpp:112

References MathLib::createZeroedMatrix(), MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::getFluidDensityAndViscosity(), NumLib::getIndices(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::mechanical_strain, ParameterLib::SpatialPosition::setElementID(), MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, and MaterialPropertyLib::VariableArray::volumetric_strain.

◆ getIntPtEpsilon()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > const & ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtEpsilon ( const double ,
std::vector< GlobalVector * > const & ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & ,
std::vector< double > & cache ) const
inlineoverridevirtual

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

Definition at line 314 of file HydroMechanicsFEM.h.

319 {
320 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
321 _ip_data, &IpData::eps, cache);
322 }
ProcessLib::getIntegrationPointKelvinVectorData
std::vector< double > const & getIntegrationPointKelvinVectorData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
Definition SetOrGetIntegrationPointData.h:68

References ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data, ProcessLib::HydroMechanics::IntegrationPointData< BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::eps, and ProcessLib::getIntegrationPointKelvinVectorData().

◆ getIntPtSigma()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > const & ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtSigma ( const double ,
std::vector< GlobalVector * > const & ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & ,
std::vector< double > & cache ) const
inlineoverridevirtual

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

Definition at line 304 of file HydroMechanicsFEM.h.

309 {
310 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
311 _ip_data, &IpData::sigma_eff, cache);
312 }

References ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data, ProcessLib::getIntegrationPointKelvinVectorData(), and ProcessLib::HydroMechanics::IntegrationPointData< BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::sigma_eff.

◆ getMaterialID()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
int ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getMaterialID ( ) const
overrideprivatevirtual

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

Definition at line 1178 of file HydroMechanicsFEM-impl.h.

1179{
1180 return _process_data.material_ids == nullptr
1181 ? 0
1182 : (*_process_data.material_ids)[_element.getID()];
1183}

◆ getMaterialStateVariablesAt()

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

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

Definition at line 1190 of file HydroMechanicsFEM-impl.h.

1192{
1193 return *_ip_data[integration_point].material_state_variables;
1194}

◆ getNumberOfIntegrationPoints()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
unsigned ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getNumberOfIntegrationPoints ( ) const
overrideprivatevirtual

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

Definition at line 1169 of file HydroMechanicsFEM-impl.h.

1170{
1171 return _integration_method.getNumberOfPoints();
1172}

◆ getShapeMatrix()

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

Provides the shape matrix at the given integration point.

Implements NumLib::ExtrapolatableElement.

Definition at line 280 of file HydroMechanicsFEM.h.

282 {
283 auto const& N_u = _secondary_data.N_u[integration_point];
284
285 // assumes N is stored contiguously in memory
286 return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
287 }

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

◆ getSigma()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getSigma ( ) const
overridevirtual

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

Definition at line 1025 of file HydroMechanicsFEM-impl.h.

1026{
1027 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
1028 _ip_data, &IpData::sigma_eff);
1029}

References ProcessLib::getIntegrationPointKelvinVectorData().

◆ getStrainRateVariable()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getStrainRateVariable ( ) const
overridevirtual

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

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

1062{
1063 unsigned const n_integration_points =
1064 _integration_method.getNumberOfPoints();
1065
1066 std::vector<double> ip_strain_rate_variables(n_integration_points);
1067
1068 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1069 {
1070 ip_strain_rate_variables[ip] = _ip_data[ip].strain_rate_variable;
1071 }
1072
1073 return ip_strain_rate_variables;
1074}

◆ initializeConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::initializeConcrete ( )
inlineoverridevirtual

Set initial stress from parameter.

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 227 of file HydroMechanicsFEM.h.

228 {
229 unsigned const n_integration_points =
230 _integration_method.getNumberOfPoints();
231
232 for (unsigned ip = 0; ip < n_integration_points; ip++)
233 {
234 auto& ip_data = _ip_data[ip];
235
236 ParameterLib::SpatialPosition const x_position{
237 std::nullopt, _element.getID(),
238 MathLib::Point3d(
239 NumLib::interpolateCoordinates<ShapeFunctionPressure,
240 ShapeMatricesTypePressure>(
241 _element, ip_data.N_p))};
242
244 if (_process_data.initial_stress.value)
245 {
246 ip_data.sigma_eff =
247 MathLib::KelvinVector::symmetricTensorToKelvinVector<
248 DisplacementDim>((*_process_data.initial_stress.value)(
249 std::numeric_limits<
250 double>::quiet_NaN() /* time independent */,
251 x_position));
252 }
253
254 double const t = 0; // TODO (naumov) pass t from top
255 ip_data.solid_material.initializeInternalStateVariables(
256 t, x_position, *ip_data.material_state_variables);
257
258 ip_data.pushBackState();
259 }
260 }
MathLib::KelvinVector::symmetricTensorToKelvinVector
Eigen::Matrix< double, Eigen::MatrixBase< Derived >::RowsAtCompileTime, 1 > symmetricTensorToKelvinVector(Eigen::MatrixBase< Derived > const &v)
Definition KelvinVector.h:205
NumLib::interpolateCoordinates
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
Definition InitShapeMatrices.h:102

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

◆ postNonLinearSolverConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::postNonLinearSolverConcrete ( Eigen::VectorXd const & local_x,
Eigen::VectorXd const & local_x_prev,
double const t,
double const dt,
int const process_id )
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 812 of file HydroMechanicsFEM-impl.h.

817{
818 // Note: local_x and local_x_prev only contain the solutions of current
819 // process in the staggered scheme. This has to be changed according to the
820 // same two arguments in postTimestepConcrete.
821
822 int const n_integration_points = _integration_method.getNumberOfPoints();
823
824 auto const staggered_scheme_ptr =
825 std::get_if<Staggered>(&_process_data.coupling_scheme);
826
827 if (staggered_scheme_ptr &&
828 process_id == _process_data.hydraulic_process_id)
829 {
830 if (!staggered_scheme_ptr->fixed_stress_over_time_step)
831 {
832 auto const p =
833 local_x.template segment<pressure_size>(pressure_index);
834 auto const p_prev =
835 local_x_prev.template segment<pressure_size>(pressure_index);
836
837 for (int ip = 0; ip < n_integration_points; ip++)
838 {
839 auto& ip_data = _ip_data[ip];
840
841 auto const& N_p = ip_data.N_p;
842
843 ip_data.strain_rate_variable = N_p.dot(p - p_prev) / dt;
844 }
845 }
846 }
847
848 if (!staggered_scheme_ptr ||
849 process_id == _process_data.mechanics_related_process_id)
850 {
851 ParameterLib::SpatialPosition x_position;
852 x_position.setElementID(_element.getID());
853
854 auto const& medium =
855 _process_data.media_map.getMedium(_element.getID());
856
857 auto const T_ref =
858 medium->property(MPL::PropertyType::reference_temperature)
859 .template value<double>(MPL::EmptyVariableArray, x_position, t,
860 dt);
861
862 auto const u =
863 local_x.template segment<displacement_size>(displacement_index);
864
865 MPL::VariableArray vars;
866 vars.temperature = T_ref;
867
868 for (int ip = 0; ip < n_integration_points; ip++)
869 {
870 auto const& N_u = _ip_data[ip].N_u;
871 auto const& dNdx_u = _ip_data[ip].dNdx_u;
872
873 auto const x_coord =
874 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
875 ShapeMatricesTypeDisplacement>(
876 _element, N_u);
877 auto const B = LinearBMatrix::computeBMatrix<
878 DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
879 typename BMatricesType::BMatrixType>(dNdx_u, N_u, x_coord,
880 _is_axially_symmetric);
881
882 auto& eps = _ip_data[ip].eps;
883 eps.noalias() = B * u;
884 vars.mechanical_strain.emplace<
885 MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(eps);
886
887 _ip_data[ip].updateConstitutiveRelation(vars, t, x_position, dt, u,
888 T_ref);
889 }
890 }
891}
MaterialPropertyLib::EmptyVariableArray
static const VariableArray EmptyVariableArray
Definition VariableType.h:208

References ProcessLib::LinearBMatrix::computeBMatrix(), MaterialPropertyLib::EmptyVariableArray, NumLib::interpolateXCoordinate(), MaterialPropertyLib::VariableArray::mechanical_strain, ParameterLib::SpatialPosition::setElementID(), and MaterialPropertyLib::VariableArray::temperature.

◆ postTimestepConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::postTimestepConcrete ( Eigen::VectorXd const & local_x,
Eigen::VectorXd const & local_x_prev,
double const t,
double const dt,
int const process_id )
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 897 of file HydroMechanicsFEM-impl.h.

901{
902 auto const staggered_scheme_ptr =
903 std::get_if<Staggered>(&_process_data.coupling_scheme);
904
905 if (staggered_scheme_ptr &&
906 process_id == _process_data.hydraulic_process_id)
907 {
908 if (staggered_scheme_ptr->fixed_stress_over_time_step)
909 {
910 auto const fixed_stress_stabilization_parameter =
911 staggered_scheme_ptr->fixed_stress_stabilization_parameter;
912
913 auto const p =
914 local_x.template segment<pressure_size>(pressure_index);
915 auto const p_prev =
916 local_x_prev.template segment<pressure_size>(pressure_index);
917
918 ParameterLib::SpatialPosition x_position;
919 x_position.setElementID(_element.getID());
920
921 auto const& solid_material =
922 MaterialLib::Solids::selectSolidConstitutiveRelation(
923 _process_data.solid_materials, _process_data.material_ids,
924 _element.getID());
925
926 auto const& medium =
927 _process_data.media_map.getMedium(_element.getID());
928 MPL::VariableArray vars;
929
930 auto const T_ref =
931 medium->property(MPL::PropertyType::reference_temperature)
932 .template value<double>(vars, x_position, t, dt);
933 vars.temperature = T_ref;
934
935 int const n_integration_points =
936 _integration_method.getNumberOfPoints();
937 for (int ip = 0; ip < n_integration_points; ip++)
938 {
939 auto& ip_data = _ip_data[ip];
940
941 auto const& N_p = ip_data.N_p;
942
943 auto const& eps = ip_data.eps;
944 auto const& eps_prev = ip_data.eps_prev;
945 const double eps_v_dot =
946 (Invariants::trace(eps) - Invariants::trace(eps_prev)) / dt;
947
948 auto const C_el = ip_data.computeElasticTangentStiffness(
949 t, x_position, dt, T_ref);
950 auto const K_S =
951 solid_material.getBulkModulus(t, x_position, &C_el);
952
953 auto const alpha_b =
954 medium->property(MPL::PropertyType::biot_coefficient)
955 .template value<double>(vars, x_position, t, dt);
956
957 ip_data.strain_rate_variable =
958 eps_v_dot - fixed_stress_stabilization_parameter * alpha_b *
959 N_p.dot(p - p_prev) / dt / K_S;
960 }
961 }
962 }
963
964 unsigned const n_integration_points =
965 _integration_method.getNumberOfPoints();
966
967 for (unsigned ip = 0; ip < n_integration_points; ip++)
968 {
969 _ip_data[ip].pushBackState();
970 }
971}

References MaterialLib::Solids::selectSolidConstitutiveRelation(), ParameterLib::SpatialPosition::setElementID(), and MaterialPropertyLib::VariableArray::temperature.

◆ setInitialConditionsConcrete()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 754 of file HydroMechanicsFEM-impl.h.

758{
759 ParameterLib::SpatialPosition x_position;
760 x_position.setElementID(_element.getID());
761
762 auto const& medium = _process_data.media_map.getMedium(_element.getID());
763
764 auto const p = local_x.template segment<pressure_size>(pressure_index);
765 auto const u =
766 local_x.template segment<displacement_size>(displacement_index);
767
768 auto const& identity2 = Invariants::identity2;
769 const double dt = 0.0;
770
771 MPL::VariableArray vars;
772
773 int const n_integration_points = _integration_method.getNumberOfPoints();
774 for (int ip = 0; ip < n_integration_points; ip++)
775 {
776 auto const& N_u = _ip_data[ip].N_u;
777 auto const& dNdx_u = _ip_data[ip].dNdx_u;
778
779 auto const x_coord =
780 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
781 ShapeMatricesTypeDisplacement>(
782 _element, N_u);
783 auto const B =
784 LinearBMatrix::computeBMatrix<DisplacementDim,
785 ShapeFunctionDisplacement::NPOINTS,
786 typename BMatricesType::BMatrixType>(
787 dNdx_u, N_u, x_coord, _is_axially_symmetric);
788
789 auto& eps = _ip_data[ip].eps;
790 eps.noalias() = B * u;
791 vars.mechanical_strain
792 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
793 eps);
794
795 if (_process_data.initial_stress.isTotalStress())
796 {
797 auto const& N_p = _ip_data[ip].N_p;
798 auto const alpha_b =
799 medium->property(MPL::PropertyType::biot_coefficient)
800 .template value<double>(vars, x_position, t, dt);
801
802 auto& sigma_eff = _ip_data[ip].sigma_eff;
803 sigma_eff.noalias() += alpha_b * N_p.dot(p) * identity2;
804 _ip_data[ip].sigma_eff_prev.noalias() = sigma_eff;
805 }
806 }
807}

References ProcessLib::LinearBMatrix::computeBMatrix(), NumLib::interpolateXCoordinate(), MaterialPropertyLib::VariableArray::mechanical_strain, and ParameterLib::SpatialPosition::setElementID().

◆ setIPDataInitialConditions()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::size_t ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::setIPDataInitialConditions ( std::string_view const name,
double const * values,
int const integration_order )
overridevirtual

Returns number of read integration points.

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

Definition at line 977 of file HydroMechanicsFEM-impl.h.

980{
981 if (integration_order !=
982 static_cast<int>(_integration_method.getIntegrationOrder()))
983 {
984 OGS_FATAL(
985 "Setting integration point initial conditions; The integration "
986 "order of the local assembler for element {:d} is different from "
987 "the integration order in the initial condition.",
988 _element.getID());
989 }
990
991 if (name == "sigma")
992 {
993 if (_process_data.initial_stress.value != nullptr)
994 {
995 OGS_FATAL(
996 "Setting initial conditions for stress from integration "
997 "point data and from a parameter '{:s}' is not possible "
998 "simultaneously.",
999 _process_data.initial_stress.value->name);
1000 }
1001
1002 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
1003 values, _ip_data, &IpData::sigma_eff);
1004 }
1005
1006 if (name == "epsilon")
1007 {
1008 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
1009 values, _ip_data, &IpData::eps);
1010 }
1011
1012 if (name == "strain_rate_variable")
1013 {
1014 return ProcessLib::setIntegrationPointScalarData(
1015 values, _ip_data, &IpData::strain_rate_variable);
1016 }
1017
1018 return 0;
1019}
ProcessLib::setIntegrationPointKelvinVectorData
std::size_t setIntegrationPointKelvinVectorData(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType IpData::*const member)
Definition SetOrGetIntegrationPointData.h:133
ProcessLib::setIntegrationPointScalarData
std::size_t setIntegrationPointScalarData(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType IpData::*const member)
Definition SetOrGetIntegrationPointData.h:215

References OGS_FATAL, ProcessLib::setIntegrationPointKelvinVectorData(), and ProcessLib::setIntegrationPointScalarData().

Member Data Documentation

◆ _element

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
MeshLib::Element const& ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_element
private

Definition at line 395 of file HydroMechanicsFEM.h.

Referenced by ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::initializeConcrete().

◆ _integration_method

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
NumLib::GenericIntegrationMethod const& ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_integration_method
private

Definition at line 394 of file HydroMechanicsFEM.h.

Referenced by ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::HydroMechanicsLocalAssembler(), and ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::initializeConcrete().

◆ _ip_data

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector<IpData, Eigen::aligned_allocator<IpData> > ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data
private

Definition at line 392 of file HydroMechanicsFEM.h.

Referenced by ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::HydroMechanicsLocalAssembler(), ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtEpsilon(), ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtSigma(), and ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::initializeConcrete().

◆ _is_axially_symmetric

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
bool const ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_is_axially_symmetric
private

Definition at line 396 of file HydroMechanicsFEM.h.

◆ _process_data

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
HydroMechanicsProcessData<DisplacementDim>& ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_process_data
private

Definition at line 384 of file HydroMechanicsFEM.h.

Referenced by ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::HydroMechanicsLocalAssembler(), and ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::initializeConcrete().

◆ _secondary_data

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

Definition at line 399 of file HydroMechanicsFEM.h.

Referenced by ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::HydroMechanicsLocalAssembler(), and ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getShapeMatrix().

◆ displacement_index

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

Definition at line 403 of file HydroMechanicsFEM.h.

◆ displacement_size

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

Definition at line 404 of file HydroMechanicsFEM.h.

◆ KelvinVectorSize

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

Definition at line 177 of file HydroMechanicsFEM.h.

◆ N_u_op

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

Definition at line 183 of file HydroMechanicsFEM.h.

◆ pressure_index

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

Definition at line 401 of file HydroMechanicsFEM.h.

◆ pressure_size

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

Definition at line 402 of file HydroMechanicsFEM.h.

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