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