OGS
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.
 
- 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:
ShapeMatricesTypePressure, DisplacementDim,
ShapeFunctionDisplacement::NPOINTS>
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
BMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > BMatricesType

Definition at line 388 of file HydroMechanicsFEM.h.

◆ ShapeMatricesTypeDisplacement

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

◆ ShapeMatricesTypePressure

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

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

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

References ProcessLib::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 }
#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 102 of file HydroMechanicsFEM-impl.h.

108{
109 assert(local_x.size() == pressure_size + displacement_size);
110
111 auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
112 pressure_size> const>(local_x.data() + pressure_index, pressure_size);
113
114 auto u =
115 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
116 displacement_size> const>(local_x.data() + displacement_index,
118
119 auto p_prev =
120 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
121 pressure_size> const>(local_x_prev.data() + pressure_index,
123 auto u_prev =
124 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
125 displacement_size> const>(local_x_prev.data() + displacement_index,
127
128 auto local_Jac = MathLib::createZeroedMatrix<
129 typename ShapeMatricesTypeDisplacement::template MatrixType<
132 local_Jac_data, displacement_size + pressure_size,
134
135 auto local_rhs = MathLib::createZeroedVector<
136 typename ShapeMatricesTypeDisplacement::template VectorType<
138 local_rhs_data, displacement_size + pressure_size);
139
141 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
143
145 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
147
148 typename ShapeMatricesTypePressure::NodalMatrixType add_p_derivative =
149 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
151
152 typename ShapeMatricesTypeDisplacement::template MatrixType<
154 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
157
158 typename ShapeMatricesTypeDisplacement::template MatrixType<
160 Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
163
164 typename ShapeMatricesTypeDisplacement::template MatrixType<
166 Kpu_k = ShapeMatricesTypeDisplacement::template MatrixType<
169
170 auto const& solid_material =
172 _process_data.solid_materials, _process_data.material_ids,
173 _element.getID());
174
176 x_position.setElementID(_element.getID());
177
178 unsigned const n_integration_points =
180
181 auto const& b = _process_data.specific_body_force;
182 auto const& medium = _process_data.media_map.getMedium(_element.getID());
183 auto const& solid = medium->phase("Solid");
184 auto const& fluid = fluidPhase(*medium);
186 auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
188
189 auto const T_ref =
190 medium->property(MPL::PropertyType::reference_temperature)
191 .template value<double>(vars, x_position, t, dt);
192 vars.temperature = T_ref;
193
194 auto const& identity2 = Invariants::identity2;
195
196 for (unsigned ip = 0; ip < n_integration_points; ip++)
197 {
198 x_position.setIntegrationPoint(ip);
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,
210 _element, N_u);
211 auto const B =
212 LinearBMatrix::computeBMatrix<DisplacementDim,
213 ShapeFunctionDisplacement::NPOINTS,
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 // Quick workaround: If fluid density is described as ideal gas, then
240 // the molar mass must be passed to the MPL::IdealGasLaw via the
241 // variable_array and the fluid must have the property
242 // MPL::PropertyType::molar_mass. For other density models (e.g.
243 // Constant), it is not mandatory to specify the molar mass.
244 if (fluid.hasProperty(MPL::PropertyType::molar_mass))
245 {
246 vars.molar_mass =
247 fluid.property(MPL::PropertyType::molar_mass)
248 .template value<double>(vars, x_position, t, dt);
249 }
250 auto const rho_fr =
251 fluid.property(MPL::PropertyType::density)
252 .template value<double>(vars, x_position, t, dt);
253 vars.density = rho_fr;
254
255 auto const mu = fluid.property(MPL::PropertyType::viscosity)
256 .template value<double>(vars, x_position, t, dt);
257
258 auto const beta_p =
259 fluid.property(MPL::PropertyType::density)
260 .template dValue<double>(vars, _process_data.phase_variable,
261 x_position, t, dt) /
262 rho_fr;
263
264 // Set mechanical variables for the intrinsic permeability model
265 // For stress dependent permeability.
266 {
267 auto const sigma_total =
268 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
269
270 vars.total_stress.emplace<SymmetricTensor>(
272 sigma_total));
273 }
274 // For strain dependent permeability
277 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
280 eps);
281
283 medium->property(MPL::PropertyType::permeability)
284 .value(vars, x_position, t, dt));
285 auto const dkde = MPL::formEigenTensor<
287 (*medium)[MPL::PropertyType::permeability].dValue(
288 vars, MPL::Variable::mechanical_strain, x_position, t, dt));
289
290 auto const K_over_mu = K / mu;
291
292 auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
293 dt, u, T_ref);
294
295 //
296 // displacement equation, displacement part
297 //
298 local_Jac
299 .template block<displacement_size, displacement_size>(
301 .noalias() += B.transpose() * C * B * w;
302
303 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
304 local_rhs.template segment<displacement_size>(displacement_index)
305 .noalias() -=
306 (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
307
308 //
309 // displacement equation, pressure part
310 //
311 Kup.noalias() += B.transpose() * alpha * identity2 * N_p * w;
312
313 //
314 // pressure equation, pressure part.
315 //
316 laplace_p.noalias() +=
317 rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
318
319 storage_p.noalias() +=
320 rho_fr * N_p.transpose() * N_p * w *
321 ((alpha - porosity) * (1.0 - alpha) / K_S + porosity * beta_p);
322
323 // density dependence on pressure evaluated for Darcy-term,
324 // for laplace and storage terms this dependence is neglected
325 add_p_derivative.noalias() += rho_fr * beta_p * dNdx_p.transpose() *
326 K_over_mu *
327 (dNdx_p * p - 2.0 * rho_fr * b) * N_p * w;
328
329 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
330 dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
331
332 //
333 // pressure equation, displacement part.
334 //
335 Kpu.noalias() +=
336 rho_fr * alpha * N_p.transpose() * identity2.transpose() * B * w;
337
338 Kpu_k.noalias() +=
339 dNdx_p.transpose() *
341 dNdx_p * p - rho_fr * b) *
342 dkde * B * rho_fr / mu * w;
343 }
344 // displacement equation, pressure part
345 local_Jac
346 .template block<displacement_size, pressure_size>(displacement_index,
348 .noalias() = -Kup;
349
350 if (_process_data.mass_lumping)
351 {
352 storage_p = storage_p.colwise().sum().eval().asDiagonal();
353
354 if constexpr (pressure_size == displacement_size)
355 {
356 Kpu = Kpu.colwise().sum().eval().asDiagonal();
357 Kpu_k = Kpu_k.colwise().sum().eval().asDiagonal();
358 }
359 }
360
361 // pressure equation, pressure part.
362 local_Jac
363 .template block<pressure_size, pressure_size>(pressure_index,
365 .noalias() += laplace_p + storage_p / dt + add_p_derivative;
366
367 // pressure equation, displacement part.
368 local_Jac
369 .template block<pressure_size, displacement_size>(pressure_index,
371 .noalias() += Kpu / dt + Kpu_k;
372
373 // pressure equation
374 local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
375 laplace_p * p + storage_p * (p - p_prev) / dt + Kpu * (u - u_prev) / dt;
376
377 // displacement equation
378 local_rhs.template segment<displacement_size>(displacement_index)
379 .noalias() += Kup * p;
380}
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:89
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Eigen::Matrix< double, KelvinVectorSize, 1 > SymmetricTensor
@ rho
density. For some models, rho substitutes p
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Phase const & fluidPhase(Medium const &medium)
Returns a gas or aqueous liquid phase of the given medium.
Definition Medium.cpp:100
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Eigen::Matrix< double, DisplacementDim, kelvin_vector_dimensions(DisplacementDim)> liftVectorToKelvin(Eigen::Matrix< double, DisplacementDim, 1 > const &v)
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
auto eval(Function &f, Tuples &... ts) -> typename detail::GetFunctionReturnType< decltype(&Function::eval)>::type
Definition Apply.h:274
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.
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
static Eigen::Matrix< double, KelvinVectorSize, 1 > const identity2
Kelvin mapping of 2nd order identity tensor.
static double trace(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Trace of the corresponding tensor.

References ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::gas_phase_pressure, NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MathLib::KelvinVector::liftVectorToKelvin(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::mechanical_strain, MaterialPropertyLib::VariableArray::molar_mass, MaterialLib::Solids::selectSolidConstitutiveRelation(), ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), 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 678 of file HydroMechanicsFEM-impl.h.

682{
683 auto const p = local_x.template segment<pressure_size>(pressure_index);
684 auto const u =
685 local_x.template segment<displacement_size>(displacement_index);
686
687 auto local_Jac = MathLib::createZeroedMatrix<
688 typename ShapeMatricesTypeDisplacement::template MatrixType<
690 local_Jac_data, displacement_size, displacement_size);
691
692 auto local_rhs =
693 MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
694 template VectorType<displacement_size>>(
695 local_b_data, displacement_size);
696
698 x_position.setElementID(_element.getID());
699
700 auto const& medium = _process_data.media_map.getMedium(_element.getID());
701 auto const& solid = medium->phase("Solid");
702 auto const& fluid = fluidPhase(*medium);
704 auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
706
707 auto const T_ref =
708 medium->property(MPL::PropertyType::reference_temperature)
709 .template value<double>(vars, x_position, t, dt);
710 vars.temperature = T_ref;
711
712 int const n_integration_points = _integration_method.getNumberOfPoints();
713 for (int ip = 0; ip < n_integration_points; ip++)
714 {
715 x_position.setIntegrationPoint(ip);
716 auto const& w = _ip_data[ip].integration_weight;
717
718 auto const& N_u = _ip_data[ip].N_u;
719 auto const& dNdx_u = _ip_data[ip].dNdx_u;
720
721 auto const& N_p = _ip_data[ip].N_p;
722
723 auto const x_coord =
724 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
726 _element, N_u);
727 auto const B =
728 LinearBMatrix::computeBMatrix<DisplacementDim,
729 ShapeFunctionDisplacement::NPOINTS,
731 dNdx_u, N_u, x_coord, _is_axially_symmetric);
732
733 auto& eps = _ip_data[ip].eps;
734 auto const& sigma_eff = _ip_data[ip].sigma_eff;
735
736 phase_pressure = N_p.dot(p);
737
738 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
739 .template value<double>(vars, x_position, t, dt);
740 auto const rho_sr =
741 solid.property(MPL::PropertyType::density)
742 .template value<double>(vars, x_position, t, dt);
743 auto const porosity =
744 medium->property(MPL::PropertyType::porosity)
745 .template value<double>(vars, x_position, t, dt);
746
747 // Quick workaround: If fluid density is described as ideal gas, then
748 // the molar mass must be passed to the MPL::IdealGasLaw via the
749 // variable_array and the fluid must have the property
750 // MPL::PropertyType::molar_mass. For other density models (e.g.
751 // Constant), it is not mandatory to specify the molar mass.
752 if (fluid.hasProperty(MPL::PropertyType::molar_mass))
753 {
754 vars.molar_mass =
755 fluid.property(MPL::PropertyType::molar_mass)
756 .template value<double>(vars, x_position, t, dt);
757 }
758 auto const rho_fr =
759 fluid.property(MPL::PropertyType::density)
760 .template value<double>(vars, x_position, t, dt);
761
762 auto const& b = _process_data.specific_body_force;
763 auto const& identity2 = MathLib::KelvinVector::Invariants<
765 DisplacementDim)>::identity2;
766
767 eps.noalias() = B * u;
770 eps);
771
772 auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
773 dt, u, T_ref);
774
775 local_Jac.noalias() += B.transpose() * C * B * w;
776
777 double p_at_xi = 0.;
778 NumLib::shapeFunctionInterpolate(p, N_p, p_at_xi);
779
780 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
781 local_rhs.noalias() -=
782 (B.transpose() * (sigma_eff - alpha * identity2 * p_at_xi) -
783 N_u_op(N_u).transpose() * rho * b) *
784 w;
785 }
786}
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)

References ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::gas_phase_pressure, NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::mechanical_strain, MaterialPropertyLib::VariableArray::molar_mass, 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 492 of file HydroMechanicsFEM-impl.h.

497{
498 auto local_rhs =
499 MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
500 template VectorType<pressure_size>>(
501 local_b_data, pressure_size);
502
504 pos.setElementID(this->_element.getID());
505
506 auto const p = local_x.template segment<pressure_size>(pressure_index);
507
508 auto const p_prev =
509 local_x_prev.template segment<pressure_size>(pressure_index);
510
511 auto local_Jac = MathLib::createZeroedMatrix<
512 typename ShapeMatricesTypeDisplacement::template MatrixType<
513 pressure_size, pressure_size>>(local_Jac_data, pressure_size,
515
517 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
519
521 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
523
524 typename ShapeMatricesTypePressure::NodalMatrixType add_p_derivative =
525 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
527
528 auto const& solid_material =
530 _process_data.solid_materials, _process_data.material_ids,
531 _element.getID());
532
534 x_position.setElementID(_element.getID());
535
536 auto const& medium = _process_data.media_map.getMedium(_element.getID());
537 auto const& fluid = fluidPhase(*medium);
539 auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
541
542 auto const T_ref =
543 medium->property(MPL::PropertyType::reference_temperature)
544 .template value<double>(vars, x_position, t, dt);
545 vars.temperature = T_ref;
546
547 auto const& identity2 = Invariants::identity2;
548
549 auto const staggered_scheme =
550 std::get<Staggered>(_process_data.coupling_scheme);
551 auto const fixed_stress_stabilization_parameter =
552 staggered_scheme.fixed_stress_stabilization_parameter;
553 auto const fixed_stress_over_time_step =
554 staggered_scheme.fixed_stress_over_time_step;
555
556 int const n_integration_points = _integration_method.getNumberOfPoints();
557 for (int ip = 0; ip < n_integration_points; ip++)
558 {
559 x_position.setIntegrationPoint(ip);
560 auto const& w = _ip_data[ip].integration_weight;
561
562 auto const& N_p = _ip_data[ip].N_p;
563 auto const& dNdx_p = _ip_data[ip].dNdx_p;
564
565 double const p_int_pt = N_p.dot(p);
566
567 phase_pressure = p_int_pt;
568
569 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
570 t, x_position, dt, T_ref);
571 auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
572
573 auto const alpha_b =
574 medium->property(MPL::PropertyType::biot_coefficient)
575 .template value<double>(vars, x_position, t, dt);
576
577 // Set mechanical variables for the intrinsic permeability model
578 // For stress dependent permeability.
579 auto const sigma_total =
580 (_ip_data[ip].sigma_eff - alpha_b * p_int_pt * identity2).eval();
581 vars.total_stress.emplace<SymmetricTensor>(
583
584 // For strain dependent permeability
587 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
590 _ip_data[ip].eps);
591
593 medium->property(MPL::PropertyType::permeability)
594 .value(vars, x_position, t, dt));
595 auto const porosity =
596 medium->property(MPL::PropertyType::porosity)
597 .template value<double>(vars, x_position, t, dt);
598
599 // Quick workaround: If fluid density is described as ideal gas, then
600 // the molar mass must be passed to the MPL::IdealGasLaw via the
601 // variable_array and the fluid must have the property
602 // MPL::PropertyType::molar_mass. For other density models (e.g.
603 // Constant), it is not mandatory to specify the molar mass.
604 if (fluid.hasProperty(MPL::PropertyType::molar_mass))
605 {
606 vars.molar_mass =
607 fluid.property(MPL::PropertyType::molar_mass)
608 .template value<double>(vars, x_position, t, dt);
609 }
610 auto const rho_fr =
611 fluid.property(MPL::PropertyType::density)
612 .template value<double>(vars, x_position, t, dt);
613 vars.density = rho_fr;
614
615 auto const mu = fluid.property(MPL::PropertyType::viscosity)
616 .template value<double>(vars, x_position, t, dt);
617 auto const beta_p =
618 fluid.property(MPL::PropertyType::density)
619 .template dValue<double>(vars, _process_data.phase_variable,
620 x_position, t, dt) /
621 rho_fr;
622
623 auto const K_over_mu = K / mu;
624
625 laplace.noalias() +=
626 rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
627
628 // Artificial compressibility from the fixed stress splitting:
629 auto const beta_FS =
630 fixed_stress_stabilization_parameter * alpha_b * alpha_b / K_S;
631
632 storage.noalias() += rho_fr * N_p.transpose() * N_p * w *
633 ((alpha_b - porosity) * (1.0 - alpha_b) / K_S +
634 porosity * beta_p + beta_FS);
635
636 auto const& b = _process_data.specific_body_force;
637
638 // bodyforce-driven Darcy flow
639 local_rhs.noalias() +=
640 dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
641
642 // density dependence on pressure evaluated for Darcy-term,
643 // for laplace and storage terms this dependence is neglected (as is
644 // done for monolithic too)
645 add_p_derivative.noalias() += rho_fr * beta_p * dNdx_p.transpose() *
646 K_over_mu *
647 (dNdx_p * p - 2.0 * rho_fr * b) * N_p * w;
648
649 if (!fixed_stress_over_time_step)
650 {
651 auto const& eps = _ip_data[ip].eps;
652 auto const& eps_prev = _ip_data[ip].eps_prev;
653 const double eps_v_dot =
654 (Invariants::trace(eps) - Invariants::trace(eps_prev)) / dt;
655
656 // Constant portion of strain rate term:
657 double const strain_rate_b =
658 alpha_b * eps_v_dot -
659 beta_FS * _ip_data[ip].strain_rate_variable;
660
661 local_rhs.noalias() -= strain_rate_b * rho_fr * N_p * w;
662 }
663 else
664 {
665 // Constant portion of strain rate term:
666 local_rhs.noalias() -=
667 alpha_b * _ip_data[ip].strain_rate_variable * rho_fr * N_p * w;
668 }
669 }
670 local_Jac.noalias() = laplace + storage / dt + add_p_derivative;
671
672 local_rhs.noalias() -= laplace * p + storage * (p - p_prev) / dt;
673}

References MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::gas_phase_pressure, MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::mechanical_strain, MaterialPropertyLib::VariableArray::molar_mass, MaterialLib::Solids::selectSolidConstitutiveRelation(), ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), 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 791 of file HydroMechanicsFEM-impl.h.

798{
799 // For the equations with pressure
800 if (process_id == _process_data.hydraulic_process_id)
801 {
802 assembleWithJacobianForPressureEquations(t, dt, local_x, local_x_prev,
803 local_b_data, local_Jac_data);
804 return;
805 }
806
807 // For the equations with deformation
808 assembleWithJacobianForDeformationEquations(t, dt, local_x, local_b_data,
809 local_Jac_data);
810}
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)
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)

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

1146{
1147 auto const p = local_x.template segment<pressure_size>(pressure_index);
1148
1150 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1151 DisplacementDim>(_element, _is_axially_symmetric, p,
1152 *_process_data.pressure_interpolated);
1153
1154 int const elem_id = _element.getID();
1156 x_position.setElementID(elem_id);
1157 unsigned const n_integration_points =
1159
1160 auto const& medium = _process_data.media_map.getMedium(elem_id);
1161 MPL::VariableArray vars;
1162 auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
1163 : vars.liquid_phase_pressure;
1164
1165 SymmetricTensor k_sum = SymmetricTensor::Zero(KelvinVectorSize);
1167 Eigen::Matrix<double, 3, 3>::Zero());
1168
1169 auto const& identity2 = Invariants::identity2;
1170
1171 for (unsigned ip = 0; ip < n_integration_points; ip++)
1172 {
1173 x_position.setIntegrationPoint(ip);
1174
1175 auto const& eps = _ip_data[ip].eps;
1176 sigma_eff_sum += _ip_data[ip].sigma_eff;
1177
1178 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
1179 .template value<double>(vars, x_position, t, dt);
1180 double const p_int_pt = _ip_data[ip].N_p.dot(p);
1181
1182 phase_pressure = p_int_pt;
1183
1184 // Set mechanical variables for the intrinsic permeability model
1185 // For stress dependent permeability.
1186 {
1187 auto const sigma_total =
1188 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
1189 vars.total_stress.emplace<SymmetricTensor>(
1191 sigma_total));
1192 }
1193 // For strain dependent permeability
1196 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
1199 eps);
1200
1202 medium->property(MPL::PropertyType::permeability)
1203 .value(vars, x_position, t, dt));
1204 }
1205
1206 Eigen::Map<Eigen::VectorXd>(
1207 &(*_process_data.permeability)[elem_id * KelvinVectorSize],
1208 KelvinVectorSize) = k_sum / n_integration_points;
1209
1210 Eigen::Matrix<double, 3, 3, 0, 3, 3> const sigma_avg =
1212 n_integration_points;
1213
1214 Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> e_s(sigma_avg);
1215
1216 Eigen::Map<Eigen::Vector3d>(
1217 &(*_process_data.principal_stress_values)[elem_id * 3], 3) =
1218 e_s.eigenvalues();
1219
1220 auto eigen_vectors = e_s.eigenvectors();
1221
1222 for (auto i = 0; i < 3; i++)
1223 {
1224 Eigen::Map<Eigen::Vector3d>(
1225 &(*_process_data.principal_stress_vector[i])[elem_id * 3], 3) =
1226 eigen_vectors.col(i);
1227 }
1228}
SymmetricTensor< GlobalDim > getSymmetricTensor(MaterialPropertyLib::PropertyDataType const &values)
Eigen::Matrix< double, 3, 3 > kelvinVectorToTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
KelvinVectorType< DisplacementDim > tensorToKelvin(Eigen::Matrix< double, 3, 3 > const &m)
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)

References MaterialPropertyLib::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(), ParameterLib::SpatialPosition::setIntegrationPoint(), 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 1098 of file HydroMechanicsFEM-impl.h.

1099{
1100 auto const kelvin_vector_size =
1102 unsigned const n_integration_points =
1104
1105 std::vector<double> ip_epsilon_values;
1106 auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
1107 double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
1108 ip_epsilon_values, n_integration_points, kelvin_vector_size);
1109
1110 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1111 {
1112 auto const& eps = _ip_data[ip].eps;
1113 cache_mat.row(ip) =
1115 }
1116
1117 return ip_epsilon_values;
1118}

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

391{
392 int const hydraulic_process_id = _process_data.hydraulic_process_id;
393 auto const indices =
394 NumLib::getIndices(_element.getID(), *dof_table[hydraulic_process_id]);
395 assert(!indices.empty());
396 auto const local_x = x[hydraulic_process_id]->get(indices);
397
398 unsigned const n_integration_points =
400 cache.clear();
401 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
402 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
403 cache, DisplacementDim, n_integration_points);
404
405 auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
406 pressure_size> const>(local_x.data() + pressure_index, pressure_size);
407
409 x_position.setElementID(_element.getID());
410
411 auto const& medium = _process_data.media_map.getMedium(_element.getID());
412 auto const& fluid = fluidPhase(*medium);
414 auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
416
417 // TODO (naumov) Temporary value not used by current material models. Need
418 // extension of secondary variables interface.
419 double const dt = std::numeric_limits<double>::quiet_NaN();
420 vars.temperature =
421 medium->property(MPL::PropertyType::reference_temperature)
422 .template value<double>(vars, x_position, t, dt);
423
424 auto const& identity2 = Invariants::identity2;
425
426 for (unsigned ip = 0; ip < n_integration_points; ip++)
427 {
428 x_position.setIntegrationPoint(ip);
429
430 double const p_int_pt = _ip_data[ip].N_p.dot(p);
431
432 phase_pressure = p_int_pt;
433
434 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
435 .template value<double>(vars, x_position, t, dt);
436
437 // Set mechanical variables for the intrinsic permeability model
438 // For stress dependent permeability.
439 auto const sigma_total =
440 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
441 vars.total_stress.emplace<SymmetricTensor>(
443
444 // For strain dependent permeability
447 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
450 _ip_data[ip].eps);
451
453 medium->property(MPL::PropertyType::permeability)
454 .value(vars, x_position, t, dt));
455
456 // Quick workaround: If fluid density is described as ideal gas, then
457 // the molar mass must be passed to the MPL::IdealGasLaw via the
458 // variable_array and the fluid must have the property
459 // MPL::PropertyType::molar_mass. For other density models (e.g.
460 // Constant), it is not mandatory to specify the molar mass.
461 if (fluid.hasProperty(MPL::PropertyType::molar_mass))
462 {
463 vars.molar_mass =
464 fluid.property(MPL::PropertyType::molar_mass)
465 .template value<double>(vars, x_position, t, dt);
466 }
467
468 auto const rho_fr =
469 fluid.property(MPL::PropertyType::density)
470 .template value<double>(vars, x_position, t, dt);
471 vars.density = rho_fr;
472
473 auto const mu = fluid.property(MPL::PropertyType::viscosity)
474 .template value<double>(vars, x_position, t, dt);
475
476 auto const K_over_mu = K / mu;
477
478 auto const& b = _process_data.specific_body_force;
479
480 // Compute the velocity
481 auto const& dNdx_p = _ip_data[ip].dNdx_p;
482 cache_matrix.col(ip).noalias() =
483 -K_over_mu * dNdx_p * p + K_over_mu * rho_fr * b;
484 }
485
486 return cache;
487}
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)

References MathLib::createZeroedMatrix(), MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::gas_phase_pressure, NumLib::getIndices(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::mechanical_strain, MaterialPropertyLib::VariableArray::molar_mass, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), 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

◆ 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

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

1244{
1245 return _process_data.material_ids == nullptr
1246 ? 0
1247 : (*_process_data.material_ids)[_element.getID()];
1248}

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

1257{
1258 return *_ip_data[integration_point].material_state_variables;
1259}

◆ getNumberOfIntegrationPoints()

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

◆ 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

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

1125{
1126 unsigned const n_integration_points =
1128
1129 std::vector<double> ip_strain_rate_variables(n_integration_points);
1130
1131 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1132 {
1133 ip_strain_rate_variables[ip] = _ip_data[ip].strain_rate_variable;
1134 }
1135
1136 return ip_strain_rate_variables;
1137}

◆ 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 =
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(), ip,
239 NumLib::interpolateCoordinates<ShapeFunctionPressure,
241 _element, ip_data.N_p))};
242
244 if (_process_data.initial_stress.value)
245 {
246 ip_data.sigma_eff =
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 }
Eigen::Matrix< double, Eigen::MatrixBase< Derived >::RowsAtCompileTime, 1 > symmetricTensorToKelvinVector(Eigen::MatrixBase< Derived > const &v)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)

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

879{
880 // Note: local_x and local_x_prev only contain the solutions of current
881 // process in the staggered scheme. This has to be changed according to the
882 // same two arguments in postTimestepConcrete.
883
884 int const n_integration_points = _integration_method.getNumberOfPoints();
885
886 auto const staggered_scheme_ptr =
887 std::get_if<Staggered>(&_process_data.coupling_scheme);
888
889 if (staggered_scheme_ptr &&
890 process_id == _process_data.hydraulic_process_id)
891 {
892 if (!staggered_scheme_ptr->fixed_stress_over_time_step)
893 {
894 auto const p =
895 local_x.template segment<pressure_size>(pressure_index);
896 auto const p_prev =
897 local_x_prev.template segment<pressure_size>(pressure_index);
898
899 for (int ip = 0; ip < n_integration_points; ip++)
900 {
901 auto& ip_data = _ip_data[ip];
902
903 auto const& N_p = ip_data.N_p;
904
905 ip_data.strain_rate_variable = N_p.dot(p - p_prev) / dt;
906 }
907 }
908 }
909
910 if (!staggered_scheme_ptr ||
911 process_id == _process_data.mechanics_related_process_id)
912 {
914 x_position.setElementID(_element.getID());
915
916 auto const& medium =
917 _process_data.media_map.getMedium(_element.getID());
918
919 auto const T_ref =
920 medium->property(MPL::PropertyType::reference_temperature)
921 .template value<double>(MPL::EmptyVariableArray, x_position, t,
922 dt);
923
924 auto const u =
925 local_x.template segment<displacement_size>(displacement_index);
926
928 vars.temperature = T_ref;
929
930 for (int ip = 0; ip < n_integration_points; ip++)
931 {
932 x_position.setIntegrationPoint(ip);
933 auto const& N_u = _ip_data[ip].N_u;
934 auto const& dNdx_u = _ip_data[ip].dNdx_u;
935
936 auto const x_coord =
937 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
939 _element, N_u);
940 auto const B = LinearBMatrix::computeBMatrix<
941 DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
942 typename BMatricesType::BMatrixType>(dNdx_u, N_u, x_coord,
944
945 auto& eps = _ip_data[ip].eps;
946 eps.noalias() = B * u;
947 vars.mechanical_strain.emplace<
949
950 _ip_data[ip].updateConstitutiveRelation(vars, t, x_position, dt, u,
951 T_ref);
952 }
953 }
954}
static const VariableArray EmptyVariableArray

References ProcessLib::LinearBMatrix::computeBMatrix(), MaterialPropertyLib::EmptyVariableArray, NumLib::interpolateXCoordinate(), MaterialPropertyLib::VariableArray::mechanical_strain, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), 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 960 of file HydroMechanicsFEM-impl.h.

964{
965 auto const staggered_scheme_ptr =
966 std::get_if<Staggered>(&_process_data.coupling_scheme);
967
968 if (staggered_scheme_ptr &&
969 process_id == _process_data.hydraulic_process_id)
970 {
971 if (staggered_scheme_ptr->fixed_stress_over_time_step)
972 {
973 auto const fixed_stress_stabilization_parameter =
974 staggered_scheme_ptr->fixed_stress_stabilization_parameter;
975
976 auto const p =
977 local_x.template segment<pressure_size>(pressure_index);
978 auto const p_prev =
979 local_x_prev.template segment<pressure_size>(pressure_index);
980
982 x_position.setElementID(_element.getID());
983
984 auto const& solid_material =
986 _process_data.solid_materials, _process_data.material_ids,
987 _element.getID());
988
989 auto const& medium =
990 _process_data.media_map.getMedium(_element.getID());
992
993 auto const T_ref =
994 medium->property(MPL::PropertyType::reference_temperature)
995 .template value<double>(vars, x_position, t, dt);
996 vars.temperature = T_ref;
997
998 int const n_integration_points =
1000 for (int ip = 0; ip < n_integration_points; ip++)
1001 {
1002 auto& ip_data = _ip_data[ip];
1003
1004 auto const& N_p = ip_data.N_p;
1005
1006 auto const& eps = ip_data.eps;
1007 auto const& eps_prev = ip_data.eps_prev;
1008 const double eps_v_dot =
1009 (Invariants::trace(eps) - Invariants::trace(eps_prev)) / dt;
1010
1011 auto const C_el = ip_data.computeElasticTangentStiffness(
1012 t, x_position, dt, T_ref);
1013 auto const K_S =
1014 solid_material.getBulkModulus(t, x_position, &C_el);
1015
1016 auto const alpha_b =
1017 medium->property(MPL::PropertyType::biot_coefficient)
1018 .template value<double>(vars, x_position, t, dt);
1019
1020 ip_data.strain_rate_variable =
1021 eps_v_dot - fixed_stress_stabilization_parameter * alpha_b *
1022 N_p.dot(p - p_prev) / dt / K_S;
1023 }
1024 }
1025 }
1026
1027 unsigned const n_integration_points =
1029
1030 for (unsigned ip = 0; ip < n_integration_points; ip++)
1031 {
1032 _ip_data[ip].pushBackState();
1033 }
1034}

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

819{
821 x_position.setElementID(_element.getID());
822
823 auto const& medium = _process_data.media_map.getMedium(_element.getID());
824
825 auto const p = local_x.template segment<pressure_size>(pressure_index);
826 auto const u =
827 local_x.template segment<displacement_size>(displacement_index);
828
829 auto const& identity2 = Invariants::identity2;
830 const double dt = 0.0;
831
833
834 int const n_integration_points = _integration_method.getNumberOfPoints();
835 for (int ip = 0; ip < n_integration_points; ip++)
836 {
837 x_position.setIntegrationPoint(ip);
838 auto const& N_u = _ip_data[ip].N_u;
839 auto const& dNdx_u = _ip_data[ip].dNdx_u;
840
841 auto const x_coord =
842 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
844 _element, N_u);
845 auto const B =
846 LinearBMatrix::computeBMatrix<DisplacementDim,
847 ShapeFunctionDisplacement::NPOINTS,
849 dNdx_u, N_u, x_coord, _is_axially_symmetric);
850
851 auto& eps = _ip_data[ip].eps;
852 eps.noalias() = B * u;
855 eps);
856
857 if (_process_data.initial_stress.isTotalStress())
858 {
859 auto const& N_p = _ip_data[ip].N_p;
860 auto const alpha_b =
861 medium->property(MPL::PropertyType::biot_coefficient)
862 .template value<double>(vars, x_position, t, dt);
863
864 auto& sigma_eff = _ip_data[ip].sigma_eff;
865 sigma_eff.noalias() += alpha_b * N_p.dot(p) * identity2;
866 _ip_data[ip].sigma_eff_prev.noalias() = sigma_eff;
867 }
868 }
869}

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

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

1043{
1044 if (integration_order !=
1045 static_cast<int>(_integration_method.getIntegrationOrder()))
1046 {
1047 OGS_FATAL(
1048 "Setting integration point initial conditions; The integration "
1049 "order of the local assembler for element {:d} is different from "
1050 "the integration order in the initial condition.",
1051 _element.getID());
1052 }
1053
1054 if (name == "sigma")
1055 {
1056 if (_process_data.initial_stress.value != nullptr)
1057 {
1058 OGS_FATAL(
1059 "Setting initial conditions for stress from integration "
1060 "point data and from a parameter '{:s}' is not possible "
1061 "simultaneously.",
1062 _process_data.initial_stress.value->name);
1063 }
1064
1066 values, _ip_data, &IpData::sigma_eff);
1067 }
1068
1069 if (name == "epsilon")
1070 {
1072 values, _ip_data, &IpData::eps);
1073 }
1074
1075 if (name == "strain_rate_variable")
1076 {
1079 }
1080
1081 return 0;
1082}
std::size_t setIntegrationPointKelvinVectorData(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType IpData::*const member)
std::size_t setIntegrationPointScalarData(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType IpData::*const member)

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

◆ _integration_method

◆ _ip_data

◆ _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

◆ _secondary_data

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

◆ 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:

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:
DisplacementDim,
constexpr Eigen::CwiseNullaryOp< EigenBlockMatrixViewFunctor< D, M >, typename EigenBlockMatrixViewFunctor< D, M >::Matrix > eigenBlockMatrixView(const Eigen::MatrixBase< M > &matrix)
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType

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: