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