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