OGS 6.1.0-1721-g6382411ad
RichardsMechanicsFEM-impl.h
Go to the documentation of this file.
1 
12 #pragma once
13 
14 #include "RichardsMechanicsFEM.h"
15 
16 #include <cassert>
17 
19 #include "MathLib/KelvinVector.h"
22 
23 namespace ProcessLib
24 {
25 namespace RichardsMechanics
26 {
27 template <int DisplacementDim>
28 Eigen::Matrix<double, DisplacementDim, DisplacementDim> intrinsicPermeability(
29  double const t, SpatialPosition const& x_position, int const material_id,
31 {
32  const Eigen::MatrixXd& permeability =
33  material.getPermeability(material_id, t, x_position, DisplacementDim);
34  if (permeability.rows() == DisplacementDim)
35  {
36  return permeability;
37  }
38  if (permeability.rows() == 1)
39  {
40  return Eigen::Matrix<double, DisplacementDim,
41  DisplacementDim>::Identity() *
42  permeability(0, 0);
43  }
44 
45  OGS_FATAL(
46  "Intrinsic permeability dimension is neither %d nor one, but %dx%d.",
47  DisplacementDim, permeability.rows(), permeability.cols());
48 }
49 
50 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
51  typename IntegrationMethod, int DisplacementDim>
52 RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
53  ShapeFunctionPressure, IntegrationMethod,
54  DisplacementDim>::
55  RichardsMechanicsLocalAssembler(
56  MeshLib::Element const& e,
57  std::size_t const /*local_matrix_size*/,
58  bool const is_axially_symmetric,
59  unsigned const integration_order,
61  : _process_data(process_data),
62  _integration_method(integration_order),
63  _element(e),
64  _is_axially_symmetric(is_axially_symmetric)
65 {
66  unsigned const n_integration_points =
67  _integration_method.getNumberOfPoints();
68 
69  _ip_data.reserve(n_integration_points);
70  _secondary_data.N_u.resize(n_integration_points);
71 
72  auto const shape_matrices_u =
73  initShapeMatrices<ShapeFunctionDisplacement,
74  ShapeMatricesTypeDisplacement, IntegrationMethod,
75  DisplacementDim>(e, is_axially_symmetric,
77 
78  auto const shape_matrices_p =
79  initShapeMatrices<ShapeFunctionPressure, ShapeMatricesTypePressure,
80  IntegrationMethod, DisplacementDim>(
81  e, is_axially_symmetric, _integration_method);
82 
83  auto const& solid_material =
85  _process_data.solid_materials, _process_data.material_ids,
86  e.getID());
87 
88  for (unsigned ip = 0; ip < n_integration_points; ip++)
89  {
90  _ip_data.emplace_back(solid_material);
91  auto& ip_data = _ip_data[ip];
92  auto const& sm_u = shape_matrices_u[ip];
93  _ip_data[ip].integration_weight =
94  _integration_method.getWeightedPoint(ip).getWeight() *
95  sm_u.integralMeasure * sm_u.detJ;
96 
97  ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
98  DisplacementDim, displacement_size>::Zero(DisplacementDim,
100  for (int i = 0; i < DisplacementDim; ++i)
101  ip_data.N_u_op
102  .template block<1, displacement_size / DisplacementDim>(
103  i, i * displacement_size / DisplacementDim)
104  .noalias() = sm_u.N;
105 
106  ip_data.N_u = sm_u.N;
107  ip_data.dNdx_u = sm_u.dNdx;
108 
109  ip_data.N_p = shape_matrices_p[ip].N;
110  ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
111 
112  _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
113  }
114 }
115 
116 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
117  typename IntegrationMethod, int DisplacementDim>
119  ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
120  DisplacementDim>::assemble(double const t,
121  std::vector<double> const& local_x,
122  std::vector<double>& local_M_data,
123  std::vector<double>& local_K_data,
124  std::vector<double>& local_rhs_data)
125 {
126  assert(local_x.size() == pressure_size + displacement_size);
127 
128  auto p_L =
129  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
130  pressure_size> const>(local_x.data() + pressure_index,
131  pressure_size);
132 
133  auto u =
134  Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
135  displacement_size> const>(local_x.data() + displacement_index,
137 
139  typename ShapeMatricesTypeDisplacement::template MatrixType<
141  displacement_size + pressure_size>>(
142  local_K_data, displacement_size + pressure_size,
144 
146  typename ShapeMatricesTypeDisplacement::template MatrixType<
148  displacement_size + pressure_size>>(
149  local_M_data, displacement_size + pressure_size,
151 
152  auto rhs = MathLib::createZeroedVector<
153  typename ShapeMatricesTypeDisplacement::template VectorType<
154  displacement_size + pressure_size>>(
155  local_rhs_data, displacement_size + pressure_size);
156 
157  double const& dt = _process_data.dt;
158  auto const material_id =
159  _process_data.flow_material->getMaterialID(_element.getID());
160 
161  SpatialPosition x_position;
162  x_position.setElementID(_element.getID());
163 
164  unsigned const n_integration_points =
165  _integration_method.getNumberOfPoints();
166  for (unsigned ip = 0; ip < n_integration_points; ip++)
167  {
168  x_position.setIntegrationPoint(ip);
169  auto const& w = _ip_data[ip].integration_weight;
170 
171  auto const& N_u_op = _ip_data[ip].N_u_op;
172 
173  auto const& N_u = _ip_data[ip].N_u;
174  auto const& dNdx_u = _ip_data[ip].dNdx_u;
175 
176  auto const& N_p = _ip_data[ip].N_p;
177  auto const& dNdx_p = _ip_data[ip].dNdx_p;
178 
179  auto const x_coord =
180  interpolateXCoordinate<ShapeFunctionDisplacement,
182  N_u);
183  auto const B =
184  LinearBMatrix::computeBMatrix<DisplacementDim,
185  ShapeFunctionDisplacement::NPOINTS,
186  typename BMatricesType::BMatrixType>(
187  dNdx_u, N_u, x_coord, _is_axially_symmetric);
188 
189  double p_cap_ip;
190  NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
191 
192  auto& eps = _ip_data[ip].eps;
193  auto& S_L = _ip_data[ip].saturation;
194 
195  auto const alpha = _process_data.biot_coefficient(t, x_position)[0];
196  auto const rho_SR = _process_data.solid_density(t, x_position)[0];
197  auto const K_SR = _process_data.solid_bulk_modulus(t, x_position)[0];
198  auto const K_LR = _process_data.fluid_bulk_modulus(t, x_position)[0];
199  auto const temperature = _process_data.temperature(t, x_position)[0];
200  auto const porosity = _process_data.flow_material->getPorosity(
201  material_id, t, x_position, -p_cap_ip, temperature, p_cap_ip);
202  auto const rho_LR = _process_data.flow_material->getFluidDensity(
203  -p_cap_ip, temperature);
204  auto const& b = _process_data.specific_body_force;
205  auto const& identity2 = MathLib::KelvinVector::Invariants<
207  DisplacementDim>::value>::identity2;
208 
209  S_L = _process_data.flow_material->getSaturation(
210  material_id, t, x_position, -p_cap_ip, temperature, p_cap_ip);
211 
212  double const dS_L_dp_cap =
213  _process_data.flow_material->getSaturationDerivative(
214  material_id, t, x_position, -p_cap_ip, temperature, S_L);
215 
216  double const k_rel =
217  _process_data.flow_material->getRelativePermeability(
218  t, x_position, -p_cap_ip, temperature, S_L);
219 
220  double const mu = _process_data.flow_material->getFluidViscosity(
221  -p_cap_ip, temperature);
222 
223  GlobalDimMatrixType const rho_K_over_mu =
224  intrinsicPermeability<DisplacementDim>(
225  t, x_position, material_id, *_process_data.flow_material) *
226  (rho_LR * k_rel / mu);
227 
228  //
229  // displacement equation, displacement part
230  //
231  eps.noalias() = B * u;
232 
233  auto C = _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
234  temperature);
235 
236  //
237  // displacement equation, displacement part
238  //
239  K.template block<displacement_size, displacement_size>(
241  .noalias() += B.transpose() * C * B * w;
242 
243  double const rho = rho_SR * (1 - porosity) + porosity * rho_LR;
244  rhs.template segment<displacement_size>(displacement_index).noalias() +=
245  N_u_op.transpose() * rho * b * w;
246 
247  //
248  // pressure equation, pressure part.
249  //
250  double const a0 = S_L * (alpha - porosity) / K_SR;
251  // Volumetric average specific storage of the solid and fluid phases.
252  double const specific_storage =
253  dS_L_dp_cap * (p_cap_ip * a0 - porosity) +
254  S_L * (porosity / K_LR + a0);
255  M.template block<pressure_size, pressure_size>(pressure_index,
257  .noalias() += N_p.transpose() * rho_LR * specific_storage * N_p * w;
258 
259  K.template block<pressure_size, pressure_size>(pressure_index,
261  .noalias() += dNdx_p.transpose() * rho_K_over_mu * dNdx_p * w;
262 
263  rhs.template segment<pressure_size>(pressure_index).noalias() +=
264  dNdx_p.transpose() * rho_LR * rho_K_over_mu * b * w;
265 
266  //
267  // displacement equation, pressure part
268  //
269  K.template block<displacement_size, pressure_size>(displacement_index,
271  .noalias() -= B.transpose() * alpha * S_L * identity2 * N_p * w;
272 
273  //
274  // pressure equation, displacement part.
275  //
276  M.template block<pressure_size, displacement_size>(pressure_index,
278  .noalias() += N_p.transpose() * S_L * rho_LR * alpha *
279  identity2.transpose() * B * w;
280  }
281 }
282 
283 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
284  typename IntegrationMethod, int DisplacementDim>
285 void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
286  ShapeFunctionPressure, IntegrationMethod,
287  DisplacementDim>::
288  assembleWithJacobian(double const t, std::vector<double> const& local_x,
289  std::vector<double> const& local_xdot,
290  const double /*dxdot_dx*/, const double /*dx_dx*/,
291  std::vector<double>& /*local_M_data*/,
292  std::vector<double>& /*local_K_data*/,
293  std::vector<double>& local_rhs_data,
294  std::vector<double>& local_Jac_data)
295 {
296  assert(local_x.size() == pressure_size + displacement_size);
297 
298  auto p_L =
299  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
300  pressure_size> const>(local_x.data() + pressure_index,
301  pressure_size);
302 
303  auto u =
304  Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
305  displacement_size> const>(local_x.data() + displacement_index,
307 
308  auto p_L_dot =
309  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
310  pressure_size> const>(local_xdot.data() + pressure_index,
311  pressure_size);
312  auto u_dot =
313  Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
314  displacement_size> const>(local_xdot.data() + displacement_index,
316 
317  auto local_Jac = MathLib::createZeroedMatrix<
318  typename ShapeMatricesTypeDisplacement::template MatrixType<
320  displacement_size + pressure_size>>(
321  local_Jac_data, displacement_size + pressure_size,
323 
324  auto local_rhs = MathLib::createZeroedVector<
325  typename ShapeMatricesTypeDisplacement::template VectorType<
326  displacement_size + pressure_size>>(
327  local_rhs_data, displacement_size + pressure_size);
328 
329  typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
330  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
331  pressure_size);
332 
333  typename ShapeMatricesTypePressure::NodalMatrixType storage_p =
334  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
335  pressure_size);
336 
337  typename ShapeMatricesTypeDisplacement::template MatrixType<
338  displacement_size, pressure_size>
339  Kup = ShapeMatricesTypeDisplacement::template MatrixType<
340  displacement_size, pressure_size>::Zero(displacement_size,
341  pressure_size);
342 
343  typename ShapeMatricesTypeDisplacement::template MatrixType<
344  pressure_size, displacement_size>
345  Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
346  pressure_size, displacement_size>::Zero(pressure_size,
347  displacement_size);
348 
349  double const& dt = _process_data.dt;
350  auto const material_id =
351  _process_data.flow_material->getMaterialID(_element.getID());
352 
353  SpatialPosition x_position;
354  x_position.setElementID(_element.getID());
355 
356  unsigned const n_integration_points =
357  _integration_method.getNumberOfPoints();
358  for (unsigned ip = 0; ip < n_integration_points; ip++)
359  {
360  x_position.setIntegrationPoint(ip);
361  auto const& w = _ip_data[ip].integration_weight;
362 
363  auto const& N_u_op = _ip_data[ip].N_u_op;
364 
365  auto const& N_u = _ip_data[ip].N_u;
366  auto const& dNdx_u = _ip_data[ip].dNdx_u;
367 
368  auto const& N_p = _ip_data[ip].N_p;
369  auto const& dNdx_p = _ip_data[ip].dNdx_p;
370 
371  auto const x_coord =
372  interpolateXCoordinate<ShapeFunctionDisplacement,
374  N_u);
375  auto const B =
376  LinearBMatrix::computeBMatrix<DisplacementDim,
377  ShapeFunctionDisplacement::NPOINTS,
378  typename BMatricesType::BMatrixType>(
379  dNdx_u, N_u, x_coord, _is_axially_symmetric);
380 
381  double p_cap_ip;
382  NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
383 
384  double p_cap_dot_ip;
385  NumLib::shapeFunctionInterpolate(-p_L_dot, N_p, p_cap_dot_ip);
386 
388  DisplacementDim);
389  for (int i = 0; i < u_ip.size(); ++i)
390  {
392  u.segment(i * ShapeFunctionDisplacement::NPOINTS,
393  ShapeFunctionDisplacement::NPOINTS),
394  N_u, u_ip.coeffRef(i));
395  }
396 
397  auto& eps = _ip_data[ip].eps;
398  auto& S_L = _ip_data[ip].saturation;
399  auto const& sigma_eff = _ip_data[ip].sigma_eff;
400 
401  auto const alpha = _process_data.biot_coefficient(t, x_position)[0];
402  auto const rho_SR = _process_data.solid_density(t, x_position)[0];
403  auto const K_SR = _process_data.solid_bulk_modulus(t, x_position)[0];
404  auto const K_LR = _process_data.fluid_bulk_modulus(t, x_position)[0];
405  auto const temperature = _process_data.temperature(t, x_position)[0];
406 
407  auto const porosity = _process_data.flow_material->getPorosity(
408  material_id, t, x_position, -p_cap_ip, temperature, p_cap_ip);
409  auto const rho_LR = _process_data.flow_material->getFluidDensity(
410  -p_cap_ip, temperature);
411  auto const& b = _process_data.specific_body_force;
412  auto const& identity2 = MathLib::KelvinVector::Invariants<
414  DisplacementDim>::value>::identity2;
415 
416  S_L = _process_data.flow_material->getSaturation(
417  material_id, t, x_position, -p_cap_ip, temperature, p_cap_ip);
418 
419  double const dS_L_dp_cap =
420  _process_data.flow_material->getSaturationDerivative(
421  material_id, t, x_position, -p_cap_ip, temperature, S_L);
422 
423  double const d2S_L_dp_cap_2 =
424  _process_data.flow_material->getSaturationDerivative2(
425  material_id, t, x_position, -p_cap_ip, temperature, S_L);
426 
427  double const k_rel =
428  _process_data.flow_material->getRelativePermeability(
429  t, x_position, -p_cap_ip, temperature, S_L);
430 
431  double const mu = _process_data.flow_material->getFluidViscosity(
432  -p_cap_ip, temperature);
433 
434  GlobalDimMatrixType const rho_Ki_over_mu =
435  intrinsicPermeability<DisplacementDim>(
436  t, x_position, material_id, *_process_data.flow_material) *
437  (rho_LR / mu);
438 
439  //
440  // displacement equation, displacement part
441  //
442  eps.noalias() = B * u;
443 
444  auto C = _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
445  temperature);
446 
447  local_Jac
448  .template block<displacement_size, displacement_size>(
450  .noalias() += B.transpose() * C * B * w;
451 
452  double const rho = rho_SR * (1 - porosity) + porosity * rho_LR;
453  local_rhs.template segment<displacement_size>(displacement_index)
454  .noalias() -=
455  (B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w;
456 
457  //
458  // displacement equation, pressure part
459  //
460  Kup.noalias() += B.transpose() * alpha * S_L * identity2 * N_p * w;
461 
462  /* For future implementation including swelling.
463  double const dsigma_eff_dp_cap = -K_intrinsic * m_swell * n *
464  std::pow(S_L, n - 1) * dS_L_dp_cap *
465  identity2;
466  local_Jac
467  .template block<displacement_size, pressure_size>(
468  displacement_index, pressure_index)
469  .noalias() -= B.transpose() * dsigma_eff_dp_cap * N_p * w;
470  */
471 
472  local_Jac
473  .template block<displacement_size, pressure_size>(
475  .noalias() -= B.transpose() * alpha *
476  (S_L + p_cap_ip * dS_L_dp_cap) * identity2 * N_p * w;
477 
478  local_Jac
479  .template block<displacement_size, pressure_size>(
481  .noalias() +=
482  N_u_op.transpose() * porosity * rho_LR * dS_L_dp_cap * b * N_p * w;
483  //
484  // pressure equation, displacement part.
485  //
486  Kpu.noalias() += N_p.transpose() * S_L * rho_LR * alpha *
487  identity2.transpose() * B * w;
488 
489  //
490  // pressure equation, pressure part.
491  //
492  laplace_p.noalias() +=
493  dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
494 
495  double const a0 = (alpha - porosity) / K_SR;
496  double const specific_storage =
497  dS_L_dp_cap * (p_cap_ip * S_L * a0 - porosity) +
498  S_L * (porosity / K_LR + S_L * a0);
499 
500  double const dspecific_storage_dp_cap =
501  d2S_L_dp_cap_2 * (p_cap_ip * S_L * a0 - porosity) +
502  dS_L_dp_cap *
503  (porosity / K_LR + a0 * 3 * S_L + dS_L_dp_cap * p_cap_ip * a0);
504 
505  storage_p.noalias() +=
506  N_p.transpose() * rho_LR * specific_storage * N_p * w;
507 
508  local_Jac
509  .template block<pressure_size, pressure_size>(pressure_index,
511  .noalias() += N_p.transpose() * rho_LR * p_cap_dot_ip *
512  dspecific_storage_dp_cap * N_p * w;
513 
514  /* In the derivation there is a div(du/dt) term in the Jacobian, but
515  * this implementation increases the total runtime by 1%. Maybe a very
516  * large step is needed to see the increase of efficiency.
517  double div_u_dot = 0;
518  for (int i = 0; i < DisplacementDim; ++i)
519  {
520  div_u_dot +=
521  (dNdx_u *
522  u_dot.template segment<ShapeFunctionDisplacement::NPOINTS>(
523  i * ShapeFunctionDisplacement::NPOINTS))[i];
524  }
525  local_Jac
526  .template block<pressure_size, pressure_size>(pressure_index,
527  pressure_index)
528  .noalias() -= N_p.transpose() * rho_LR * dS_L_dp_cap * alpha *
529  div_u_dot * N_p * w;
530  */
531 
532  double const dk_rel_dS_l =
533  _process_data.flow_material->getRelativePermeabilityDerivative(
534  t, x_position, -p_cap_ip, temperature, S_L);
536  grad_p_cap = -dNdx_p * p_L;
537  local_Jac
538  .template block<pressure_size, pressure_size>(pressure_index,
540  .noalias() += dNdx_p.transpose() * rho_Ki_over_mu * grad_p_cap *
541  dk_rel_dS_l * dS_L_dp_cap * N_p * w;
542 
543  local_Jac
544  .template block<pressure_size, pressure_size>(pressure_index,
546  .noalias() += dNdx_p.transpose() * rho_LR * rho_Ki_over_mu * b *
547  dk_rel_dS_l * dS_L_dp_cap * N_p * w;
548 
549  local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
550  dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
551  }
552 
553  // pressure equation, pressure part.
554  local_Jac
555  .template block<pressure_size, pressure_size>(pressure_index,
557  .noalias() += laplace_p + storage_p / dt;
558 
559  // pressure equation, displacement part.
560  local_Jac
561  .template block<pressure_size, displacement_size>(pressure_index,
563  .noalias() = Kpu / dt;
564 
565  // pressure equation
566  local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
567  laplace_p * p_L + storage_p * p_L_dot + Kpu * u_dot;
568 
569  // displacement equation
570  local_rhs.template segment<displacement_size>(displacement_index)
571  .noalias() += Kup * p_L;
572 }
573 
574 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
575  typename IntegrationMethod, int DisplacementDim>
576 std::vector<double> const& RichardsMechanicsLocalAssembler<
577  ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
578  DisplacementDim>::
579  getIntPtSigma(const double /*t*/,
580  GlobalVector const& /*current_solution*/,
581  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
582  std::vector<double>& cache) const
583 {
584  static const int kelvin_vector_size =
586  auto const num_intpts = _ip_data.size();
587 
588  cache.clear();
589  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
590  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
591  cache, kelvin_vector_size, num_intpts);
592 
593  for (unsigned ip = 0; ip < num_intpts; ++ip)
594  {
595  auto const& sigma = _ip_data[ip].sigma_eff;
596  cache_mat.col(ip) =
598  }
599 
600  return cache;
601 }
602 
603 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
604  typename IntegrationMethod, int DisplacementDim>
605 std::vector<double> const& RichardsMechanicsLocalAssembler<
606  ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
607  DisplacementDim>::
608  getIntPtEpsilon(const double /*t*/,
609  GlobalVector const& /*current_solution*/,
610  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
611  std::vector<double>& cache) const
612 {
613  auto const kelvin_vector_size =
615  auto const num_intpts = _ip_data.size();
616 
617  cache.clear();
618  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
619  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
620  cache, kelvin_vector_size, num_intpts);
621 
622  for (unsigned ip = 0; ip < num_intpts; ++ip)
623  {
624  auto const& eps = _ip_data[ip].eps;
625  cache_mat.col(ip) =
627  }
628 
629  return cache;
630 }
631 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
632  typename IntegrationMethod, int DisplacementDim>
633 std::vector<double> const& RichardsMechanicsLocalAssembler<
634  ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
635  DisplacementDim>::getIntPtDarcyVelocity(const double t,
636  GlobalVector const&
637  current_solution,
639  dof_table,
640  std::vector<double>& cache) const
641 {
642  auto const num_intpts = _ip_data.size();
643 
644  auto const indices = NumLib::getIndices(_element.getID(), dof_table);
645  assert(!indices.empty());
646  auto const local_x = current_solution.get(indices);
647 
648  cache.clear();
649  auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
650  double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
651  cache, DisplacementDim, num_intpts);
652 
653  auto p_L =
654  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
655  pressure_size> const>(local_x.data() + pressure_index,
656  pressure_size);
657 
658  auto const material_id =
659  _process_data.flow_material->getMaterialID(_element.getID());
660 
661  unsigned const n_integration_points =
662  _integration_method.getNumberOfPoints();
663 
664  SpatialPosition x_position;
665  x_position.setElementID(_element.getID());
666  for (unsigned ip = 0; ip < n_integration_points; ip++)
667  {
668  x_position.setIntegrationPoint(ip);
669  auto const& N_p = _ip_data[ip].N_p;
670 
671  double p_cap_ip;
672  NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
673 
674  auto const temperature = _process_data.temperature(t, x_position)[0];
675  GlobalDimMatrixType const K_over_mu =
676  intrinsicPermeability<DisplacementDim>(
677  t, x_position, material_id, *_process_data.flow_material) /
678  _process_data.flow_material->getFluidViscosity(-p_cap_ip,
679  temperature);
680  auto const rho_LR = _process_data.flow_material->getFluidDensity(
681  -p_cap_ip, temperature);
682  auto const& b = _process_data.specific_body_force;
683 
684  // Compute the velocity
685  auto const& dNdx_p = _ip_data[ip].dNdx_p;
686  cache_matrix.col(ip).noalias() =
687  -K_over_mu * dNdx_p * p_L - rho_LR * K_over_mu * b;
688  }
689 
690  return cache;
691 }
692 
693 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
694  typename IntegrationMethod, int DisplacementDim>
695 std::vector<double> const& RichardsMechanicsLocalAssembler<
696  ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
697  DisplacementDim>::
698  getIntPtSaturation(const double /*t*/,
699  GlobalVector const& /*current_solution*/,
700  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
701  std::vector<double>& cache) const
702 {
703  auto const num_intpts = _ip_data.size();
704 
705  cache.clear();
706  auto cache_mat = MathLib::createZeroedMatrix<
707  Eigen::Matrix<double, 1, Eigen::Dynamic, Eigen::RowMajor>>(cache, 1,
708  num_intpts);
709 
710  for (unsigned ip = 0; ip < num_intpts; ++ip)
711  {
712  cache_mat[ip] = _ip_data[ip].saturation;
713  }
714 
715  return cache;
716 }
717 
718 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
719  typename IntegrationMethod, int DisplacementDim>
720 void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
721  ShapeFunctionPressure, IntegrationMethod,
722  DisplacementDim>::
724  const double /*t*/, const std::vector<double>& /*local_xdot*/,
725  const double /*dxdot_dx*/, const double /*dx_dx*/,
726  std::vector<double>& /*local_M_data*/,
727  std::vector<double>& /*local_K_data*/,
728  std::vector<double>& /*local_b_data*/,
729  std::vector<double>& /*local_Jac_data*/,
730  const LocalCoupledSolutions& /*local_coupled_solutions*/)
731 {
732  OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
733 }
734 
735 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
736  typename IntegrationMethod, int DisplacementDim>
737 void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
738  ShapeFunctionPressure, IntegrationMethod,
739  DisplacementDim>::
741  const double /*t*/, const std::vector<double>& /*local_xdot*/,
742  const double /*dxdot_dx*/, const double /*dx_dx*/,
743  std::vector<double>& /*local_M_data*/,
744  std::vector<double>& /*local_K_data*/,
745  std::vector<double>& /*local_b_data*/,
746  std::vector<double>& /*local_Jac_data*/,
747  const LocalCoupledSolutions& /*local_coupled_solutions*/)
748 {
749  OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
750 }
751 
752 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
753  typename IntegrationMethod, int DisplacementDim>
754 void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
755  ShapeFunctionPressure, IntegrationMethod,
756  DisplacementDim>::
758  const double t,
759  const std::vector<double>& local_xdot,
760  const double dxdot_dx,
761  const double dx_dx,
762  std::vector<double>& local_M_data,
763  std::vector<double>& local_K_data,
764  std::vector<double>& local_b_data,
765  std::vector<double>& local_Jac_data,
766  const LocalCoupledSolutions& local_coupled_solutions)
767 {
768  // For the equations with pressure
769  if (local_coupled_solutions.process_id == 0)
770  {
772  t, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
773  local_b_data, local_Jac_data, local_coupled_solutions);
774  return;
775  }
776 
777  // For the equations with deformation
779  t, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
780  local_b_data, local_Jac_data, local_coupled_solutions);
781 }
782 
783 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
784  typename IntegrationMethod, int DisplacementDim>
785 void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
786  ShapeFunctionPressure, IntegrationMethod,
787  DisplacementDim>::
788  postNonLinearSolverConcrete(std::vector<double> const& local_x,
789  double const t,
790  bool const use_monolithic_scheme)
791 {
792  const int displacement_offset =
793  use_monolithic_scheme ? displacement_index : 0;
794 
795  auto u =
796  Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
797  displacement_size> const>(local_x.data() + displacement_offset,
799  double const& dt = _process_data.dt;
800  SpatialPosition x_position;
801  x_position.setElementID(_element.getID());
802 
803  int const n_integration_points = _integration_method.getNumberOfPoints();
804  for (int ip = 0; ip < n_integration_points; ip++)
805  {
806  x_position.setIntegrationPoint(ip);
807  auto const& N_u = _ip_data[ip].N_u;
808  auto const& dNdx_u = _ip_data[ip].dNdx_u;
809  auto const temperature = _process_data.temperature(t, x_position)[0];
810 
811  auto const x_coord =
812  interpolateXCoordinate<ShapeFunctionDisplacement,
814  N_u);
815  auto const B =
816  LinearBMatrix::computeBMatrix<DisplacementDim,
817  ShapeFunctionDisplacement::NPOINTS,
818  typename BMatricesType::BMatrixType>(
819  dNdx_u, N_u, x_coord, _is_axially_symmetric);
820 
821  auto& eps = _ip_data[ip].eps;
822  eps.noalias() = B * u;
823 
824  _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
825  temperature);
826  }
827 }
828 
829 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
830  typename IntegrationMethod, int DisplacementDim>
831 void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
832  ShapeFunctionPressure, IntegrationMethod,
833  DisplacementDim>::
835  std::vector<double> const& local_x)
836 {
837  auto p_L =
838  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
839  pressure_size> const>(local_x.data() + pressure_index,
840  pressure_size);
841  auto const material_id =
842  _process_data.flow_material->getMaterialID(_element.getID());
843 
844  SpatialPosition x_position;
845  x_position.setElementID(_element.getID());
846 
847  unsigned const n_integration_points =
848  _integration_method.getNumberOfPoints();
849 
850  double saturation_avg = 0;
851 
852  for (unsigned ip = 0; ip < n_integration_points; ip++)
853  {
854  x_position.setIntegrationPoint(ip);
855  auto const& N_p = _ip_data[ip].N_p;
856  auto const temperature = _process_data.temperature(t, x_position)[0];
857 
858  double p_cap_ip;
859  NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
860  auto& S_L = _ip_data[ip].saturation;
861  S_L = _process_data.flow_material->getSaturation(
862  material_id, t, x_position, -p_cap_ip, temperature, p_cap_ip);
863  saturation_avg += S_L;
864  }
865  saturation_avg /= n_integration_points;
866 
867  (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
868 
870  ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
871  DisplacementDim>(_element, _is_axially_symmetric, p_L,
872  *_process_data.pressure_interpolated);
873 }
874 
875 } // namespace RichardsMechanics
876 } // namespace ProcessLib
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
void assembleWithJacobianForDeformationEquations(double const t, std::vector< double > const &local_xdot, const double dxdot_dx, const double dx_dx, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data, LocalCoupledSolutions const &local_coupled_solutions)
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:23
std::vector< double > const & getIntPtSigma(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
void assembleWithJacobian(double const t, std::vector< double > const &local_x, std::vector< double > const &local_xdot, const double, const double, std::vector< double > &, std::vector< double > &, std::vector< double > &local_rhs_data, std::vector< double > &local_Jac_data) override
Eigen::Matrix< double, DisplacementDim, DisplacementDim > intrinsicPermeability(double const t, SpatialPosition const &x_position, int const material_id, RichardsFlow::RichardsFlowMaterialProperties const &material)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool is_axially_symmetric, IntegrationMethod const &integration_method)
void computeSecondaryVariableConcrete(double const t, std::vector< double > const &local_x) override
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)
void setElementID(std::size_t element_id)
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
void postNonLinearSolverConcrete(std::vector< double > const &local_x, double const t, bool const use_monolithic_scheme) override
void shapeFunctionInterpolate(const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
Definition: Interpolation.h:73
Eigen::MatrixXd const & getPermeability(const int material_id, const double t, const ProcessLib::SpatialPosition &pos, const int dim) const
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:31
std::vector< double > const & getIntPtDarcyVelocity(const double t, GlobalVector const &current_solution, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< double > &cache) const override
void assemble(double const t, std::vector< double > const &local_x, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_rhs_data) override
std::vector< double > const & getIntPtSaturation(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &) const override
void assembleWithJacobianForPressureEquations(double const t, std::vector< double > const &local_xdot, const double dxdot_dx, const double dx_dx, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data, LocalCoupledSolutions const &local_coupled_solutions)
void assembleWithJacobianForStaggeredScheme(double const t, std::vector< double > const &local_xdot, const double dxdot_dx, const double dx_dx, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data, LocalCoupledSolutions const &local_coupled_solutions) override
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:41
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
typename ShapeMatricesTypePressure::GlobalDimMatrixType GlobalDimMatrixType
VectorType< GlobalDim > GlobalDimVectorType
std::vector< double > const & getIntPtEpsilon(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:90
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
#define OGS_FATAL(fmt,...)
Definition: Error.h:71
RichardsMechanicsProcessData< DisplacementDim > & _process_data
void setIntegrationPoint(unsigned integration_point)
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition: BMatrixPolicy.h:54