OGS
HydroMechanicsFEM.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <memory>
14#include <vector>
15
32
33namespace ProcessLib
34{
35namespace HydroMechanics
36{
37namespace MPL = MaterialPropertyLib;
38
39template <typename BMatricesType, typename ShapeMatricesTypeDisplacement,
40 typename ShapeMatricesTypePressure, int DisplacementDim, int NPoints>
42{
51
52 typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev;
53 typename BMatricesType::KelvinVectorType eps, eps_prev;
54
55 typename ShapeMatricesTypeDisplacement::NodalRowVectorType N_u;
56 typename ShapeMatricesTypeDisplacement::GlobalDimNodalMatrixType dNdx_u;
57
58 typename ShapeMatricesTypePressure::NodalRowVectorType N_p;
59 typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p;
60
62 std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
63 DisplacementDim>::MaterialStateVariables>
66
67 // previous pressure rate for the fixed stress splitting
68 // approach in the staggered scheme.
69 // TODO: disable in monolithic scheme to save memory.
71
73 {
74 eps_prev = eps;
76 material_state_variables->pushBackState();
77 }
78
81 double const t,
82 ParameterLib::SpatialPosition const& x_position,
83 double const dt,
84 double const temperature)
85 {
86 namespace MPL = MaterialPropertyLib;
87
88 MPL::VariableArray variable_array;
89 MPL::VariableArray variable_array_prev;
90
91 auto const null_state = solid_material.createMaterialStateVariables();
92 solid_material.initializeInternalStateVariables(t, x_position,
93 *null_state);
94
96
97 variable_array.stress.emplace<KV>(KV::Zero());
98 variable_array.mechanical_strain.emplace<KV>(KV::Zero());
99 variable_array.temperature = temperature;
100
101 variable_array_prev.stress.emplace<KV>(KV::Zero());
102 variable_array_prev.mechanical_strain.emplace<KV>(KV::Zero());
103 variable_array_prev.temperature = temperature;
104
105 auto&& solution =
106 solid_material.integrateStress(variable_array_prev, variable_array,
107 t, x_position, dt, *null_state);
108
109 if (!solution)
110 {
111 OGS_FATAL("Computation of elastic tangent stiffness failed.");
112 }
113
115 std::move(std::get<2>(*solution));
116
117 return C;
118 }
119
120 template <typename DisplacementVectorType>
121 typename BMatricesType::KelvinMatrixType updateConstitutiveRelation(
122 MPL::VariableArray const& variable_array,
123 double const t,
124 ParameterLib::SpatialPosition const& x_position,
125 double const dt,
126 DisplacementVectorType const& /*u*/,
127 double const T)
128 {
129 MaterialPropertyLib::VariableArray variable_array_prev;
130 variable_array_prev.stress
133 variable_array_prev.mechanical_strain
135 eps_prev);
136 variable_array_prev.temperature = T;
137
138 auto&& solution = solid_material.integrateStress(
139 variable_array_prev, variable_array, t, x_position, dt,
141
142 if (!solution)
143 {
144 OGS_FATAL("Computation of local constitutive relation failed.");
145 }
146
148 std::tie(sigma_eff, material_state_variables, C) = std::move(*solution);
149
150 return C;
151 }
152
154};
155
158template <typename ShapeMatrixType>
160{
161 std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N_u;
162};
163
164template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
165 int DisplacementDim>
167 : public LocalAssemblerInterface<DisplacementDim>
168{
169public:
172
173 // Types for pressure.
176
177 static int const KelvinVectorSize =
180
181 using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>;
182
183 static constexpr auto& N_u_op = MathLib::eigenBlockMatrixView<
184 DisplacementDim,
186
189
191 MeshLib::Element const& e,
192 std::size_t const /*local_matrix_size*/,
193 NumLib::GenericIntegrationMethod const& integration_method,
194 bool const is_axially_symmetric,
196
198 std::size_t setIPDataInitialConditions(
199 std::string_view const name,
200 double const* values,
201 int const integration_order) override;
202
203 void assemble(double const /*t*/, double const /*dt*/,
204 std::vector<double> const& /*local_x*/,
205 std::vector<double> const& /*local_x_prev*/,
206 std::vector<double>& /*local_M_data*/,
207 std::vector<double>& /*local_K_data*/,
208 std::vector<double>& /*local_rhs_data*/) override
209 {
210 OGS_FATAL(
211 "HydroMechanicsLocalAssembler: assembly without jacobian is not "
212 "implemented.");
213 }
214
215 void assembleWithJacobian(double const t, double const dt,
216 std::vector<double> const& local_x,
217 std::vector<double> const& local_x_prev,
218 std::vector<double>& /*local_M_data*/,
219 std::vector<double>& /*local_K_data*/,
220 std::vector<double>& local_rhs_data,
221 std::vector<double>& local_Jac_data) override;
222
224 const double t, double const dt, Eigen::VectorXd const& local_x,
225 Eigen::VectorXd const& local_x_prev, int const process_id,
226 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
227 std::vector<double>& local_b_data,
228 std::vector<double>& local_Jac_data) override;
229
230 void initializeConcrete() override
231 {
232 unsigned const n_integration_points =
234
235 for (unsigned ip = 0; ip < n_integration_points; ip++)
236 {
237 auto& ip_data = _ip_data[ip];
238
239 ParameterLib::SpatialPosition const x_position{
240 std::nullopt, _element.getID(), ip,
242 NumLib::interpolateCoordinates<ShapeFunctionPressure,
244 _element, ip_data.N_p))};
245
247 if (_process_data.initial_stress.value)
248 {
249 ip_data.sigma_eff =
251 DisplacementDim>((*_process_data.initial_stress.value)(
252 std::numeric_limits<
253 double>::quiet_NaN() /* time independent */,
254 x_position));
255 }
256
257 double const t = 0; // TODO (naumov) pass t from top
258 ip_data.solid_material.initializeInternalStateVariables(
259 t, x_position, *ip_data.material_state_variables);
260
261 ip_data.pushBackState();
262 }
263 }
264
265 void postTimestepConcrete(Eigen::VectorXd const& local_x,
266 Eigen::VectorXd const& local_x_prev,
267 double const t, double const dt,
268 int const process_id) override;
269
271 double const t, double const dt, Eigen::VectorXd const& local_xs,
272 Eigen::VectorXd const& local_x_prev) override;
273
274 void postNonLinearSolverConcrete(std::vector<double> const& local_x,
275 std::vector<double> const& local_x_prev,
276 double const t, double const dt,
277 int const process_id) override;
278
279 void setInitialConditionsConcrete(Eigen::VectorXd const local_x,
280 double const t,
281 int const process_id) override;
282
283 Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
284 const unsigned integration_point) const override
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 }
291
292 // TODO (naumov) This method is same as getIntPtSigma but for arguments and
293 // the ordering of the cache_mat.
294 // There should be only one.
295 std::vector<double> getSigma() const override;
296
297 std::vector<double> getEpsilon() const override;
298
299 std::vector<double> getStrainRateVariable() const override;
300
301 std::vector<double> const& getIntPtDarcyVelocity(
302 const double t,
303 std::vector<GlobalVector*> const& x,
304 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
305 std::vector<double>& cache) const override;
306
307 std::vector<double> const& getIntPtSigma(
308 const double /*t*/,
309 std::vector<GlobalVector*> const& /*x*/,
310 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
311 std::vector<double>& cache) const override
312 {
313 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
314 _ip_data, &IpData::sigma_eff, cache);
315 }
316
317 std::vector<double> const& getIntPtEpsilon(
318 const double /*t*/,
319 std::vector<GlobalVector*> const& /*x*/,
320 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
321 std::vector<double>& cache) const override
322 {
323 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
324 _ip_data, &IpData::eps, cache);
325 }
326
327private:
347 const double t, double const dt, Eigen::VectorXd const& local_x,
348 std::vector<double>& local_b_data, std::vector<double>& local_Jac_data);
349
374 const double t, double const dt, Eigen::VectorXd const& local_x,
375 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_b_data,
376 std::vector<double>& local_Jac_data);
377
378 unsigned getNumberOfIntegrationPoints() const override;
379
380 int getMaterialID() const override;
381
383 DisplacementDim>::MaterialStateVariables const&
384 getMaterialStateVariablesAt(unsigned integration_point) const override;
385
386private:
388
391 using IpData =
393 ShapeMatricesTypePressure, DisplacementDim,
394 ShapeFunctionDisplacement::NPOINTS>;
395 std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
396
401 typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType>
403
404 static const int pressure_index = 0;
405 static const int pressure_size = ShapeFunctionPressure::NPOINTS;
406 static const int displacement_index = ShapeFunctionPressure::NPOINTS;
407 static const int displacement_size =
408 ShapeFunctionDisplacement::NPOINTS * DisplacementDim;
409};
410
411} // namespace HydroMechanics
412} // namespace ProcessLib
413
#define OGS_FATAL(...)
Definition Error.h:26
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > mechanical_strain
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > stress
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:89
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
NumLib::GenericIntegrationMethod const & _integration_method
void assemble(double const, double const, std::vector< double > const &, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
HydroMechanicsLocalAssembler(HydroMechanicsLocalAssembler &&)=delete
Eigen::Matrix< double, KelvinVectorSize, 1 > SymmetricTensor
void postNonLinearSolverConcrete(std::vector< double > const &local_x, std::vector< double > const &local_x_prev, double const t, double const dt, int const process_id) 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 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 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)
std::size_t setIPDataInitialConditions(std::string_view const name, double const *values, int const integration_order) override
Returns number of read integration points.
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
void computeSecondaryVariableConcrete(double const t, double const dt, Eigen::VectorXd const &local_xs, Eigen::VectorXd const &local_x_prev) override
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt(unsigned integration_point) const override
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)
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
void setInitialConditionsConcrete(Eigen::VectorXd const local_x, double const t, int const process_id) override
HydroMechanicsLocalAssembler(HydroMechanicsLocalAssembler const &)=delete
std::vector< double > const & getIntPtEpsilon(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const 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
HydroMechanicsProcessData< DisplacementDim > & _process_data
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
BMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > BMatricesType
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
std::vector< double > const & getIntPtSigma(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Eigen::Matrix< double, Eigen::MatrixBase< Derived >::RowsAtCompileTime, 1 > symmetricTensorToKelvinVector(Eigen::MatrixBase< Derived > const &v)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType
constexpr Eigen::CwiseNullaryOp< EigenBlockMatrixViewFunctor< D, M >, typename EigenBlockMatrixViewFunctor< D, M >::Matrix > eigenBlockMatrixView(const Eigen::MatrixBase< M > &matrix)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
BMatricesType::KelvinMatrixType updateConstitutiveRelation(MPL::VariableArray const &variable_array, double const t, ParameterLib::SpatialPosition const &x_position, double const dt, DisplacementVectorType const &, double const T)
ShapeMatricesTypeDisplacement::GlobalDimNodalMatrixType dNdx_u
std::unique_ptr< typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables > material_state_variables
ShapeMatricesTypeDisplacement::NodalRowVectorType N_u
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > computeElasticTangentStiffness(double const t, ParameterLib::SpatialPosition const &x_position, double const dt, double const temperature)
ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p
ShapeMatricesTypePressure::NodalRowVectorType N_p
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material
IntegrationPointData(MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material)
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u