OGS 6.2.1-97-g73d1aeda3
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) + S_L * 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  if (_process_data.apply_mass_lumping)
286  {
287  auto Mpp = M.template block<pressure_size, pressure_size>(
289  Mpp = Mpp.colwise().sum().eval().asDiagonal();
290  }
291 }
292 
293 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
294  typename IntegrationMethod, int DisplacementDim>
295 void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
296  ShapeFunctionPressure, IntegrationMethod,
297  DisplacementDim>::
298  assembleWithJacobian(double const t, std::vector<double> const& local_x,
299  std::vector<double> const& local_xdot,
300  const double /*dxdot_dx*/, const double /*dx_dx*/,
301  std::vector<double>& /*local_M_data*/,
302  std::vector<double>& /*local_K_data*/,
303  std::vector<double>& local_rhs_data,
304  std::vector<double>& local_Jac_data)
305 {
306  assert(local_x.size() == pressure_size + displacement_size);
307 
308  auto p_L =
309  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
310  pressure_size> const>(local_x.data() + pressure_index,
311  pressure_size);
312 
313  auto u =
314  Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
315  displacement_size> const>(local_x.data() + displacement_index,
317 
318  auto p_L_dot =
319  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
320  pressure_size> const>(local_xdot.data() + pressure_index,
321  pressure_size);
322  auto u_dot =
323  Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
324  displacement_size> const>(local_xdot.data() + displacement_index,
326 
327  auto local_Jac = MathLib::createZeroedMatrix<
328  typename ShapeMatricesTypeDisplacement::template MatrixType<
330  displacement_size + pressure_size>>(
331  local_Jac_data, displacement_size + pressure_size,
333 
334  auto local_rhs = MathLib::createZeroedVector<
335  typename ShapeMatricesTypeDisplacement::template VectorType<
336  displacement_size + pressure_size>>(
337  local_rhs_data, displacement_size + pressure_size);
338 
339  typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
340  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
341  pressure_size);
342 
343  typename ShapeMatricesTypePressure::NodalMatrixType storage_p =
344  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
345  pressure_size);
346 
347  typename ShapeMatricesTypeDisplacement::template MatrixType<
348  displacement_size, pressure_size>
349  Kup = ShapeMatricesTypeDisplacement::template MatrixType<
350  displacement_size, pressure_size>::Zero(displacement_size,
351  pressure_size);
352 
353  typename ShapeMatricesTypeDisplacement::template MatrixType<
354  pressure_size, displacement_size>
355  Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
356  pressure_size, displacement_size>::Zero(pressure_size,
357  displacement_size);
358 
359  double const& dt = _process_data.dt;
360  auto const material_id =
361  _process_data.flow_material->getMaterialID(_element.getID());
362 
364  x_position.setElementID(_element.getID());
365 
366  unsigned const n_integration_points =
367  _integration_method.getNumberOfPoints();
368  for (unsigned ip = 0; ip < n_integration_points; ip++)
369  {
370  x_position.setIntegrationPoint(ip);
371  auto const& w = _ip_data[ip].integration_weight;
372 
373  auto const& N_u_op = _ip_data[ip].N_u_op;
374 
375  auto const& N_u = _ip_data[ip].N_u;
376  auto const& dNdx_u = _ip_data[ip].dNdx_u;
377 
378  auto const& N_p = _ip_data[ip].N_p;
379  auto const& dNdx_p = _ip_data[ip].dNdx_p;
380 
381  auto const x_coord =
382  interpolateXCoordinate<ShapeFunctionDisplacement,
384  N_u);
385  auto const B =
386  LinearBMatrix::computeBMatrix<DisplacementDim,
387  ShapeFunctionDisplacement::NPOINTS,
388  typename BMatricesType::BMatrixType>(
389  dNdx_u, N_u, x_coord, _is_axially_symmetric);
390 
391  double p_cap_ip;
392  NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
393 
394  double p_cap_dot_ip;
395  NumLib::shapeFunctionInterpolate(-p_L_dot, N_p, p_cap_dot_ip);
396 
398  DisplacementDim);
399  for (int i = 0; i < u_ip.size(); ++i)
400  {
402  u.segment(i * ShapeFunctionDisplacement::NPOINTS,
403  ShapeFunctionDisplacement::NPOINTS),
404  N_u, u_ip.coeffRef(i));
405  }
406 
407  auto& eps = _ip_data[ip].eps;
408  auto& S_L = _ip_data[ip].saturation;
409  auto const& sigma_eff = _ip_data[ip].sigma_eff;
410 
411  auto const alpha = _process_data.biot_coefficient(t, x_position)[0];
412  auto const rho_SR = _process_data.solid_density(t, x_position)[0];
413  auto const K_SR = _process_data.solid_bulk_modulus(t, x_position)[0];
414  auto const K_LR = _process_data.fluid_bulk_modulus(t, x_position)[0];
415  auto const temperature = _process_data.temperature(t, x_position)[0];
416 
417  auto const porosity = _process_data.flow_material->getPorosity(
418  material_id, t, x_position, -p_cap_ip, temperature, p_cap_ip);
419  auto const rho_LR = _process_data.flow_material->getFluidDensity(
420  -p_cap_ip, temperature);
421  auto const& b = _process_data.specific_body_force;
422  auto const& identity2 = MathLib::KelvinVector::Invariants<
424  DisplacementDim>::value>::identity2;
425 
426  S_L = _process_data.flow_material->getSaturation(
427  material_id, t, x_position, -p_cap_ip, temperature, p_cap_ip);
428 
429  double const dS_L_dp_cap =
430  _process_data.flow_material->getSaturationDerivative(
431  material_id, t, x_position, -p_cap_ip, temperature, S_L);
432 
433  double const d2S_L_dp_cap_2 =
434  _process_data.flow_material->getSaturationDerivative2(
435  material_id, t, x_position, -p_cap_ip, temperature, S_L);
436 
437  double const k_rel =
438  _process_data.flow_material->getRelativePermeability(
439  t, x_position, -p_cap_ip, temperature, S_L);
440 
441  double const mu = _process_data.flow_material->getFluidViscosity(
442  -p_cap_ip, temperature);
443 
444  GlobalDimMatrixType const rho_Ki_over_mu =
445  intrinsicPermeability<DisplacementDim>(
446  t, x_position, material_id, *_process_data.flow_material) *
447  (rho_LR / mu);
448 
449  //
450  // displacement equation, displacement part
451  //
452  eps.noalias() = B * u;
453 
454  auto C = _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
455  temperature);
456 
457  local_Jac
458  .template block<displacement_size, displacement_size>(
460  .noalias() += B.transpose() * C * B * w;
461 
462  double const rho = rho_SR * (1 - porosity) + S_L * porosity * rho_LR;
463  local_rhs.template segment<displacement_size>(displacement_index)
464  .noalias() -=
465  (B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w;
466 
467  //
468  // displacement equation, pressure part
469  //
470  Kup.noalias() += B.transpose() * alpha * S_L * identity2 * N_p * w;
471 
472  /* For future implementation including swelling.
473  double const dsigma_eff_dp_cap = -K_intrinsic * m_swell * n *
474  std::pow(S_L, n - 1) * dS_L_dp_cap *
475  identity2;
476  local_Jac
477  .template block<displacement_size, pressure_size>(
478  displacement_index, pressure_index)
479  .noalias() -= B.transpose() * dsigma_eff_dp_cap * N_p * w;
480  */
481 
482  local_Jac
483  .template block<displacement_size, pressure_size>(
485  .noalias() -= B.transpose() * alpha *
486  (S_L + p_cap_ip * dS_L_dp_cap) * identity2 * N_p * w;
487 
488  local_Jac
489  .template block<displacement_size, pressure_size>(
491  .noalias() +=
492  N_u_op.transpose() * porosity * rho_LR * dS_L_dp_cap * b * N_p * w;
493  //
494  // pressure equation, displacement part.
495  //
496  Kpu.noalias() += N_p.transpose() * S_L * rho_LR * alpha *
497  identity2.transpose() * B * w;
498 
499  //
500  // pressure equation, pressure part.
501  //
502  laplace_p.noalias() +=
503  dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
504 
505  double const a0 = (alpha - porosity) / K_SR;
506  double const specific_storage =
507  dS_L_dp_cap * (p_cap_ip * S_L * a0 - porosity) +
508  S_L * (porosity / K_LR + S_L * a0);
509 
510  double const dspecific_storage_dp_cap =
511  d2S_L_dp_cap_2 * (p_cap_ip * S_L * a0 - porosity) +
512  dS_L_dp_cap *
513  (porosity / K_LR + a0 * 3 * S_L + dS_L_dp_cap * p_cap_ip * a0);
514 
515  storage_p.noalias() +=
516  N_p.transpose() * rho_LR * specific_storage * N_p * w;
517 
518  local_Jac
519  .template block<pressure_size, pressure_size>(pressure_index,
521  .noalias() += N_p.transpose() * rho_LR * p_cap_dot_ip *
522  dspecific_storage_dp_cap * N_p * w;
523 
524  /* In the derivation there is a div(du/dt) term in the Jacobian, but
525  * this implementation increases the total runtime by 1%. Maybe a very
526  * large step is needed to see the increase of efficiency.
527  double div_u_dot = 0;
528  for (int i = 0; i < DisplacementDim; ++i)
529  {
530  div_u_dot +=
531  (dNdx_u *
532  u_dot.template segment<ShapeFunctionDisplacement::NPOINTS>(
533  i * ShapeFunctionDisplacement::NPOINTS))[i];
534  }
535  local_Jac
536  .template block<pressure_size, pressure_size>(pressure_index,
537  pressure_index)
538  .noalias() -= N_p.transpose() * rho_LR * dS_L_dp_cap * alpha *
539  div_u_dot * N_p * w;
540  */
541 
542  double const dk_rel_dS_l =
543  _process_data.flow_material->getRelativePermeabilityDerivative(
544  t, x_position, -p_cap_ip, temperature, S_L);
546  grad_p_cap = -dNdx_p * p_L;
547  local_Jac
548  .template block<pressure_size, pressure_size>(pressure_index,
550  .noalias() += dNdx_p.transpose() * rho_Ki_over_mu * grad_p_cap *
551  dk_rel_dS_l * dS_L_dp_cap * N_p * w;
552 
553  local_Jac
554  .template block<pressure_size, pressure_size>(pressure_index,
556  .noalias() += dNdx_p.transpose() * rho_LR * rho_Ki_over_mu * b *
557  dk_rel_dS_l * dS_L_dp_cap * N_p * w;
558 
559  local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
560  dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
561  }
562 
563  if (_process_data.apply_mass_lumping)
564  {
565  storage_p = storage_p.colwise().sum().eval().asDiagonal();
566  }
567 
568  // pressure equation, pressure part.
569  local_Jac
570  .template block<pressure_size, pressure_size>(pressure_index,
572  .noalias() += laplace_p + storage_p / dt;
573 
574  // pressure equation, displacement part.
575  local_Jac
576  .template block<pressure_size, displacement_size>(pressure_index,
578  .noalias() = Kpu / dt;
579 
580  // pressure equation
581  local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
582  laplace_p * p_L + storage_p * p_L_dot + Kpu * u_dot;
583 
584  // displacement equation
585  local_rhs.template segment<displacement_size>(displacement_index)
586  .noalias() += Kup * p_L;
587 }
588 
589 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
590  typename IntegrationMethod, int DisplacementDim>
591 std::vector<double> const& RichardsMechanicsLocalAssembler<
592  ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
593  DisplacementDim>::
594  getIntPtSigma(const double /*t*/,
595  GlobalVector const& /*current_solution*/,
596  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
597  std::vector<double>& cache) const
598 {
599  static const int kelvin_vector_size =
601  auto const num_intpts = _ip_data.size();
602 
603  cache.clear();
604  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
605  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
606  cache, kelvin_vector_size, num_intpts);
607 
608  for (unsigned ip = 0; ip < num_intpts; ++ip)
609  {
610  auto const& sigma = _ip_data[ip].sigma_eff;
611  cache_mat.col(ip) =
613  }
614 
615  return cache;
616 }
617 
618 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
619  typename IntegrationMethod, int DisplacementDim>
620 std::vector<double> const& RichardsMechanicsLocalAssembler<
621  ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
622  DisplacementDim>::
623  getIntPtEpsilon(const double /*t*/,
624  GlobalVector const& /*current_solution*/,
625  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
626  std::vector<double>& cache) const
627 {
628  auto const kelvin_vector_size =
630  auto const num_intpts = _ip_data.size();
631 
632  cache.clear();
633  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
634  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
635  cache, kelvin_vector_size, num_intpts);
636 
637  for (unsigned ip = 0; ip < num_intpts; ++ip)
638  {
639  auto const& eps = _ip_data[ip].eps;
640  cache_mat.col(ip) =
642  }
643 
644  return cache;
645 }
646 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
647  typename IntegrationMethod, int DisplacementDim>
648 std::vector<double> const& RichardsMechanicsLocalAssembler<
649  ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
650  DisplacementDim>::getIntPtDarcyVelocity(const double t,
651  GlobalVector const&
652  current_solution,
654  dof_table,
655  std::vector<double>& cache) const
656 {
657  auto const num_intpts = _ip_data.size();
658 
659  auto const indices = NumLib::getIndices(_element.getID(), dof_table);
660  assert(!indices.empty());
661  auto const local_x = current_solution.get(indices);
662 
663  cache.clear();
664  auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
665  double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
666  cache, DisplacementDim, num_intpts);
667 
668  auto p_L =
669  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
670  pressure_size> const>(local_x.data() + pressure_index,
671  pressure_size);
672 
673  auto const material_id =
674  _process_data.flow_material->getMaterialID(_element.getID());
675 
676  unsigned const n_integration_points =
677  _integration_method.getNumberOfPoints();
678 
680  x_position.setElementID(_element.getID());
681  for (unsigned ip = 0; ip < n_integration_points; ip++)
682  {
683  x_position.setIntegrationPoint(ip);
684  auto const& N_p = _ip_data[ip].N_p;
685 
686  double p_cap_ip;
687  NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
688 
689  auto const temperature = _process_data.temperature(t, x_position)[0];
690  GlobalDimMatrixType const K_over_mu =
691  intrinsicPermeability<DisplacementDim>(
692  t, x_position, material_id, *_process_data.flow_material) /
693  _process_data.flow_material->getFluidViscosity(-p_cap_ip,
694  temperature);
695  auto const rho_LR = _process_data.flow_material->getFluidDensity(
696  -p_cap_ip, temperature);
697  auto const& b = _process_data.specific_body_force;
698 
699  // Compute the velocity
700  auto const& dNdx_p = _ip_data[ip].dNdx_p;
701  cache_matrix.col(ip).noalias() =
702  -K_over_mu * dNdx_p * p_L - rho_LR * K_over_mu * b;
703  }
704 
705  return cache;
706 }
707 
708 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
709  typename IntegrationMethod, int DisplacementDim>
710 std::vector<double> const& RichardsMechanicsLocalAssembler<
711  ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
712  DisplacementDim>::
713  getIntPtSaturation(const double /*t*/,
714  GlobalVector const& /*current_solution*/,
715  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
716  std::vector<double>& cache) const
717 {
718  auto const num_intpts = _ip_data.size();
719 
720  cache.clear();
721  auto cache_mat = MathLib::createZeroedMatrix<
722  Eigen::Matrix<double, 1, Eigen::Dynamic, Eigen::RowMajor>>(cache, 1,
723  num_intpts);
724 
725  for (unsigned ip = 0; ip < num_intpts; ++ip)
726  {
727  cache_mat[ip] = _ip_data[ip].saturation;
728  }
729 
730  return cache;
731 }
732 
733 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
734  typename IntegrationMethod, int DisplacementDim>
735 void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
736  ShapeFunctionPressure, IntegrationMethod,
737  DisplacementDim>::
739  const double /*t*/, const std::vector<double>& /*local_xdot*/,
740  const double /*dxdot_dx*/, const double /*dx_dx*/,
741  std::vector<double>& /*local_M_data*/,
742  std::vector<double>& /*local_K_data*/,
743  std::vector<double>& /*local_b_data*/,
744  std::vector<double>& /*local_Jac_data*/,
745  const LocalCoupledSolutions& /*local_coupled_solutions*/)
746 {
747  OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
748 }
749 
750 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
751  typename IntegrationMethod, int DisplacementDim>
752 void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
753  ShapeFunctionPressure, IntegrationMethod,
754  DisplacementDim>::
756  const double /*t*/, const std::vector<double>& /*local_xdot*/,
757  const double /*dxdot_dx*/, const double /*dx_dx*/,
758  std::vector<double>& /*local_M_data*/,
759  std::vector<double>& /*local_K_data*/,
760  std::vector<double>& /*local_b_data*/,
761  std::vector<double>& /*local_Jac_data*/,
762  const LocalCoupledSolutions& /*local_coupled_solutions*/)
763 {
764  OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
765 }
766 
767 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
768  typename IntegrationMethod, int DisplacementDim>
769 void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
770  ShapeFunctionPressure, IntegrationMethod,
771  DisplacementDim>::
773  const double t,
774  const std::vector<double>& local_xdot,
775  const double dxdot_dx,
776  const double dx_dx,
777  std::vector<double>& local_M_data,
778  std::vector<double>& local_K_data,
779  std::vector<double>& local_b_data,
780  std::vector<double>& local_Jac_data,
781  const LocalCoupledSolutions& local_coupled_solutions)
782 {
783  // For the equations with pressure
784  if (local_coupled_solutions.process_id == 0)
785  {
787  t, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
788  local_b_data, local_Jac_data, local_coupled_solutions);
789  return;
790  }
791 
792  // For the equations with deformation
794  t, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
795  local_b_data, local_Jac_data, local_coupled_solutions);
796 }
797 
798 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
799  typename IntegrationMethod, int DisplacementDim>
800 void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
801  ShapeFunctionPressure, IntegrationMethod,
802  DisplacementDim>::
803  postNonLinearSolverConcrete(std::vector<double> const& local_x,
804  double const t,
805  bool const use_monolithic_scheme)
806 {
807  const int displacement_offset =
808  use_monolithic_scheme ? displacement_index : 0;
809 
810  auto u =
811  Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
812  displacement_size> const>(local_x.data() + displacement_offset,
814  double const& dt = _process_data.dt;
816  x_position.setElementID(_element.getID());
817 
818  int const n_integration_points = _integration_method.getNumberOfPoints();
819  for (int ip = 0; ip < n_integration_points; ip++)
820  {
821  x_position.setIntegrationPoint(ip);
822  auto const& N_u = _ip_data[ip].N_u;
823  auto const& dNdx_u = _ip_data[ip].dNdx_u;
824  auto const temperature = _process_data.temperature(t, x_position)[0];
825 
826  auto const x_coord =
827  interpolateXCoordinate<ShapeFunctionDisplacement,
829  N_u);
830  auto const B =
831  LinearBMatrix::computeBMatrix<DisplacementDim,
832  ShapeFunctionDisplacement::NPOINTS,
833  typename BMatricesType::BMatrixType>(
834  dNdx_u, N_u, x_coord, _is_axially_symmetric);
835 
836  auto& eps = _ip_data[ip].eps;
837  eps.noalias() = B * u;
838 
839  _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
840  temperature);
841  }
842 }
843 
844 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
845  typename IntegrationMethod, int DisplacementDim>
846 void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
847  ShapeFunctionPressure, IntegrationMethod,
848  DisplacementDim>::
850  std::vector<double> const& local_x)
851 {
852  auto p_L =
853  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
854  pressure_size> const>(local_x.data() + pressure_index,
855  pressure_size);
856  auto const material_id =
857  _process_data.flow_material->getMaterialID(_element.getID());
858 
860  x_position.setElementID(_element.getID());
861 
862  unsigned const n_integration_points =
863  _integration_method.getNumberOfPoints();
864 
865  double saturation_avg = 0;
866 
867  for (unsigned ip = 0; ip < n_integration_points; ip++)
868  {
869  x_position.setIntegrationPoint(ip);
870  auto const& N_p = _ip_data[ip].N_p;
871  auto const temperature = _process_data.temperature(t, x_position)[0];
872 
873  double p_cap_ip;
874  NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
875  auto& S_L = _ip_data[ip].saturation;
876  S_L = _process_data.flow_material->getSaturation(
877  material_id, t, x_position, -p_cap_ip, temperature, p_cap_ip);
878  saturation_avg += S_L;
879  }
880  saturation_avg /= n_integration_points;
881 
882  (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
883 
885  ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
886  DisplacementDim>(_element, _is_axially_symmetric, p_L,
887  *_process_data.pressure_interpolated);
888 }
889 
890 } // namespace RichardsMechanics
891 } // 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
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)
Eigen::MatrixXd getPermeability(const int material_id, const double t, const ParameterLib::SpatialPosition &pos, const int dim) const
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