OGS
HydroMechanicsFEM-impl.h
Go to the documentation of this file.
1
11
12#pragma once
13
14#include <Eigen/Eigenvalues>
15
16#include "HydroMechanicsFEM.h"
17#include "HydroMechanicsProcessData.h"
18#include "MaterialLib/MPL/Medium.h"
19#include "MaterialLib/MPL/Property.h"
20#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
21#include "MaterialLib/MPL/Utils/GetFluidDensityAndViscosity.h"
22#include "MaterialLib/MPL/Utils/GetSymmetricTensor.h"
23#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
24#include "MathLib/KelvinVector.h"
25#include "NumLib/Fem/Interpolation.h"
26#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
27
28namespace ProcessLib
29{
30namespace HydroMechanics
31{
32namespace MPL = MaterialPropertyLib;
33
34template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
35 int DisplacementDim>
36HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
37 DisplacementDim>::
38 HydroMechanicsLocalAssembler(
39 MeshLib::Element const& e,
40 std::size_t const /*local_matrix_size*/,
41 NumLib::GenericIntegrationMethod const& integration_method,
42 bool const is_axially_symmetric,
43 HydroMechanicsProcessData<DisplacementDim>& process_data)
44 : _process_data(process_data),
45 _integration_method(integration_method),
46 _element(e),
47 _is_axially_symmetric(is_axially_symmetric)
48{
49 unsigned const n_integration_points =
50 _integration_method.getNumberOfPoints();
51
52 _ip_data.reserve(n_integration_points);
53 _secondary_data.N_u.resize(n_integration_points);
54
55 auto const shape_matrices_u =
56 NumLib::initShapeMatrices<ShapeFunctionDisplacement,
57 ShapeMatricesTypeDisplacement,
58 DisplacementDim>(e, is_axially_symmetric,
59 _integration_method);
60
61 auto const shape_matrices_p =
62 NumLib::initShapeMatrices<ShapeFunctionPressure,
63 ShapeMatricesTypePressure, DisplacementDim>(
64 e, is_axially_symmetric, _integration_method);
65
66 auto const& solid_material =
67 MaterialLib::Solids::selectSolidConstitutiveRelation(
68 _process_data.solid_materials, _process_data.material_ids,
69 e.getID());
70
71 for (unsigned ip = 0; ip < n_integration_points; ip++)
72 {
73 _ip_data.emplace_back(solid_material);
74 auto& ip_data = _ip_data[ip];
75 auto const& sm_u = shape_matrices_u[ip];
76 ip_data.integration_weight =
77 _integration_method.getWeightedPoint(ip).getWeight() *
78 sm_u.integralMeasure * sm_u.detJ;
79
80 // Initialize current time step values
81 static const int kelvin_vector_size =
82 MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
83 ip_data.sigma_eff.setZero(kelvin_vector_size);
84 ip_data.eps.setZero(kelvin_vector_size);
85
86 // Previous time step values are not initialized and are set later.
87 ip_data.eps_prev.resize(kelvin_vector_size);
88 ip_data.sigma_eff_prev.resize(kelvin_vector_size);
89
90 ip_data.N_u = sm_u.N;
91 ip_data.dNdx_u = sm_u.dNdx;
92
93 ip_data.N_p = shape_matrices_p[ip].N;
94 ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
95
96 _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
97 }
98}
99
100template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
101 int DisplacementDim>
102void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
103 ShapeFunctionPressure, DisplacementDim>::
104 assembleWithJacobian(double const t, double const dt,
105 std::vector<double> const& local_x,
106 std::vector<double> const& local_x_prev,
107 std::vector<double>& local_rhs_data,
108 std::vector<double>& local_Jac_data)
109{
110 assert(local_x.size() == pressure_size + displacement_size);
111
112 auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
113 pressure_size> const>(local_x.data() + pressure_index, pressure_size);
114
115 auto u =
116 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
117 displacement_size> const>(local_x.data() + displacement_index,
118 displacement_size);
119
120 auto p_prev =
121 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
122 pressure_size> const>(local_x_prev.data() + pressure_index,
123 pressure_size);
124 auto u_prev =
125 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
126 displacement_size> const>(local_x_prev.data() + displacement_index,
127 displacement_size);
128
129 auto local_Jac = MathLib::createZeroedMatrix<
130 typename ShapeMatricesTypeDisplacement::template MatrixType<
131 displacement_size + pressure_size,
132 displacement_size + pressure_size>>(
133 local_Jac_data, displacement_size + pressure_size,
134 displacement_size + pressure_size);
135
136 auto local_rhs = MathLib::createZeroedVector<
137 typename ShapeMatricesTypeDisplacement::template VectorType<
138 displacement_size + pressure_size>>(
139 local_rhs_data, displacement_size + pressure_size);
140
141 typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
142 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
143 pressure_size);
144
145 typename ShapeMatricesTypePressure::NodalMatrixType storage_p =
146 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
147 pressure_size);
148
149 typename ShapeMatricesTypePressure::NodalMatrixType add_p_derivative =
150 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
151 pressure_size);
152
153 typename ShapeMatricesTypeDisplacement::template MatrixType<
154 displacement_size, pressure_size>
155 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
156 displacement_size, pressure_size>::Zero(displacement_size,
157 pressure_size);
158
159 typename ShapeMatricesTypeDisplacement::template MatrixType<
160 pressure_size, displacement_size>
161 Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
162 pressure_size, displacement_size>::Zero(pressure_size,
163 displacement_size);
164
165 typename ShapeMatricesTypeDisplacement::template MatrixType<
166 pressure_size, displacement_size>
167 Kpu_k = ShapeMatricesTypeDisplacement::template MatrixType<
168 pressure_size, displacement_size>::Zero(pressure_size,
169 displacement_size);
170
171 auto const& solid_material =
172 MaterialLib::Solids::selectSolidConstitutiveRelation(
173 _process_data.solid_materials, _process_data.material_ids,
174 _element.getID());
175
176 ParameterLib::SpatialPosition x_position;
177 x_position.setElementID(_element.getID());
178
179 unsigned const n_integration_points =
180 _integration_method.getNumberOfPoints();
181
182 auto const& b = _process_data.specific_body_force;
183 auto const& medium = _process_data.media_map.getMedium(_element.getID());
184 auto const& solid = medium->phase("Solid");
185 auto const& fluid = fluidPhase(*medium);
186 MPL::VariableArray vars;
187 auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
188 : vars.liquid_phase_pressure;
189
190 auto const T_ref =
191 medium->property(MPL::PropertyType::reference_temperature)
192 .template value<double>(vars, x_position, t, dt);
193 vars.temperature = T_ref;
194
195 auto const& identity2 = Invariants::identity2;
196
197 for (unsigned ip = 0; ip < n_integration_points; ip++)
198 {
199 auto const& w = _ip_data[ip].integration_weight;
200
201 auto const& N_u = _ip_data[ip].N_u;
202 auto const& dNdx_u = _ip_data[ip].dNdx_u;
203
204 auto const& N_p = _ip_data[ip].N_p;
205 auto const& dNdx_p = _ip_data[ip].dNdx_p;
206
207 x_position = {
208 std::nullopt, _element.getID(),
209 MathLib::Point3d(
210 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
211 ShapeMatricesTypeDisplacement>(
212 _element, N_u))};
213
214 auto const x_coord = x_position.getCoordinates().value()[0];
215
216 auto const B =
217 LinearBMatrix::computeBMatrix<DisplacementDim,
218 ShapeFunctionDisplacement::NPOINTS,
219 typename BMatricesType::BMatrixType>(
220 dNdx_u, N_u, x_coord, _is_axially_symmetric);
221
222 auto& eps = _ip_data[ip].eps;
223 eps.noalias() = B * u;
224 auto const& sigma_eff = _ip_data[ip].sigma_eff;
225
226 double const p_int_pt = N_p.dot(p);
227
228 phase_pressure = p_int_pt;
229
230 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
231 t, x_position, dt, T_ref);
232 auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
233
234 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
235 .template value<double>(vars, x_position, t, dt);
236
237 auto const rho_sr =
238 solid.property(MPL::PropertyType::density)
239 .template value<double>(vars, x_position, t, dt);
240 auto const porosity =
241 medium->property(MPL::PropertyType::porosity)
242 .template value<double>(vars, x_position, t, dt);
243
244 auto const [rho_fr, mu] =
245 MaterialPropertyLib::getFluidDensityAndViscosity(t, dt, x_position,
246 fluid, vars);
247
248 auto const beta_p =
249 fluid.property(MPL::PropertyType::density)
250 .template dValue<double>(vars, _process_data.phase_variable,
251 x_position, t, dt) /
252 rho_fr;
253
254 // Set mechanical variables for the intrinsic permeability model
255 // For stress dependent permeability.
256 {
257 auto const sigma_total =
258 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
259
260 vars.total_stress.emplace<SymmetricTensor>(
261 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
262 sigma_total));
263 }
264 // For strain dependent permeability
265 vars.volumetric_strain = Invariants::trace(eps);
266 vars.equivalent_plastic_strain =
267 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
268 vars.mechanical_strain
269 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
270 eps);
271
272 auto const K = MPL::formEigenTensor<DisplacementDim>(
273 medium->property(MPL::PropertyType::permeability)
274 .value(vars, x_position, t, dt));
275 auto const dkde = MPL::formEigenTensor<
276 MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim)>(
277 (*medium)[MPL::PropertyType::permeability].dValue(
278 vars, MPL::Variable::mechanical_strain, x_position, t, dt));
279
280 auto const K_over_mu = K / mu;
281
282 auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
283 dt, u, T_ref);
284
285 //
286 // displacement equation, displacement part
287 //
288 local_Jac
289 .template block<displacement_size, displacement_size>(
290 displacement_index, displacement_index)
291 .noalias() += B.transpose() * C * B * w;
292
293 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
294 local_rhs.template segment<displacement_size>(displacement_index)
295 .noalias() -=
296 (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
297
298 //
299 // displacement equation, pressure part
300 //
301 Kup.noalias() += B.transpose() * alpha * identity2 * N_p * w;
302
303 //
304 // pressure equation, pressure part.
305 //
306 laplace_p.noalias() +=
307 rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
308
309 storage_p.noalias() +=
310 rho_fr * N_p.transpose() * N_p * w *
311 ((alpha - porosity) * (1.0 - alpha) / K_S + porosity * beta_p);
312
313 // density dependence on pressure evaluated for Darcy-term,
314 // for laplace and storage terms this dependence is neglected
315 add_p_derivative.noalias() += rho_fr * beta_p * dNdx_p.transpose() *
316 K_over_mu *
317 (dNdx_p * p - 2.0 * rho_fr * b) * N_p * w;
318
319 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
320 dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
321
322 //
323 // pressure equation, displacement part.
324 //
325 Kpu.noalias() +=
326 rho_fr * alpha * N_p.transpose() * identity2.transpose() * B * w;
327
328 Kpu_k.noalias() +=
329 dNdx_p.transpose() *
330 MathLib::KelvinVector::liftVectorToKelvin<DisplacementDim>(
331 dNdx_p * p - rho_fr * b) *
332 dkde * B * rho_fr / mu * w;
333 }
334 // displacement equation, pressure part
335 local_Jac
336 .template block<displacement_size, pressure_size>(displacement_index,
337 pressure_index)
338 .noalias() = -Kup;
339
340 if (_process_data.mass_lumping)
341 {
342 storage_p = storage_p.colwise().sum().eval().asDiagonal();
343
344 if constexpr (pressure_size == displacement_size)
345 {
346 Kpu = Kpu.colwise().sum().eval().asDiagonal();
347 Kpu_k = Kpu_k.colwise().sum().eval().asDiagonal();
348 }
349 }
350
351 // pressure equation, pressure part.
352 local_Jac
353 .template block<pressure_size, pressure_size>(pressure_index,
354 pressure_index)
355 .noalias() += laplace_p + storage_p / dt + add_p_derivative;
356
357 // pressure equation, displacement part.
358 local_Jac
359 .template block<pressure_size, displacement_size>(pressure_index,
360 displacement_index)
361 .noalias() += Kpu / dt + Kpu_k;
362
363 // pressure equation
364 local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
365 laplace_p * p + storage_p * (p - p_prev) / dt + Kpu * (u - u_prev) / dt;
366
367 // displacement equation
368 local_rhs.template segment<displacement_size>(displacement_index)
369 .noalias() += Kup * p;
370}
371
372template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
373 int DisplacementDim>
374std::vector<double> const& HydroMechanicsLocalAssembler<
375 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
376 getIntPtDarcyVelocity(
377 const double t,
378 std::vector<GlobalVector*> const& x,
379 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
380 std::vector<double>& cache) const
381{
382 int const hydraulic_process_id = _process_data.hydraulic_process_id;
383 auto const indices =
384 NumLib::getIndices(_element.getID(), *dof_table[hydraulic_process_id]);
385 assert(!indices.empty());
386 auto const local_x = x[hydraulic_process_id]->get(indices);
387
388 unsigned const n_integration_points =
389 _integration_method.getNumberOfPoints();
390 cache.clear();
391 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
392 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
393 cache, DisplacementDim, n_integration_points);
394
395 auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
396 pressure_size> const>(local_x.data() + pressure_index, pressure_size);
397
398 ParameterLib::SpatialPosition x_position;
399 x_position.setElementID(_element.getID());
400
401 auto const& medium = _process_data.media_map.getMedium(_element.getID());
402 auto const& fluid = fluidPhase(*medium);
403 MPL::VariableArray vars;
404 auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
405 : vars.liquid_phase_pressure;
406
407 // TODO (naumov) Temporary value not used by current material models. Need
408 // extension of secondary variables interface.
409 double const dt = std::numeric_limits<double>::quiet_NaN();
410 vars.temperature =
411 medium->property(MPL::PropertyType::reference_temperature)
412 .template value<double>(vars, x_position, t, dt);
413
414 auto const& identity2 = Invariants::identity2;
415
416 for (unsigned ip = 0; ip < n_integration_points; ip++)
417 {
418 x_position = {
419 std::nullopt, _element.getID(),
420 MathLib::Point3d(
421 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
422 ShapeMatricesTypeDisplacement>(
423 _element, _ip_data[ip].N_u))};
424
425 double const p_int_pt = _ip_data[ip].N_p.dot(p);
426
427 phase_pressure = p_int_pt;
428
429 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
430 .template value<double>(vars, x_position, t, dt);
431
432 // Set mechanical variables for the intrinsic permeability model
433 // For stress dependent permeability.
434 auto const sigma_total =
435 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
436 vars.total_stress.emplace<SymmetricTensor>(
437 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_total));
438
439 // For strain dependent permeability
440 vars.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
441 vars.equivalent_plastic_strain =
442 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
443 vars.mechanical_strain
444 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
445 _ip_data[ip].eps);
446
447 auto const K = MPL::formEigenTensor<DisplacementDim>(
448 medium->property(MPL::PropertyType::permeability)
449 .value(vars, x_position, t, dt));
450
451 auto const [rho_fr, mu] =
452 MaterialPropertyLib::getFluidDensityAndViscosity(t, dt, x_position,
453 fluid, vars);
454
455 auto const K_over_mu = K / mu;
456
457 auto const& b = _process_data.specific_body_force;
458
459 // Compute the velocity
460 auto const& dNdx_p = _ip_data[ip].dNdx_p;
461 cache_matrix.col(ip).noalias() =
462 -K_over_mu * dNdx_p * p + K_over_mu * rho_fr * b;
463 }
464
465 return cache;
466}
467
468template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
469 int DisplacementDim>
470void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
471 ShapeFunctionPressure, DisplacementDim>::
472 assembleWithJacobianForPressureEquations(
473 const double t, double const dt, Eigen::VectorXd const& local_x,
474 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_b_data,
475 std::vector<double>& local_Jac_data)
476{
477 auto local_rhs =
478 MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
479 template VectorType<pressure_size>>(
480 local_b_data, pressure_size);
481
482 ParameterLib::SpatialPosition pos;
483 pos.setElementID(this->_element.getID());
484
485 auto const p = local_x.template segment<pressure_size>(pressure_index);
486
487 auto const p_prev =
488 local_x_prev.template segment<pressure_size>(pressure_index);
489
490 auto local_Jac = MathLib::createZeroedMatrix<
491 typename ShapeMatricesTypeDisplacement::template MatrixType<
492 pressure_size, pressure_size>>(local_Jac_data, pressure_size,
493 pressure_size);
494
495 typename ShapeMatricesTypePressure::NodalMatrixType storage =
496 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
497 pressure_size);
498
499 typename ShapeMatricesTypePressure::NodalMatrixType laplace =
500 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
501 pressure_size);
502
503 typename ShapeMatricesTypePressure::NodalMatrixType add_p_derivative =
504 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
505 pressure_size);
506
507 auto const& solid_material =
508 MaterialLib::Solids::selectSolidConstitutiveRelation(
509 _process_data.solid_materials, _process_data.material_ids,
510 _element.getID());
511
512 ParameterLib::SpatialPosition x_position;
513 x_position.setElementID(_element.getID());
514
515 auto const& medium = _process_data.media_map.getMedium(_element.getID());
516 auto const& fluid = fluidPhase(*medium);
517 MPL::VariableArray vars;
518 auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
519 : vars.liquid_phase_pressure;
520
521 auto const T_ref =
522 medium->property(MPL::PropertyType::reference_temperature)
523 .template value<double>(vars, x_position, t, dt);
524 vars.temperature = T_ref;
525
526 auto const& identity2 = Invariants::identity2;
527
528 auto const staggered_scheme =
529 std::get<Staggered>(_process_data.coupling_scheme);
530 auto const fixed_stress_stabilization_parameter =
531 staggered_scheme.fixed_stress_stabilization_parameter;
532 auto const fixed_stress_over_time_step =
533 staggered_scheme.fixed_stress_over_time_step;
534
535 int const n_integration_points = _integration_method.getNumberOfPoints();
536 for (int ip = 0; ip < n_integration_points; ip++)
537 {
538 auto const& w = _ip_data[ip].integration_weight;
539
540 auto const& N_p = _ip_data[ip].N_p;
541 auto const& dNdx_p = _ip_data[ip].dNdx_p;
542
543 x_position = {
544 std::nullopt, _element.getID(),
545 MathLib::Point3d(
546 NumLib::interpolateCoordinates<ShapeFunctionPressure,
547 ShapeMatricesTypePressure>(
548 _element, _ip_data[ip].N_p))};
549
550 double const p_int_pt = N_p.dot(p);
551
552 phase_pressure = p_int_pt;
553
554 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
555 t, x_position, dt, T_ref);
556 auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
557
558 auto const alpha_b =
559 medium->property(MPL::PropertyType::biot_coefficient)
560 .template value<double>(vars, x_position, t, dt);
561
562 // Set mechanical variables for the intrinsic permeability model
563 // For stress dependent permeability.
564 auto const sigma_total =
565 (_ip_data[ip].sigma_eff - alpha_b * p_int_pt * identity2).eval();
566 vars.total_stress.emplace<SymmetricTensor>(
567 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_total));
568
569 // For strain dependent permeability
570 vars.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
571 vars.equivalent_plastic_strain =
572 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
573 vars.mechanical_strain
574 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
575 _ip_data[ip].eps);
576
577 auto const K = MPL::formEigenTensor<DisplacementDim>(
578 medium->property(MPL::PropertyType::permeability)
579 .value(vars, x_position, t, dt));
580 auto const porosity =
581 medium->property(MPL::PropertyType::porosity)
582 .template value<double>(vars, x_position, t, dt);
583
584 auto const [rho_fr, mu] =
585 MaterialPropertyLib::getFluidDensityAndViscosity(t, dt, x_position,
586 fluid, vars);
587
588 auto const beta_p =
589 fluid.property(MPL::PropertyType::density)
590 .template dValue<double>(vars, _process_data.phase_variable,
591 x_position, t, dt) /
592 rho_fr;
593
594 auto const K_over_mu = K / mu;
595
596 laplace.noalias() +=
597 rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
598
599 // Artificial compressibility from the fixed stress splitting:
600 auto const beta_FS =
601 fixed_stress_stabilization_parameter * alpha_b * alpha_b / K_S;
602
603 storage.noalias() += rho_fr * N_p.transpose() * N_p * w *
604 ((alpha_b - porosity) * (1.0 - alpha_b) / K_S +
605 porosity * beta_p + beta_FS);
606
607 auto const& b = _process_data.specific_body_force;
608
609 // bodyforce-driven Darcy flow
610 local_rhs.noalias() +=
611 dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
612
613 // density dependence on pressure evaluated for Darcy-term,
614 // for laplace and storage terms this dependence is neglected (as is
615 // done for monolithic too)
616 add_p_derivative.noalias() += rho_fr * beta_p * dNdx_p.transpose() *
617 K_over_mu *
618 (dNdx_p * p - 2.0 * rho_fr * b) * N_p * w;
619
620 if (!fixed_stress_over_time_step)
621 {
622 auto const& eps = _ip_data[ip].eps;
623 auto const& eps_prev = _ip_data[ip].eps_prev;
624 const double eps_v_dot =
625 (Invariants::trace(eps) - Invariants::trace(eps_prev)) / dt;
626
627 // Constant portion of strain rate term:
628 double const strain_rate_b =
629 alpha_b * eps_v_dot -
630 beta_FS * _ip_data[ip].strain_rate_variable;
631
632 local_rhs.noalias() -= strain_rate_b * rho_fr * N_p * w;
633 }
634 else
635 {
636 // Constant portion of strain rate term:
637 local_rhs.noalias() -=
638 alpha_b * _ip_data[ip].strain_rate_variable * rho_fr * N_p * w;
639 }
640 }
641 local_Jac.noalias() = laplace + storage / dt + add_p_derivative;
642
643 local_rhs.noalias() -= laplace * p + storage * (p - p_prev) / dt;
644}
645
646template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
647 int DisplacementDim>
648void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
649 ShapeFunctionPressure, DisplacementDim>::
650 assembleWithJacobianForDeformationEquations(
651 const double t, double const dt, Eigen::VectorXd const& local_x,
652 std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
653{
654 auto const p = local_x.template segment<pressure_size>(pressure_index);
655 auto const u =
656 local_x.template segment<displacement_size>(displacement_index);
657
658 auto local_Jac = MathLib::createZeroedMatrix<
659 typename ShapeMatricesTypeDisplacement::template MatrixType<
660 displacement_size, displacement_size>>(
661 local_Jac_data, displacement_size, displacement_size);
662
663 auto local_rhs =
664 MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
665 template VectorType<displacement_size>>(
666 local_b_data, displacement_size);
667
668 ParameterLib::SpatialPosition x_position;
669 x_position.setElementID(_element.getID());
670
671 auto const& medium = _process_data.media_map.getMedium(_element.getID());
672 auto const& solid = medium->phase("Solid");
673 auto const& fluid = fluidPhase(*medium);
674 MPL::VariableArray vars;
675 auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
676 : vars.liquid_phase_pressure;
677
678 auto const T_ref =
679 medium->property(MPL::PropertyType::reference_temperature)
680 .template value<double>(vars, x_position, t, dt);
681 vars.temperature = T_ref;
682
683 int const n_integration_points = _integration_method.getNumberOfPoints();
684 for (int ip = 0; ip < n_integration_points; ip++)
685 {
686 auto const& w = _ip_data[ip].integration_weight;
687
688 auto const& N_u = _ip_data[ip].N_u;
689 auto const& dNdx_u = _ip_data[ip].dNdx_u;
690
691 auto const& N_p = _ip_data[ip].N_p;
692
693 x_position = {
694 std::nullopt, _element.getID(),
695 MathLib::Point3d(
696 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
697 ShapeMatricesTypeDisplacement>(
698 _element, N_u))};
699
700 auto const x_coord = x_position.getCoordinates().value()[0];
701 auto const B =
702 LinearBMatrix::computeBMatrix<DisplacementDim,
703 ShapeFunctionDisplacement::NPOINTS,
704 typename BMatricesType::BMatrixType>(
705 dNdx_u, N_u, x_coord, _is_axially_symmetric);
706
707 auto& eps = _ip_data[ip].eps;
708 auto const& sigma_eff = _ip_data[ip].sigma_eff;
709
710 phase_pressure = N_p.dot(p);
711
712 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
713 .template value<double>(vars, x_position, t, dt);
714 auto const rho_sr =
715 solid.property(MPL::PropertyType::density)
716 .template value<double>(vars, x_position, t, dt);
717 auto const porosity =
718 medium->property(MPL::PropertyType::porosity)
719 .template value<double>(vars, x_position, t, dt);
720
721 auto const rho_fr = MaterialPropertyLib::getFluidDensity(
722 t, dt, x_position, fluid, vars);
723
724 auto const& b = _process_data.specific_body_force;
725 auto const& identity2 = MathLib::KelvinVector::Invariants<
726 MathLib::KelvinVector::kelvin_vector_dimensions(
727 DisplacementDim)>::identity2;
728
729 eps.noalias() = B * u;
730 vars.mechanical_strain
731 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
732 eps);
733
734 auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
735 dt, u, T_ref);
736
737 local_Jac.noalias() += B.transpose() * C * B * w;
738
739 double p_at_xi = 0.;
740 NumLib::shapeFunctionInterpolate(p, N_p, p_at_xi);
741
742 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
743 local_rhs.noalias() -=
744 (B.transpose() * (sigma_eff - alpha * identity2 * p_at_xi) -
745 N_u_op(N_u).transpose() * rho * b) *
746 w;
747 }
748}
749
750template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
751 int DisplacementDim>
752void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
753 ShapeFunctionPressure, DisplacementDim>::
754 assembleWithJacobianForStaggeredScheme(const double t, double const dt,
755 Eigen::VectorXd const& local_x,
756 Eigen::VectorXd const& local_x_prev,
757 int const process_id,
758 std::vector<double>& local_b_data,
759 std::vector<double>& local_Jac_data)
760{
761 // For the equations with pressure
762 if (process_id == _process_data.hydraulic_process_id)
763 {
764 assembleWithJacobianForPressureEquations(t, dt, local_x, local_x_prev,
765 local_b_data, local_Jac_data);
766 return;
767 }
768
769 // For the equations with deformation
770 assembleWithJacobianForDeformationEquations(t, dt, local_x, local_b_data,
771 local_Jac_data);
772}
773
774template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
775 int DisplacementDim>
776void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
777 ShapeFunctionPressure, DisplacementDim>::
778 setInitialConditionsConcrete(Eigen::VectorXd const local_x,
779 double const t,
780 int const /*process_id*/)
781{
782 ParameterLib::SpatialPosition x_position;
783 x_position.setElementID(_element.getID());
784
785 auto const& medium = _process_data.media_map.getMedium(_element.getID());
786
787 auto const p = local_x.template segment<pressure_size>(pressure_index);
788 auto const u =
789 local_x.template segment<displacement_size>(displacement_index);
790
791 auto const& identity2 = Invariants::identity2;
792 const double dt = 0.0;
793
794 MPL::VariableArray vars;
795
796 int const n_integration_points = _integration_method.getNumberOfPoints();
797 for (int ip = 0; ip < n_integration_points; ip++)
798 {
799 auto const& N_u = _ip_data[ip].N_u;
800 auto const& dNdx_u = _ip_data[ip].dNdx_u;
801
802 x_position = {
803 std::nullopt, _element.getID(),
804 MathLib::Point3d(
805 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
806 ShapeMatricesTypeDisplacement>(
807 _element, N_u))};
808
809 auto const x_coord = x_position.getCoordinates().value()[0];
810 auto const B =
811 LinearBMatrix::computeBMatrix<DisplacementDim,
812 ShapeFunctionDisplacement::NPOINTS,
813 typename BMatricesType::BMatrixType>(
814 dNdx_u, N_u, x_coord, _is_axially_symmetric);
815
816 auto& eps = _ip_data[ip].eps;
817 eps.noalias() = B * u;
818 vars.mechanical_strain
819 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
820 eps);
821
822 if (_process_data.initial_stress.isTotalStress())
823 {
824 auto const& N_p = _ip_data[ip].N_p;
825 auto const alpha_b =
826 medium->property(MPL::PropertyType::biot_coefficient)
827 .template value<double>(vars, x_position, t, dt);
828
829 auto& sigma_eff = _ip_data[ip].sigma_eff;
830 sigma_eff.noalias() += alpha_b * N_p.dot(p) * identity2;
831 _ip_data[ip].sigma_eff_prev.noalias() = sigma_eff;
832 }
833 }
834}
835
836template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
837 int DisplacementDim>
838void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
839 ShapeFunctionPressure, DisplacementDim>::
840 postNonLinearSolverConcrete(Eigen::VectorXd const& local_x,
841 Eigen::VectorXd const& local_x_prev,
842 double const t, double const dt,
843 int const process_id)
844{
845 // Note: local_x and local_x_prev only contain the solutions of current
846 // process in the staggered scheme. This has to be changed according to the
847 // same two arguments in postTimestepConcrete.
848
849 int const n_integration_points = _integration_method.getNumberOfPoints();
850
851 auto const staggered_scheme_ptr =
852 std::get_if<Staggered>(&_process_data.coupling_scheme);
853
854 if (staggered_scheme_ptr &&
855 process_id == _process_data.hydraulic_process_id)
856 {
857 if (!staggered_scheme_ptr->fixed_stress_over_time_step)
858 {
859 auto const p =
860 local_x.template segment<pressure_size>(pressure_index);
861 auto const p_prev =
862 local_x_prev.template segment<pressure_size>(pressure_index);
863
864 for (int ip = 0; ip < n_integration_points; ip++)
865 {
866 auto& ip_data = _ip_data[ip];
867
868 auto const& N_p = ip_data.N_p;
869
870 ip_data.strain_rate_variable = N_p.dot(p - p_prev) / dt;
871 }
872 }
873 }
874
875 if (!staggered_scheme_ptr ||
876 process_id == _process_data.mechanics_related_process_id)
877 {
878 ParameterLib::SpatialPosition x_position;
879 x_position.setElementID(_element.getID());
880
881 auto const& medium =
882 _process_data.media_map.getMedium(_element.getID());
883
884 auto const T_ref =
885 medium->property(MPL::PropertyType::reference_temperature)
886 .template value<double>(MPL::EmptyVariableArray, x_position, t,
887 dt);
888
889 auto const u =
890 local_x.template segment<displacement_size>(displacement_index);
891
892 MPL::VariableArray vars;
893 vars.temperature = T_ref;
894
895 for (int ip = 0; ip < n_integration_points; ip++)
896 {
897 auto const& N_u = _ip_data[ip].N_u;
898 auto const& dNdx_u = _ip_data[ip].dNdx_u;
899
900 x_position = {std::nullopt, _element.getID(),
901 MathLib::Point3d(NumLib::interpolateCoordinates<
902 ShapeFunctionDisplacement,
903 ShapeMatricesTypeDisplacement>(
904 _element, N_u))};
905
906 auto const x_coord = x_position.getCoordinates().value()[0];
907 auto const B = LinearBMatrix::computeBMatrix<
908 DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
909 typename BMatricesType::BMatrixType>(dNdx_u, N_u, x_coord,
910 _is_axially_symmetric);
911
912 auto& eps = _ip_data[ip].eps;
913 eps.noalias() = B * u;
914 vars.mechanical_strain.emplace<
915 MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(eps);
916
917 _ip_data[ip].updateConstitutiveRelation(vars, t, x_position, dt, u,
918 T_ref);
919 }
920 }
921}
922
923template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
924 int DisplacementDim>
925void HydroMechanicsLocalAssembler<
926 ShapeFunctionDisplacement, ShapeFunctionPressure,
927 DisplacementDim>::postTimestepConcrete(Eigen::VectorXd const& local_x,
928 Eigen::VectorXd const& local_x_prev,
929 double const t, double const dt,
930 int const process_id)
931{
932 auto const staggered_scheme_ptr =
933 std::get_if<Staggered>(&_process_data.coupling_scheme);
934
935 if (staggered_scheme_ptr &&
936 process_id == _process_data.hydraulic_process_id)
937 {
938 if (staggered_scheme_ptr->fixed_stress_over_time_step)
939 {
940 auto const fixed_stress_stabilization_parameter =
941 staggered_scheme_ptr->fixed_stress_stabilization_parameter;
942
943 auto const p =
944 local_x.template segment<pressure_size>(pressure_index);
945 auto const p_prev =
946 local_x_prev.template segment<pressure_size>(pressure_index);
947
948 ParameterLib::SpatialPosition x_position;
949 x_position.setElementID(_element.getID());
950
951 auto const& solid_material =
952 MaterialLib::Solids::selectSolidConstitutiveRelation(
953 _process_data.solid_materials, _process_data.material_ids,
954 _element.getID());
955
956 auto const& medium =
957 _process_data.media_map.getMedium(_element.getID());
958 MPL::VariableArray vars;
959
960 auto const T_ref =
961 medium->property(MPL::PropertyType::reference_temperature)
962 .template value<double>(vars, x_position, t, dt);
963 vars.temperature = T_ref;
964
965 int const n_integration_points =
966 _integration_method.getNumberOfPoints();
967 for (int ip = 0; ip < n_integration_points; ip++)
968 {
969 auto& ip_data = _ip_data[ip];
970
971 auto const& N_p = ip_data.N_p;
972
973 x_position = {
974 std::nullopt, _element.getID(),
975 MathLib::Point3d(
976 NumLib::interpolateCoordinates<
977 ShapeFunctionPressure, ShapeMatricesTypePressure>(
978 _element, N_p))};
979
980 auto const& eps = ip_data.eps;
981 auto const& eps_prev = ip_data.eps_prev;
982 const double eps_v_dot =
983 (Invariants::trace(eps) - Invariants::trace(eps_prev)) / dt;
984
985 auto const C_el = ip_data.computeElasticTangentStiffness(
986 t, x_position, dt, T_ref);
987 auto const K_S =
988 solid_material.getBulkModulus(t, x_position, &C_el);
989
990 auto const alpha_b =
991 medium->property(MPL::PropertyType::biot_coefficient)
992 .template value<double>(vars, x_position, t, dt);
993
994 ip_data.strain_rate_variable =
995 eps_v_dot - fixed_stress_stabilization_parameter * alpha_b *
996 N_p.dot(p - p_prev) / dt / K_S;
997 }
998 }
999 }
1000
1001 unsigned const n_integration_points =
1002 _integration_method.getNumberOfPoints();
1003
1004 for (unsigned ip = 0; ip < n_integration_points; ip++)
1005 {
1006 _ip_data[ip].pushBackState();
1007 }
1008}
1009
1010template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1011 int DisplacementDim>
1012std::size_t HydroMechanicsLocalAssembler<
1013 ShapeFunctionDisplacement, ShapeFunctionPressure,
1014 DisplacementDim>::setIPDataInitialConditions(std::string_view const name,
1015 double const* values,
1016 int const integration_order)
1017{
1018 if (integration_order !=
1019 static_cast<int>(_integration_method.getIntegrationOrder()))
1020 {
1021 OGS_FATAL(
1022 "Setting integration point initial conditions; The integration "
1023 "order of the local assembler for element {:d} is different from "
1024 "the integration order in the initial condition.",
1025 _element.getID());
1026 }
1027
1028 if (name == "sigma")
1029 {
1030 if (_process_data.initial_stress.value != nullptr)
1031 {
1032 OGS_FATAL(
1033 "Setting initial conditions for stress from integration "
1034 "point data and from a parameter '{:s}' is not possible "
1035 "simultaneously.",
1036 _process_data.initial_stress.value->name);
1037 }
1038
1039 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
1040 values, _ip_data, &IpData::sigma_eff);
1041 }
1042
1043 if (name == "epsilon")
1044 {
1045 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
1046 values, _ip_data, &IpData::eps);
1047 }
1048
1049 if (name == "strain_rate_variable")
1050 {
1051 return ProcessLib::setIntegrationPointScalarData(
1052 values, _ip_data, &IpData::strain_rate_variable);
1053 }
1054
1055 return 0;
1056}
1057
1058template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1059 int DisplacementDim>
1060std::vector<double>
1061HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
1067
1068template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1069 int DisplacementDim>
1070std::vector<double>
1071HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
1072 DisplacementDim>::getEpsilon() const
1073{
1074 auto const kelvin_vector_size =
1075 MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
1076 unsigned const n_integration_points =
1077 _integration_method.getNumberOfPoints();
1078
1079 std::vector<double> ip_epsilon_values;
1080 auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
1081 double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
1082 ip_epsilon_values, n_integration_points, kelvin_vector_size);
1083
1084 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1085 {
1086 auto const& eps = _ip_data[ip].eps;
1087 cache_mat.row(ip) =
1088 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(eps);
1089 }
1090
1091 return ip_epsilon_values;
1092}
1093
1094template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1095 int DisplacementDim>
1096std::vector<double>
1097HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
1098 DisplacementDim>::getStrainRateVariable() const
1099{
1100 unsigned const n_integration_points =
1101 _integration_method.getNumberOfPoints();
1102
1103 std::vector<double> ip_strain_rate_variables(n_integration_points);
1104
1105 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1106 {
1107 ip_strain_rate_variables[ip] = _ip_data[ip].strain_rate_variable;
1108 }
1109
1110 return ip_strain_rate_variables;
1111}
1112
1113template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1114 int DisplacementDim>
1115void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
1116 ShapeFunctionPressure, DisplacementDim>::
1117 computeSecondaryVariableConcrete(double const t, double const dt,
1118 Eigen::VectorXd const& local_x,
1119 Eigen::VectorXd const& /*local_x_prev*/)
1120{
1121 auto const p = local_x.template segment<pressure_size>(pressure_index);
1122
1123 NumLib::interpolateToHigherOrderNodes<
1124 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1125 DisplacementDim>(_element, _is_axially_symmetric, p,
1126 *_process_data.pressure_interpolated);
1127
1128 int const elem_id = _element.getID();
1129
1130 unsigned const n_integration_points =
1131 _integration_method.getNumberOfPoints();
1132
1133 auto const& medium = _process_data.media_map.getMedium(elem_id);
1134 MPL::VariableArray vars;
1135 auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
1136 : vars.liquid_phase_pressure;
1137
1138 SymmetricTensor k_sum = SymmetricTensor::Zero(KelvinVectorSize);
1139 auto sigma_eff_sum = MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
1140 Eigen::Matrix<double, 3, 3>::Zero());
1141
1142 auto const& identity2 = Invariants::identity2;
1143
1144 for (unsigned ip = 0; ip < n_integration_points; ip++)
1145 {
1146 auto const& eps = _ip_data[ip].eps;
1147 sigma_eff_sum += _ip_data[ip].sigma_eff;
1148
1149 ParameterLib::SpatialPosition x_position = {
1150 std::nullopt, _element.getID(),
1151 MathLib::Point3d(
1152 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
1153 ShapeMatricesTypeDisplacement>(
1154 _element, _ip_data[ip].N_u))};
1155
1156 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
1157 .template value<double>(vars, x_position, t, dt);
1158 double const p_int_pt = _ip_data[ip].N_p.dot(p);
1159
1160 phase_pressure = p_int_pt;
1161
1162 // Set mechanical variables for the intrinsic permeability model
1163 // For stress dependent permeability.
1164 {
1165 auto const sigma_total =
1166 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
1167 vars.total_stress.emplace<SymmetricTensor>(
1168 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
1169 sigma_total));
1170 }
1171 // For strain dependent permeability
1172 vars.volumetric_strain = Invariants::trace(eps);
1173 vars.equivalent_plastic_strain =
1174 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
1175 vars.mechanical_strain
1176 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
1177 eps);
1178
1179 k_sum += MPL::getSymmetricTensor<DisplacementDim>(
1180 medium->property(MPL::PropertyType::permeability)
1181 .value(vars, x_position, t, dt));
1182 }
1183
1184 Eigen::Map<Eigen::VectorXd>(
1185 &(*_process_data.permeability)[elem_id * KelvinVectorSize],
1186 KelvinVectorSize) = k_sum / n_integration_points;
1187
1188 Eigen::Matrix<double, 3, 3, 0, 3, 3> const sigma_avg =
1189 MathLib::KelvinVector::kelvinVectorToTensor(sigma_eff_sum) /
1190 n_integration_points;
1191
1192 Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> e_s(sigma_avg);
1193
1194 Eigen::Map<Eigen::Vector3d>(
1195 &(*_process_data.principal_stress_values)[elem_id * 3], 3) =
1196 e_s.eigenvalues();
1197
1198 auto eigen_vectors = e_s.eigenvectors();
1199
1200 for (auto i = 0; i < 3; i++)
1201 {
1202 Eigen::Map<Eigen::Vector3d>(
1203 &(*_process_data.principal_stress_vector[i])[elem_id * 3], 3) =
1204 eigen_vectors.col(i);
1205 }
1206}
1207
1208template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1209 int DisplacementDim>
1210unsigned HydroMechanicsLocalAssembler<
1211 ShapeFunctionDisplacement, ShapeFunctionPressure,
1212 DisplacementDim>::getNumberOfIntegrationPoints() const
1213{
1214 return _integration_method.getNumberOfPoints();
1215}
1216
1217template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1218 int DisplacementDim>
1219int HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
1220 ShapeFunctionPressure,
1221 DisplacementDim>::getMaterialID() const
1222{
1223 return _process_data.material_ids == nullptr
1224 ? 0
1225 : (*_process_data.material_ids)[_element.getID()];
1226}
1227
1228template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1229 int DisplacementDim>
1230typename MaterialLib::Solids::MechanicsBase<
1231 DisplacementDim>::MaterialStateVariables const&
1232HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
1233 DisplacementDim>::
1234 getMaterialStateVariablesAt(unsigned integration_point) const
1235{
1236 return *_ip_data[integration_point].material_state_variables;
1237}
1238} // namespace HydroMechanics
1239} // namespace ProcessLib
OGS_FATAL
#define OGS_FATAL(...)
Definition Error.h:26
MeshLib::Element::getID
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:91
ParameterLib::SpatialPosition::setElementID
void setElementID(std::size_t element_id)
Definition SpatialPosition.h:76
ParameterLib::SpatialPosition::getCoordinates
std::optional< MathLib::Point3d > const getCoordinates() const
Definition SpatialPosition.h:63
ProcessLib::BMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim >::BMatrixType
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition BMatrixPolicy.h:55
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_integration_method
NumLib::GenericIntegrationMethod const & _integration_method
Definition HydroMechanicsFEM.h:394
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::SymmetricTensor
Eigen::Matrix< double, KelvinVectorSize, 1 > SymmetricTensor
Definition HydroMechanicsFEM.h:181
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::postTimestepConcrete
void postTimestepConcrete(Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, int const process_id) override
Definition HydroMechanicsFEM-impl.h:927
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::assembleWithJacobianForPressureEquations
void assembleWithJacobianForPressureEquations(const double t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
Definition HydroMechanicsFEM-impl.h:472
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::ShapeMatricesTypeDisplacement
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
Definition HydroMechanicsFEM.h:170
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::ShapeMatricesTypePressure
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
Definition HydroMechanicsFEM.h:174
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::assembleWithJacobianForDeformationEquations
void assembleWithJacobianForDeformationEquations(const double t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
Definition HydroMechanicsFEM-impl.h:650
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_process_data
HydroMechanicsProcessData< DisplacementDim > & _process_data
Definition HydroMechanicsFEM.h:384
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_secondary_data
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
Definition HydroMechanicsFEM.h:399
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_ip_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
Definition HydroMechanicsFEM.h:392
MaterialLib::Solids::selectSolidConstitutiveRelation
auto & selectSolidConstitutiveRelation(SolidMaterialsMap const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
Definition SelectSolidConstitutiveRelation.h:36
MaterialPropertyLib::getFluidDensityAndViscosity
std::tuple< double, double > getFluidDensityAndViscosity(double const t, double const dt, ParameterLib::SpatialPosition const &pos, MaterialPropertyLib::Phase const &fluid_phase, MaterialPropertyLib::VariableArray &vars)
It computes fluid density and viscosity for single phase flow model.
Definition GetFluidDensityAndViscosity.cpp:42
MaterialPropertyLib::getFluidDensity
double getFluidDensity(double const t, double const dt, ParameterLib::SpatialPosition const &pos, Phase const &fluid_phase, VariableArray &vars)
It computes fluid density for single phase flow model.
Definition GetFluidDensityAndViscosity.cpp:21
MaterialPropertyLib::formEigenTensor
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Definition FormEigenTensor.cpp:126
MaterialPropertyLib::EmptyVariableArray
static const VariableArray EmptyVariableArray
Definition VariableType.h:208
MaterialPropertyLib::getSymmetricTensor
SymmetricTensor< GlobalDim > getSymmetricTensor(MaterialPropertyLib::PropertyDataType const &values)
Definition GetSymmetricTensor.cpp:107
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::kelvin_vector_dimensions
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition KelvinVector.h:24
MathLib::KelvinVector::KelvinVectorType
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition KelvinVector.h:47
MathLib::KelvinVector::kelvinVectorToTensor
Eigen::Matrix< double, 3, 3 > kelvinVectorToTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Definition KelvinVector.cpp:64
MathLib::KelvinVector::liftVectorToKelvin
Eigen::Matrix< double, DisplacementDim, kelvin_vector_dimensions(DisplacementDim)> liftVectorToKelvin(Eigen::Matrix< double, DisplacementDim, 1 > const &v)
MathLib::KelvinVector::tensorToKelvin
KelvinVectorType< DisplacementDim > tensorToKelvin(Eigen::Matrix< double, 3, 3 > const &m)
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::detail::shapeFunctionInterpolate
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
Definition Interpolation.h:27
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
NumLib::getIndices
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
Definition DOFTableUtil.cpp:112
NumLib::initShapeMatrices
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
Definition InitShapeMatrices.h:54
NumLib::interpolateCoordinates
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
Definition InitShapeMatrices.h:102
ProcessLib::LinearBMatrix::computeBMatrix
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
Definition LinearBMatrix.h:44
ProcessLib::getIntegrationPointKelvinVectorData
std::vector< double > const & getIntegrationPointKelvinVectorData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
Definition SetOrGetIntegrationPointData.h:68
ProcessLib::setIntegrationPointKelvinVectorData
std::size_t setIntegrationPointKelvinVectorData(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType IpData::*const member)
Definition SetOrGetIntegrationPointData.h:133
ProcessLib::setIPDataInitialConditions
void setIPDataInitialConditions(std::vector< std::unique_ptr< MeshLib::IntegrationPointWriter > > const &_integration_point_writer, MeshLib::Properties const &mesh_properties, LocalAssemblersVector &local_assemblers)
Definition SetIPDataInitialConditions.h:36
ProcessLib::setIntegrationPointScalarData
std::size_t setIntegrationPointScalarData(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType IpData::*const member)
Definition SetOrGetIntegrationPointData.h:215
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