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