OGS
HydroMechanicsFEM-impl.h
Go to the documentation of this file.
1 
12 #pragma once
13 
14 #include <Eigen/Eigenvalues>
15 
16 #include "MaterialLib/MPL/Medium.h"
17 #include "MaterialLib/MPL/Property.h"
18 #include "MaterialLib/MPL/Utils/FormEigenTensor.h"
19 #include "MaterialLib/MPL/Utils/GetSymmetricTensor.h"
20 #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
21 #include "MathLib/KelvinVector.h"
22 #include "NumLib/Function/Interpolation.h"
23 #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
24 #include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
25 
26 namespace ProcessLib
27 {
28 namespace HydroMechanics
29 {
30 namespace MPL = MaterialPropertyLib;
31 
32 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
33  int DisplacementDim>
34 HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
35  DisplacementDim>::
36  HydroMechanicsLocalAssembler(
37  MeshLib::Element const& e,
38  std::size_t const /*local_matrix_size*/,
39  NumLib::GenericIntegrationMethod const& integration_method,
40  bool const is_axially_symmetric,
41  HydroMechanicsProcessData<DisplacementDim>& process_data)
42  : _process_data(process_data),
43  _integration_method(integration_method),
44  _element(e),
45  _is_axially_symmetric(is_axially_symmetric)
46 {
47  unsigned const n_integration_points =
48  _integration_method.getNumberOfPoints();
49 
50  _ip_data.reserve(n_integration_points);
51  _secondary_data.N_u.resize(n_integration_points);
52 
53  auto const shape_matrices_u =
54  NumLib::initShapeMatrices<ShapeFunctionDisplacement,
55  ShapeMatricesTypeDisplacement,
56  DisplacementDim>(e, is_axially_symmetric,
57  _integration_method);
58 
59  auto const shape_matrices_p =
60  NumLib::initShapeMatrices<ShapeFunctionPressure,
61  ShapeMatricesTypePressure, DisplacementDim>(
62  e, is_axially_symmetric, _integration_method);
63 
64  auto const& solid_material =
65  MaterialLib::Solids::selectSolidConstitutiveRelation(
66  _process_data.solid_materials, _process_data.material_ids,
67  e.getID());
68 
69  for (unsigned ip = 0; ip < n_integration_points; ip++)
70  {
71  _ip_data.emplace_back(solid_material);
72  auto& ip_data = _ip_data[ip];
73  auto const& sm_u = shape_matrices_u[ip];
74  ip_data.integration_weight =
75  _integration_method.getWeightedPoint(ip).getWeight() *
76  sm_u.integralMeasure * sm_u.detJ;
77 
78  // Initialize current time step values
79  static const int kelvin_vector_size =
80  MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
81  ip_data.sigma_eff.setZero(kelvin_vector_size);
82  ip_data.eps.setZero(kelvin_vector_size);
83 
84  // Previous time step values are not initialized and are set later.
85  ip_data.eps_prev.resize(kelvin_vector_size);
86  ip_data.sigma_eff_prev.resize(kelvin_vector_size);
87 
88  ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
89  DisplacementDim, displacement_size>::Zero(DisplacementDim,
90  displacement_size);
91  for (int i = 0; i < DisplacementDim; ++i)
92  {
93  ip_data.N_u_op
94  .template block<1, displacement_size / DisplacementDim>(
95  i, i * displacement_size / DisplacementDim)
96  .noalias() = sm_u.N;
97  }
98 
99  ip_data.N_u = sm_u.N;
100  ip_data.dNdx_u = sm_u.dNdx;
101 
102  ip_data.N_p = shape_matrices_p[ip].N;
103  ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
104 
105  _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
106  }
107 }
108 
109 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
110  int DisplacementDim>
111 void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
112  ShapeFunctionPressure, DisplacementDim>::
113  assembleWithJacobian(double const t, double const dt,
114  std::vector<double> const& local_x,
115  std::vector<double> const& local_xdot,
116  std::vector<double>& /*local_M_data*/,
117  std::vector<double>& /*local_K_data*/,
118  std::vector<double>& local_rhs_data,
119  std::vector<double>& local_Jac_data)
120 {
121  assert(local_x.size() == pressure_size + displacement_size);
122 
123  auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
124  pressure_size> const>(local_x.data() + pressure_index, pressure_size);
125 
126  auto u =
127  Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
128  displacement_size> const>(local_x.data() + displacement_index,
129  displacement_size);
130 
131  auto p_dot =
132  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
133  pressure_size> const>(local_xdot.data() + pressure_index,
134  pressure_size);
135  auto u_dot =
136  Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
137  displacement_size> const>(local_xdot.data() + displacement_index,
138  displacement_size);
139 
140  auto local_Jac = MathLib::createZeroedMatrix<
141  typename ShapeMatricesTypeDisplacement::template MatrixType<
142  displacement_size + pressure_size,
143  displacement_size + pressure_size>>(
144  local_Jac_data, displacement_size + pressure_size,
145  displacement_size + pressure_size);
146 
147  auto local_rhs = MathLib::createZeroedVector<
148  typename ShapeMatricesTypeDisplacement::template VectorType<
149  displacement_size + pressure_size>>(
150  local_rhs_data, displacement_size + pressure_size);
151 
152  typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
153  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
154  pressure_size);
155 
156  typename ShapeMatricesTypePressure::NodalMatrixType storage_p =
157  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
158  pressure_size);
159 
160  typename ShapeMatricesTypePressure::NodalMatrixType add_p_derivative =
161  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
162  pressure_size);
163 
164  typename ShapeMatricesTypeDisplacement::template MatrixType<
165  displacement_size, pressure_size>
166  Kup = ShapeMatricesTypeDisplacement::template MatrixType<
167  displacement_size, pressure_size>::Zero(displacement_size,
168  pressure_size);
169 
170  typename ShapeMatricesTypeDisplacement::template MatrixType<
171  pressure_size, displacement_size>
172  Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
173  pressure_size, displacement_size>::Zero(pressure_size,
174  displacement_size);
175 
176  typename ShapeMatricesTypeDisplacement::template MatrixType<
177  pressure_size, displacement_size>
178  Kpu_k = ShapeMatricesTypeDisplacement::template MatrixType<
179  pressure_size, displacement_size>::Zero(pressure_size,
180  displacement_size);
181 
182  auto const& solid_material =
183  MaterialLib::Solids::selectSolidConstitutiveRelation(
184  _process_data.solid_materials, _process_data.material_ids,
185  _element.getID());
186 
187  ParameterLib::SpatialPosition x_position;
188  x_position.setElementID(_element.getID());
189 
190  unsigned const n_integration_points =
191  _integration_method.getNumberOfPoints();
192 
193  auto const& b = _process_data.specific_body_force;
194  auto const& medium = _process_data.media_map->getMedium(_element.getID());
195  auto const& solid = medium->phase("Solid");
196  auto const& fluid = fluidPhase(*medium);
197  MPL::VariableArray vars;
198 
199  auto const T_ref =
200  medium->property(MPL::PropertyType::reference_temperature)
201  .template value<double>(vars, x_position, t, dt);
202  vars.temperature = T_ref;
203 
204  auto const& identity2 = Invariants::identity2;
205 
206  for (unsigned ip = 0; ip < n_integration_points; ip++)
207  {
208  x_position.setIntegrationPoint(ip);
209  auto const& w = _ip_data[ip].integration_weight;
210 
211  auto const& N_u_op = _ip_data[ip].N_u_op;
212 
213  auto const& N_u = _ip_data[ip].N_u;
214  auto const& dNdx_u = _ip_data[ip].dNdx_u;
215 
216  auto const& N_p = _ip_data[ip].N_p;
217  auto const& dNdx_p = _ip_data[ip].dNdx_p;
218 
219  auto const x_coord =
220  NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
221  ShapeMatricesTypeDisplacement>(
222  _element, N_u);
223  auto const B =
224  LinearBMatrix::computeBMatrix<DisplacementDim,
225  ShapeFunctionDisplacement::NPOINTS,
226  typename BMatricesType::BMatrixType>(
227  dNdx_u, N_u, x_coord, _is_axially_symmetric);
228 
229  auto& eps = _ip_data[ip].eps;
230  eps.noalias() = B * u;
231  auto const& sigma_eff = _ip_data[ip].sigma_eff;
232 
233  double const p_int_pt = N_p.dot(p);
234  vars.phase_pressure = p_int_pt;
235 
236  auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
237  t, x_position, dt, T_ref);
238  auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
239 
240  auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
241  .template value<double>(vars, x_position, t, dt);
242 
243  auto const rho_sr =
244  solid.property(MPL::PropertyType::density)
245  .template value<double>(vars, x_position, t, dt);
246  auto const porosity =
247  medium->property(MPL::PropertyType::porosity)
248  .template value<double>(vars, x_position, t, dt);
249 
250  auto const mu = fluid.property(MPL::PropertyType::viscosity)
251  .template value<double>(vars, x_position, t, dt);
252 
253  // Quick workaround: If fluid density is described as ideal gas, then
254  // the molar mass must be passed to the MPL::IdealGasLaw via the
255  // variable_array and the fluid must have the property
256  // MPL::PropertyType::molar_mass. For other density models (e.g.
257  // Constant), it is not mandatory to specify the molar mass.
258  if (fluid.hasProperty(MPL::PropertyType::molar_mass))
259  {
260  vars.molar_mass =
261  fluid.property(MPL::PropertyType::molar_mass)
262  .template value<double>(vars, x_position, t, dt);
263  }
264  auto const rho_fr =
265  fluid.property(MPL::PropertyType::density)
266  .template value<double>(vars, x_position, t, dt);
267  auto const beta_p =
268  fluid.property(MPL::PropertyType::density)
269  .template dValue<double>(vars, MPL::Variable::phase_pressure,
270  x_position, t, dt) /
271  rho_fr;
272 
273  // Set mechanical variables for the intrinsic permeability model
274  // For stress dependent permeability.
275  {
276  auto const sigma_total =
277  (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
278 
279  vars.total_stress.emplace<SymmetricTensor>(
280  MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
281  sigma_total));
282  }
283  // For strain dependent permeability
284  vars.volumetric_strain = Invariants::trace(eps);
285  vars.equivalent_plastic_strain =
286  _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
287  vars.mechanical_strain
288  .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
289  eps);
290 
291  auto const K = MPL::formEigenTensor<DisplacementDim>(
292  medium->property(MPL::PropertyType::permeability)
293  .value(vars, x_position, t, dt));
294  auto const dkde = MPL::formEigenTensor<
295  MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim)>(
296  (*medium)[MPL::PropertyType::permeability].dValue(
297  vars, MPL::Variable::mechanical_strain, x_position, t, dt));
298 
299  auto const K_over_mu = K / mu;
300 
301  auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
302  dt, u, T_ref);
303 
304  //
305  // displacement equation, displacement part
306  //
307  local_Jac
308  .template block<displacement_size, displacement_size>(
309  displacement_index, displacement_index)
310  .noalias() += B.transpose() * C * B * w;
311 
312  double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
313  local_rhs.template segment<displacement_size>(displacement_index)
314  .noalias() -=
315  (B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w;
316 
317  //
318  // displacement equation, pressure part
319  //
320  Kup.noalias() += B.transpose() * alpha * identity2 * N_p * w;
321 
322  //
323  // pressure equation, pressure part.
324  //
325  laplace_p.noalias() +=
326  rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
327 
328  storage_p.noalias() +=
329  rho_fr * N_p.transpose() * N_p * w *
330  ((alpha - porosity) * (1.0 - alpha) / K_S + porosity * beta_p);
331 
332  // density dependence on pressure evaluated for Darcy-term,
333  // for laplace and storage terms this dependence is neglected
334  add_p_derivative.noalias() += rho_fr * beta_p * dNdx_p.transpose() *
335  K_over_mu *
336  (dNdx_p * p - 2.0 * rho_fr * b) * N_p * w;
337 
338  local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
339  dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
340 
341  //
342  // pressure equation, displacement part.
343  //
344  Kpu.noalias() +=
345  rho_fr * alpha * N_p.transpose() * identity2.transpose() * B * w;
346 
347  Kpu_k.noalias() +=
348  dNdx_p.transpose() *
349  MathLib::KelvinVector::liftVectorToKelvin<DisplacementDim>(
350  dNdx_p * p - rho_fr * b) *
351  dkde * B * rho_fr / mu * w;
352  }
353  // displacement equation, pressure part
354  local_Jac
355  .template block<displacement_size, pressure_size>(displacement_index,
356  pressure_index)
357  .noalias() = -Kup;
358 
359  if (_process_data.mass_lumping)
360  {
361  storage_p = storage_p.colwise().sum().eval().asDiagonal();
362 
363  if constexpr (pressure_size == displacement_size)
364  {
365  Kpu = Kpu.colwise().sum().eval().asDiagonal();
366  Kpu_k = Kpu_k.colwise().sum().eval().asDiagonal();
367  }
368  }
369 
370  // pressure equation, pressure part.
371  local_Jac
372  .template block<pressure_size, pressure_size>(pressure_index,
373  pressure_index)
374  .noalias() += laplace_p + storage_p / dt + add_p_derivative;
375 
376  // pressure equation, displacement part.
377  local_Jac
378  .template block<pressure_size, displacement_size>(pressure_index,
379  displacement_index)
380  .noalias() += Kpu / dt + Kpu_k;
381 
382  // pressure equation
383  local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
384  laplace_p * p + storage_p * p_dot + Kpu * u_dot;
385 
386  // displacement equation
387  local_rhs.template segment<displacement_size>(displacement_index)
388  .noalias() += Kup * p;
389 }
390 
391 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
392  int DisplacementDim>
393 std::vector<double> const& HydroMechanicsLocalAssembler<
394  ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
395  getIntPtDarcyVelocity(
396  const double t,
397  std::vector<GlobalVector*> const& x,
398  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
399  std::vector<double>& cache) const
400 {
401  int const hydraulic_process_id = _process_data.hydraulic_process_id;
402  auto const indices =
403  NumLib::getIndices(_element.getID(), *dof_table[hydraulic_process_id]);
404  assert(!indices.empty());
405  auto const local_x = x[hydraulic_process_id]->get(indices);
406 
407  unsigned const n_integration_points =
408  _integration_method.getNumberOfPoints();
409  cache.clear();
410  auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
411  double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
412  cache, DisplacementDim, n_integration_points);
413 
414  auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
415  pressure_size> const>(local_x.data() + pressure_index, pressure_size);
416 
417  ParameterLib::SpatialPosition x_position;
418  x_position.setElementID(_element.getID());
419 
420  auto const& medium = _process_data.media_map->getMedium(_element.getID());
421  auto const& fluid = fluidPhase(*medium);
422  MPL::VariableArray vars;
423 
424  // TODO (naumov) Temporary value not used by current material models. Need
425  // extension of secondary variables interface.
426  double const dt = std::numeric_limits<double>::quiet_NaN();
427  vars.temperature =
428  medium->property(MPL::PropertyType::reference_temperature)
429  .template value<double>(vars, x_position, t, dt);
430 
431  auto const& identity2 = Invariants::identity2;
432 
433  for (unsigned ip = 0; ip < n_integration_points; ip++)
434  {
435  x_position.setIntegrationPoint(ip);
436 
437  double const p_int_pt = _ip_data[ip].N_p.dot(p);
438  vars.phase_pressure = p_int_pt;
439 
440  auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
441  .template value<double>(vars, x_position, t, dt);
442 
443  // Set mechanical variables for the intrinsic permeability model
444  // For stress dependent permeability.
445  auto const sigma_total =
446  (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
447  vars.total_stress.emplace<SymmetricTensor>(
448  MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_total));
449 
450  // For strain dependent permeability
451  vars.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
452  vars.equivalent_plastic_strain =
453  _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
454  vars.mechanical_strain
455  .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
456  _ip_data[ip].eps);
457 
458  auto const K = MPL::formEigenTensor<DisplacementDim>(
459  medium->property(MPL::PropertyType::permeability)
460  .value(vars, x_position, t, dt));
461 
462  auto const mu = fluid.property(MPL::PropertyType::viscosity)
463  .template value<double>(vars, x_position, t, dt);
464  // Quick workaround: If fluid density is described as ideal gas, then
465  // the molar mass must be passed to the MPL::IdealGasLaw via the
466  // variable_array and the fluid must have the property
467  // MPL::PropertyType::molar_mass. For other density models (e.g.
468  // Constant), it is not mandatory to specify the molar mass.
469  if (fluid.hasProperty(MPL::PropertyType::molar_mass))
470  {
471  vars.molar_mass =
472  fluid.property(MPL::PropertyType::molar_mass)
473  .template value<double>(vars, x_position, t, dt);
474  }
475  auto const rho_fr =
476  fluid.property(MPL::PropertyType::density)
477  .template value<double>(vars, x_position, t, dt);
478 
479  auto const K_over_mu = K / mu;
480 
481  auto const& b = _process_data.specific_body_force;
482 
483  // Compute the velocity
484  auto const& dNdx_p = _ip_data[ip].dNdx_p;
485  cache_matrix.col(ip).noalias() =
486  -K_over_mu * dNdx_p * p + K_over_mu * rho_fr * b;
487  }
488 
489  return cache;
490 }
491 
492 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
493  int DisplacementDim>
494 void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
495  ShapeFunctionPressure, DisplacementDim>::
496  assembleWithJacobianForPressureEquations(
497  const double t, double const dt, Eigen::VectorXd const& local_x,
498  Eigen::VectorXd const& local_xdot, std::vector<double>& local_b_data,
499  std::vector<double>& local_Jac_data)
500 {
501  auto local_rhs =
502  MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
503  template VectorType<pressure_size>>(
504  local_b_data, pressure_size);
505 
506  ParameterLib::SpatialPosition pos;
507  pos.setElementID(this->_element.getID());
508 
509  auto const p = local_x.template segment<pressure_size>(pressure_index);
510 
511  auto const p_dot =
512  local_xdot.template segment<pressure_size>(pressure_index);
513 
514  auto local_Jac = MathLib::createZeroedMatrix<
515  typename ShapeMatricesTypeDisplacement::template MatrixType<
516  pressure_size, pressure_size>>(local_Jac_data, pressure_size,
517  pressure_size);
518 
519  typename ShapeMatricesTypePressure::NodalMatrixType storage =
520  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
521  pressure_size);
522 
523  typename ShapeMatricesTypePressure::NodalMatrixType laplace =
524  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
525  pressure_size);
526 
527  typename ShapeMatricesTypePressure::NodalMatrixType coupling =
528  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
529  pressure_size);
530 
531  typename ShapeMatricesTypePressure::NodalMatrixType add_p_derivative =
532  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
533  pressure_size);
534 
535  auto const& solid_material =
536  MaterialLib::Solids::selectSolidConstitutiveRelation(
537  _process_data.solid_materials, _process_data.material_ids,
538  _element.getID());
539 
540  ParameterLib::SpatialPosition x_position;
541  x_position.setElementID(_element.getID());
542 
543  auto const& medium = _process_data.media_map->getMedium(_element.getID());
544  auto const& fluid = fluidPhase(*medium);
545  MPL::VariableArray vars;
546 
547  auto const T_ref =
548  medium->property(MPL::PropertyType::reference_temperature)
549  .template value<double>(vars, x_position, t, dt);
550  vars.temperature = T_ref;
551 
552  auto const& identity2 = Invariants::identity2;
553 
554  auto const fixed_stress_stabilization_parameter =
555  _process_data.fixed_stress_stabilization_parameter;
556 
557  int const n_integration_points = _integration_method.getNumberOfPoints();
558  for (int ip = 0; ip < n_integration_points; ip++)
559  {
560  x_position.setIntegrationPoint(ip);
561  auto const& w = _ip_data[ip].integration_weight;
562 
563  auto const& N_p = _ip_data[ip].N_p;
564  auto const& dNdx_p = _ip_data[ip].dNdx_p;
565 
566  double const p_int_pt = N_p.dot(p);
567  vars.phase_pressure = p_int_pt;
568 
569  auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
570  t, x_position, dt, T_ref);
571  auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
572 
573  auto const alpha_b =
574  medium->property(MPL::PropertyType::biot_coefficient)
575  .template value<double>(vars, x_position, t, dt);
576 
577  // Set mechanical variables for the intrinsic permeability model
578  // For stress dependent permeability.
579  auto const sigma_total =
580  (_ip_data[ip].sigma_eff - alpha_b * p_int_pt * identity2).eval();
581  vars.total_stress.emplace<SymmetricTensor>(
582  MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_total));
583 
584  // For strain dependent permeability
585  vars.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
586  vars.equivalent_plastic_strain =
587  _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
588  vars.mechanical_strain
589  .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
590  _ip_data[ip].eps);
591 
592  auto const K = MPL::formEigenTensor<DisplacementDim>(
593  medium->property(MPL::PropertyType::permeability)
594  .value(vars, x_position, t, dt));
595  auto const porosity =
596  medium->property(MPL::PropertyType::porosity)
597  .template value<double>(vars, x_position, t, dt);
598 
599  auto const mu = fluid.property(MPL::PropertyType::viscosity)
600  .template value<double>(vars, x_position, t, dt);
601  // Quick workaround: If fluid density is described as ideal gas, then
602  // the molar mass must be passed to the MPL::IdealGasLaw via the
603  // variable_array and the fluid must have the property
604  // MPL::PropertyType::molar_mass. For other density models (e.g.
605  // Constant), it is not mandatory to specify the molar mass.
606  if (fluid.hasProperty(MPL::PropertyType::molar_mass))
607  {
608  vars.molar_mass =
609  fluid.property(MPL::PropertyType::molar_mass)
610  .template value<double>(vars, x_position, t, dt);
611  }
612  auto const rho_fr =
613  fluid.property(MPL::PropertyType::density)
614  .template value<double>(vars, x_position, t, dt);
615  auto const beta_p =
616  fluid.property(MPL::PropertyType::density)
617  .template dValue<double>(vars, MPL::Variable::phase_pressure,
618  x_position, t, dt) /
619  rho_fr;
620 
621  auto const K_over_mu = K / mu;
622 
623  // artificial compressibility to account for coupling
624  auto const beta_FS =
625  fixed_stress_stabilization_parameter * alpha_b * alpha_b / K_S;
626 
627  laplace.noalias() +=
628  rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
629 
630  storage.noalias() +=
631  rho_fr * N_p.transpose() * N_p * w *
632  ((alpha_b - porosity) * (1.0 - alpha_b) / K_S + porosity * beta_p);
633 
634  auto const& b = _process_data.specific_body_force;
635 
636  // bodyforce-driven Darcy flow
637  local_rhs.noalias() +=
638  dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
639 
640  // density dependence on pressure evaluated for Darcy-term,
641  // for laplace and storage terms this dependence is neglected (as is
642  // done for monolithic too)
643  add_p_derivative.noalias() += rho_fr * beta_p * dNdx_p.transpose() *
644  K_over_mu *
645  (dNdx_p * p - 2.0 * rho_fr * b) * N_p * w;
646 
647  auto& eps = _ip_data[ip].eps;
648  auto& eps_prev = _ip_data[ip].eps_prev;
649  const double eps_v_dot =
650  (Invariants::trace(eps) - Invariants::trace(eps_prev)) / dt;
651 
652  const double p_c = _ip_data[ip].coupling_pressure;
653 
654  // pressure-dependent part of fixed-stress coupling
655  coupling.noalias() += rho_fr * N_p.transpose() * N_p * w * beta_FS / dt;
656 
657  // constant part of fixed-stress coupling
658  local_rhs.noalias() -=
659  (alpha_b * eps_v_dot - beta_FS * p_c / dt) * rho_fr * N_p * w;
660  }
661  local_Jac.noalias() = laplace + storage / dt + add_p_derivative + coupling;
662 
663  local_rhs.noalias() -= laplace * p + storage * p_dot + coupling * p;
664 }
665 
666 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
667  int DisplacementDim>
668 void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
669  ShapeFunctionPressure, DisplacementDim>::
670  assembleWithJacobianForDeformationEquations(
671  const double t, double const dt, Eigen::VectorXd const& local_x,
672  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
673 {
674  auto const p = local_x.template segment<pressure_size>(pressure_index);
675  auto const u =
676  local_x.template segment<displacement_size>(displacement_index);
677 
678  auto local_Jac = MathLib::createZeroedMatrix<
679  typename ShapeMatricesTypeDisplacement::template MatrixType<
680  displacement_size, displacement_size>>(
681  local_Jac_data, displacement_size, displacement_size);
682 
683  auto local_rhs =
684  MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
685  template VectorType<displacement_size>>(
686  local_b_data, displacement_size);
687 
688  ParameterLib::SpatialPosition x_position;
689  x_position.setElementID(_element.getID());
690 
691  auto const& medium = _process_data.media_map->getMedium(_element.getID());
692  auto const& solid = medium->phase("Solid");
693  auto const& fluid = fluidPhase(*medium);
694  MPL::VariableArray vars;
695 
696  auto const T_ref =
697  medium->property(MPL::PropertyType::reference_temperature)
698  .template value<double>(vars, x_position, t, dt);
699  vars.temperature = T_ref;
700 
701  int const n_integration_points = _integration_method.getNumberOfPoints();
702  for (int ip = 0; ip < n_integration_points; ip++)
703  {
704  x_position.setIntegrationPoint(ip);
705  auto const& w = _ip_data[ip].integration_weight;
706 
707  auto const& N_u_op = _ip_data[ip].N_u_op;
708 
709  auto const& N_u = _ip_data[ip].N_u;
710  auto const& dNdx_u = _ip_data[ip].dNdx_u;
711 
712  auto const& N_p = _ip_data[ip].N_p;
713 
714  auto const x_coord =
715  NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
716  ShapeMatricesTypeDisplacement>(
717  _element, N_u);
718  auto const B =
719  LinearBMatrix::computeBMatrix<DisplacementDim,
720  ShapeFunctionDisplacement::NPOINTS,
721  typename BMatricesType::BMatrixType>(
722  dNdx_u, N_u, x_coord, _is_axially_symmetric);
723 
724  auto& eps = _ip_data[ip].eps;
725  auto const& sigma_eff = _ip_data[ip].sigma_eff;
726 
727  vars.phase_pressure = N_p.dot(p);
728 
729  auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
730  .template value<double>(vars, x_position, t, dt);
731  auto const rho_sr =
732  solid.property(MPL::PropertyType::density)
733  .template value<double>(vars, x_position, t, dt);
734  auto const porosity =
735  medium->property(MPL::PropertyType::porosity)
736  .template value<double>(vars, x_position, t, dt);
737 
738  // Quick workaround: If fluid density is described as ideal gas, then
739  // the molar mass must be passed to the MPL::IdealGasLaw via the
740  // variable_array and the fluid must have the property
741  // MPL::PropertyType::molar_mass. For other density models (e.g.
742  // Constant), it is not mandatory to specify the molar mass.
743  if (fluid.hasProperty(MPL::PropertyType::molar_mass))
744  {
745  vars.molar_mass =
746  fluid.property(MPL::PropertyType::molar_mass)
747  .template value<double>(vars, x_position, t, dt);
748  }
749  auto const rho_fr =
750  fluid.property(MPL::PropertyType::density)
751  .template value<double>(vars, x_position, t, dt);
752 
753  auto const& b = _process_data.specific_body_force;
754  auto const& identity2 = MathLib::KelvinVector::Invariants<
755  MathLib::KelvinVector::kelvin_vector_dimensions(
756  DisplacementDim)>::identity2;
757 
758  eps.noalias() = B * u;
759  vars.mechanical_strain
760  .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
761  eps);
762 
763  auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
764  dt, u, T_ref);
765 
766  local_Jac.noalias() += B.transpose() * C * B * w;
767 
768  double p_at_xi = 0.;
769  NumLib::shapeFunctionInterpolate(p, N_p, p_at_xi);
770 
771  double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
772  local_rhs.noalias() -=
773  (B.transpose() * (sigma_eff - alpha * identity2 * p_at_xi) -
774  N_u_op.transpose() * rho * b) *
775  w;
776  }
777 }
778 
779 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
780  int DisplacementDim>
781 void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
782  ShapeFunctionPressure, DisplacementDim>::
783  assembleWithJacobianForStaggeredScheme(
784  const double t, double const dt, Eigen::VectorXd const& local_x,
785  Eigen::VectorXd const& local_xdot, int const process_id,
786  std::vector<double>& /*local_M_data*/,
787  std::vector<double>& /*local_K_data*/,
788  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
789 {
790  // For the equations with pressure
791  if (process_id == _process_data.hydraulic_process_id)
792  {
793  assembleWithJacobianForPressureEquations(t, dt, local_x, local_xdot,
794  local_b_data, local_Jac_data);
795  return;
796  }
797 
798  // For the equations with deformation
799  assembleWithJacobianForDeformationEquations(t, dt, local_x, local_b_data,
800  local_Jac_data);
801 }
802 
803 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
804  int DisplacementDim>
805 void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
806  ShapeFunctionPressure, DisplacementDim>::
807  setInitialConditionsConcrete(std::vector<double> const& local_x,
808  double const /*t*/,
809  bool const use_monolithic_scheme,
810  int const process_id)
811 {
812  int const n_integration_points = _integration_method.getNumberOfPoints();
813 
814  if (use_monolithic_scheme || process_id == 0)
815  {
816  auto p =
817  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
818  pressure_size> const>(local_x.data(), pressure_size);
819 
820  for (int ip = 0; ip < n_integration_points; ip++)
821  {
822  auto const& N_p = _ip_data[ip].N_p;
823  _ip_data[ip].coupling_pressure = N_p.dot(p);
824  }
825  }
826 
827  if (use_monolithic_scheme || process_id == 1)
828  {
829  ParameterLib::SpatialPosition x_position;
830  x_position.setElementID(_element.getID());
831 
832  const int displacement_offset =
833  use_monolithic_scheme ? displacement_index : 0;
834 
835  auto u = Eigen::Map<typename ShapeMatricesTypeDisplacement::
836  template VectorType<displacement_size> const>(
837  local_x.data() + displacement_offset, displacement_size);
838 
839  MPL::VariableArray vars;
840 
841  for (int ip = 0; ip < n_integration_points; ip++)
842  {
843  x_position.setIntegrationPoint(ip);
844  auto const& N_u = _ip_data[ip].N_u;
845  auto const& dNdx_u = _ip_data[ip].dNdx_u;
846 
847  auto const x_coord =
848  NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
849  ShapeMatricesTypeDisplacement>(
850  _element, N_u);
851  auto const B = LinearBMatrix::computeBMatrix<
852  DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
853  typename BMatricesType::BMatrixType>(dNdx_u, N_u, x_coord,
854  _is_axially_symmetric);
855 
856  auto& eps = _ip_data[ip].eps;
857  eps.noalias() = B * u;
858  vars.mechanical_strain.emplace<
859  MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(eps);
860  }
861  }
862 }
863 
864 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
865  int DisplacementDim>
866 void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
867  ShapeFunctionPressure, DisplacementDim>::
868  postNonLinearSolverConcrete(std::vector<double> const& local_x,
869  std::vector<double> const& /*local_xdot*/,
870  double const t, double const dt,
871  bool const use_monolithic_scheme,
872  int const process_id)
873 {
874  int const n_integration_points = _integration_method.getNumberOfPoints();
875 
876  if (use_monolithic_scheme || process_id == 0)
877  {
878  auto p =
879  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
880  pressure_size> const>(local_x.data(), pressure_size);
881 
882  for (int ip = 0; ip < n_integration_points; ip++)
883  {
884  auto const& N_p = _ip_data[ip].N_p;
885  _ip_data[ip].coupling_pressure = N_p.dot(p);
886  }
887  }
888 
889  if (use_monolithic_scheme || process_id == 1)
890  {
891  ParameterLib::SpatialPosition x_position;
892  x_position.setElementID(_element.getID());
893 
894  auto const& medium =
895  _process_data.media_map->getMedium(_element.getID());
896 
897  auto const T_ref =
898  medium->property(MPL::PropertyType::reference_temperature)
899  .template value<double>(MPL::VariableArray(), x_position, t,
900  dt);
901 
902  const int displacement_offset =
903  use_monolithic_scheme ? displacement_index : 0;
904 
905  auto u = Eigen::Map<typename ShapeMatricesTypeDisplacement::
906  template VectorType<displacement_size> const>(
907  local_x.data() + displacement_offset, displacement_size);
908 
909  MPL::VariableArray vars;
910  vars.temperature = T_ref;
911 
912  for (int ip = 0; ip < n_integration_points; ip++)
913  {
914  x_position.setIntegrationPoint(ip);
915  auto const& N_u = _ip_data[ip].N_u;
916  auto const& dNdx_u = _ip_data[ip].dNdx_u;
917 
918  auto const x_coord =
919  NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
920  ShapeMatricesTypeDisplacement>(
921  _element, N_u);
922  auto const B = LinearBMatrix::computeBMatrix<
923  DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
924  typename BMatricesType::BMatrixType>(dNdx_u, N_u, x_coord,
925  _is_axially_symmetric);
926 
927  auto& eps = _ip_data[ip].eps;
928  eps.noalias() = B * u;
929  vars.mechanical_strain.emplace<
930  MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(eps);
931 
932  _ip_data[ip].updateConstitutiveRelation(vars, t, x_position, dt, u,
933  T_ref);
934  }
935  }
936 }
937 
938 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
939  int DisplacementDim>
940 std::size_t HydroMechanicsLocalAssembler<
941  ShapeFunctionDisplacement, ShapeFunctionPressure,
942  DisplacementDim>::setIPDataInitialConditions(std::string const& name,
943  double const* values,
944  int const integration_order)
945 {
946  if (integration_order !=
947  static_cast<int>(_integration_method.getIntegrationOrder()))
948  {
949  OGS_FATAL(
950  "Setting integration point initial conditions; The integration "
951  "order of the local assembler for element {:d} is different from "
952  "the integration order in the initial condition.",
953  _element.getID());
954  }
955 
956  if (name == "sigma_ip")
957  {
958  if (_process_data.initial_stress != nullptr)
959  {
960  OGS_FATAL(
961  "Setting initial conditions for stress from integration "
962  "point data and from a parameter '{:s}' is not possible "
963  "simultaneously.",
964  _process_data.initial_stress->name);
965  }
966 
967  return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
968  values, _ip_data, &IpData::sigma_eff);
969  }
970 
971  if (name == "epsilon_ip")
972  {
973  return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
974  values, _ip_data, &IpData::eps);
975  }
976 
977  return 0;
978 }
979 
980 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
981  int DisplacementDim>
982 std::vector<double>
983 HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
984  DisplacementDim>::getSigma() const
985 {
986  return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
987  _ip_data, &IpData::sigma_eff);
988 }
989 
990 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
991  int DisplacementDim>
992 std::vector<double>
993 HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
994  DisplacementDim>::getEpsilon() const
995 {
996  auto const kelvin_vector_size =
997  MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
998  unsigned const n_integration_points =
999  _integration_method.getNumberOfPoints();
1000 
1001  std::vector<double> ip_epsilon_values;
1002  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
1003  double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
1004  ip_epsilon_values, n_integration_points, kelvin_vector_size);
1005 
1006  for (unsigned ip = 0; ip < n_integration_points; ++ip)
1007  {
1008  auto const& eps = _ip_data[ip].eps;
1009  cache_mat.row(ip) =
1010  MathLib::KelvinVector::kelvinVectorToSymmetricTensor(eps);
1011  }
1012 
1013  return ip_epsilon_values;
1014 }
1015 
1016 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1017  int DisplacementDim>
1018 void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
1019  ShapeFunctionPressure, DisplacementDim>::
1020  computeSecondaryVariableConcrete(double const t, double const dt,
1021  Eigen::VectorXd const& local_x,
1022  Eigen::VectorXd const& /*local_x_dot*/)
1023 {
1024  auto const p = local_x.template segment<pressure_size>(pressure_index);
1025 
1026  NumLib::interpolateToHigherOrderNodes<
1027  ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1028  DisplacementDim>(_element, _is_axially_symmetric, p,
1029  *_process_data.pressure_interpolated);
1030 
1031  int const elem_id = _element.getID();
1032  ParameterLib::SpatialPosition x_position;
1033  x_position.setElementID(elem_id);
1034  unsigned const n_integration_points =
1035  _integration_method.getNumberOfPoints();
1036 
1037  auto const& medium = _process_data.media_map->getMedium(elem_id);
1038  MPL::VariableArray vars;
1039 
1040  SymmetricTensor k_sum = SymmetricTensor::Zero(KelvinVectorSize);
1041  auto sigma_eff_sum = MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
1042  Eigen::Matrix<double, 3, 3>::Zero());
1043 
1044  auto const& identity2 = Invariants::identity2;
1045 
1046  for (unsigned ip = 0; ip < n_integration_points; ip++)
1047  {
1048  x_position.setIntegrationPoint(ip);
1049 
1050  auto const& eps = _ip_data[ip].eps;
1051  sigma_eff_sum += _ip_data[ip].sigma_eff;
1052 
1053  auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
1054  .template value<double>(vars, x_position, t, dt);
1055  double const p_int_pt = _ip_data[ip].N_p.dot(p);
1056  vars.phase_pressure = p_int_pt;
1057 
1058  // Set mechanical variables for the intrinsic permeability model
1059  // For stress dependent permeability.
1060  {
1061  auto const sigma_total =
1062  (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
1063  vars.total_stress.emplace<SymmetricTensor>(
1064  MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
1065  sigma_total));
1066  }
1067  // For strain dependent permeability
1068  vars.volumetric_strain = Invariants::trace(eps);
1069  vars.equivalent_plastic_strain =
1070  _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
1071  vars.mechanical_strain
1072  .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
1073  eps);
1074 
1075  k_sum += MPL::getSymmetricTensor<DisplacementDim>(
1076  medium->property(MPL::PropertyType::permeability)
1077  .value(vars, x_position, t, dt));
1078  }
1079 
1080  Eigen::Map<Eigen::VectorXd>(
1081  &(*_process_data.permeability)[elem_id * KelvinVectorSize],
1082  KelvinVectorSize) = k_sum / n_integration_points;
1083 
1084  Eigen::Matrix<double, 3, 3, 0, 3, 3> const sigma_avg =
1085  MathLib::KelvinVector::kelvinVectorToTensor(sigma_eff_sum) /
1086  n_integration_points;
1087 
1088  Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> e_s(sigma_avg);
1089 
1090  Eigen::Map<Eigen::Vector3d>(
1091  &(*_process_data.principal_stress_values)[elem_id * 3], 3) =
1092  e_s.eigenvalues();
1093 
1094  auto eigen_vectors = e_s.eigenvectors();
1095 
1096  for (auto i = 0; i < 3; i++)
1097  {
1098  Eigen::Map<Eigen::Vector3d>(
1099  &(*_process_data.principal_stress_vector[i])[elem_id * 3], 3) =
1100  eigen_vectors.col(i);
1101  }
1102 }
1103 
1104 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1105  int DisplacementDim>
1106 unsigned HydroMechanicsLocalAssembler<
1107  ShapeFunctionDisplacement, ShapeFunctionPressure,
1108  DisplacementDim>::getNumberOfIntegrationPoints() const
1109 {
1110  return _integration_method.getNumberOfPoints();
1111 }
1112 
1113 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1114  int DisplacementDim>
1115 typename MaterialLib::Solids::MechanicsBase<
1116  DisplacementDim>::MaterialStateVariables const&
1117 HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
1118  DisplacementDim>::
1119  getMaterialStateVariablesAt(unsigned integration_point) const
1120 {
1121  return *_ip_data[integration_point].material_state_variables;
1122 }
1123 } // namespace HydroMechanics
1124 } // namespace ProcessLib
OGS_FATAL
#define OGS_FATAL(...)
Definition: Error.h:26
MaterialPropertyLib::VariableArray::mechanical_strain
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > mechanical_strain
Definition: VariableType.h:165
MaterialPropertyLib::VariableArray::total_stress
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > total_stress
Definition: VariableType.h:184
MathLib::WeightedPoint::getWeight
double getWeight() const
Definition: WeightedPoint.h:80
MeshLib::Element::getID
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:89
NumLib::GenericIntegrationMethod::getWeightedPoint
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
Definition: GenericIntegrationMethod.h:36
ParameterLib::SpatialPosition::setElementID
void setElementID(std::size_t element_id)
Definition: SpatialPosition.h:59
ParameterLib::SpatialPosition::setIntegrationPoint
void setIntegrationPoint(unsigned integration_point)
Definition: SpatialPosition.h:65
ProcessLib::BMatrixPolicyType::BMatrixType
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition: BMatrixPolicy.h:55
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_integration_method
NumLib::GenericIntegrationMethod const & _integration_method
Definition: HydroMechanicsFEM.h:395
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::SymmetricTensor
Eigen::Matrix< double, KelvinVectorSize, 1 > SymmetricTensor
Definition: HydroMechanicsFEM.h:181
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::ShapeMatricesTypePressure
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
Definition: HydroMechanicsFEM.h:175
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::ShapeMatricesTypeDisplacement
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
Definition: HydroMechanicsFEM.h:171
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_process_data
HydroMechanicsProcessData< DisplacementDim > & _process_data
Definition: HydroMechanicsFEM.h:385
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_secondary_data
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
Definition: HydroMechanicsFEM.h:400
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler::_ip_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
Definition: HydroMechanicsFEM.h:393
MaterialLib::Solids::selectSolidConstitutiveRelation
MechanicsBase< DisplacementDim > & selectSolidConstitutiveRelation(std::map< int, std::unique_ptr< MechanicsBase< DisplacementDim >>> const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
Definition: SelectSolidConstitutiveRelation.h:33
MaterialPropertyLib::formEigenTensor
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Definition: FormEigenTensor.cpp:115
MaterialPropertyLib::fluidPhase
Phase const & fluidPhase(Medium const &medium)
Returns a gas or aqueous liquid phase of the given medium.
Definition: Medium.cpp:84
MathLib::KelvinVector::kelvinVectorToSymmetricTensor
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Definition: KelvinVector.cpp:142
MathLib::KelvinVector::KelvinVectorType
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:49
MathLib::KelvinVector::kelvin_vector_dimensions
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:24
MathLib::KelvinVector::kelvinVectorToTensor
Eigen::Matrix< double, 3, 3 > kelvinVectorToTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Definition: KelvinVector.cpp:64
MathLib::createZeroedMatrix
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:32
MathLib::p
static const double p
Definition: GaussLegendreTet.cpp:62
MathLib::u
static const double u
Definition: GaussLegendreTet.cpp:104
MathLib::createZeroedVector
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Definition: EigenMapTools.h:153
MathLib::t
static const double t
Definition: GaussLegendreTet.cpp:103
NumLib::interpolateXCoordinate
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
Definition: InitShapeMatrices.h:71
NumLib::shapeFunctionInterpolate
void shapeFunctionInterpolate(const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
Definition: Interpolation.h:80
NumLib::interpolateToHigherOrderNodes
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
Definition: Interpolation.h:107
NumLib::getIndices
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
Definition: DOFTableUtil.cpp:124
NumLib::initShapeMatrices
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
Definition: InitShapeMatrices.h:53
ProcessLib::LinearBMatrix::computeBMatrix
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
Definition: LinearBMatrix.h:42
ProcessLib::setIPDataInitialConditions
void setIPDataInitialConditions(std::vector< std::unique_ptr< IntegrationPointWriter >> const &_integration_point_writer, MeshLib::Properties const &mesh_properties, LocalAssemblersVector &local_assemblers)
Definition: SetIPDataInitialConditions.h:20
EigenFixedShapeMatrixPolicy::NodalMatrixType
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
Definition: ShapeMatrixPolicy.h:70
ProcessLib::HydroMechanics::SecondaryData::N_u
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u
Definition: HydroMechanicsFEM.h:161