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 = ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim >
 
using ShapeMatricesTypePressure = ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim >
 
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 > &, std::vector< double > &, 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_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
 
void initializeConcrete () override
 
void postTimestepConcrete (Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, bool const use_monolithic_scheme, 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 (std::vector< double > const &local_x, std::vector< double > const &local_x_prev, double const t, double const dt, bool const use_monolithic_scheme, int const process_id) override
 
void setInitialConditionsConcrete (std::vector< double > const &local_x, double const t, bool const use_monolithic_scheme, 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
 
virtual std::size_t setIPDataInitialConditions (std::string_view const name, double const *values, int const integration_order)=0
 
virtual std::vector< double > getSigma () const =0
 
virtual std::vector< double > getEpsilon () const =0
 
virtual std::vector< double > getStrainRateVariable () const =0
 
virtual std::vector< double > const & getIntPtSigma (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual std::vector< double > const & getIntPtEpsilon (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
 
virtual 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 =0
 
virtual unsigned getNumberOfIntegrationPoints () const =0
 
virtual int getMaterialID () const =0
 
virtual MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt (unsigned) const =0
 
- Public Member Functions inherited from ProcessLib::LocalAssemblerInterface
virtual ~LocalAssemblerInterface ()=default
 
void setInitialConditions (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, bool const use_monolithic_scheme, int const process_id)
 
virtual void initialize (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
 
virtual void preAssemble (double const, double const, std::vector< double > const &)
 
virtual void assemble (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
 
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 assembleWithJacobian (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
 
virtual void assembleWithJacobianForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_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, bool const use_monolithic_scheme, int const process_id)
 
void postNonLinearSolver (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, bool const use_monolithic_scheme, 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 Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix (const unsigned integration_point) const =0
 Provides the shape matrix at the given integration point.
 
virtual ~ExtrapolatableElement ()=default
 

Static Public Attributes

static int const KelvinVectorSize
 
static constexpr auto & N_u_op
 

Private Types

using BMatricesType = BMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim >
 
using IpData = IntegrationPointData< BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, ShapeFunctionDisplacement::NPOINTS >
 

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 = BMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>
private

Definition at line 392 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 = IntegrationPointData<BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, ShapeFunctionDisplacement::NPOINTS>
private

Definition at line 394 of file HydroMechanicsFEM.h.

◆ ShapeMatricesTypeDisplacement

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ShapeMatricesTypeDisplacement = 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 = 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 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
Definition: WeightedPoint.h:80
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
NumLib::GenericIntegrationMethod const & _integration_method
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
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.
Definition: KelvinVector.h:24
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 > &  ,
std::vector< double > &  ,
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.

110{
111 assert(local_x.size() == pressure_size + displacement_size);
112
113 auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
114 pressure_size> const>(local_x.data() + pressure_index, pressure_size);
115
116 auto u =
117 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
118 displacement_size> const>(local_x.data() + displacement_index,
120
121 auto p_prev =
122 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
123 pressure_size> const>(local_x_prev.data() + pressure_index,
125 auto u_prev =
126 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
127 displacement_size> const>(local_x_prev.data() + displacement_index,
129
130 auto local_Jac = MathLib::createZeroedMatrix<
131 typename ShapeMatricesTypeDisplacement::template MatrixType<
134 local_Jac_data, displacement_size + pressure_size,
136
137 auto local_rhs = MathLib::createZeroedVector<
138 typename ShapeMatricesTypeDisplacement::template VectorType<
140 local_rhs_data, displacement_size + pressure_size);
141
143 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
145
147 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
149
150 typename ShapeMatricesTypePressure::NodalMatrixType add_p_derivative =
151 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
153
154 typename ShapeMatricesTypeDisplacement::template MatrixType<
156 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
159
160 typename ShapeMatricesTypeDisplacement::template MatrixType<
162 Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
165
166 typename ShapeMatricesTypeDisplacement::template MatrixType<
168 Kpu_k = ShapeMatricesTypeDisplacement::template MatrixType<
171
172 auto const& solid_material =
174 _process_data.solid_materials, _process_data.material_ids,
175 _element.getID());
176
178 x_position.setElementID(_element.getID());
179
180 unsigned const n_integration_points =
182
183 auto const& b = _process_data.specific_body_force;
184 auto const& medium = _process_data.media_map.getMedium(_element.getID());
185 auto const& solid = medium->phase("Solid");
186 auto const& fluid = fluidPhase(*medium);
187 auto const& phase_pressure = _process_data.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 x_position.setIntegrationPoint(ip);
200 auto const& w = _ip_data[ip].integration_weight;
201
202 auto const& N_u = _ip_data[ip].N_u;
203 auto const& dNdx_u = _ip_data[ip].dNdx_u;
204
205 auto const& N_p = _ip_data[ip].N_p;
206 auto const& dNdx_p = _ip_data[ip].dNdx_p;
207
208 auto const x_coord =
209 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
211 _element, N_u);
212 auto const B =
213 LinearBMatrix::computeBMatrix<DisplacementDim,
214 ShapeFunctionDisplacement::NPOINTS,
216 dNdx_u, N_u, x_coord, _is_axially_symmetric);
217
218 auto& eps = _ip_data[ip].eps;
219 eps.noalias() = B * u;
220 auto const& sigma_eff = _ip_data[ip].sigma_eff;
221
222 double const p_int_pt = N_p.dot(p);
223
224 // setting both vars equal to enable MPL access for gas and liquid
225 // properties
226 vars.liquid_phase_pressure = vars.gas_phase_pressure = p_int_pt;
227
228 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
229 t, x_position, dt, T_ref);
230 auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
231
232 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
233 .template value<double>(vars, x_position, t, dt);
234
235 auto const rho_sr =
236 solid.property(MPL::PropertyType::density)
237 .template value<double>(vars, x_position, t, dt);
238 auto const porosity =
239 medium->property(MPL::PropertyType::porosity)
240 .template value<double>(vars, x_position, t, dt);
241
242 // Quick workaround: If fluid density is described as ideal gas, then
243 // the molar mass must be passed to the MPL::IdealGasLaw via the
244 // variable_array and the fluid must have the property
245 // MPL::PropertyType::molar_mass. For other density models (e.g.
246 // Constant), it is not mandatory to specify the molar mass.
247 if (fluid.hasProperty(MPL::PropertyType::molar_mass))
248 {
249 vars.molar_mass =
250 fluid.property(MPL::PropertyType::molar_mass)
251 .template value<double>(vars, x_position, t, dt);
252 }
253 auto const rho_fr =
254 fluid.property(MPL::PropertyType::density)
255 .template value<double>(vars, x_position, t, dt);
256 vars.density = rho_fr;
257
258 auto const mu = fluid.property(MPL::PropertyType::viscosity)
259 .template value<double>(vars, x_position, t, dt);
260
261 auto const beta_p = fluid.property(MPL::PropertyType::density)
262 .template dValue<double>(vars, phase_pressure,
263 x_position, t, dt) /
264 rho_fr;
265
266 // Set mechanical variables for the intrinsic permeability model
267 // For stress dependent permeability.
268 {
269 auto const sigma_total =
270 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
271
272 vars.total_stress.emplace<SymmetricTensor>(
274 sigma_total));
275 }
276 // For strain dependent permeability
279 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
282 eps);
283
284 auto const K = MPL::formEigenTensor<DisplacementDim>(
285 medium->property(MPL::PropertyType::permeability)
286 .value(vars, x_position, t, dt));
287 auto const dkde = MPL::formEigenTensor<
289 (*medium)[MPL::PropertyType::permeability].dValue(
290 vars, MPL::Variable::mechanical_strain, x_position, t, dt));
291
292 auto const K_over_mu = K / mu;
293
294 auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
295 dt, u, T_ref);
296
297 //
298 // displacement equation, displacement part
299 //
300 local_Jac
301 .template block<displacement_size, displacement_size>(
303 .noalias() += B.transpose() * C * B * w;
304
305 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
306 local_rhs.template segment<displacement_size>(displacement_index)
307 .noalias() -=
308 (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
309
310 //
311 // displacement equation, pressure part
312 //
313 Kup.noalias() += B.transpose() * alpha * identity2 * N_p * w;
314
315 //
316 // pressure equation, pressure part.
317 //
318 laplace_p.noalias() +=
319 rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
320
321 storage_p.noalias() +=
322 rho_fr * N_p.transpose() * N_p * w *
323 ((alpha - porosity) * (1.0 - alpha) / K_S + porosity * beta_p);
324
325 // density dependence on pressure evaluated for Darcy-term,
326 // for laplace and storage terms this dependence is neglected
327 add_p_derivative.noalias() += rho_fr * beta_p * dNdx_p.transpose() *
328 K_over_mu *
329 (dNdx_p * p - 2.0 * rho_fr * b) * N_p * w;
330
331 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
332 dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
333
334 //
335 // pressure equation, displacement part.
336 //
337 Kpu.noalias() +=
338 rho_fr * alpha * N_p.transpose() * identity2.transpose() * B * w;
339
340 Kpu_k.noalias() +=
341 dNdx_p.transpose() *
342 MathLib::KelvinVector::liftVectorToKelvin<DisplacementDim>(
343 dNdx_p * p - rho_fr * b) *
344 dkde * B * rho_fr / mu * w;
345 }
346 // displacement equation, pressure part
347 local_Jac
348 .template block<displacement_size, pressure_size>(displacement_index,
350 .noalias() = -Kup;
351
352 if (_process_data.mass_lumping)
353 {
354 storage_p = storage_p.colwise().sum().eval().asDiagonal();
355
356 if constexpr (pressure_size == displacement_size)
357 {
358 Kpu = Kpu.colwise().sum().eval().asDiagonal();
359 Kpu_k = Kpu_k.colwise().sum().eval().asDiagonal();
360 }
361 }
362
363 // pressure equation, pressure part.
364 local_Jac
365 .template block<pressure_size, pressure_size>(pressure_index,
367 .noalias() += laplace_p + storage_p / dt + add_p_derivative;
368
369 // pressure equation, displacement part.
370 local_Jac
371 .template block<pressure_size, displacement_size>(pressure_index,
373 .noalias() += Kpu / dt + Kpu_k;
374
375 // pressure equation
376 local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
377 laplace_p * p + storage_p * (p - p_prev) / dt + Kpu * (u - u_prev) / dt;
378
379 // displacement equation
380 local_rhs.template segment<displacement_size>(displacement_index)
381 .noalias() += Kup * p;
382}
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > mechanical_strain
Definition: VariableType.h:182
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > total_stress
Definition: VariableType.h:201
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
Definition: BMatrixPolicy.h:55
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, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:49
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 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)
Definition: EigenMapTools.h:32
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:267
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
Definition: LinearBMatrix.h:42
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
static Eigen::Matrix< double, KelvinVectorSize, 1 > const identity2
Kelvin mapping of 2nd order identity tensor.
Definition: KelvinVector.h:107
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(), 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 680 of file HydroMechanicsFEM-impl.h.

684{
685 auto const p = local_x.template segment<pressure_size>(pressure_index);
686 auto const u =
687 local_x.template segment<displacement_size>(displacement_index);
688
689 auto local_Jac = MathLib::createZeroedMatrix<
690 typename ShapeMatricesTypeDisplacement::template MatrixType<
692 local_Jac_data, displacement_size, displacement_size);
693
694 auto local_rhs =
695 MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
696 template VectorType<displacement_size>>(
697 local_b_data, displacement_size);
698
700 x_position.setElementID(_element.getID());
701
702 auto const& medium = _process_data.media_map.getMedium(_element.getID());
703 auto const& solid = medium->phase("Solid");
704 auto const& fluid = fluidPhase(*medium);
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 // setting both vars equal to enable MPL access for gas and liquid
737 // properties
738 vars.liquid_phase_pressure = vars.gas_phase_pressure = N_p.dot(p);
739
740 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
741 .template value<double>(vars, x_position, t, dt);
742 auto const rho_sr =
743 solid.property(MPL::PropertyType::density)
744 .template value<double>(vars, x_position, t, dt);
745 auto const porosity =
746 medium->property(MPL::PropertyType::porosity)
747 .template value<double>(vars, x_position, t, dt);
748
749 // Quick workaround: If fluid density is described as ideal gas, then
750 // the molar mass must be passed to the MPL::IdealGasLaw via the
751 // variable_array and the fluid must have the property
752 // MPL::PropertyType::molar_mass. For other density models (e.g.
753 // Constant), it is not mandatory to specify the molar mass.
754 if (fluid.hasProperty(MPL::PropertyType::molar_mass))
755 {
756 vars.molar_mass =
757 fluid.property(MPL::PropertyType::molar_mass)
758 .template value<double>(vars, x_position, t, dt);
759 }
760 auto const rho_fr =
761 fluid.property(MPL::PropertyType::density)
762 .template value<double>(vars, x_position, t, dt);
763
764 auto const& b = _process_data.specific_body_force;
765 auto const& identity2 = MathLib::KelvinVector::Invariants<
767 DisplacementDim)>::identity2;
768
769 eps.noalias() = B * u;
772 eps);
773
774 auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
775 dt, u, T_ref);
776
777 local_Jac.noalias() += B.transpose() * C * B * w;
778
779 double p_at_xi = 0.;
780 NumLib::shapeFunctionInterpolate(p, N_p, p_at_xi);
781
782 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
783 local_rhs.noalias() -=
784 (B.transpose() * (sigma_eff - alpha * identity2 * p_at_xi) -
785 N_u_op(N_u).transpose() * rho * b) *
786 w;
787 }
788}
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
Definition: Interpolation.h:27

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

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

References MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, 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_M_data,
std::vector< double > &  local_K_data,
std::vector< double > &  local_b_data,
std::vector< double > &  local_Jac_data 
)
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

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

1180{
1181 auto const p = local_x.template segment<pressure_size>(pressure_index);
1182
1184 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1185 DisplacementDim>(_element, _is_axially_symmetric, p,
1186 *_process_data.pressure_interpolated);
1187
1188 int const elem_id = _element.getID();
1190 x_position.setElementID(elem_id);
1191 unsigned const n_integration_points =
1193
1194 auto const& medium = _process_data.media_map.getMedium(elem_id);
1195 MPL::VariableArray vars;
1196
1197 SymmetricTensor k_sum = SymmetricTensor::Zero(KelvinVectorSize);
1198 auto sigma_eff_sum = MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
1199 Eigen::Matrix<double, 3, 3>::Zero());
1200
1201 auto const& identity2 = Invariants::identity2;
1202
1203 for (unsigned ip = 0; ip < n_integration_points; ip++)
1204 {
1205 x_position.setIntegrationPoint(ip);
1206
1207 auto const& eps = _ip_data[ip].eps;
1208 sigma_eff_sum += _ip_data[ip].sigma_eff;
1209
1210 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
1211 .template value<double>(vars, x_position, t, dt);
1212 double const p_int_pt = _ip_data[ip].N_p.dot(p);
1213
1214 // setting both vars equal to enable MPL access for gas and liquid
1215 // properties
1216 vars.liquid_phase_pressure = vars.gas_phase_pressure = p_int_pt;
1217
1218 // Set mechanical variables for the intrinsic permeability model
1219 // For stress dependent permeability.
1220 {
1221 auto const sigma_total =
1222 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
1223 vars.total_stress.emplace<SymmetricTensor>(
1225 sigma_total));
1226 }
1227 // For strain dependent permeability
1230 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
1233 eps);
1234
1235 k_sum += MPL::getSymmetricTensor<DisplacementDim>(
1236 medium->property(MPL::PropertyType::permeability)
1237 .value(vars, x_position, t, dt));
1238 }
1239
1240 Eigen::Map<Eigen::VectorXd>(
1241 &(*_process_data.permeability)[elem_id * KelvinVectorSize],
1242 KelvinVectorSize) = k_sum / n_integration_points;
1243
1244 Eigen::Matrix<double, 3, 3, 0, 3, 3> const sigma_avg =
1246 n_integration_points;
1247
1248 Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> e_s(sigma_avg);
1249
1250 Eigen::Map<Eigen::Vector3d>(
1251 &(*_process_data.principal_stress_values)[elem_id * 3], 3) =
1252 e_s.eigenvalues();
1253
1254 auto eigen_vectors = e_s.eigenvectors();
1255
1256 for (auto i = 0; i < 3; i++)
1257 {
1258 Eigen::Map<Eigen::Vector3d>(
1259 &(*_process_data.principal_stress_vector[i])[elem_id * 3], 3) =
1260 eigen_vectors.col(i);
1261 }
1262}
Eigen::Matrix< double, 3, 3 > kelvinVectorToTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
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, NumLib::interpolateToHigherOrderNodes(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MathLib::KelvinVector::kelvinVectorToTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::mechanical_strain, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), 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
overridevirtual

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

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

1133{
1134 auto const kelvin_vector_size =
1136 unsigned const n_integration_points =
1138
1139 std::vector<double> ip_epsilon_values;
1140 auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
1141 double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
1142 ip_epsilon_values, n_integration_points, kelvin_vector_size);
1143
1144 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1145 {
1146 auto const& eps = _ip_data[ip].eps;
1147 cache_mat.row(ip) =
1149 }
1150
1151 return ip_epsilon_values;
1152}

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

393{
394 int const hydraulic_process_id = _process_data.hydraulic_process_id;
395 auto const indices =
396 NumLib::getIndices(_element.getID(), *dof_table[hydraulic_process_id]);
397 assert(!indices.empty());
398 auto const local_x = x[hydraulic_process_id]->get(indices);
399
400 unsigned const n_integration_points =
402 cache.clear();
403 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
404 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
405 cache, DisplacementDim, n_integration_points);
406
407 auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
408 pressure_size> const>(local_x.data() + pressure_index, pressure_size);
409
411 x_position.setElementID(_element.getID());
412
413 auto const& medium = _process_data.media_map.getMedium(_element.getID());
414 auto const& fluid = fluidPhase(*medium);
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 // setting both vars equal to enable MPL access for gas and liquid
433 // properties
434 vars.liquid_phase_pressure = vars.gas_phase_pressure = p_int_pt;
435
436 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
437 .template value<double>(vars, x_position, t, dt);
438
439 // Set mechanical variables for the intrinsic permeability model
440 // For stress dependent permeability.
441 auto const sigma_total =
442 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
443 vars.total_stress.emplace<SymmetricTensor>(
445
446 // For strain dependent permeability
449 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
452 _ip_data[ip].eps);
453
454 auto const K = MPL::formEigenTensor<DisplacementDim>(
455 medium->property(MPL::PropertyType::permeability)
456 .value(vars, x_position, t, dt));
457
458 // Quick workaround: If fluid density is described as ideal gas, then
459 // the molar mass must be passed to the MPL::IdealGasLaw via the
460 // variable_array and the fluid must have the property
461 // MPL::PropertyType::molar_mass. For other density models (e.g.
462 // Constant), it is not mandatory to specify the molar mass.
463 if (fluid.hasProperty(MPL::PropertyType::molar_mass))
464 {
465 vars.molar_mass =
466 fluid.property(MPL::PropertyType::molar_mass)
467 .template value<double>(vars, x_position, t, dt);
468 }
469
470 auto const rho_fr =
471 fluid.property(MPL::PropertyType::density)
472 .template value<double>(vars, x_position, t, dt);
473 vars.density = rho_fr;
474
475 auto const mu = fluid.property(MPL::PropertyType::viscosity)
476 .template value<double>(vars, x_position, t, dt);
477
478 auto const K_over_mu = K / mu;
479
480 auto const& b = _process_data.specific_body_force;
481
482 // Compute the velocity
483 auto const& dNdx_p = _ip_data[ip].dNdx_p;
484 cache_matrix.col(ip).noalias() =
485 -K_over_mu * dNdx_p * p + K_over_mu * rho_fr * b;
486 }
487
488 return cache;
489}
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::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(), 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
overrideprivatevirtual

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

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

1278{
1279 return _process_data.material_ids == nullptr
1280 ? 0
1281 : (*_process_data.material_ids)[_element.getID()];
1282}

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

1291{
1292 return *_ip_data[integration_point].material_state_variables;
1293}

◆ getNumberOfIntegrationPoints()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
unsigned ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getNumberOfIntegrationPoints
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 286 of file HydroMechanicsFEM.h.

288 {
289 auto const& N_u = _secondary_data.N_u[integration_point];
290
291 // assumes N is stored contiguously in memory
292 return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
293 }

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
overridevirtual

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

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

1123{
1124 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
1126}

◆ getStrainRateVariable()

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

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

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

1159{
1160 unsigned const n_integration_points =
1162
1163 std::vector<double> ip_strain_rate_variables(n_integration_points);
1164
1165 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1166 {
1167 ip_strain_rate_variables[ip] = _ip_data[ip].strain_rate_variable;
1168 }
1169
1170 return ip_strain_rate_variables;
1171}

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

231 {
232 unsigned const n_integration_points =
234
235 for (unsigned ip = 0; ip < n_integration_points; ip++)
236 {
237 auto& ip_data = _ip_data[ip];
238
239 ParameterLib::SpatialPosition const x_position{
240 std::nullopt, _element.getID(), ip,
242 NumLib::interpolateCoordinates<ShapeFunctionPressure,
244 _element, ip_data.N_p))};
245
247 if (_process_data.initial_stress.value)
248 {
249 ip_data.sigma_eff =
251 DisplacementDim>((*_process_data.initial_stress.value)(
252 std::numeric_limits<
253 double>::quiet_NaN() /* time independent */,
254 x_position));
255 }
256
257 double const t = 0; // TODO (naumov) pass t from top
258 ip_data.solid_material.initializeInternalStateVariables(
259 t, x_position, *ip_data.material_state_variables);
260
261 ip_data.pushBackState();
262 }
263 }
Eigen::Matrix< double, Eigen::MatrixBase< Derived >::RowsAtCompileTime, 1 > symmetricTensorToKelvinVector(Eigen::MatrixBase< Derived > const &v)
Definition: KelvinVector.h:201
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 ( std::vector< double > const &  local_x,
std::vector< double > const &  local_x_prev,
double const  t,
double const  dt,
bool const  use_monolithic_scheme,
int const  process_id 
)
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

903{
904 // Note: local_x and local_x_prev only contain the solutions of current
905 // process in the staggered scheme. This has to be changed according to the
906 // same two arguments in postTimestepConcrete.
907
908 int const n_integration_points = _integration_method.getNumberOfPoints();
909
910 auto const staggered_scheme_ptr =
911 std::get_if<Staggered>(&_process_data.coupling_scheme);
912
913 if (staggered_scheme_ptr &&
914 process_id == _process_data.hydraulic_process_id)
915 {
916 if (!staggered_scheme_ptr->fixed_stress_over_time_step)
917 {
918 auto const p =
919 Eigen::Map<typename ShapeMatricesTypePressure::
920 template VectorType<pressure_size> const>(
921 local_x.data(), pressure_size);
922
923 auto const p_prev =
924 Eigen::Map<typename ShapeMatricesTypePressure::
925 template VectorType<pressure_size> const>(
926 local_x_prev.data(), pressure_size);
927
928 for (int ip = 0; ip < n_integration_points; ip++)
929 {
930 auto& ip_data = _ip_data[ip];
931
932 auto const& N_p = ip_data.N_p;
933
934 ip_data.strain_rate_variable = N_p.dot(p - p_prev) / dt;
935 }
936 }
937 }
938
939 if (!staggered_scheme_ptr ||
940 process_id == _process_data.mechanics_related_process_id)
941 {
943 x_position.setElementID(_element.getID());
944
945 auto const& medium =
946 _process_data.media_map.getMedium(_element.getID());
947
948 auto const T_ref =
949 medium->property(MPL::PropertyType::reference_temperature)
950 .template value<double>(MPL::EmptyVariableArray, x_position, t,
951 dt);
952
953 const int displacement_offset =
954 (!staggered_scheme_ptr) ? displacement_index : 0;
955
956 auto u = Eigen::Map<typename ShapeMatricesTypeDisplacement::
957 template VectorType<displacement_size> const>(
958 local_x.data() + displacement_offset, displacement_size);
959
961 vars.temperature = T_ref;
962
963 for (int ip = 0; ip < n_integration_points; ip++)
964 {
965 x_position.setIntegrationPoint(ip);
966 auto const& N_u = _ip_data[ip].N_u;
967 auto const& dNdx_u = _ip_data[ip].dNdx_u;
968
969 auto const x_coord =
970 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
972 _element, N_u);
973 auto const B = LinearBMatrix::computeBMatrix<
974 DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
975 typename BMatricesType::BMatrixType>(dNdx_u, N_u, x_coord,
977
978 auto& eps = _ip_data[ip].eps;
979 eps.noalias() = B * u;
980 vars.mechanical_strain.emplace<
982
983 _ip_data[ip].updateConstitutiveRelation(vars, t, x_position, dt, u,
984 T_ref);
985 }
986 }
987}
static const VariableArray EmptyVariableArray
Definition: VariableType.h:210

References ProcessLib::LinearBMatrix::computeBMatrix(), MaterialPropertyLib::EmptyVariableArray, NumLib::interpolateXCoordinate(), 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,
bool const  use_monolithic_scheme,
int const  process_id 
)
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

998{
999 auto const staggered_scheme_ptr =
1000 std::get_if<Staggered>(&_process_data.coupling_scheme);
1001
1002 if (staggered_scheme_ptr &&
1003 process_id == _process_data.hydraulic_process_id)
1004 {
1005 if (staggered_scheme_ptr->fixed_stress_over_time_step)
1006 {
1007 auto const fixed_stress_stabilization_parameter =
1008 staggered_scheme_ptr->fixed_stress_stabilization_parameter;
1009
1010 auto const p =
1011 local_x.template segment<pressure_size>(pressure_index);
1012 auto const p_prev =
1013 local_x_prev.template segment<pressure_size>(pressure_index);
1014
1016 x_position.setElementID(_element.getID());
1017
1018 auto const& solid_material =
1020 _process_data.solid_materials, _process_data.material_ids,
1021 _element.getID());
1022
1023 auto const& medium =
1024 _process_data.media_map.getMedium(_element.getID());
1025 MPL::VariableArray vars;
1026
1027 auto const T_ref =
1028 medium->property(MPL::PropertyType::reference_temperature)
1029 .template value<double>(vars, x_position, t, dt);
1030 vars.temperature = T_ref;
1031
1032 int const n_integration_points =
1034 for (int ip = 0; ip < n_integration_points; ip++)
1035 {
1036 auto& ip_data = _ip_data[ip];
1037
1038 auto const& N_p = ip_data.N_p;
1039
1040 auto const& eps = ip_data.eps;
1041 auto const& eps_prev = ip_data.eps_prev;
1042 const double eps_v_dot =
1043 (Invariants::trace(eps) - Invariants::trace(eps_prev)) / dt;
1044
1045 auto const C_el = ip_data.computeElasticTangentStiffness(
1046 t, x_position, dt, T_ref);
1047 auto const K_S =
1048 solid_material.getBulkModulus(t, x_position, &C_el);
1049
1050 auto const alpha_b =
1051 medium->property(MPL::PropertyType::biot_coefficient)
1052 .template value<double>(vars, x_position, t, dt);
1053
1054 ip_data.strain_rate_variable =
1055 eps_v_dot - fixed_stress_stabilization_parameter * alpha_b *
1056 N_p.dot(p - p_prev) / dt / K_S;
1057 }
1058 }
1059 }
1060
1061 unsigned const n_integration_points =
1063
1064 for (unsigned ip = 0; ip < n_integration_points; ip++)
1065 {
1066 _ip_data[ip].pushBackState();
1067 }
1068}

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 ( std::vector< double > const &  local_x,
double const  t,
bool const  use_monolithic_scheme,
int const  process_id 
)
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

822{
823 if (!use_monolithic_scheme &&
824 process_id == _process_data.hydraulic_process_id)
825 {
826 return;
827 }
828
829 if (use_monolithic_scheme ||
830 process_id == _process_data.mechanics_related_process_id)
831 {
833 x_position.setElementID(_element.getID());
834
835 auto const& medium =
836 _process_data.media_map.getMedium(_element.getID());
837
838 int const displacement_offset =
839 use_monolithic_scheme ? displacement_index : 0;
840
841 auto const u =
842 Eigen::Map<typename ShapeMatricesTypeDisplacement::
843 template VectorType<displacement_size> const>(
844 local_x.data() + displacement_offset, displacement_size);
845
846 int const pressure_offset = use_monolithic_scheme ? pressure_index : 0;
847
848 auto const p =
849 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
850 pressure_size> const>(local_x.data() + pressure_offset,
852 auto const& identity2 = Invariants::identity2;
853 const double dt = 0.0;
854
856
857 int const n_integration_points =
859 for (int ip = 0; ip < n_integration_points; ip++)
860 {
861 x_position.setIntegrationPoint(ip);
862 auto const& N_u = _ip_data[ip].N_u;
863 auto const& dNdx_u = _ip_data[ip].dNdx_u;
864
865 auto const x_coord =
866 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
868 _element, N_u);
869 auto const B = LinearBMatrix::computeBMatrix<
870 DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
871 typename BMatricesType::BMatrixType>(dNdx_u, N_u, x_coord,
873
874 auto& eps = _ip_data[ip].eps;
875 eps.noalias() = B * u;
876 vars.mechanical_strain.emplace<
878
879 if (_process_data.initial_stress.isTotalStress())
880 {
881 auto const& N_p = _ip_data[ip].N_p;
882 auto const alpha_b =
883 medium->property(MPL::PropertyType::biot_coefficient)
884 .template value<double>(vars, x_position, t, dt);
885
886 auto& sigma_eff = _ip_data[ip].sigma_eff;
887 sigma_eff.noalias() += alpha_b * N_p.dot(p) * identity2;
888 _ip_data[ip].sigma_eff_prev.noalias() = sigma_eff;
889 }
890 }
891 }
892}

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

1077{
1078 if (integration_order !=
1079 static_cast<int>(_integration_method.getIntegrationOrder()))
1080 {
1081 OGS_FATAL(
1082 "Setting integration point initial conditions; The integration "
1083 "order of the local assembler for element {:d} is different from "
1084 "the integration order in the initial condition.",
1085 _element.getID());
1086 }
1087
1088 if (name == "sigma")
1089 {
1090 if (_process_data.initial_stress.value != nullptr)
1091 {
1092 OGS_FATAL(
1093 "Setting initial conditions for stress from integration "
1094 "point data and from a parameter '{:s}' is not possible "
1095 "simultaneously.",
1096 _process_data.initial_stress.value->name);
1097 }
1098
1099 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
1100 values, _ip_data, &IpData::sigma_eff);
1101 }
1102
1103 if (name == "epsilon")
1104 {
1105 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
1106 values, _ip_data, &IpData::eps);
1107 }
1108
1109 if (name == "strain_rate_variable")
1110 {
1113 }
1114
1115 return 0;
1116}
std::size_t setIntegrationPointScalarData(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType IpData::*const member)

References OGS_FATAL, 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 402 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 409 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 410 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>
constexpr 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 407 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 408 of file HydroMechanicsFEM.h.


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