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