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