Processing math: 100%
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 x_position = {
208 std::nullopt, _element.getID(),
209 MathLib::Point3d(
210 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
211 ShapeMatricesTypeDisplacement>(
212 _element, N_u))};
213
214 auto const x_coord = x_position.getCoordinates().value()[0];
215
216 auto const B =
217 LinearBMatrix::computeBMatrix<DisplacementDim,
218 ShapeFunctionDisplacement::NPOINTS,
219 typename BMatricesType::BMatrixType>(
220 dNdx_u, N_u, x_coord, _is_axially_symmetric);
221
222 auto& eps = _ip_data[ip].eps;
223 eps.noalias() = B * u;
224 auto const& sigma_eff = _ip_data[ip].sigma_eff;
225
226 double const p_int_pt = N_p.dot(p);
227
228 phase_pressure = p_int_pt;
229
230 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
231 t, x_position, dt, T_ref);
232 auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
233
234 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
235 .template value<double>(vars, x_position, t, dt);
236
237 auto const rho_sr =
238 solid.property(MPL::PropertyType::density)
239 .template value<double>(vars, x_position, t, dt);
240 auto const porosity =
241 medium->property(MPL::PropertyType::porosity)
242 .template value<double>(vars, x_position, t, dt);
243
244 auto const [rho_fr, mu] =
245 MaterialPropertyLib::getFluidDensityAndViscosity(t, dt, x_position,
246 fluid, vars);
247
248 auto const beta_p =
249 fluid.property(MPL::PropertyType::density)
250 .template dValue<double>(vars, _process_data.phase_variable,
251 x_position, t, dt) /
252 rho_fr;
253
254 // Set mechanical variables for the intrinsic permeability model
255 // For stress dependent permeability.
256 {
257 auto const sigma_total =
258 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
259
260 vars.total_stress.emplace<SymmetricTensor>(
261 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
262 sigma_total));
263 }
264 // For strain dependent permeability
265 vars.volumetric_strain = Invariants::trace(eps);
266 vars.equivalent_plastic_strain =
267 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
268 vars.mechanical_strain
269 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
270 eps);
271
272 auto const K = MPL::formEigenTensor<DisplacementDim>(
273 medium->property(MPL::PropertyType::permeability)
274 .value(vars, x_position, t, dt));
275 auto const dkde = MPL::formEigenTensor<
276 MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim)>(
277 (*medium)[MPL::PropertyType::permeability].dValue(
278 vars, MPL::Variable::mechanical_strain, x_position, t, dt));
279
280 auto const K_over_mu = K / mu;
281
282 auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
283 dt, u, T_ref);
284
285 //
286 // displacement equation, displacement part
287 //
288 local_Jac
289 .template block<displacement_size, displacement_size>(
290 displacement_index, displacement_index)
291 .noalias() += B.transpose() * C * B * w;
292
293 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
294 local_rhs.template segment<displacement_size>(displacement_index)
295 .noalias() -=
296 (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
297
298 //
299 // displacement equation, pressure part
300 //
301 Kup.noalias() += B.transpose() * alpha * identity2 * N_p * w;
302
303 //
304 // pressure equation, pressure part.
305 //
306 laplace_p.noalias() +=
307 rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
308
309 storage_p.noalias() +=
310 rho_fr * N_p.transpose() * N_p * w *
311 ((alpha - porosity) * (1.0 - alpha) / K_S + porosity * beta_p);
312
313 // density dependence on pressure evaluated for Darcy-term,
314 // for laplace and storage terms this dependence is neglected
315 add_p_derivative.noalias() += rho_fr * beta_p * dNdx_p.transpose() *
316 K_over_mu *
317 (dNdx_p * p - 2.0 * rho_fr * b) * N_p * w;
318
319 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
320 dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
321
322 //
323 // pressure equation, displacement part.
324 //
325 Kpu.noalias() +=
326 rho_fr * alpha * N_p.transpose() * identity2.transpose() * B * w;
327
328 Kpu_k.noalias() +=
329 dNdx_p.transpose() *
330 MathLib::KelvinVector::liftVectorToKelvin<DisplacementDim>(
331 dNdx_p * p - rho_fr * b) *
332 dkde * B * rho_fr / mu * w;
333 }
334 // displacement equation, pressure part
335 local_Jac
336 .template block<displacement_size, pressure_size>(displacement_index,
337 pressure_index)
338 .noalias() = -Kup;
339
340 if (_process_data.mass_lumping)
341 {
342 storage_p = storage_p.colwise().sum().eval().asDiagonal();
343
344 if constexpr (pressure_size == displacement_size)
345 {
346 Kpu = Kpu.colwise().sum().eval().asDiagonal();
347 Kpu_k = Kpu_k.colwise().sum().eval().asDiagonal();
348 }
349 }
350
351 // pressure equation, pressure part.
352 local_Jac
353 .template block<pressure_size, pressure_size>(pressure_index,
354 pressure_index)
355 .noalias() += laplace_p + storage_p / dt + add_p_derivative;
356
357 // pressure equation, displacement part.
358 local_Jac
359 .template block<pressure_size, displacement_size>(pressure_index,
360 displacement_index)
361 .noalias() += Kpu / dt + Kpu_k;
362
363 // pressure equation
364 local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
365 laplace_p * p + storage_p * (p - p_prev) / dt + Kpu * (u - u_prev) / dt;
366
367 // displacement equation
368 local_rhs.template segment<displacement_size>(displacement_index)
369 .noalias() += Kup * p;
370}
MeshLib::Element::getID
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:91
ParameterLib::SpatialPosition::setElementID
void setElementID(std::size_t element_id)
Definition SpatialPosition.h:76
ParameterLib::SpatialPosition::getCoordinates
std::optional< MathLib::Point3d > const getCoordinates() const
Definition SpatialPosition.h:63
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
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::interpolateCoordinates
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
Definition InitShapeMatrices.h:102
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, ParameterLib::SpatialPosition::getCoordinates(), MaterialPropertyLib::getFluidDensityAndViscosity(), NumLib::interpolateCoordinates(), 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 649 of file HydroMechanicsFEM-impl.h.

653{
654 auto const p = local_x.template segment<pressure_size>(pressure_index);
655 auto const u =
656 local_x.template segment<displacement_size>(displacement_index);
657
658 auto local_Jac = MathLib::createZeroedMatrix<
659 typename ShapeMatricesTypeDisplacement::template MatrixType<
660 displacement_size, displacement_size>>(
661 local_Jac_data, displacement_size, displacement_size);
662
663 auto local_rhs =
664 MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
665 template VectorType<displacement_size>>(
666 local_b_data, displacement_size);
667
668 ParameterLib::SpatialPosition x_position;
669 x_position.setElementID(_element.getID());
670
671 auto const& medium = _process_data.media_map.getMedium(_element.getID());
672 auto const& solid = medium->phase("Solid");
673 auto const& fluid = fluidPhase(*medium);
674 MPL::VariableArray vars;
675 auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
676 : vars.liquid_phase_pressure;
677
678 auto const T_ref =
679 medium->property(MPL::PropertyType::reference_temperature)
680 .template value<double>(vars, x_position, t, dt);
681 vars.temperature = T_ref;
682
683 int const n_integration_points = _integration_method.getNumberOfPoints();
684 for (int ip = 0; ip < n_integration_points; ip++)
685 {
686 auto const& w = _ip_data[ip].integration_weight;
687
688 auto const& N_u = _ip_data[ip].N_u;
689 auto const& dNdx_u = _ip_data[ip].dNdx_u;
690
691 auto const& N_p = _ip_data[ip].N_p;
692
693 x_position = {
694 std::nullopt, _element.getID(),
695 MathLib::Point3d(
696 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
697 ShapeMatricesTypeDisplacement>(
698 _element, N_u))};
699
700 auto const x_coord = x_position.getCoordinates().value()[0];
701 auto const B =
702 LinearBMatrix::computeBMatrix<DisplacementDim,
703 ShapeFunctionDisplacement::NPOINTS,
704 typename BMatricesType::BMatrixType>(
705 dNdx_u, N_u, x_coord, _is_axially_symmetric);
706
707 auto& eps = _ip_data[ip].eps;
708 auto const& sigma_eff = _ip_data[ip].sigma_eff;
709
710 phase_pressure = N_p.dot(p);
711
712 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
713 .template value<double>(vars, x_position, t, dt);
714 auto const rho_sr =
715 solid.property(MPL::PropertyType::density)
716 .template value<double>(vars, x_position, t, dt);
717 auto const porosity =
718 medium->property(MPL::PropertyType::porosity)
719 .template value<double>(vars, x_position, t, dt);
720
721 auto const rho_fr = MaterialPropertyLib::getFluidDensity(
722 t, dt, x_position, fluid, vars);
723
724 auto const& b = _process_data.specific_body_force;
725 auto const& identity2 = MathLib::KelvinVector::Invariants<
726 MathLib::KelvinVector::kelvin_vector_dimensions(
727 DisplacementDim)>::identity2;
728
729 eps.noalias() = B * u;
730 vars.mechanical_strain
731 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
732 eps);
733
734 auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
735 dt, u, T_ref);
736
737 local_Jac.noalias() += B.transpose() * C * B * w;
738
739 double p_at_xi = 0.;
740 NumLib::shapeFunctionInterpolate(p, N_p, p_at_xi);
741
742 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
743 local_rhs.noalias() -=
744 (B.transpose() * (sigma_eff - alpha * identity2 * p_at_xi) -
745 N_u_op(N_u).transpose() * rho * b) *
746 w;
747 }
748}
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::interpolateCoordinates(), 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 471 of file HydroMechanicsFEM-impl.h.

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

References MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::getFluidDensityAndViscosity(), NumLib::interpolateCoordinates(), 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 753 of file HydroMechanicsFEM-impl.h.

760{
761 // For the equations with pressure
762 if (process_id == _process_data.hydraulic_process_id)
763 {
764 assembleWithJacobianForPressureEquations(t, dt, local_x, local_x_prev,
765 local_b_data, local_Jac_data);
766 return;
767 }
768
769 // For the equations with deformation
770 assembleWithJacobianForDeformationEquations(t, dt, local_x, local_b_data,
771 local_Jac_data);
772}
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:472
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:650

◆ 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 1116 of file HydroMechanicsFEM-impl.h.

1120{
1121 auto const p = local_x.template segment<pressure_size>(pressure_index);
1122
1123 NumLib::interpolateToHigherOrderNodes<
1124 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1125 DisplacementDim>(_element, _is_axially_symmetric, p,
1126 *_process_data.pressure_interpolated);
1127
1128 int const elem_id = _element.getID();
1129
1130 unsigned const n_integration_points =
1131 _integration_method.getNumberOfPoints();
1132
1133 auto const& medium = _process_data.media_map.getMedium(elem_id);
1134 MPL::VariableArray vars;
1135 auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
1136 : vars.liquid_phase_pressure;
1137
1138 SymmetricTensor k_sum = SymmetricTensor::Zero(KelvinVectorSize);
1139 auto sigma_eff_sum = MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
1140 Eigen::Matrix<double, 3, 3>::Zero());
1141
1142 auto const& identity2 = Invariants::identity2;
1143
1144 for (unsigned ip = 0; ip < n_integration_points; ip++)
1145 {
1146 auto const& eps = _ip_data[ip].eps;
1147 sigma_eff_sum += _ip_data[ip].sigma_eff;
1148
1149 ParameterLib::SpatialPosition x_position = {
1150 std::nullopt, _element.getID(),
1151 MathLib::Point3d(
1152 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
1153 ShapeMatricesTypeDisplacement>(
1154 _element, _ip_data[ip].N_u))};
1155
1156 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
1157 .template value<double>(vars, x_position, t, dt);
1158 double const p_int_pt = _ip_data[ip].N_p.dot(p);
1159
1160 phase_pressure = p_int_pt;
1161
1162 // Set mechanical variables for the intrinsic permeability model
1163 // For stress dependent permeability.
1164 {
1165 auto const sigma_total =
1166 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
1167 vars.total_stress.emplace<SymmetricTensor>(
1168 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
1169 sigma_total));
1170 }
1171 // For strain dependent permeability
1172 vars.volumetric_strain = Invariants::trace(eps);
1173 vars.equivalent_plastic_strain =
1174 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
1175 vars.mechanical_strain
1176 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
1177 eps);
1178
1179 k_sum += MPL::getSymmetricTensor<DisplacementDim>(
1180 medium->property(MPL::PropertyType::permeability)
1181 .value(vars, x_position, t, dt));
1182 }
1183
1184 Eigen::Map<Eigen::VectorXd>(
1185 &(*_process_data.permeability)[elem_id * KelvinVectorSize],
1186 KelvinVectorSize) = k_sum / n_integration_points;
1187
1188 Eigen::Matrix<double, 3, 3, 0, 3, 3> const sigma_avg =
1189 MathLib::KelvinVector::kelvinVectorToTensor(sigma_eff_sum) /
1190 n_integration_points;
1191
1192 Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> e_s(sigma_avg);
1193
1194 Eigen::Map<Eigen::Vector3d>(
1195 &(*_process_data.principal_stress_values)[elem_id * 3], 3) =
1196 e_s.eigenvalues();
1197
1198 auto eigen_vectors = e_s.eigenvectors();
1199
1200 for (auto i = 0; i < 3; i++)
1201 {
1202 Eigen::Map<Eigen::Vector3d>(
1203 &(*_process_data.principal_stress_vector[i])[elem_id * 3], 3) =
1204 eigen_vectors.col(i);
1205 }
1206}
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::interpolateCoordinates(), NumLib::interpolateToHigherOrderNodes(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MathLib::KelvinVector::kelvinVectorToTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::mechanical_strain, 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 1072 of file HydroMechanicsFEM-impl.h.

1073{
1074 auto const kelvin_vector_size =
1075 MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
1076 unsigned const n_integration_points =
1077 _integration_method.getNumberOfPoints();
1078
1079 std::vector<double> ip_epsilon_values;
1080 auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
1081 double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
1082 ip_epsilon_values, n_integration_points, kelvin_vector_size);
1083
1084 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1085 {
1086 auto const& eps = _ip_data[ip].eps;
1087 cache_mat.row(ip) =
1088 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(eps);
1089 }
1090
1091 return ip_epsilon_values;
1092}

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 375 of file HydroMechanicsFEM-impl.h.

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

1222{
1223 return _process_data.material_ids == nullptr
1224 ? 0
1225 : (*_process_data.material_ids)[_element.getID()];
1226}

◆ 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 1233 of file HydroMechanicsFEM-impl.h.

1235{
1236 return *_ip_data[integration_point].material_state_variables;
1237}

◆ 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 1212 of file HydroMechanicsFEM-impl.h.

1213{
1214 return _integration_method.getNumberOfPoints();
1215}

◆ 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 1062 of file HydroMechanicsFEM-impl.h.

1063{
1064 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
1065 _ip_data, &IpData::sigma_eff);
1066}

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 1098 of file HydroMechanicsFEM-impl.h.

1099{
1100 unsigned const n_integration_points =
1101 _integration_method.getNumberOfPoints();
1102
1103 std::vector<double> ip_strain_rate_variables(n_integration_points);
1104
1105 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1106 {
1107 ip_strain_rate_variables[ip] = _ip_data[ip].strain_rate_variable;
1108 }
1109
1110 return ip_strain_rate_variables;
1111}

◆ 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

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 839 of file HydroMechanicsFEM-impl.h.

844{
845 // Note: local_x and local_x_prev only contain the solutions of current
846 // process in the staggered scheme. This has to be changed according to the
847 // same two arguments in postTimestepConcrete.
848
849 int const n_integration_points = _integration_method.getNumberOfPoints();
850
851 auto const staggered_scheme_ptr =
852 std::get_if<Staggered>(&_process_data.coupling_scheme);
853
854 if (staggered_scheme_ptr &&
855 process_id == _process_data.hydraulic_process_id)
856 {
857 if (!staggered_scheme_ptr->fixed_stress_over_time_step)
858 {
859 auto const p =
860 local_x.template segment<pressure_size>(pressure_index);
861 auto const p_prev =
862 local_x_prev.template segment<pressure_size>(pressure_index);
863
864 for (int ip = 0; ip < n_integration_points; ip++)
865 {
866 auto& ip_data = _ip_data[ip];
867
868 auto const& N_p = ip_data.N_p;
869
870 ip_data.strain_rate_variable = N_p.dot(p - p_prev) / dt;
871 }
872 }
873 }
874
875 if (!staggered_scheme_ptr ||
876 process_id == _process_data.mechanics_related_process_id)
877 {
878 ParameterLib::SpatialPosition x_position;
879 x_position.setElementID(_element.getID());
880
881 auto const& medium =
882 _process_data.media_map.getMedium(_element.getID());
883
884 auto const T_ref =
885 medium->property(MPL::PropertyType::reference_temperature)
886 .template value<double>(MPL::EmptyVariableArray, x_position, t,
887 dt);
888
889 auto const u =
890 local_x.template segment<displacement_size>(displacement_index);
891
892 MPL::VariableArray vars;
893 vars.temperature = T_ref;
894
895 for (int ip = 0; ip < n_integration_points; ip++)
896 {
897 auto const& N_u = _ip_data[ip].N_u;
898 auto const& dNdx_u = _ip_data[ip].dNdx_u;
899
900 x_position = {std::nullopt, _element.getID(),
901 MathLib::Point3d(NumLib::interpolateCoordinates<
902 ShapeFunctionDisplacement,
903 ShapeMatricesTypeDisplacement>(
904 _element, N_u))};
905
906 auto const x_coord = x_position.getCoordinates().value()[0];
907 auto const B = LinearBMatrix::computeBMatrix<
908 DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
909 typename BMatricesType::BMatrixType>(dNdx_u, N_u, x_coord,
910 _is_axially_symmetric);
911
912 auto& eps = _ip_data[ip].eps;
913 eps.noalias() = B * u;
914 vars.mechanical_strain.emplace<
915 MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(eps);
916
917 _ip_data[ip].updateConstitutiveRelation(vars, t, x_position, dt, u,
918 T_ref);
919 }
920 }
921}
MaterialPropertyLib::EmptyVariableArray
static const VariableArray EmptyVariableArray
Definition VariableType.h:208

References ProcessLib::LinearBMatrix::computeBMatrix(), MaterialPropertyLib::EmptyVariableArray, ParameterLib::SpatialPosition::getCoordinates(), NumLib::interpolateCoordinates(), 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 927 of file HydroMechanicsFEM-impl.h.

931{
932 auto const staggered_scheme_ptr =
933 std::get_if<Staggered>(&_process_data.coupling_scheme);
934
935 if (staggered_scheme_ptr &&
936 process_id == _process_data.hydraulic_process_id)
937 {
938 if (staggered_scheme_ptr->fixed_stress_over_time_step)
939 {
940 auto const fixed_stress_stabilization_parameter =
941 staggered_scheme_ptr->fixed_stress_stabilization_parameter;
942
943 auto const p =
944 local_x.template segment<pressure_size>(pressure_index);
945 auto const p_prev =
946 local_x_prev.template segment<pressure_size>(pressure_index);
947
948 ParameterLib::SpatialPosition x_position;
949 x_position.setElementID(_element.getID());
950
951 auto const& solid_material =
952 MaterialLib::Solids::selectSolidConstitutiveRelation(
953 _process_data.solid_materials, _process_data.material_ids,
954 _element.getID());
955
956 auto const& medium =
957 _process_data.media_map.getMedium(_element.getID());
958 MPL::VariableArray vars;
959
960 auto const T_ref =
961 medium->property(MPL::PropertyType::reference_temperature)
962 .template value<double>(vars, x_position, t, dt);
963 vars.temperature = T_ref;
964
965 int const n_integration_points =
966 _integration_method.getNumberOfPoints();
967 for (int ip = 0; ip < n_integration_points; ip++)
968 {
969 auto& ip_data = _ip_data[ip];
970
971 auto const& N_p = ip_data.N_p;
972
973 x_position = {
974 std::nullopt, _element.getID(),
975 MathLib::Point3d(
976 NumLib::interpolateCoordinates<
977 ShapeFunctionPressure, ShapeMatricesTypePressure>(
978 _element, N_p))};
979
980 auto const& eps = ip_data.eps;
981 auto const& eps_prev = ip_data.eps_prev;
982 const double eps_v_dot =
983 (Invariants::trace(eps) - Invariants::trace(eps_prev)) / dt;
984
985 auto const C_el = ip_data.computeElasticTangentStiffness(
986 t, x_position, dt, T_ref);
987 auto const K_S =
988 solid_material.getBulkModulus(t, x_position, &C_el);
989
990 auto const alpha_b =
991 medium->property(MPL::PropertyType::biot_coefficient)
992 .template value<double>(vars, x_position, t, dt);
993
994 ip_data.strain_rate_variable =
995 eps_v_dot - fixed_stress_stabilization_parameter * alpha_b *
996 N_p.dot(p - p_prev) / dt / K_S;
997 }
998 }
999 }
1000
1001 unsigned const n_integration_points =
1002 _integration_method.getNumberOfPoints();
1003
1004 for (unsigned ip = 0; ip < n_integration_points; ip++)
1005 {
1006 _ip_data[ip].pushBackState();
1007 }
1008}

References NumLib::interpolateCoordinates(), 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 777 of file HydroMechanicsFEM-impl.h.

781{
782 ParameterLib::SpatialPosition x_position;
783 x_position.setElementID(_element.getID());
784
785 auto const& medium = _process_data.media_map.getMedium(_element.getID());
786
787 auto const p = local_x.template segment<pressure_size>(pressure_index);
788 auto const u =
789 local_x.template segment<displacement_size>(displacement_index);
790
791 auto const& identity2 = Invariants::identity2;
792 const double dt = 0.0;
793
794 MPL::VariableArray vars;
795
796 int const n_integration_points = _integration_method.getNumberOfPoints();
797 for (int ip = 0; ip < n_integration_points; ip++)
798 {
799 auto const& N_u = _ip_data[ip].N_u;
800 auto const& dNdx_u = _ip_data[ip].dNdx_u;
801
802 x_position = {
803 std::nullopt, _element.getID(),
804 MathLib::Point3d(
805 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
806 ShapeMatricesTypeDisplacement>(
807 _element, N_u))};
808
809 auto const x_coord = x_position.getCoordinates().value()[0];
810 auto const B =
811 LinearBMatrix::computeBMatrix<DisplacementDim,
812 ShapeFunctionDisplacement::NPOINTS,
813 typename BMatricesType::BMatrixType>(
814 dNdx_u, N_u, x_coord, _is_axially_symmetric);
815
816 auto& eps = _ip_data[ip].eps;
817 eps.noalias() = B * u;
818 vars.mechanical_strain
819 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
820 eps);
821
822 if (_process_data.initial_stress.isTotalStress())
823 {
824 auto const& N_p = _ip_data[ip].N_p;
825 auto const alpha_b =
826 medium->property(MPL::PropertyType::biot_coefficient)
827 .template value<double>(vars, x_position, t, dt);
828
829 auto& sigma_eff = _ip_data[ip].sigma_eff;
830 sigma_eff.noalias() += alpha_b * N_p.dot(p) * identity2;
831 _ip_data[ip].sigma_eff_prev.noalias() = sigma_eff;
832 }
833 }
834}

References ProcessLib::LinearBMatrix::computeBMatrix(), ParameterLib::SpatialPosition::getCoordinates(), NumLib::interpolateCoordinates(), 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 1014 of file HydroMechanicsFEM-impl.h.

1017{
1018 if (integration_order !=
1019 static_cast<int>(_integration_method.getIntegrationOrder()))
1020 {
1021 OGS_FATAL(
1022 "Setting integration point initial conditions; The integration "
1023 "order of the local assembler for element {:d} is different from "
1024 "the integration order in the initial condition.",
1025 _element.getID());
1026 }
1027
1028 if (name == "sigma")
1029 {
1030 if (_process_data.initial_stress.value != nullptr)
1031 {
1032 OGS_FATAL(
1033 "Setting initial conditions for stress from integration "
1034 "point data and from a parameter '{:s}' is not possible "
1035 "simultaneously.",
1036 _process_data.initial_stress.value->name);
1037 }
1038
1039 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
1040 values, _ip_data, &IpData::sigma_eff);
1041 }
1042
1043 if (name == "epsilon")
1044 {
1045 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
1046 values, _ip_data, &IpData::eps);
1047 }
1048
1049 if (name == "strain_rate_variable")
1050 {
1051 return ProcessLib::setIntegrationPointScalarData(
1052 values, _ip_data, &IpData::strain_rate_variable);
1053 }
1054
1055 return 0;
1056}
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: