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