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 > &, std::vector< double > &, std::vector< double > &local_rhs_data, std::vector< double > &local_Jac_data) override
 
void assembleWithJacobianForStaggeredScheme (const double t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
 
void initializeConcrete () override
 
void postTimestepConcrete (Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, 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 389 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:389

Definition at line 391 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:397
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_process_data
HydroMechanicsProcessData< DisplacementDim > & _process_data
Definition HydroMechanicsFEM.h:387
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_secondary_data
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
Definition HydroMechanicsFEM.h:402
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_ip_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
Definition HydroMechanicsFEM.h:395
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 > & ,
std::vector< double > & ,
std::vector< double > & local_rhs_data,
std::vector< double > & local_Jac_data )
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

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

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

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

References MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::VariableArray::gas_phase_pressure, MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::mechanical_strain, MaterialPropertyLib::VariableArray::molar_mass, MaterialLib::Solids::selectSolidConstitutiveRelation(), ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, and MaterialPropertyLib::VariableArray::volumetric_strain.

◆ assembleWithJacobianForStaggeredScheme()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

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

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

1148{
1149 auto const p = local_x.template segment<pressure_size>(pressure_index);
1150
1151 NumLib::interpolateToHigherOrderNodes<
1152 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1153 DisplacementDim>(_element, _is_axially_symmetric, p,
1154 *_process_data.pressure_interpolated);
1155
1156 int const elem_id = _element.getID();
1157 ParameterLib::SpatialPosition x_position;
1158 x_position.setElementID(elem_id);
1159 unsigned const n_integration_points =
1160 _integration_method.getNumberOfPoints();
1161
1162 auto const& medium = _process_data.media_map.getMedium(elem_id);
1163 MPL::VariableArray vars;
1164
1165 SymmetricTensor k_sum = SymmetricTensor::Zero(KelvinVectorSize);
1166 auto sigma_eff_sum = MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
1167 Eigen::Matrix<double, 3, 3>::Zero());
1168
1169 auto const& identity2 = Invariants::identity2;
1170
1171 for (unsigned ip = 0; ip < n_integration_points; ip++)
1172 {
1173 x_position.setIntegrationPoint(ip);
1174
1175 auto const& eps = _ip_data[ip].eps;
1176 sigma_eff_sum += _ip_data[ip].sigma_eff;
1177
1178 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
1179 .template value<double>(vars, x_position, t, dt);
1180 double const p_int_pt = _ip_data[ip].N_p.dot(p);
1181
1182 // setting both vars equal to enable MPL access for gas and liquid
1183 // properties
1184 vars.liquid_phase_pressure = vars.gas_phase_pressure = p_int_pt;
1185
1186 // Set mechanical variables for the intrinsic permeability model
1187 // For stress dependent permeability.
1188 {
1189 auto const sigma_total =
1190 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
1191 vars.total_stress.emplace<SymmetricTensor>(
1192 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
1193 sigma_total));
1194 }
1195 // For strain dependent permeability
1196 vars.volumetric_strain = Invariants::trace(eps);
1197 vars.equivalent_plastic_strain =
1198 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
1199 vars.mechanical_strain
1200 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
1201 eps);
1202
1203 k_sum += MPL::getSymmetricTensor<DisplacementDim>(
1204 medium->property(MPL::PropertyType::permeability)
1205 .value(vars, x_position, t, dt));
1206 }
1207
1208 Eigen::Map<Eigen::VectorXd>(
1209 &(*_process_data.permeability)[elem_id * KelvinVectorSize],
1210 KelvinVectorSize) = k_sum / n_integration_points;
1211
1212 Eigen::Matrix<double, 3, 3, 0, 3, 3> const sigma_avg =
1213 MathLib::KelvinVector::kelvinVectorToTensor(sigma_eff_sum) /
1214 n_integration_points;
1215
1216 Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> e_s(sigma_avg);
1217
1218 Eigen::Map<Eigen::Vector3d>(
1219 &(*_process_data.principal_stress_values)[elem_id * 3], 3) =
1220 e_s.eigenvalues();
1221
1222 auto eigen_vectors = e_s.eigenvectors();
1223
1224 for (auto i = 0; i < 3; i++)
1225 {
1226 Eigen::Map<Eigen::Vector3d>(
1227 &(*_process_data.principal_stress_vector[i])[elem_id * 3], 3) =
1228 eigen_vectors.col(i);
1229 }
1230}
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 1100 of file HydroMechanicsFEM-impl.h.

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

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

◆ getIntPtDarcyVelocity()

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

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

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

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

322 {
323 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
324 _ip_data, &IpData::eps, cache);
325 }

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

312 {
313 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
314 _ip_data, &IpData::sigma_eff, cache);
315 }

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

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

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

1259{
1260 return *_ip_data[integration_point].material_state_variables;
1261}

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

1237{
1238 return _integration_method.getNumberOfPoints();
1239}

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

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

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

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

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

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

◆ initializeConcrete()

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

Set initial stress from parameter.

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 230 of file HydroMechanicsFEM.h.

231 {
232 unsigned const n_integration_points =
233 _integration_method.getNumberOfPoints();
234
235 for (unsigned ip = 0; ip < n_integration_points; ip++)
236 {
237 auto& ip_data = _ip_data[ip];
238
239 ParameterLib::SpatialPosition const x_position{
240 std::nullopt, _element.getID(), ip,
241 MathLib::Point3d(
242 NumLib::interpolateCoordinates<ShapeFunctionPressure,
243 ShapeMatricesTypePressure>(
244 _element, ip_data.N_p))};
245
247 if (_process_data.initial_stress.value)
248 {
249 ip_data.sigma_eff =
250 MathLib::KelvinVector::symmetricTensorToKelvinVector<
251 DisplacementDim>((*_process_data.initial_stress.value)(
252 std::numeric_limits<
253 double>::quiet_NaN() /* time independent */,
254 x_position));
255 }
256
257 double const t = 0; // TODO (naumov) pass t from top
258 ip_data.solid_material.initializeInternalStateVariables(
259 t, x_position, *ip_data.material_state_variables);
260
261 ip_data.pushBackState();
262 }
263 }
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 876 of file HydroMechanicsFEM-impl.h.

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

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

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

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

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

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

1045{
1046 if (integration_order !=
1047 static_cast<int>(_integration_method.getIntegrationOrder()))
1048 {
1049 OGS_FATAL(
1050 "Setting integration point initial conditions; The integration "
1051 "order of the local assembler for element {:d} is different from "
1052 "the integration order in the initial condition.",
1053 _element.getID());
1054 }
1055
1056 if (name == "sigma")
1057 {
1058 if (_process_data.initial_stress.value != nullptr)
1059 {
1060 OGS_FATAL(
1061 "Setting initial conditions for stress from integration "
1062 "point data and from a parameter '{:s}' is not possible "
1063 "simultaneously.",
1064 _process_data.initial_stress.value->name);
1065 }
1066
1067 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
1068 values, _ip_data, &IpData::sigma_eff);
1069 }
1070
1071 if (name == "epsilon")
1072 {
1073 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
1074 values, _ip_data, &IpData::eps);
1075 }
1076
1077 if (name == "strain_rate_variable")
1078 {
1079 return ProcessLib::setIntegrationPointScalarData(
1080 values, _ip_data, &IpData::strain_rate_variable);
1081 }
1082
1083 return 0;
1084}
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 398 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 397 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 395 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 399 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 387 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 402 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 406 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 407 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 404 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 405 of file HydroMechanicsFEM.h.

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