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