OGS
RichardsMechanicsFEM-impl.h
Go to the documentation of this file.
1
12#pragma once
13
14#include <Eigen/LU>
15#include <cassert>
16
17#include "ComputeMicroPorosity.h"
18#include "IntegrationPointData.h"
19#include "MaterialLib/MPL/Medium.h"
20#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
21#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
22#include "MathLib/EigenBlockMatrixView.h"
23#include "MathLib/KelvinVector.h"
24#include "NumLib/Function/Interpolation.h"
25#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
26#include "ProcessLib/Utils/TransposeInPlace.h"
27
28namespace ProcessLib
29{
30namespace RichardsMechanics
31{
32template <int DisplacementDim, typename IPData>
33void updateSwellingStressAndVolumetricStrain(
34 IPData& ip_data, MaterialPropertyLib::Medium const& medium,
35 MaterialPropertyLib::Phase const& solid_phase,
36 MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> const& C_el,
37 double const rho_LR, double const mu,
38 std::optional<MicroPorosityParameters> micro_porosity_parameters,
39 double const alpha, double const phi, double const p_cap_ip,
40 MPL::VariableArray& variables, MPL::VariableArray& variables_prev,
41 ParameterLib::SpatialPosition const& x_position, double const t,
42 double const dt)
43{
44 auto const& identity2 = MathLib::KelvinVector::Invariants<
45 MathLib::KelvinVector::kelvin_vector_dimensions(
46 DisplacementDim)>::identity2;
47
48 auto& sigma_sw = ip_data.sigma_sw;
49 auto const& sigma_sw_prev = ip_data.sigma_sw_prev;
50 if (!medium.hasProperty(MPL::PropertyType::saturation_micro))
51 {
52 // If there is swelling, compute it. Update volumetric strain rate,
53 // s.t. it corresponds to the mechanical part only.
54 sigma_sw = sigma_sw_prev;
55 if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
56 {
57 auto const sigma_sw_dot =
58 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
59 MPL::formEigenTensor<3>(
60 solid_phase[MPL::PropertyType::swelling_stress_rate]
61 .value(variables, variables_prev, x_position, t,
62 dt)));
63 sigma_sw += sigma_sw_dot * dt;
64
65 // !!! Misusing volumetric strain for mechanical volumetric
66 // strain just to update the transport porosity !!!
67 variables.volumetric_strain +=
68 identity2.transpose() * C_el.inverse() * sigma_sw;
69 variables_prev.volumetric_strain +=
70 identity2.transpose() * C_el.inverse() * sigma_sw_prev;
71 }
72 }
73
74 // TODO (naumov) saturation_micro must be always defined together with
75 // the micro_porosity_parameters.
76 if (medium.hasProperty(MPL::PropertyType::saturation_micro))
77 {
78 double const phi_M_prev = ip_data.transport_porosity_prev;
79 double const phi_prev = ip_data.porosity_prev;
80 double const phi_m_prev = phi_prev - phi_M_prev;
81 double const p_L_m_prev = ip_data.liquid_pressure_m_prev;
82
83 auto const S_L_m_prev = ip_data.saturation_m_prev;
84
85 auto const [delta_phi_m, delta_e_sw, delta_p_L_m, delta_sigma_sw] =
86 computeMicroPorosity<DisplacementDim>(
87 identity2.transpose() * C_el.inverse(), rho_LR, mu,
88 *micro_porosity_parameters, alpha, phi, -p_cap_ip, p_L_m_prev,
89 variables_prev, S_L_m_prev, phi_m_prev, x_position, t, dt,
90 medium.property(MPL::PropertyType::saturation_micro),
91 solid_phase.property(MPL::PropertyType::swelling_stress_rate));
92
93 auto& phi_M = ip_data.transport_porosity;
94 phi_M = phi - (phi_m_prev + delta_phi_m);
95 variables_prev.transport_porosity = phi_M_prev;
96 variables.transport_porosity = phi_M;
97
98 auto& p_L_m = ip_data.liquid_pressure_m;
99 p_L_m = p_L_m_prev + delta_p_L_m;
100 { // Update micro saturation.
101 MPL::VariableArray variables_prev;
102 variables_prev.capillary_pressure = -p_L_m_prev;
103 MPL::VariableArray variables;
104 variables.capillary_pressure = -p_L_m;
105
106 ip_data.saturation_m =
107 medium.property(MPL::PropertyType::saturation_micro)
108 .template value<double>(variables, x_position, t, dt);
109 }
110 sigma_sw = sigma_sw_prev + delta_sigma_sw;
111 }
112}
113
114template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
115 int DisplacementDim>
116RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
117 ShapeFunctionPressure, DisplacementDim>::
118 RichardsMechanicsLocalAssembler(
119 MeshLib::Element const& e,
120 std::size_t const /*local_matrix_size*/,
121 NumLib::GenericIntegrationMethod const& integration_method,
122 bool const is_axially_symmetric,
123 RichardsMechanicsProcessData<DisplacementDim>& process_data)
124 : _process_data(process_data),
125 _integration_method(integration_method),
126 _element(e),
127 _is_axially_symmetric(is_axially_symmetric)
128{
129 unsigned const n_integration_points =
130 _integration_method.getNumberOfPoints();
131
132 _ip_data.reserve(n_integration_points);
133 _secondary_data.N_u.resize(n_integration_points);
134
135 auto const shape_matrices_u =
136 NumLib::initShapeMatrices<ShapeFunctionDisplacement,
137 ShapeMatricesTypeDisplacement,
138 DisplacementDim>(e, is_axially_symmetric,
139 _integration_method);
140
141 auto const shape_matrices_p =
142 NumLib::initShapeMatrices<ShapeFunctionPressure,
143 ShapeMatricesTypePressure, DisplacementDim>(
144 e, is_axially_symmetric, _integration_method);
145
146 auto const& solid_material =
147 MaterialLib::Solids::selectSolidConstitutiveRelation(
148 _process_data.solid_materials, _process_data.material_ids,
149 e.getID());
150
151 auto const& medium = _process_data.media_map.getMedium(_element.getID());
152
153 ParameterLib::SpatialPosition x_position;
154 x_position.setElementID(_element.getID());
155 for (unsigned ip = 0; ip < n_integration_points; ip++)
156 {
157 x_position.setIntegrationPoint(ip);
158 _ip_data.emplace_back(solid_material);
159 auto& ip_data = _ip_data[ip];
160 auto const& sm_u = shape_matrices_u[ip];
161 _ip_data[ip].integration_weight =
162 _integration_method.getWeightedPoint(ip).getWeight() *
163 sm_u.integralMeasure * sm_u.detJ;
164
165 ip_data.N_u = sm_u.N;
166 ip_data.dNdx_u = sm_u.dNdx;
167
168 ip_data.N_p = shape_matrices_p[ip].N;
169 ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
170
171 // Initial porosity. Could be read from integration point data or mesh.
172 ip_data.porosity =
173 medium->property(MPL::porosity)
174 .template initialValue<double>(
175 x_position,
176 std::numeric_limits<
177 double>::quiet_NaN() /* t independent */);
178
179 ip_data.transport_porosity = ip_data.porosity;
180 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
181 {
182 ip_data.transport_porosity =
183 medium->property(MPL::transport_porosity)
184 .template initialValue<double>(
185 x_position,
186 std::numeric_limits<
187 double>::quiet_NaN() /* t independent */);
188 }
189
190 _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
191 }
192}
193
194template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
195 int DisplacementDim>
196std::size_t RichardsMechanicsLocalAssembler<
197 ShapeFunctionDisplacement, ShapeFunctionPressure,
198 DisplacementDim>::setIPDataInitialConditions(std::string_view name,
199 double const* values,
200 int const integration_order)
201{
202 if (integration_order !=
203 static_cast<int>(_integration_method.getIntegrationOrder()))
204 {
205 OGS_FATAL(
206 "Setting integration point initial conditions; The integration "
207 "order of the local assembler for element {:d} is different "
208 "from the integration order in the initial condition.",
209 _element.getID());
210 }
211
212 if (name == "sigma")
213 {
214 if (_process_data.initial_stress != nullptr)
215 {
216 OGS_FATAL(
217 "Setting initial conditions for stress from integration "
218 "point data and from a parameter '{:s}' is not possible "
219 "simultaneously.",
220 _process_data.initial_stress->name);
221 }
222 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
223 values, _ip_data, &IpData::sigma_eff);
224 }
225
226 if (name == "saturation")
227 {
228 return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
229 &IpData::saturation);
230 }
231 if (name == "porosity")
232 {
233 return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
234 &IpData::porosity);
235 }
236 if (name == "transport_porosity")
237 {
238 return ProcessLib::setIntegrationPointScalarData(
239 values, _ip_data, &IpData::transport_porosity);
240 }
241 if (name == "swelling_stress")
242 {
243 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
244 values, _ip_data, &IpData::sigma_sw);
245 }
246 if (name == "epsilon")
247 {
248 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
249 values, _ip_data, &IpData::eps);
250 }
251 if (name.starts_with("material_state_variable_"))
252 {
253 name.remove_prefix(24);
254
255 // Using first ip data for solid material. TODO (naumov) move solid
256 // material into element, store only material state in IPs.
257 auto const& internal_variables =
258 _ip_data[0].solid_material.getInternalVariables();
259 if (auto const iv = std::find_if(
260 begin(internal_variables), end(internal_variables),
261 [&name](auto const& iv) { return iv.name == name; });
262 iv != end(internal_variables))
263 {
264 DBUG("Setting material state variable '{:s}'", name);
265 return ProcessLib::setIntegrationPointDataMaterialStateVariables(
266 values, _ip_data, &IpData::material_state_variables,
267 iv->reference);
268 }
269
270 ERR("Could not find variable {:s} in solid material model's internal "
271 "variables.",
272 name);
273 }
274 return 0;
275}
276
277template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
278 int DisplacementDim>
279void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
280 ShapeFunctionPressure, DisplacementDim>::
281 setInitialConditionsConcrete(Eigen::VectorXd const local_x,
282 double const t,
283 int const /*process_id*/)
284{
285 assert(local_x.size() == pressure_size + displacement_size);
286
287 auto const p_L = local_x.template segment<pressure_size>(pressure_index);
288
289 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
290 auto const& medium = _process_data.media_map.getMedium(_element.getID());
291 MPL::VariableArray variables;
292
293 ParameterLib::SpatialPosition x_position;
294 x_position.setElementID(_element.getID());
295
296 auto const& solid_phase = medium->phase("Solid");
297
298 unsigned const n_integration_points =
299 _integration_method.getNumberOfPoints();
300 for (unsigned ip = 0; ip < n_integration_points; ip++)
301 {
302 x_position.setIntegrationPoint(ip);
303
304 auto const& N_p = _ip_data[ip].N_p;
305
306 double p_cap_ip;
307 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
308
309 variables.capillary_pressure = p_cap_ip;
310 variables.liquid_phase_pressure = -p_cap_ip;
311 // setting pG to 1 atm
312 // TODO : rewrite equations s.t. p_L = pG-p_cap
313 variables.gas_phase_pressure = 1.0e5;
314 _ip_data[ip].liquid_pressure_m_prev = -p_cap_ip;
315 _ip_data[ip].liquid_pressure_m = -p_cap_ip;
316
317 auto const temperature =
318 medium->property(MPL::PropertyType::reference_temperature)
319 .template value<double>(variables, x_position, t, dt);
320 variables.temperature = temperature;
321
322 _ip_data[ip].saturation_prev =
323 medium->property(MPL::PropertyType::saturation)
324 .template value<double>(variables, x_position, t, dt);
325
326 if (medium->hasProperty(MPL::PropertyType::saturation_micro))
327 {
328 MPL::VariableArray vars;
329 vars.capillary_pressure = p_cap_ip;
330 double const S_L_m =
331 medium->property(MPL::PropertyType::saturation_micro)
332 .template value<double>(vars, x_position, t, dt);
333 _ip_data[ip].saturation_m_prev = S_L_m;
334 }
335
336 // Set eps_m_prev from potentially non-zero eps and sigma_sw from
337 // restart.
338 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
339 t, x_position, dt, temperature);
340 auto& eps = _ip_data[ip].eps;
341 auto& sigma_sw = _ip_data[ip].sigma_sw;
342
343 _ip_data[ip].eps_m_prev.noalias() =
344 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
345 ? eps + C_el.inverse() * sigma_sw
346 : eps;
347 }
348}
349
350template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
351 int DisplacementDim>
352void RichardsMechanicsLocalAssembler<
353 ShapeFunctionDisplacement, ShapeFunctionPressure,
354 DisplacementDim>::assemble(double const t, double const dt,
355 std::vector<double> const& local_x,
356 std::vector<double> const& local_x_prev,
357 std::vector<double>& local_M_data,
358 std::vector<double>& local_K_data,
359 std::vector<double>& local_rhs_data)
360{
361 assert(local_x.size() == pressure_size + displacement_size);
362
363 auto p_L =
364 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
365 pressure_size> const>(local_x.data() + pressure_index,
366 pressure_size);
367
368 auto u =
369 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
370 displacement_size> const>(local_x.data() + displacement_index,
371 displacement_size);
372
373 auto p_L_prev =
374 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
375 pressure_size> const>(local_x_prev.data() + pressure_index,
376 pressure_size);
377
378 auto u_prev =
379 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
380 displacement_size> const>(local_x_prev.data() + displacement_index,
381 displacement_size);
382
383 auto K = MathLib::createZeroedMatrix<
384 typename ShapeMatricesTypeDisplacement::template MatrixType<
385 displacement_size + pressure_size,
386 displacement_size + pressure_size>>(
387 local_K_data, displacement_size + pressure_size,
388 displacement_size + pressure_size);
389
390 auto M = MathLib::createZeroedMatrix<
391 typename ShapeMatricesTypeDisplacement::template MatrixType<
392 displacement_size + pressure_size,
393 displacement_size + pressure_size>>(
394 local_M_data, displacement_size + pressure_size,
395 displacement_size + pressure_size);
396
397 auto rhs = MathLib::createZeroedVector<
398 typename ShapeMatricesTypeDisplacement::template VectorType<
399 displacement_size + pressure_size>>(
400 local_rhs_data, displacement_size + pressure_size);
401
402 auto const& identity2 = MathLib::KelvinVector::Invariants<
403 MathLib::KelvinVector::kelvin_vector_dimensions(
404 DisplacementDim)>::identity2;
405
406 auto const& medium = _process_data.media_map.getMedium(_element.getID());
407 auto const& liquid_phase = medium->phase("AqueousLiquid");
408 auto const& solid_phase = medium->phase("Solid");
409 MPL::VariableArray variables;
410 MPL::VariableArray variables_prev;
411
412 ParameterLib::SpatialPosition x_position;
413 x_position.setElementID(_element.getID());
414
415 unsigned const n_integration_points =
416 _integration_method.getNumberOfPoints();
417 for (unsigned ip = 0; ip < n_integration_points; ip++)
418 {
419 x_position.setIntegrationPoint(ip);
420 auto const& w = _ip_data[ip].integration_weight;
421
422 auto const& N_u = _ip_data[ip].N_u;
423 auto const& dNdx_u = _ip_data[ip].dNdx_u;
424
425 auto const& N_p = _ip_data[ip].N_p;
426 auto const& dNdx_p = _ip_data[ip].dNdx_p;
427
428 auto const x_coord =
429 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
430 ShapeMatricesTypeDisplacement>(
431 _element, N_u);
432 auto const B =
433 LinearBMatrix::computeBMatrix<DisplacementDim,
434 ShapeFunctionDisplacement::NPOINTS,
435 typename BMatricesType::BMatrixType>(
436 dNdx_u, N_u, x_coord, _is_axially_symmetric);
437
438 auto& eps = _ip_data[ip].eps;
439 auto& eps_m = _ip_data[ip].eps_m;
440 eps.noalias() = B * u;
441
442 auto& S_L = _ip_data[ip].saturation;
443 auto const S_L_prev = _ip_data[ip].saturation_prev;
444
445 double p_cap_ip;
446 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
447
448 double p_cap_prev_ip;
449 NumLib::shapeFunctionInterpolate(-p_L_prev, N_p, p_cap_prev_ip);
450
451 variables.capillary_pressure = p_cap_ip;
452 variables.liquid_phase_pressure = -p_cap_ip;
453 // setting pG to 1 atm
454 // TODO : rewrite equations s.t. p_L = pG-p_cap
455 variables.gas_phase_pressure = 1.0e5;
456
457 auto const temperature =
458 medium->property(MPL::PropertyType::reference_temperature)
459 .template value<double>(variables, x_position, t, dt);
460 variables.temperature = temperature;
461
462 auto const alpha =
463 medium->property(MPL::PropertyType::biot_coefficient)
464 .template value<double>(variables, x_position, t, dt);
465 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
466 t, x_position, dt, temperature);
467
468 auto const beta_SR =
469 (1 - alpha) /
470 _ip_data[ip].solid_material.getBulkModulus(t, x_position, &C_el);
471 variables.grain_compressibility = beta_SR;
472
473 auto const rho_LR =
474 liquid_phase.property(MPL::PropertyType::density)
475 .template value<double>(variables, x_position, t, dt);
476 variables.density = rho_LR;
477
478 auto const& b = _process_data.specific_body_force;
479
480 S_L = medium->property(MPL::PropertyType::saturation)
481 .template value<double>(variables, x_position, t, dt);
482 variables.liquid_saturation = S_L;
483 variables_prev.liquid_saturation = S_L_prev;
484
485 // tangent derivative for Jacobian
486 double const dS_L_dp_cap =
487 medium->property(MPL::PropertyType::saturation)
488 .template dValue<double>(variables,
489 MPL::Variable::capillary_pressure,
490 x_position, t, dt);
491 // secant derivative from time discretization for storage
492 // use tangent, if secant is not available
493 double const DeltaS_L_Deltap_cap =
494 (p_cap_ip == p_cap_prev_ip)
495 ? dS_L_dp_cap
496 : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
497
498 auto const chi = [medium, x_position, t, dt](double const S_L)
499 {
500 MPL::VariableArray vs;
501 vs.liquid_saturation = S_L;
502 return medium->property(MPL::PropertyType::bishops_effective_stress)
503 .template value<double>(vs, x_position, t, dt);
504 };
505 double const chi_S_L = chi(S_L);
506 double const chi_S_L_prev = chi(S_L_prev);
507
508 double const p_FR = -chi_S_L * p_cap_ip;
509 variables.effective_pore_pressure = p_FR;
510 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
511
512 // Set volumetric strain rate for the general case without swelling.
513 variables.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
514 variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
515
516 auto& phi = _ip_data[ip].porosity;
517 { // Porosity update
518 variables_prev.porosity = _ip_data[ip].porosity_prev;
519 phi = medium->property(MPL::PropertyType::porosity)
520 .template value<double>(variables, variables_prev,
521 x_position, t, dt);
522 variables.porosity = phi;
523 }
524
525 if (alpha < phi)
526 {
527 OGS_FATAL(
528 "RichardsMechanics: Biot-coefficient {} is smaller than "
529 "porosity {} in element/integration point {}/{}.",
530 alpha, phi, _element.getID(), ip);
531 }
532
533 // Swelling and possibly volumetric strain rate update.
534 {
535 auto& sigma_sw = _ip_data[ip].sigma_sw;
536 auto const& sigma_sw_prev = _ip_data[ip].sigma_sw_prev;
537
538 // If there is swelling, compute it. Update volumetric strain rate,
539 // s.t. it corresponds to the mechanical part only.
540 sigma_sw = sigma_sw_prev;
541 if (solid_phase.hasProperty(
542 MPL::PropertyType::swelling_stress_rate))
543 {
544 auto const sigma_sw_dot =
545 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
546 MPL::formEigenTensor<3>(
547 solid_phase[MPL::PropertyType::swelling_stress_rate]
548 .value(variables, variables_prev, x_position, t,
549 dt)));
550 sigma_sw += sigma_sw_dot * dt;
551
552 // !!! Misusing volumetric strain for mechanical volumetric
553 // strain just to update the transport porosity !!!
554 variables.volumetric_strain +=
555 identity2.transpose() * C_el.inverse() * sigma_sw;
556 variables_prev.volumetric_strain +=
557 identity2.transpose() * C_el.inverse() * sigma_sw_prev;
558 }
559
560 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
561 {
562 variables_prev.transport_porosity =
563 _ip_data[ip].transport_porosity_prev;
564
565 _ip_data[ip].transport_porosity =
566 medium->property(MPL::PropertyType::transport_porosity)
567 .template value<double>(variables, variables_prev,
568 x_position, t, dt);
569 variables.transport_porosity = _ip_data[ip].transport_porosity;
570 }
571 else
572 {
573 variables.transport_porosity = phi;
574 }
575 }
576
577 double const k_rel =
578 medium->property(MPL::PropertyType::relative_permeability)
579 .template value<double>(variables, x_position, t, dt);
580 auto const mu =
581 liquid_phase.property(MPL::PropertyType::viscosity)
582 .template value<double>(variables, x_position, t, dt);
583
584 auto const& sigma_sw = _ip_data[ip].sigma_sw;
585 auto const& sigma_eff = _ip_data[ip].sigma_eff;
586
587 // Set mechanical variables for the intrinsic permeability model
588 // For stress dependent permeability.
589 {
590 auto const sigma_total =
591 (sigma_eff - alpha * p_FR * identity2).eval();
592
593 // For stress dependent permeability.
594 variables.total_stress.emplace<SymmetricTensor>(
595 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
596 sigma_total));
597 }
598
599 variables.equivalent_plastic_strain =
600 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
601
602 auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
603 medium->property(MPL::PropertyType::permeability)
604 .value(variables, x_position, t, dt));
605
606 GlobalDimMatrixType const rho_K_over_mu =
607 K_intrinsic * rho_LR * k_rel / mu;
608
609 //
610 // displacement equation, displacement part
611 //
612 eps_m.noalias() =
613 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
614 ? eps + C_el.inverse() * sigma_sw
615 : eps;
616 variables.mechanical_strain
617 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
618 eps_m);
619
620 _ip_data[ip].updateConstitutiveRelation(variables, t, x_position, dt,
621 temperature);
622
623 // p_SR
624 variables.solid_grain_pressure =
625 p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
626 auto const rho_SR =
627 solid_phase.property(MPL::PropertyType::density)
628 .template value<double>(variables, x_position, t, dt);
629
630 //
631 // displacement equation, displacement part
632 //
633 double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
634 rhs.template segment<displacement_size>(displacement_index).noalias() -=
635 (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
636
637 //
638 // pressure equation, pressure part.
639 //
640 auto const beta_LR =
641 1 / rho_LR *
642 liquid_phase.property(MPL::PropertyType::density)
643 .template dValue<double>(variables,
644 MPL::Variable::liquid_phase_pressure,
645 x_position, t, dt);
646
647 double const a0 = S_L * (alpha - phi) * beta_SR;
648 // Volumetric average specific storage of the solid and fluid phases.
649 double const specific_storage =
650 DeltaS_L_Deltap_cap * (p_cap_ip * a0 - phi) +
651 S_L * (phi * beta_LR + a0);
652 M.template block<pressure_size, pressure_size>(pressure_index,
653 pressure_index)
654 .noalias() += N_p.transpose() * rho_LR * specific_storage * N_p * w;
655
656 K.template block<pressure_size, pressure_size>(pressure_index,
657 pressure_index)
658 .noalias() += dNdx_p.transpose() * rho_K_over_mu * dNdx_p * w;
659
660 rhs.template segment<pressure_size>(pressure_index).noalias() +=
661 dNdx_p.transpose() * rho_LR * rho_K_over_mu * b * w;
662
663 //
664 // displacement equation, pressure part
665 //
666 K.template block<displacement_size, pressure_size>(displacement_index,
667 pressure_index)
668 .noalias() -= B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
669
670 //
671 // pressure equation, displacement part.
672 //
673 M.template block<pressure_size, displacement_size>(pressure_index,
674 displacement_index)
675 .noalias() += N_p.transpose() * S_L * rho_LR * alpha *
676 identity2.transpose() * B * w;
677 }
678
679 if (_process_data.apply_mass_lumping)
680 {
681 auto Mpp = M.template block<pressure_size, pressure_size>(
682 pressure_index, pressure_index);
683 Mpp = Mpp.colwise().sum().eval().asDiagonal();
684 }
685}
686
687template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
688 int DisplacementDim>
689void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
690 ShapeFunctionPressure, DisplacementDim>::
691 assembleWithJacobian(double const t, double const dt,
692 std::vector<double> const& local_x,
693 std::vector<double> const& local_x_prev,
694 std::vector<double>& /*local_M_data*/,
695 std::vector<double>& /*local_K_data*/,
696 std::vector<double>& local_rhs_data,
697 std::vector<double>& local_Jac_data)
698{
699 assert(local_x.size() == pressure_size + displacement_size);
700
701 auto p_L =
702 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
703 pressure_size> const>(local_x.data() + pressure_index,
704 pressure_size);
705
706 auto u =
707 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
708 displacement_size> const>(local_x.data() + displacement_index,
709 displacement_size);
710
711 auto p_L_prev =
712 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
713 pressure_size> const>(local_x_prev.data() + pressure_index,
714 pressure_size);
715 auto u_prev =
716 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
717 displacement_size> const>(local_x_prev.data() + displacement_index,
718 displacement_size);
719
720 auto local_Jac = MathLib::createZeroedMatrix<
721 typename ShapeMatricesTypeDisplacement::template MatrixType<
722 displacement_size + pressure_size,
723 displacement_size + pressure_size>>(
724 local_Jac_data, displacement_size + pressure_size,
725 displacement_size + pressure_size);
726
727 auto local_rhs = MathLib::createZeroedVector<
728 typename ShapeMatricesTypeDisplacement::template VectorType<
729 displacement_size + pressure_size>>(
730 local_rhs_data, displacement_size + pressure_size);
731
732 auto const& identity2 = MathLib::KelvinVector::Invariants<
733 MathLib::KelvinVector::kelvin_vector_dimensions(
734 DisplacementDim)>::identity2;
735
736 typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
737 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
738 pressure_size);
739
740 typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_p =
741 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
742 pressure_size);
743
744 typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_S_Jpp =
745 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
746 pressure_size);
747
748 typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_S =
749 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
750 pressure_size);
751
752 typename ShapeMatricesTypeDisplacement::template MatrixType<
753 displacement_size, pressure_size>
754 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
755 displacement_size, pressure_size>::Zero(displacement_size,
756 pressure_size);
757
758 typename ShapeMatricesTypeDisplacement::template MatrixType<
759 pressure_size, displacement_size>
760 Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
761 pressure_size, displacement_size>::Zero(pressure_size,
762 displacement_size);
763
764 auto const& medium = _process_data.media_map.getMedium(_element.getID());
765 auto const& liquid_phase = medium->phase("AqueousLiquid");
766 auto const& solid_phase = medium->phase("Solid");
767 MPL::VariableArray variables;
768 MPL::VariableArray variables_prev;
769
770 ParameterLib::SpatialPosition x_position;
771 x_position.setElementID(_element.getID());
772
773 unsigned const n_integration_points =
774 _integration_method.getNumberOfPoints();
775 for (unsigned ip = 0; ip < n_integration_points; ip++)
776 {
777 x_position.setIntegrationPoint(ip);
778 auto const& w = _ip_data[ip].integration_weight;
779
780 auto const& N_u = _ip_data[ip].N_u;
781 auto const& dNdx_u = _ip_data[ip].dNdx_u;
782
783 auto const& N_p = _ip_data[ip].N_p;
784 auto const& dNdx_p = _ip_data[ip].dNdx_p;
785
786 auto const x_coord =
787 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
788 ShapeMatricesTypeDisplacement>(
789 _element, N_u);
790 auto const B =
791 LinearBMatrix::computeBMatrix<DisplacementDim,
792 ShapeFunctionDisplacement::NPOINTS,
793 typename BMatricesType::BMatrixType>(
794 dNdx_u, N_u, x_coord, _is_axially_symmetric);
795
796 double p_cap_ip;
797 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
798
799 double p_cap_prev_ip;
800 NumLib::shapeFunctionInterpolate(-p_L_prev, N_p, p_cap_prev_ip);
801
802 variables.capillary_pressure = p_cap_ip;
803 variables.liquid_phase_pressure = -p_cap_ip;
804 // setting pG to 1 atm
805 // TODO : rewrite equations s.t. p_L = pG-p_cap
806 variables.gas_phase_pressure = 1.0e5;
807
808 auto const temperature =
809 medium->property(MPL::PropertyType::reference_temperature)
810 .template value<double>(variables, x_position, t, dt);
811 variables.temperature = temperature;
812
813 auto& eps = _ip_data[ip].eps;
814 auto& eps_m = _ip_data[ip].eps_m;
815 eps.noalias() = B * u;
816 auto const& sigma_eff = _ip_data[ip].sigma_eff;
817 auto& S_L = _ip_data[ip].saturation;
818 auto const S_L_prev = _ip_data[ip].saturation_prev;
819 auto const alpha =
820 medium->property(MPL::PropertyType::biot_coefficient)
821 .template value<double>(variables, x_position, t, dt);
822
823 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
824 t, x_position, dt, temperature);
825
826 auto const beta_SR =
827 (1 - alpha) /
828 _ip_data[ip].solid_material.getBulkModulus(t, x_position, &C_el);
829 variables.grain_compressibility = beta_SR;
830
831 auto const rho_LR =
832 liquid_phase.property(MPL::PropertyType::density)
833 .template value<double>(variables, x_position, t, dt);
834 variables.density = rho_LR;
835
836 auto const& b = _process_data.specific_body_force;
837
838 S_L = medium->property(MPL::PropertyType::saturation)
839 .template value<double>(variables, x_position, t, dt);
840 variables.liquid_saturation = S_L;
841 variables_prev.liquid_saturation = S_L_prev;
842
843 // tangent derivative for Jacobian
844 double const dS_L_dp_cap =
845 medium->property(MPL::PropertyType::saturation)
846 .template dValue<double>(variables,
847 MPL::Variable::capillary_pressure,
848 x_position, t, dt);
849 // secant derivative from time discretization for storage
850 // use tangent, if secant is not available
851 double const DeltaS_L_Deltap_cap =
852 (p_cap_ip == p_cap_prev_ip)
853 ? dS_L_dp_cap
854 : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
855
856 auto const chi = [medium, x_position, t, dt](double const S_L)
857 {
858 MPL::VariableArray vs;
859 vs.liquid_saturation = S_L;
860 return medium->property(MPL::PropertyType::bishops_effective_stress)
861 .template value<double>(vs, x_position, t, dt);
862 };
863 double const chi_S_L = chi(S_L);
864 double const chi_S_L_prev = chi(S_L_prev);
865
866 double const p_FR = -chi_S_L * p_cap_ip;
867 variables.effective_pore_pressure = p_FR;
868 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
869
870 // Set volumetric strain rate for the general case without swelling.
871 variables.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
872 variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
873
874 auto& phi = _ip_data[ip].porosity;
875 { // Porosity update
876
877 variables_prev.porosity = _ip_data[ip].porosity_prev;
878 phi = medium->property(MPL::PropertyType::porosity)
879 .template value<double>(variables, variables_prev,
880 x_position, t, dt);
881 variables.porosity = phi;
882 }
883
884 if (alpha < phi)
885 {
886 OGS_FATAL(
887 "RichardsMechanics: Biot-coefficient {} is smaller than "
888 "porosity {} in element/integration point {}/{}.",
889 alpha, phi, _element.getID(), ip);
890 }
891
892 auto const mu =
893 liquid_phase.property(MPL::PropertyType::viscosity)
894 .template value<double>(variables, x_position, t, dt);
895
896 // Swelling and possibly volumetric strain rate update.
897 updateSwellingStressAndVolumetricStrain<DisplacementDim>(
898 _ip_data[ip], *medium, solid_phase, C_el, rho_LR, mu,
899 _process_data.micro_porosity_parameters, alpha, phi, p_cap_ip,
900 variables, variables_prev, x_position, t, dt);
901
902 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
903 {
904 if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
905 {
906 variables_prev.transport_porosity =
907 _ip_data[ip].transport_porosity_prev;
908
909 _ip_data[ip].transport_porosity =
910 medium->property(MPL::PropertyType::transport_porosity)
911 .template value<double>(variables, variables_prev,
912 x_position, t, dt);
913 variables.transport_porosity = _ip_data[ip].transport_porosity;
914 }
915 }
916 else
917 {
918 variables.transport_porosity = phi;
919 }
920
921 double const k_rel =
922 medium->property(MPL::PropertyType::relative_permeability)
923 .template value<double>(variables, x_position, t, dt);
924
925 // Set mechanical variables for the intrinsic permeability model
926 // For stress dependent permeability.
927 {
928 auto const sigma_total =
929 (_ip_data[ip].sigma_eff + alpha * p_FR * identity2).eval();
930 // For stress dependent permeability.
931 variables.total_stress.emplace<SymmetricTensor>(
932 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
933 sigma_total));
934 }
935
936 variables.equivalent_plastic_strain =
937 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
938
939 auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
940 medium->property(MPL::PropertyType::permeability)
941 .value(variables, x_position, t, dt));
942
943 GlobalDimMatrixType const rho_Ki_over_mu = K_intrinsic * rho_LR / mu;
944
945 //
946 // displacement equation, displacement part
947 //
948 auto& sigma_sw = _ip_data[ip].sigma_sw;
949
950 eps_m.noalias() =
951 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
952 ? eps + C_el.inverse() * sigma_sw
953 : eps;
954 variables.mechanical_strain
955 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
956 eps_m);
957
958 auto C = _ip_data[ip].updateConstitutiveRelation(
959 variables, t, x_position, dt, temperature);
960
961 local_Jac
962 .template block<displacement_size, displacement_size>(
963 displacement_index, displacement_index)
964 .noalias() += B.transpose() * C * B * w;
965
966 // p_SR
967 variables.solid_grain_pressure =
968 p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
969 auto const rho_SR =
970 solid_phase.property(MPL::PropertyType::density)
971 .template value<double>(variables, x_position, t, dt);
972
973 double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
974 local_rhs.template segment<displacement_size>(displacement_index)
975 .noalias() -=
976 (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
977
978 //
979 // displacement equation, pressure part
980 //
981 Kup.noalias() += B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
982
983 auto const dchi_dS_L =
984 medium->property(MPL::PropertyType::bishops_effective_stress)
985 .template dValue<double>(variables,
986 MPL::Variable::liquid_saturation,
987 x_position, t, dt);
988 local_Jac
989 .template block<displacement_size, pressure_size>(
990 displacement_index, pressure_index)
991 .noalias() -= B.transpose() * alpha *
992 (chi_S_L + dchi_dS_L * p_cap_ip * dS_L_dp_cap) *
993 identity2 * N_p * w;
994
995 local_Jac
996 .template block<displacement_size, pressure_size>(
997 displacement_index, pressure_index)
998 .noalias() +=
999 N_u_op(N_u).transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w;
1000
1001 // For the swelling stress with double structure model the corresponding
1002 // Jacobian u-p entry would be required, but it does not improve
1003 // convergence and sometimes worsens it:
1004 // if (medium->hasProperty(MPL::PropertyType::saturation_micro))
1005 // {
1006 // -B.transpose() *
1007 // dsigma_sw_dS_L_m* dS_L_m_dp_cap_m*(p_L_m - p_L_m_prev) /
1008 // (p_cap_ip - p_cap_prev_ip) * N_p* w;
1009 // }
1010 if (!medium->hasProperty(MPL::PropertyType::saturation_micro) &&
1011 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
1012 {
1013 using DimMatrix = Eigen::Matrix<double, 3, 3>;
1014 auto const dsigma_sw_dS_L =
1015 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
1016 solid_phase
1017 .property(MPL::PropertyType::swelling_stress_rate)
1018 .template dValue<DimMatrix>(
1019 variables, variables_prev,
1020 MPL::Variable::liquid_saturation, x_position, t,
1021 dt));
1022 local_Jac
1023 .template block<displacement_size, pressure_size>(
1024 displacement_index, pressure_index)
1025 .noalias() +=
1026 B.transpose() * dsigma_sw_dS_L * dS_L_dp_cap * N_p * w;
1027 }
1028 //
1029 // pressure equation, displacement part.
1030 //
1031 if (_process_data.explicit_hm_coupling_in_unsaturated_zone)
1032 {
1033 Kpu.noalias() += N_p.transpose() * chi_S_L_prev * rho_LR * alpha *
1034 identity2.transpose() * B * w;
1035 }
1036 else
1037 {
1038 Kpu.noalias() += N_p.transpose() * S_L * rho_LR * alpha *
1039 identity2.transpose() * B * w;
1040 }
1041
1042 //
1043 // pressure equation, pressure part.
1044 //
1045 laplace_p.noalias() +=
1046 dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
1047
1048 auto const beta_LR =
1049 1 / rho_LR *
1050 liquid_phase.property(MPL::PropertyType::density)
1051 .template dValue<double>(variables,
1052 MPL::Variable::liquid_phase_pressure,
1053 x_position, t, dt);
1054
1055 double const a0 = (alpha - phi) * beta_SR;
1056 double const specific_storage_a_p = S_L * (phi * beta_LR + S_L * a0);
1057 double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
1058
1059 double const dspecific_storage_a_p_dp_cap =
1060 dS_L_dp_cap * (phi * beta_LR + 2 * S_L * a0);
1061 double const dspecific_storage_a_S_dp_cap =
1062 -a0 * (S_L + p_cap_ip * dS_L_dp_cap);
1063
1064 storage_p_a_p.noalias() +=
1065 N_p.transpose() * rho_LR * specific_storage_a_p * N_p * w;
1066
1067 storage_p_a_S.noalias() -= N_p.transpose() * rho_LR *
1068 specific_storage_a_S * DeltaS_L_Deltap_cap *
1069 N_p * w;
1070
1071 local_Jac
1072 .template block<pressure_size, pressure_size>(pressure_index,
1073 pressure_index)
1074 .noalias() += N_p.transpose() * (p_cap_ip - p_cap_prev_ip) / dt *
1075 rho_LR * dspecific_storage_a_p_dp_cap * N_p * w;
1076
1077 storage_p_a_S_Jpp.noalias() -=
1078 N_p.transpose() * rho_LR *
1079 ((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap +
1080 specific_storage_a_S * dS_L_dp_cap) /
1081 dt * N_p * w;
1082
1083 if (!_process_data.explicit_hm_coupling_in_unsaturated_zone)
1084 {
1085 local_Jac
1086 .template block<pressure_size, pressure_size>(pressure_index,
1087 pressure_index)
1088 .noalias() -= N_p.transpose() * rho_LR * dS_L_dp_cap * alpha *
1089 identity2.transpose() * B * (u - u_prev) / dt *
1090 N_p * w;
1091 }
1092
1093 double const dk_rel_dS_l =
1094 medium->property(MPL::PropertyType::relative_permeability)
1095 .template dValue<double>(variables,
1096 MPL::Variable::liquid_saturation,
1097 x_position, t, dt);
1098 typename ShapeMatricesTypeDisplacement::GlobalDimVectorType const
1099 grad_p_cap = -dNdx_p * p_L;
1100 local_Jac
1101 .template block<pressure_size, pressure_size>(pressure_index,
1102 pressure_index)
1103 .noalias() += dNdx_p.transpose() * rho_Ki_over_mu * grad_p_cap *
1104 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
1105
1106 local_Jac
1107 .template block<pressure_size, pressure_size>(pressure_index,
1108 pressure_index)
1109 .noalias() += dNdx_p.transpose() * rho_LR * rho_Ki_over_mu * b *
1110 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
1111
1112 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
1113 dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
1114
1115 if (medium->hasProperty(MPL::PropertyType::saturation_micro))
1116 {
1117 double const alpha_bar = _process_data.micro_porosity_parameters
1118 ->mass_exchange_coefficient;
1119 auto const p_L_m = _ip_data[ip].liquid_pressure_m;
1120 local_rhs.template segment<pressure_size>(pressure_index)
1121 .noalias() -=
1122 N_p.transpose() * alpha_bar / mu * (-p_cap_ip - p_L_m) * w;
1123
1124 local_Jac
1125 .template block<pressure_size, pressure_size>(pressure_index,
1126 pressure_index)
1127 .noalias() += N_p.transpose() * alpha_bar / mu * N_p * w;
1128 if (p_cap_ip != p_cap_prev_ip)
1129 {
1130 double const p_L_m_prev = _ip_data[ip].liquid_pressure_m_prev;
1131 local_Jac
1132 .template block<pressure_size, pressure_size>(
1133 pressure_index, pressure_index)
1134 .noalias() += N_p.transpose() * alpha_bar / mu *
1135 (p_L_m - p_L_m_prev) /
1136 (p_cap_ip - p_cap_prev_ip) * N_p * w;
1137 }
1138 }
1139 }
1140
1141 if (_process_data.apply_mass_lumping)
1142 {
1143 storage_p_a_p = storage_p_a_p.colwise().sum().eval().asDiagonal();
1144 storage_p_a_S = storage_p_a_S.colwise().sum().eval().asDiagonal();
1145 storage_p_a_S_Jpp =
1146 storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
1147 }
1148
1149 // pressure equation, pressure part.
1150 local_Jac
1151 .template block<pressure_size, pressure_size>(pressure_index,
1152 pressure_index)
1153 .noalias() += laplace_p + storage_p_a_p / dt + storage_p_a_S_Jpp;
1154
1155 // pressure equation, displacement part.
1156 local_Jac
1157 .template block<pressure_size, displacement_size>(pressure_index,
1158 displacement_index)
1159 .noalias() = Kpu / dt;
1160
1161 // pressure equation
1162 local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
1163 laplace_p * p_L +
1164 (storage_p_a_p + storage_p_a_S) * (p_L - p_L_prev) / dt +
1165 Kpu * (u - u_prev) / dt;
1166
1167 // displacement equation
1168 local_rhs.template segment<displacement_size>(displacement_index)
1169 .noalias() += Kup * p_L;
1170}
1171
1172template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1173 int DisplacementDim>
1174std::vector<double> RichardsMechanicsLocalAssembler<
1175 ShapeFunctionDisplacement, ShapeFunctionPressure,
1176 DisplacementDim>::getSigma() const
1177{
1178 constexpr int kelvin_vector_size =
1179 MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
1180
1181 return transposeInPlace<kelvin_vector_size>(
1182 [this](std::vector<double>& values)
1183 { return getIntPtSigma(0, {}, {}, values); });
1184}
1185
1186template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1187 int DisplacementDim>
1188std::vector<double> const& RichardsMechanicsLocalAssembler<
1189 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
1190 getIntPtSigma(
1191 const double /*t*/,
1192 std::vector<GlobalVector*> const& /*x*/,
1193 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
1194 std::vector<double>& cache) const
1195{
1196 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
1197 _ip_data, &IpData::sigma_eff, cache);
1198}
1199
1200template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1201 int DisplacementDim>
1202std::vector<double> RichardsMechanicsLocalAssembler<
1203 ShapeFunctionDisplacement, ShapeFunctionPressure,
1204 DisplacementDim>::getSwellingStress() const
1205{
1206 constexpr int kelvin_vector_size =
1207 MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
1208
1209 return transposeInPlace<kelvin_vector_size>(
1210 [this](std::vector<double>& values)
1211 { return getIntPtSwellingStress(0, {}, {}, values); });
1212}
1213
1214template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1215 int DisplacementDim>
1216std::vector<double> const& RichardsMechanicsLocalAssembler<
1217 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
1218 getIntPtSwellingStress(
1219 const double /*t*/,
1220 std::vector<GlobalVector*> const& /*x*/,
1221 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
1222 std::vector<double>& cache) const
1223{
1224 constexpr int kelvin_vector_size =
1225 MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
1226 auto const n_integration_points = _ip_data.size();
1227
1228 cache.clear();
1229 auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
1230 double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
1231 cache, kelvin_vector_size, n_integration_points);
1232
1233 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1234 {
1235 auto const& sigma_sw = _ip_data[ip].sigma_sw;
1236 cache_mat.col(ip) =
1237 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_sw);
1238 }
1239
1240 return cache;
1241}
1242
1243template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1244 int DisplacementDim>
1245std::vector<double> RichardsMechanicsLocalAssembler<
1246 ShapeFunctionDisplacement, ShapeFunctionPressure,
1247 DisplacementDim>::getEpsilon() const
1248{
1249 constexpr int kelvin_vector_size =
1250 MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
1251
1252 return transposeInPlace<kelvin_vector_size>(
1253 [this](std::vector<double>& values)
1254 { return getIntPtEpsilon(0, {}, {}, values); });
1255}
1256
1257template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1258 int DisplacementDim>
1259std::vector<double> const& RichardsMechanicsLocalAssembler<
1260 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
1261 getIntPtEpsilon(
1262 const double /*t*/,
1263 std::vector<GlobalVector*> const& /*x*/,
1264 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
1265 std::vector<double>& cache) const
1266{
1267 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
1268 _ip_data, &IpData::eps, cache);
1269}
1270
1271template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1272 int DisplacementDim>
1273int RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
1274 ShapeFunctionPressure,
1275 DisplacementDim>::getMaterialID() const
1276{
1277 return _process_data.material_ids == nullptr
1278 ? 0
1279 : (*_process_data.material_ids)[_element.getID()];
1280}
1281
1282template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1283 int DisplacementDim>
1284std::vector<double> RichardsMechanicsLocalAssembler<
1285 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
1286 getMaterialStateVariableInternalState(
1287 std::function<std::span<double>(
1288 typename MaterialLib::Solids::MechanicsBase<DisplacementDim>::
1289 MaterialStateVariables&)> const& get_values_span,
1290 int const& n_components) const
1291{
1292 return ProcessLib::getIntegrationPointDataMaterialStateVariables(
1293 _ip_data, &IpData::material_state_variables, get_values_span,
1294 n_components);
1295}
1296
1297template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1298 int DisplacementDim>
1299std::vector<double> const& RichardsMechanicsLocalAssembler<
1300 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
1301 getIntPtDarcyVelocity(
1302 const double /*t*/,
1303 std::vector<GlobalVector*> const& /*x*/,
1304 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
1305 std::vector<double>& cache) const
1306{
1307 unsigned const n_integration_points =
1308 _integration_method.getNumberOfPoints();
1309
1310 cache.clear();
1311 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
1312 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
1313 cache, DisplacementDim, n_integration_points);
1314
1315 for (unsigned ip = 0; ip < n_integration_points; ip++)
1316 {
1317 cache_matrix.col(ip).noalias() = _ip_data[ip].v_darcy;
1318 }
1319
1320 return cache;
1321}
1322
1323template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1324 int DisplacementDim>
1325std::vector<double> RichardsMechanicsLocalAssembler<
1326 ShapeFunctionDisplacement, ShapeFunctionPressure,
1327 DisplacementDim>::getSaturation() const
1328{
1329 std::vector<double> result;
1330 getIntPtSaturation(0, {}, {}, result);
1331 return result;
1332}
1333
1334template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1335 int DisplacementDim>
1336std::vector<double> const& RichardsMechanicsLocalAssembler<
1337 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
1338 getIntPtSaturation(
1339 const double /*t*/,
1340 std::vector<GlobalVector*> const& /*x*/,
1341 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
1342 std::vector<double>& cache) const
1343{
1344 return ProcessLib::getIntegrationPointScalarData(
1345 _ip_data, &IpData::saturation, cache);
1346}
1347
1348template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1349 int DisplacementDim>
1350std::vector<double> RichardsMechanicsLocalAssembler<
1351 ShapeFunctionDisplacement, ShapeFunctionPressure,
1352 DisplacementDim>::getMicroSaturation() const
1353{
1354 std::vector<double> result;
1355 getIntPtMicroSaturation(0, {}, {}, result);
1356 return result;
1357}
1358
1359template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1360 int DisplacementDim>
1361std::vector<double> const& RichardsMechanicsLocalAssembler<
1362 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
1363 getIntPtMicroSaturation(
1364 const double /*t*/,
1365 std::vector<GlobalVector*> const& /*x*/,
1366 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
1367 std::vector<double>& cache) const
1368{
1369 return ProcessLib::getIntegrationPointScalarData(
1370 _ip_data, &IpData::saturation_m, cache);
1371}
1372
1373template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1374 int DisplacementDim>
1375std::vector<double> RichardsMechanicsLocalAssembler<
1376 ShapeFunctionDisplacement, ShapeFunctionPressure,
1377 DisplacementDim>::getMicroPressure() const
1378{
1379 std::vector<double> result;
1380 getIntPtMicroPressure(0, {}, {}, result);
1381 return result;
1382}
1383
1384template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1385 int DisplacementDim>
1386std::vector<double> const& RichardsMechanicsLocalAssembler<
1387 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
1388 getIntPtMicroPressure(
1389 const double /*t*/,
1390 std::vector<GlobalVector*> const& /*x*/,
1391 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
1392 std::vector<double>& cache) const
1393{
1394 return ProcessLib::getIntegrationPointScalarData(
1395 _ip_data, &IpData::liquid_pressure_m, cache);
1396}
1397
1398template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1399 int DisplacementDim>
1400std::vector<double> RichardsMechanicsLocalAssembler<
1401 ShapeFunctionDisplacement, ShapeFunctionPressure,
1402 DisplacementDim>::getPorosity() const
1403{
1404 std::vector<double> result;
1405 getIntPtPorosity(0, {}, {}, result);
1406 return result;
1407}
1408
1409template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1410 int DisplacementDim>
1411std::vector<double> const& RichardsMechanicsLocalAssembler<
1412 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
1413 getIntPtPorosity(
1414 const double /*t*/,
1415 std::vector<GlobalVector*> const& /*x*/,
1416 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
1417 std::vector<double>& cache) const
1418{
1419 return ProcessLib::getIntegrationPointScalarData(_ip_data,
1420 &IpData::porosity, cache);
1421}
1422
1423template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1424 int DisplacementDim>
1425std::vector<double> RichardsMechanicsLocalAssembler<
1426 ShapeFunctionDisplacement, ShapeFunctionPressure,
1427 DisplacementDim>::getTransportPorosity() const
1428{
1429 std::vector<double> result;
1430 getIntPtTransportPorosity(0, {}, {}, result);
1431 return result;
1432}
1433
1434template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1435 int DisplacementDim>
1436std::vector<double> const& RichardsMechanicsLocalAssembler<
1437 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
1438 getIntPtTransportPorosity(
1439 const double /*t*/,
1440 std::vector<GlobalVector*> const& /*x*/,
1441 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
1442 std::vector<double>& cache) const
1443{
1444 return ProcessLib::getIntegrationPointScalarData(
1445 _ip_data, &IpData::transport_porosity, cache);
1446}
1447
1448template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1449 int DisplacementDim>
1450std::vector<double> const& RichardsMechanicsLocalAssembler<
1451 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
1452 getIntPtDryDensitySolid(
1453 const double /*t*/,
1454 std::vector<GlobalVector*> const& /*x*/,
1455 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
1456 std::vector<double>& cache) const
1457{
1458 return ProcessLib::getIntegrationPointScalarData(
1459 _ip_data, &IpData::dry_density_solid, cache);
1460}
1461
1462template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1463 int DisplacementDim>
1464void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
1465 ShapeFunctionPressure, DisplacementDim>::
1466 assembleWithJacobianForPressureEquations(
1467 const double /*t*/, double const /*dt*/,
1468 Eigen::VectorXd const& /*local_x*/,
1469 Eigen::VectorXd const& /*local_x_prev*/,
1470 std::vector<double>& /*local_M_data*/,
1471 std::vector<double>& /*local_K_data*/,
1472 std::vector<double>& /*local_b_data*/,
1473 std::vector<double>& /*local_Jac_data*/)
1474{
1475 OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
1476}
1477
1478template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1479 int DisplacementDim>
1480void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
1481 ShapeFunctionPressure, DisplacementDim>::
1482 assembleWithJacobianForDeformationEquations(
1483 const double /*t*/, double const /*dt*/,
1484 Eigen::VectorXd const& /*local_x*/,
1485 Eigen::VectorXd const& /*local_x_prev*/,
1486 std::vector<double>& /*local_M_data*/,
1487 std::vector<double>& /*local_K_data*/,
1488 std::vector<double>& /*local_b_data*/,
1489 std::vector<double>& /*local_Jac_data*/)
1490{
1491 OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
1492}
1493
1494template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1495 int DisplacementDim>
1496void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
1497 ShapeFunctionPressure, DisplacementDim>::
1498 assembleWithJacobianForStaggeredScheme(
1499 double const t, double const dt, Eigen::VectorXd const& local_x,
1500 Eigen::VectorXd const& local_x_prev, int const process_id,
1501 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
1502 std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
1503{
1504 // For the equations with pressure
1505 if (process_id == 0)
1506 {
1507 assembleWithJacobianForPressureEquations(t, dt, local_x, local_x_prev,
1508 local_M_data, local_K_data,
1509 local_b_data, local_Jac_data);
1510 return;
1511 }
1512
1513 // For the equations with deformation
1514 assembleWithJacobianForDeformationEquations(t, dt, local_x, local_x_prev,
1515 local_M_data, local_K_data,
1516 local_b_data, local_Jac_data);
1517}
1518
1519template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1520 int DisplacementDim>
1521void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
1522 ShapeFunctionPressure, DisplacementDim>::
1523 computeSecondaryVariableConcrete(double const t, double const dt,
1524 Eigen::VectorXd const& local_x,
1525 Eigen::VectorXd const& local_x_prev)
1526{
1527 auto p_L = local_x.template segment<pressure_size>(pressure_index);
1528 auto u = local_x.template segment<displacement_size>(displacement_index);
1529
1530 auto p_L_prev =
1531 local_x_prev.template segment<pressure_size>(pressure_index);
1532 auto u_prev =
1533 local_x_prev.template segment<displacement_size>(displacement_index);
1534
1535 auto const& identity2 = MathLib::KelvinVector::Invariants<
1536 MathLib::KelvinVector::kelvin_vector_dimensions(
1537 DisplacementDim)>::identity2;
1538
1539 auto const& medium = _process_data.media_map.getMedium(_element.getID());
1540 auto const& liquid_phase = medium->phase("AqueousLiquid");
1541 auto const& solid_phase = medium->phase("Solid");
1542 MPL::VariableArray variables;
1543 MPL::VariableArray variables_prev;
1544
1545 ParameterLib::SpatialPosition x_position;
1546 x_position.setElementID(_element.getID());
1547
1548 unsigned const n_integration_points =
1549 _integration_method.getNumberOfPoints();
1550
1551 double saturation_avg = 0;
1552 double porosity_avg = 0;
1553
1554 using KV = MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
1555 KV sigma_avg = KV::Zero();
1556
1557 for (unsigned ip = 0; ip < n_integration_points; ip++)
1558 {
1559 x_position.setIntegrationPoint(ip);
1560 auto const& N_p = _ip_data[ip].N_p;
1561 auto const& N_u = _ip_data[ip].N_u;
1562 auto const& dNdx_u = _ip_data[ip].dNdx_u;
1563
1564 auto const x_coord =
1565 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
1566 ShapeMatricesTypeDisplacement>(
1567 _element, N_u);
1568 auto const B =
1569 LinearBMatrix::computeBMatrix<DisplacementDim,
1570 ShapeFunctionDisplacement::NPOINTS,
1571 typename BMatricesType::BMatrixType>(
1572 dNdx_u, N_u, x_coord, _is_axially_symmetric);
1573
1574 double p_cap_ip;
1575 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
1576
1577 double p_cap_prev_ip;
1578 NumLib::shapeFunctionInterpolate(-p_L_prev, N_p, p_cap_prev_ip);
1579
1580 variables.capillary_pressure = p_cap_ip;
1581 variables.liquid_phase_pressure = -p_cap_ip;
1582 // setting pG to 1 atm
1583 // TODO : rewrite equations s.t. p_L = pG-p_cap
1584 variables.gas_phase_pressure = 1.0e5;
1585
1586 auto const temperature =
1587 medium->property(MPL::PropertyType::reference_temperature)
1588 .template value<double>(variables, x_position, t, dt);
1589 variables.temperature = temperature;
1590
1591 auto& eps = _ip_data[ip].eps;
1592 eps.noalias() = B * u;
1593 auto& eps_m = _ip_data[ip].eps_m;
1594 auto& S_L = _ip_data[ip].saturation;
1595 auto const S_L_prev = _ip_data[ip].saturation_prev;
1596 S_L = medium->property(MPL::PropertyType::saturation)
1597 .template value<double>(variables, x_position, t, dt);
1598 variables.liquid_saturation = S_L;
1599 variables_prev.liquid_saturation = S_L_prev;
1600
1601 auto const chi = [medium, x_position, t, dt](double const S_L)
1602 {
1603 MPL::VariableArray vs;
1604 vs.liquid_saturation = S_L;
1605 return medium->property(MPL::PropertyType::bishops_effective_stress)
1606 .template value<double>(vs, x_position, t, dt);
1607 };
1608 double const chi_S_L = chi(S_L);
1609 double const chi_S_L_prev = chi(S_L_prev);
1610
1611 auto const alpha =
1612 medium->property(MPL::PropertyType::biot_coefficient)
1613 .template value<double>(variables, x_position, t, dt);
1614
1615 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
1616 t, x_position, dt, temperature);
1617
1618 auto const beta_SR =
1619 (1 - alpha) /
1620 _ip_data[ip].solid_material.getBulkModulus(t, x_position, &C_el);
1621 variables.grain_compressibility = beta_SR;
1622
1623 variables.effective_pore_pressure = -chi_S_L * p_cap_ip;
1624 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
1625
1626 // Set volumetric strain rate for the general case without swelling.
1627 variables.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
1628 variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
1629
1630 auto& phi = _ip_data[ip].porosity;
1631 { // Porosity update
1632 variables_prev.porosity = _ip_data[ip].porosity_prev;
1633 phi = medium->property(MPL::PropertyType::porosity)
1634 .template value<double>(variables, variables_prev,
1635 x_position, t, dt);
1636 variables.porosity = phi;
1637 }
1638
1639 auto const rho_LR =
1640 liquid_phase.property(MPL::PropertyType::density)
1641 .template value<double>(variables, x_position, t, dt);
1642 variables.density = rho_LR;
1643 auto const mu =
1644 liquid_phase.property(MPL::PropertyType::viscosity)
1645 .template value<double>(variables, x_position, t, dt);
1646
1647 // Swelling and possibly volumetric strain rate update.
1648 updateSwellingStressAndVolumetricStrain<DisplacementDim>(
1649 _ip_data[ip], *medium, solid_phase, C_el, rho_LR, mu,
1650 _process_data.micro_porosity_parameters, alpha, phi, p_cap_ip,
1651 variables, variables_prev, x_position, t, dt);
1652
1653 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
1654 {
1655 if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
1656 {
1657 variables_prev.transport_porosity =
1658 _ip_data[ip].transport_porosity_prev;
1659
1660 _ip_data[ip].transport_porosity =
1661 medium->property(MPL::PropertyType::transport_porosity)
1662 .template value<double>(variables, variables_prev,
1663 x_position, t, dt);
1664 variables.transport_porosity = _ip_data[ip].transport_porosity;
1665 }
1666 }
1667 else
1668 {
1669 variables.transport_porosity = phi;
1670 }
1671
1672 // Set mechanical variables for the intrinsic permeability model
1673 // For stress dependent permeability.
1674 {
1675 auto const sigma_total = (_ip_data[ip].sigma_eff +
1676 alpha * chi_S_L * identity2 * p_cap_ip)
1677 .eval();
1678 // For stress dependent permeability.
1679 variables.total_stress.emplace<SymmetricTensor>(
1680 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
1681 sigma_total));
1682 }
1683
1684 variables.equivalent_plastic_strain =
1685 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
1686
1687 auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
1688 medium->property(MPL::PropertyType::permeability)
1689 .value(variables, x_position, t, dt));
1690
1691 double const k_rel =
1692 medium->property(MPL::PropertyType::relative_permeability)
1693 .template value<double>(variables, x_position, t, dt);
1694
1695 GlobalDimMatrixType const K_over_mu = k_rel * K_intrinsic / mu;
1696
1697 auto const& sigma_eff = _ip_data[ip].sigma_eff;
1698 double const p_FR = -chi_S_L * p_cap_ip;
1699 // p_SR
1700 variables.solid_grain_pressure =
1701 p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
1702 auto const rho_SR =
1703 solid_phase.property(MPL::PropertyType::density)
1704 .template value<double>(variables, x_position, t, dt);
1705 _ip_data[ip].dry_density_solid = (1 - phi) * rho_SR;
1706
1707 auto& sigma_sw = _ip_data[ip].sigma_sw;
1708
1709 eps_m.noalias() =
1710 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
1711 ? eps + C_el.inverse() * sigma_sw
1712 : eps;
1713 variables.mechanical_strain
1714 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
1715 eps_m);
1716
1717 _ip_data[ip].updateConstitutiveRelation(variables, t, x_position, dt,
1718 temperature);
1719
1720 auto const& b = _process_data.specific_body_force;
1721
1722 // Compute the velocity
1723 auto const& dNdx_p = _ip_data[ip].dNdx_p;
1724 _ip_data[ip].v_darcy.noalias() =
1725 -K_over_mu * dNdx_p * p_L + rho_LR * K_over_mu * b;
1726
1727 saturation_avg += S_L;
1728 porosity_avg += phi;
1729 sigma_avg += sigma_eff;
1730 }
1731 saturation_avg /= n_integration_points;
1732 porosity_avg /= n_integration_points;
1733 sigma_avg /= n_integration_points;
1734
1735 (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
1736 (*_process_data.element_porosity)[_element.getID()] = porosity_avg;
1737
1738 Eigen::Map<KV>(&(*_process_data.element_stresses)[_element.getID() *
1739 KV::RowsAtCompileTime]) =
1740 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_avg);
1741
1742 NumLib::interpolateToHigherOrderNodes<
1743 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1744 DisplacementDim>(_element, _is_axially_symmetric, p_L,
1745 *_process_data.pressure_interpolated);
1746}
1747
1748template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1749 int DisplacementDim>
1750unsigned RichardsMechanicsLocalAssembler<
1751 ShapeFunctionDisplacement, ShapeFunctionPressure,
1752 DisplacementDim>::getNumberOfIntegrationPoints() const
1753{
1754 return _integration_method.getNumberOfPoints();
1755}
1756
1757template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1758 int DisplacementDim>
1759typename MaterialLib::Solids::MechanicsBase<
1760 DisplacementDim>::MaterialStateVariables const&
1761RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
1762 ShapeFunctionPressure, DisplacementDim>::
1763 getMaterialStateVariablesAt(unsigned integration_point) const
1764{
1765 return *_ip_data[integration_point].material_state_variables;
1766}
1767} // namespace RichardsMechanics
1768} // namespace ProcessLib
OGS_FATAL
#define OGS_FATAL(...)
Definition Error.h:26
DBUG
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
ERR
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
MaterialPropertyLib::Medium::property
Property const & property(PropertyType const &p) const
Definition Medium.cpp:54
MaterialPropertyLib::Medium::hasProperty
bool hasProperty(PropertyType const &p) const
Definition Medium.cpp:70
MaterialPropertyLib::Phase::property
Property const & property(PropertyType const &p) const
Definition Phase.cpp:53
MaterialPropertyLib::Phase::hasProperty
bool hasProperty(PropertyType const &p) const
Definition Phase.cpp:69
MaterialPropertyLib::VariableArray::mechanical_strain
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > mechanical_strain
Definition VariableType.h:178
MaterialPropertyLib::VariableArray::total_stress
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > total_stress
Definition VariableType.h:197
MathLib::WeightedPoint::getWeight
double getWeight() const
Definition WeightedPoint.h:80
MeshLib::Element::getID
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:89
NumLib::GenericIntegrationMethod::getWeightedPoint
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
Definition GenericIntegrationMethod.h:36
ParameterLib::SpatialPosition::setElementID
void setElementID(std::size_t element_id)
Definition SpatialPosition.h:59
ParameterLib::SpatialPosition::setIntegrationPoint
void setIntegrationPoint(unsigned integration_point)
Definition SpatialPosition.h:65
ProcessLib::BMatrixPolicyType::BMatrixType
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition BMatrixPolicy.h:55
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler::ShapeMatricesTypeDisplacement
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
Definition RichardsMechanicsFEM.h:53
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler::_secondary_data
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
Definition RichardsMechanicsFEM.h:340
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler::SymmetricTensor
Eigen::Matrix< double, KelvinVectorSize, 1 > SymmetricTensor
Definition RichardsMechanicsFEM.h:74
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler::GlobalDimMatrixType
typename ShapeMatricesTypePressure::GlobalDimMatrixType GlobalDimMatrixType
Definition RichardsMechanicsFEM.h:58
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler::ShapeMatricesTypePressure
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
Definition RichardsMechanicsFEM.h:55
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler::_process_data
RichardsMechanicsProcessData< DisplacementDim > & _process_data
Definition RichardsMechanicsFEM.h:331
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler::_ip_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
Definition RichardsMechanicsFEM.h:333
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< 3 >
template Eigen::Matrix< double, 3, 3 > formEigenTensor< 3 >(MaterialPropertyLib::PropertyDataType const &values)
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::KelvinMatrixType
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType
Definition KelvinVector.h:55
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::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::RichardsMechanics::updateSwellingStressAndVolumetricStrain
void updateSwellingStressAndVolumetricStrain(IPData &ip_data, MaterialPropertyLib::Medium const &medium, MaterialPropertyLib::Phase const &solid_phase, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > const &C_el, double const rho_LR, double const mu, std::optional< MicroPorosityParameters > micro_porosity_parameters, double const alpha, double const phi, double const p_cap_ip, MPL::VariableArray &variables, MPL::VariableArray &variables_prev, ParameterLib::SpatialPosition const &x_position, double const t, double const dt)
Definition RichardsMechanicsFEM-impl.h:33
ProcessLib::getIntegrationPointScalarData
std::vector< double > const & getIntegrationPointScalarData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
Definition SetOrGetIntegrationPointData.h:176
ProcessLib::setIntegrationPointDataMaterialStateVariables
std::size_t setIntegrationPointDataMaterialStateVariables(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType member, std::function< std::span< double >(MaterialStateVariables &)> get_values_span)
Definition SetOrGetIntegrationPointData.h:270
ProcessLib::getIntegrationPointDataMaterialStateVariables
std::vector< double > getIntegrationPointDataMaterialStateVariables(IntegrationPointDataVector const &ip_data_vector, MemberType member, std::function< std::span< double >(MaterialStateVariables &)> get_values_span, int const n_components)
Definition SetOrGetIntegrationPointData.h:247
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
EigenFixedShapeMatrixPolicy::GlobalDimVectorType
VectorType< GlobalDim > GlobalDimVectorType
Definition ShapeMatrixPolicy.h:85
ProcessLib::HydroMechanics::SecondaryData::N_u
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u
Definition HydroMechanicsFEM.h:161
ProcessLib::RichardsMechanics::IntegrationPointData::material_state_variables
std::unique_ptr< typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables > material_state_variables
Definition IntegrationPointData.h:78