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