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

Detailed Description

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

Definition at line 166 of file HydroMechanicsFEM.h.

#include <HydroMechanicsFEM.h>

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

Public Types

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

Public Member Functions

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

Static Public Attributes

static int const KelvinVectorSize
 
static constexpr auto & N_u_op
 

Private Types

using BMatricesType
 
using IpData
 

Private Member Functions

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

Private Attributes

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

Static Private Attributes

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

Member Typedef Documentation

◆ BMatricesType

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

Definition at line 386 of file HydroMechanicsFEM.h.

◆ Invariants

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

Definition at line 179 of file HydroMechanicsFEM.h.

◆ IpData

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

Definition at line 388 of file HydroMechanicsFEM.h.

◆ ShapeMatricesTypeDisplacement

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

Definition at line 170 of file HydroMechanicsFEM.h.

◆ ShapeMatricesTypePressure

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

Definition at line 174 of file HydroMechanicsFEM.h.

◆ SymmetricTensor

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

Definition at line 181 of file HydroMechanicsFEM.h.

Constructor & Destructor Documentation

◆ HydroMechanicsLocalAssembler() [1/3]

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

◆ HydroMechanicsLocalAssembler() [2/3]

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

◆ HydroMechanicsLocalAssembler() [3/3]

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

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

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

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

Member Function Documentation

◆ assemble()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 203 of file HydroMechanicsFEM.h.

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

References OGS_FATAL.

◆ assembleWithJacobian()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

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

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

◆ assembleWithJacobianForDeformationEquations()

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

Assemble local matrices and vectors arise from the linearized discretized weak form of the residual of the momentum balance equation,

\[ \nabla (\sigma - \alpha_b p \mathrm{I}) = f \]

where \( \sigma\) is the effective stress tensor, \(p\) is the pore pressure, \(\alpha_b\) is the Biot constant, \(\mathrm{I}\) is the identity tensor, and \(f\) is the body force.

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

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

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

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

◆ assembleWithJacobianForPressureEquations()

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

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

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

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

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

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

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

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

◆ assembleWithJacobianForStaggeredScheme()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

760{
761 // For the equations with pressure
762 if (process_id == _process_data.hydraulic_process_id)
763 {
764 assembleWithJacobianForPressureEquations(t, dt, local_x, local_x_prev,
765 local_b_data, local_Jac_data);
766 return;
767 }
768
769 // For the equations with deformation
770 assembleWithJacobianForDeformationEquations(t, dt, local_x, local_b_data,
771 local_Jac_data);
772}
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::assembleWithJacobianForPressureEquations
void assembleWithJacobianForPressureEquations(const double t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
Definition HydroMechanicsFEM-impl.h:472
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::assembleWithJacobianForDeformationEquations
void assembleWithJacobianForDeformationEquations(const double t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
Definition HydroMechanicsFEM-impl.h:650

◆ computeSecondaryVariableConcrete()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

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

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

◆ getEpsilon()

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

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

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

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

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

◆ getIntPtDarcyVelocity()

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

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

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

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

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

◆ getIntPtEpsilon()

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

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

Definition at line 314 of file HydroMechanicsFEM.h.

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

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

◆ getIntPtSigma()

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

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

Definition at line 304 of file HydroMechanicsFEM.h.

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

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

◆ getMaterialID()

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

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

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

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

◆ getMaterialStateVariablesAt()

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

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

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

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

◆ getNumberOfIntegrationPoints()

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

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

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

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

◆ getShapeMatrix()

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

Provides the shape matrix at the given integration point.

Implements NumLib::ExtrapolatableElement.

Definition at line 280 of file HydroMechanicsFEM.h.

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

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

◆ getSigma()

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

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

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

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

References ProcessLib::getIntegrationPointKelvinVectorData().

◆ getStrainRateVariable()

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

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

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

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

◆ initializeConcrete()

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

Set initial stress from parameter.

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 227 of file HydroMechanicsFEM.h.

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

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

◆ postNonLinearSolverConcrete()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

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

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

◆ postTimestepConcrete()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

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

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

◆ setInitialConditionsConcrete()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

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

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

◆ setIPDataInitialConditions()

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

Returns number of read integration points.

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

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

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

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

Member Data Documentation

◆ _element

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

Definition at line 395 of file HydroMechanicsFEM.h.

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

◆ _integration_method

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

Definition at line 394 of file HydroMechanicsFEM.h.

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

◆ _ip_data

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

Definition at line 392 of file HydroMechanicsFEM.h.

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

◆ _is_axially_symmetric

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

Definition at line 396 of file HydroMechanicsFEM.h.

◆ _process_data

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

Definition at line 384 of file HydroMechanicsFEM.h.

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

◆ _secondary_data

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

Definition at line 399 of file HydroMechanicsFEM.h.

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

◆ displacement_index

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

Definition at line 403 of file HydroMechanicsFEM.h.

◆ displacement_size

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

Definition at line 404 of file HydroMechanicsFEM.h.

◆ KelvinVectorSize

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

Definition at line 177 of file HydroMechanicsFEM.h.

◆ N_u_op

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

Definition at line 183 of file HydroMechanicsFEM.h.

◆ pressure_index

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

Definition at line 401 of file HydroMechanicsFEM.h.

◆ pressure_size

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

Definition at line 402 of file HydroMechanicsFEM.h.

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