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