OGS
List of all members | Public Types | Public Member Functions | Static Public Attributes | Private Types | Private Member Functions | Private Attributes | Static Private Attributes
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > Class Template Reference

Detailed Description

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

Definition at line 166 of file HydroMechanicsFEM.h.

#include <HydroMechanicsFEM.h>

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

Public Types

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

Public Member Functions

 HydroMechanicsLocalAssembler (HydroMechanicsLocalAssembler const &)=delete
 
 HydroMechanicsLocalAssembler (HydroMechanicsLocalAssembler &&)=delete
 
 HydroMechanicsLocalAssembler (MeshLib::Element const &e, std::size_t const, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, HydroMechanicsProcessData< DisplacementDim > &process_data)
 
std::size_t setIPDataInitialConditions (std::string_view const name, double const *values, int const integration_order) override
 Returns number of read integration points.
 
void assemble (double const, double const, std::vector< double > const &, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
 
void assembleWithJacobian (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_rhs_data, std::vector< double > &local_Jac_data) override
 
void assembleWithJacobianForStaggeredScheme (const double t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
 
void initializeConcrete () override
 
void postTimestepConcrete (Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, int const process_id) override
 
void computeSecondaryVariableConcrete (double const t, double const dt, Eigen::VectorXd const &local_xs, Eigen::VectorXd const &local_x_prev) override
 
void postNonLinearSolverConcrete (Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, int const process_id) override
 
void setInitialConditionsConcrete (Eigen::VectorXd const local_x, double const t, int const process_id) override
 
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix (const unsigned integration_point) const override
 Provides the shape matrix at the given integration point.
 
std::vector< double > getSigma () const override
 
std::vector< double > getEpsilon () const override
 
std::vector< double > getStrainRateVariable () const override
 
std::vector< double > const & getIntPtDarcyVelocity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
 
std::vector< double > const & getIntPtSigma (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
std::vector< double > const & getIntPtEpsilon (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
 
- Public Member Functions inherited from ProcessLib::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:
IntegrationPointData<BMatricesType, ShapeMatricesTypeDisplacement,
ShapeMatricesTypePressure, DisplacementDim,
ShapeFunctionDisplacement::NPOINTS>
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::ShapeMatricesTypeDisplacement
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
Definition HydroMechanicsFEM.h:170
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::ShapeMatricesTypePressure
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
Definition HydroMechanicsFEM.h:174
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::BMatricesType
BMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > BMatricesType
Definition HydroMechanicsFEM.h:386

Definition at line 388 of file HydroMechanicsFEM.h.

◆ ShapeMatricesTypeDisplacement

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

Definition at line 170 of file HydroMechanicsFEM.h.

◆ ShapeMatricesTypePressure

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

Definition at line 174 of file HydroMechanicsFEM.h.

◆ SymmetricTensor

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

Definition at line 181 of file HydroMechanicsFEM.h.

Constructor & Destructor Documentation

◆ HydroMechanicsLocalAssembler() [1/3]

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

◆ HydroMechanicsLocalAssembler() [2/3]

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

◆ HydroMechanicsLocalAssembler() [3/3]

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

Definition at line 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 =
49 _integration_method.getNumberOfPoints();
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,
56 ShapeMatricesTypeDisplacement,
57 DisplacementDim>(e, is_axially_symmetric,
58 _integration_method);
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 =
66 MaterialLib::Solids::selectSolidConstitutiveRelation(
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 =
76 _integration_method.getWeightedPoint(ip).getWeight() *
77 sm_u.integralMeasure * sm_u.detJ;
78
79 // Initialize current time step values
80 static const int kelvin_vector_size =
81 MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
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}
MathLib::WeightedPoint::getWeight
double getWeight() const
Definition WeightedPoint.h:80
NumLib::GenericIntegrationMethod::getWeightedPoint
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
Definition GenericIntegrationMethod.h:36
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_integration_method
NumLib::GenericIntegrationMethod const & _integration_method
Definition HydroMechanicsFEM.h:394
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_process_data
HydroMechanicsProcessData< DisplacementDim > & _process_data
Definition HydroMechanicsFEM.h:384
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_secondary_data
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
Definition HydroMechanicsFEM.h:399
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_ip_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
Definition HydroMechanicsFEM.h:392
MaterialLib::Solids::selectSolidConstitutiveRelation
auto & selectSolidConstitutiveRelation(SolidMaterialsMap const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
Definition SelectSolidConstitutiveRelation.h:36
MathLib::KelvinVector::kelvin_vector_dimensions
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition KelvinVector.h:24
NumLib::initShapeMatrices
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
Definition InitShapeMatrices.h:53
ProcessLib::HydroMechanics::SecondaryData::N_u
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u
Definition HydroMechanicsFEM.h:161

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

Member Function Documentation

◆ assemble()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 203 of file HydroMechanicsFEM.h.

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

References OGS_FATAL.

◆ assembleWithJacobian()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 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,
117 displacement_size);
118
119 auto p_prev =
120 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
121 pressure_size> const>(local_x_prev.data() + pressure_index,
122 pressure_size);
123 auto u_prev =
124 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
125 displacement_size> const>(local_x_prev.data() + displacement_index,
126 displacement_size);
127
128 auto local_Jac = MathLib::createZeroedMatrix<
129 typename ShapeMatricesTypeDisplacement::template MatrixType<
130 displacement_size + pressure_size,
131 displacement_size + pressure_size>>(
132 local_Jac_data, displacement_size + pressure_size,
133 displacement_size + pressure_size);
134
135 auto local_rhs = MathLib::createZeroedVector<
136 typename ShapeMatricesTypeDisplacement::template VectorType<
137 displacement_size + pressure_size>>(
138 local_rhs_data, displacement_size + pressure_size);
139
140 typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
141 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
142 pressure_size);
143
144 typename ShapeMatricesTypePressure::NodalMatrixType storage_p =
145 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
146 pressure_size);
147
148 typename ShapeMatricesTypePressure::NodalMatrixType add_p_derivative =
149 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
150 pressure_size);
151
152 typename ShapeMatricesTypeDisplacement::template MatrixType<
153 displacement_size, pressure_size>
154 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
155 displacement_size, pressure_size>::Zero(displacement_size,
156 pressure_size);
157
158 typename ShapeMatricesTypeDisplacement::template MatrixType<
159 pressure_size, displacement_size>
160 Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
161 pressure_size, displacement_size>::Zero(pressure_size,
162 displacement_size);
163
164 typename ShapeMatricesTypeDisplacement::template MatrixType<
165 pressure_size, displacement_size>
166 Kpu_k = ShapeMatricesTypeDisplacement::template MatrixType<
167 pressure_size, displacement_size>::Zero(pressure_size,
168 displacement_size);
169
170 auto const& solid_material =
171 MaterialLib::Solids::selectSolidConstitutiveRelation(
172 _process_data.solid_materials, _process_data.material_ids,
173 _element.getID());
174
175 ParameterLib::SpatialPosition x_position;
176 x_position.setElementID(_element.getID());
177
178 unsigned const n_integration_points =
179 _integration_method.getNumberOfPoints();
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);
185 auto const& phase_pressure = _process_data.phase_pressure;
186 MPL::VariableArray vars;
187
188 auto const T_ref =
189 medium->property(MPL::PropertyType::reference_temperature)
190 .template value<double>(vars, x_position, t, dt);
191 vars.temperature = T_ref;
192
193 auto const& identity2 = Invariants::identity2;
194
195 for (unsigned ip = 0; ip < n_integration_points; ip++)
196 {
197 x_position.setIntegrationPoint(ip);
198 auto const& w = _ip_data[ip].integration_weight;
199
200 auto const& N_u = _ip_data[ip].N_u;
201 auto const& dNdx_u = _ip_data[ip].dNdx_u;
202
203 auto const& N_p = _ip_data[ip].N_p;
204 auto const& dNdx_p = _ip_data[ip].dNdx_p;
205
206 auto const x_coord =
207 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
208 ShapeMatricesTypeDisplacement>(
209 _element, N_u);
210 auto const B =
211 LinearBMatrix::computeBMatrix<DisplacementDim,
212 ShapeFunctionDisplacement::NPOINTS,
213 typename BMatricesType::BMatrixType>(
214 dNdx_u, N_u, x_coord, _is_axially_symmetric);
215
216 auto& eps = _ip_data[ip].eps;
217 eps.noalias() = B * u;
218 auto const& sigma_eff = _ip_data[ip].sigma_eff;
219
220 double const p_int_pt = N_p.dot(p);
221
222 // setting both vars equal to enable MPL access for gas and liquid
223 // properties
224 vars.liquid_phase_pressure = vars.gas_phase_pressure = p_int_pt;
225
226 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
227 t, x_position, dt, T_ref);
228 auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
229
230 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
231 .template value<double>(vars, x_position, t, dt);
232
233 auto const rho_sr =
234 solid.property(MPL::PropertyType::density)
235 .template value<double>(vars, x_position, t, dt);
236 auto const porosity =
237 medium->property(MPL::PropertyType::porosity)
238 .template value<double>(vars, x_position, t, dt);
239
240 // Quick workaround: If fluid density is described as ideal gas, then
241 // the molar mass must be passed to the MPL::IdealGasLaw via the
242 // variable_array and the fluid must have the property
243 // MPL::PropertyType::molar_mass. For other density models (e.g.
244 // Constant), it is not mandatory to specify the molar mass.
245 if (fluid.hasProperty(MPL::PropertyType::molar_mass))
246 {
247 vars.molar_mass =
248 fluid.property(MPL::PropertyType::molar_mass)
249 .template value<double>(vars, x_position, t, dt);
250 }
251 auto const rho_fr =
252 fluid.property(MPL::PropertyType::density)
253 .template value<double>(vars, x_position, t, dt);
254 vars.density = rho_fr;
255
256 auto const mu = fluid.property(MPL::PropertyType::viscosity)
257 .template value<double>(vars, x_position, t, dt);
258
259 auto const beta_p = fluid.property(MPL::PropertyType::density)
260 .template dValue<double>(vars, phase_pressure,
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>(
271 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
272 sigma_total));
273 }
274 // For strain dependent permeability
275 vars.volumetric_strain = Invariants::trace(eps);
276 vars.equivalent_plastic_strain =
277 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
278 vars.mechanical_strain
279 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
280 eps);
281
282 auto const K = MPL::formEigenTensor<DisplacementDim>(
283 medium->property(MPL::PropertyType::permeability)
284 .value(vars, x_position, t, dt));
285 auto const dkde = MPL::formEigenTensor<
286 MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim)>(
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>(
300 displacement_index, displacement_index)
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() *
340 MathLib::KelvinVector::liftVectorToKelvin<DisplacementDim>(
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,
347 pressure_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,
364 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,
370 displacement_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}
MeshLib::Element::getID
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:89
ParameterLib::SpatialPosition::setElementID
void setElementID(std::size_t element_id)
Definition SpatialPosition.h:59
ParameterLib::SpatialPosition::setIntegrationPoint
void setIntegrationPoint(unsigned integration_point)
Definition SpatialPosition.h:65
ProcessLib::BMatrixPolicyType::BMatrixType
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition BMatrixPolicy.h:55
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::SymmetricTensor
Eigen::Matrix< double, KelvinVectorSize, 1 > SymmetricTensor
Definition HydroMechanicsFEM.h:181
MaterialLib::Fluid::PropertyVariableType::rho
@ rho
density. For some models, rho substitutes p
MaterialPropertyLib::formEigenTensor
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Definition FormEigenTensor.cpp:115
MaterialPropertyLib::fluidPhase
Phase const & fluidPhase(Medium const &medium)
Returns a gas or aqueous liquid phase of the given medium.
Definition Medium.cpp:100
MathLib::KelvinVector::kelvinVectorToSymmetricTensor
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Definition KelvinVector.cpp:142
MathLib::KelvinVector::KelvinVectorType
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition KelvinVector.h:47
MathLib::createZeroedVector
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Definition EigenMapTools.h:153
MathLib::createZeroedMatrix
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition EigenMapTools.h:32
NumLib::interpolateXCoordinate
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
Definition InitShapeMatrices.h:71
ProcessLib::Graph::eval
auto eval(Function &f, Tuples &... ts) -> typename detail::GetFunctionReturnType< decltype(&Function::eval)>::type
Definition Apply.h:274
ProcessLib::LinearBMatrix::computeBMatrix
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
Definition LinearBMatrix.h:46
EigenFixedShapeMatrixPolicy::NodalMatrixType
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
Definition ShapeMatrixPolicy.h:74
MathLib::KelvinVector::Invariants::identity2
static Eigen::Matrix< double, KelvinVectorSize, 1 > const identity2
Kelvin mapping of 2nd order identity tensor.
Definition KelvinVector.h:107
MathLib::KelvinVector::Invariants::trace
static double trace(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Trace of the corresponding tensor.
Definition KelvinVector-impl.h:62

References ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::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 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<
689 displacement_size, displacement_size>>(
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
697 ParameterLib::SpatialPosition x_position;
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);
703 MPL::VariableArray vars;
704
705 auto const T_ref =
706 medium->property(MPL::PropertyType::reference_temperature)
707 .template value<double>(vars, x_position, t, dt);
708 vars.temperature = T_ref;
709
710 int const n_integration_points = _integration_method.getNumberOfPoints();
711 for (int ip = 0; ip < n_integration_points; ip++)
712 {
713 x_position.setIntegrationPoint(ip);
714 auto const& w = _ip_data[ip].integration_weight;
715
716 auto const& N_u = _ip_data[ip].N_u;
717 auto const& dNdx_u = _ip_data[ip].dNdx_u;
718
719 auto const& N_p = _ip_data[ip].N_p;
720
721 auto const x_coord =
722 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
723 ShapeMatricesTypeDisplacement>(
724 _element, N_u);
725 auto const B =
726 LinearBMatrix::computeBMatrix<DisplacementDim,
727 ShapeFunctionDisplacement::NPOINTS,
728 typename BMatricesType::BMatrixType>(
729 dNdx_u, N_u, x_coord, _is_axially_symmetric);
730
731 auto& eps = _ip_data[ip].eps;
732 auto const& sigma_eff = _ip_data[ip].sigma_eff;
733
734 // setting both vars equal to enable MPL access for gas and liquid
735 // properties
736 vars.liquid_phase_pressure = vars.gas_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<
764 MathLib::KelvinVector::kelvin_vector_dimensions(
765 DisplacementDim)>::identity2;
766
767 eps.noalias() = B * u;
768 vars.mechanical_strain
769 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
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}
NumLib::detail::shapeFunctionInterpolate
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
Definition Interpolation.h:27

References ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::gas_phase_pressure, 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
503 ParameterLib::SpatialPosition pos;
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,
514 pressure_size);
515
516 typename ShapeMatricesTypePressure::NodalMatrixType storage =
517 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
518 pressure_size);
519
520 typename ShapeMatricesTypePressure::NodalMatrixType laplace =
521 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
522 pressure_size);
523
524 typename ShapeMatricesTypePressure::NodalMatrixType add_p_derivative =
525 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
526 pressure_size);
527
528 auto const& solid_material =
529 MaterialLib::Solids::selectSolidConstitutiveRelation(
530 _process_data.solid_materials, _process_data.material_ids,
531 _element.getID());
532
533 ParameterLib::SpatialPosition x_position;
534 x_position.setElementID(_element.getID());
535
536 auto const& medium = _process_data.media_map.getMedium(_element.getID());
537 auto const& fluid = fluidPhase(*medium);
538 auto const& phase_pressure = _process_data.phase_pressure;
539 MPL::VariableArray vars;
540
541 auto const T_ref =
542 medium->property(MPL::PropertyType::reference_temperature)
543 .template value<double>(vars, x_position, t, dt);
544 vars.temperature = T_ref;
545
546 auto const& identity2 = Invariants::identity2;
547
548 auto const staggered_scheme =
549 std::get<Staggered>(_process_data.coupling_scheme);
550 auto const fixed_stress_stabilization_parameter =
551 staggered_scheme.fixed_stress_stabilization_parameter;
552 auto const fixed_stress_over_time_step =
553 staggered_scheme.fixed_stress_over_time_step;
554
555 int const n_integration_points = _integration_method.getNumberOfPoints();
556 for (int ip = 0; ip < n_integration_points; ip++)
557 {
558 x_position.setIntegrationPoint(ip);
559 auto const& w = _ip_data[ip].integration_weight;
560
561 auto const& N_p = _ip_data[ip].N_p;
562 auto const& dNdx_p = _ip_data[ip].dNdx_p;
563
564 double const p_int_pt = N_p.dot(p);
565
566 // setting both vars equal to enable MPL access for gas and liquid
567 // properties
568 vars.liquid_phase_pressure = vars.gas_phase_pressure = p_int_pt;
569
570 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
571 t, x_position, dt, T_ref);
572 auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
573
574 auto const alpha_b =
575 medium->property(MPL::PropertyType::biot_coefficient)
576 .template value<double>(vars, x_position, t, dt);
577
578 // Set mechanical variables for the intrinsic permeability model
579 // For stress dependent permeability.
580 auto const sigma_total =
581 (_ip_data[ip].sigma_eff - alpha_b * p_int_pt * identity2).eval();
582 vars.total_stress.emplace<SymmetricTensor>(
583 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_total));
584
585 // For strain dependent permeability
586 vars.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
587 vars.equivalent_plastic_strain =
588 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
589 vars.mechanical_strain
590 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
591 _ip_data[ip].eps);
592
593 auto const K = MPL::formEigenTensor<DisplacementDim>(
594 medium->property(MPL::PropertyType::permeability)
595 .value(vars, x_position, t, dt));
596 auto const porosity =
597 medium->property(MPL::PropertyType::porosity)
598 .template value<double>(vars, x_position, t, dt);
599
600 // Quick workaround: If fluid density is described as ideal gas, then
601 // the molar mass must be passed to the MPL::IdealGasLaw via the
602 // variable_array and the fluid must have the property
603 // MPL::PropertyType::molar_mass. For other density models (e.g.
604 // Constant), it is not mandatory to specify the molar mass.
605 if (fluid.hasProperty(MPL::PropertyType::molar_mass))
606 {
607 vars.molar_mass =
608 fluid.property(MPL::PropertyType::molar_mass)
609 .template value<double>(vars, x_position, t, dt);
610 }
611 auto const rho_fr =
612 fluid.property(MPL::PropertyType::density)
613 .template value<double>(vars, x_position, t, dt);
614 vars.density = rho_fr;
615
616 auto const mu = fluid.property(MPL::PropertyType::viscosity)
617 .template value<double>(vars, x_position, t, dt);
618 auto const beta_p = fluid.property(MPL::PropertyType::density)
619 .template dValue<double>(vars, phase_pressure,
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::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}
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::assembleWithJacobianForPressureEquations
void assembleWithJacobianForPressureEquations(const double t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
Definition HydroMechanicsFEM-impl.h:493
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::assembleWithJacobianForDeformationEquations
void assembleWithJacobianForDeformationEquations(const double t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
Definition HydroMechanicsFEM-impl.h:679

◆ 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
1149 NumLib::interpolateToHigherOrderNodes<
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();
1155 ParameterLib::SpatialPosition x_position;
1156 x_position.setElementID(elem_id);
1157 unsigned const n_integration_points =
1158 _integration_method.getNumberOfPoints();
1159
1160 auto const& medium = _process_data.media_map.getMedium(elem_id);
1161 MPL::VariableArray vars;
1162
1163 SymmetricTensor k_sum = SymmetricTensor::Zero(KelvinVectorSize);
1164 auto sigma_eff_sum = MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
1165 Eigen::Matrix<double, 3, 3>::Zero());
1166
1167 auto const& identity2 = Invariants::identity2;
1168
1169 for (unsigned ip = 0; ip < n_integration_points; ip++)
1170 {
1171 x_position.setIntegrationPoint(ip);
1172
1173 auto const& eps = _ip_data[ip].eps;
1174 sigma_eff_sum += _ip_data[ip].sigma_eff;
1175
1176 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
1177 .template value<double>(vars, x_position, t, dt);
1178 double const p_int_pt = _ip_data[ip].N_p.dot(p);
1179
1180 // setting both vars equal to enable MPL access for gas and liquid
1181 // properties
1182 vars.liquid_phase_pressure = vars.gas_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>(
1190 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
1191 sigma_total));
1192 }
1193 // For strain dependent permeability
1194 vars.volumetric_strain = Invariants::trace(eps);
1195 vars.equivalent_plastic_strain =
1196 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
1197 vars.mechanical_strain
1198 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
1199 eps);
1200
1201 k_sum += MPL::getSymmetricTensor<DisplacementDim>(
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 =
1211 MathLib::KelvinVector::kelvinVectorToTensor(sigma_eff_sum) /
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}
MathLib::KelvinVector::kelvinVectorToTensor
Eigen::Matrix< double, 3, 3 > kelvinVectorToTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Definition KelvinVector.cpp:64
NumLib::interpolateToHigherOrderNodes
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
Definition Interpolation.h:107

References MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::VariableArray::gas_phase_pressure, 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 ( ) const
overridevirtual

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

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

1099{
1100 auto const kelvin_vector_size =
1101 MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
1102 unsigned const n_integration_points =
1103 _integration_method.getNumberOfPoints();
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) =
1114 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(eps);
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 =
399 _integration_method.getNumberOfPoints();
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
408 ParameterLib::SpatialPosition x_position;
409 x_position.setElementID(_element.getID());
410
411 auto const& medium = _process_data.media_map.getMedium(_element.getID());
412 auto const& fluid = fluidPhase(*medium);
413 MPL::VariableArray vars;
414
415 // TODO (naumov) Temporary value not used by current material models. Need
416 // extension of secondary variables interface.
417 double const dt = std::numeric_limits<double>::quiet_NaN();
418 vars.temperature =
419 medium->property(MPL::PropertyType::reference_temperature)
420 .template value<double>(vars, x_position, t, dt);
421
422 auto const& identity2 = Invariants::identity2;
423
424 for (unsigned ip = 0; ip < n_integration_points; ip++)
425 {
426 x_position.setIntegrationPoint(ip);
427
428 double const p_int_pt = _ip_data[ip].N_p.dot(p);
429
430 // setting both vars equal to enable MPL access for gas and liquid
431 // properties
432 vars.liquid_phase_pressure = vars.gas_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>(
442 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_total));
443
444 // For strain dependent permeability
445 vars.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
446 vars.equivalent_plastic_strain =
447 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
448 vars.mechanical_strain
449 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
450 _ip_data[ip].eps);
451
452 auto const K = MPL::formEigenTensor<DisplacementDim>(
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}
NumLib::getIndices
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
Definition DOFTableUtil.cpp:127

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

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

Definition at line 314 of file HydroMechanicsFEM.h.

319 {
320 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
321 _ip_data, &IpData::eps, cache);
322 }

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

◆ getIntPtSigma()

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

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

Definition at line 304 of file HydroMechanicsFEM.h.

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

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

◆ getMaterialID()

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

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

Definition at line 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

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

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

1235{
1236 return _integration_method.getNumberOfPoints();
1237}

◆ getShapeMatrix()

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

Provides the shape matrix at the given integration point.

Implements NumLib::ExtrapolatableElement.

Definition at line 280 of file HydroMechanicsFEM.h.

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

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

◆ getSigma()

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

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

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

1089{
1090 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
1091 _ip_data, &IpData::sigma_eff);
1092}

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

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 {
913 ParameterLib::SpatialPosition x_position;
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
927 MPL::VariableArray vars;
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,
938 ShapeMatricesTypeDisplacement>(
939 _element, N_u);
940 auto const B = LinearBMatrix::computeBMatrix<
941 DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
942 typename BMatricesType::BMatrixType>(dNdx_u, N_u, x_coord,
943 _is_axially_symmetric);
944
945 auto& eps = _ip_data[ip].eps;
946 eps.noalias() = B * u;
947 vars.mechanical_strain.emplace<
948 MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(eps);
949
950 _ip_data[ip].updateConstitutiveRelation(vars, t, x_position, dt, u,
951 T_ref);
952 }
953 }
954}
MaterialPropertyLib::EmptyVariableArray
static const VariableArray EmptyVariableArray
Definition VariableType.h:202

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
981 ParameterLib::SpatialPosition x_position;
982 x_position.setElementID(_element.getID());
983
984 auto const& solid_material =
985 MaterialLib::Solids::selectSolidConstitutiveRelation(
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());
991 MPL::VariableArray vars;
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 =
999 _integration_method.getNumberOfPoints();
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 =
1028 _integration_method.getNumberOfPoints();
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{
820 ParameterLib::SpatialPosition x_position;
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
832 MPL::VariableArray vars;
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,
843 ShapeMatricesTypeDisplacement>(
844 _element, N_u);
845 auto const B =
846 LinearBMatrix::computeBMatrix<DisplacementDim,
847 ShapeFunctionDisplacement::NPOINTS,
848 typename BMatricesType::BMatrixType>(
849 dNdx_u, N_u, x_coord, _is_axially_symmetric);
850
851 auto& eps = _ip_data[ip].eps;
852 eps.noalias() = B * u;
853 vars.mechanical_strain
854 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
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
1065 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
1066 values, _ip_data, &IpData::sigma_eff);
1067 }
1068
1069 if (name == "epsilon")
1070 {
1071 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
1072 values, _ip_data, &IpData::eps);
1073 }
1074
1075 if (name == "strain_rate_variable")
1076 {
1077 return ProcessLib::setIntegrationPointScalarData(
1078 values, _ip_data, &IpData::strain_rate_variable);
1079 }
1080
1081 return 0;
1082}
ProcessLib::setIntegrationPointScalarData
std::size_t setIntegrationPointScalarData(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType IpData::*const member)
Definition SetOrGetIntegrationPointData.h:215

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

Definition at line 395 of file HydroMechanicsFEM.h.

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

◆ _integration_method

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

Definition at line 394 of file HydroMechanicsFEM.h.

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

◆ _ip_data

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

Definition at line 392 of file HydroMechanicsFEM.h.

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

◆ _is_axially_symmetric

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

Definition at line 396 of file HydroMechanicsFEM.h.

◆ _process_data

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

Definition at line 384 of file HydroMechanicsFEM.h.

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

◆ _secondary_data

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

Definition at line 399 of file HydroMechanicsFEM.h.

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

◆ displacement_index

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

Definition at line 403 of file HydroMechanicsFEM.h.

◆ displacement_size

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

Definition at line 404 of file HydroMechanicsFEM.h.

◆ KelvinVectorSize

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

Definition at line 177 of file HydroMechanicsFEM.h.

◆ N_u_op

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

Definition at line 183 of file HydroMechanicsFEM.h.

◆ pressure_index

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

Definition at line 401 of file HydroMechanicsFEM.h.

◆ pressure_size

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

Definition at line 402 of file HydroMechanicsFEM.h.

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