OGS 6.1.0-1721-g6382411ad
HydroMechanicsLocalAssemblerMatrix-impl.h
Go to the documentation of this file.
1 
10 #pragma once
11 
13 
16 #include "MathLib/KelvinVector.h"
17 #include "MeshLib/ElementStatus.h"
21 
22 namespace ProcessLib
23 {
24 namespace LIE
25 {
26 namespace HydroMechanics
27 {
28 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
29  typename IntegrationMethod, int GlobalDim>
30 HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
31  ShapeFunctionPressure, IntegrationMethod,
32  GlobalDim>::
33  HydroMechanicsLocalAssemblerMatrix(
34  MeshLib::Element const& e,
35  std::size_t const n_variables,
36  std::size_t const /*local_matrix_size*/,
37  std::vector<unsigned> const& dofIndex_to_localIndex,
38  bool const is_axially_symmetric,
39  unsigned const integration_order,
42  e, is_axially_symmetric,
43  (n_variables - 1) * ShapeFunctionDisplacement::NPOINTS * GlobalDim +
44  ShapeFunctionPressure::NPOINTS,
45  dofIndex_to_localIndex),
46  _process_data(process_data)
47 {
48  IntegrationMethod integration_method(integration_order);
49  unsigned const n_integration_points =
50  integration_method.getNumberOfPoints();
51 
52  _ip_data.reserve(n_integration_points);
53 
54  auto const shape_matrices_u =
55  initShapeMatrices<ShapeFunctionDisplacement,
56  ShapeMatricesTypeDisplacement, IntegrationMethod,
57  GlobalDim>(e, is_axially_symmetric,
58  integration_method);
59 
60  auto const shape_matrices_p =
61  initShapeMatrices<ShapeFunctionPressure, ShapeMatricesTypePressure,
62  IntegrationMethod, GlobalDim>(e, is_axially_symmetric,
63  integration_method);
64 
66  _process_data.solid_materials, _process_data.material_ids, e.getID());
67 
68  SpatialPosition x_position;
69  x_position.setElementID(e.getID());
70  for (unsigned ip = 0; ip < n_integration_points; ip++)
71  {
72  x_position.setIntegrationPoint(ip);
73 
74  _ip_data.emplace_back(solid_material);
75  auto& ip_data = _ip_data[ip];
76  auto const& sm_u = shape_matrices_u[ip];
77  auto const& sm_p = shape_matrices_p[ip];
78  ip_data.integration_weight =
79  sm_u.detJ * sm_u.integralMeasure *
80  integration_method.getWeightedPoint(ip).getWeight();
81 
82  ip_data.N_u = sm_u.N;
83  ip_data.dNdx_u = sm_u.dNdx;
84  ip_data.H_u.setZero(GlobalDim, displacement_size);
85  for (int i = 0; i < GlobalDim; ++i)
86  {
87  ip_data.H_u
88  .template block<1, displacement_size / GlobalDim>(
89  i, i * displacement_size / GlobalDim)
90  .noalias() = ip_data.N_u;
91  }
92 
93  ip_data.N_p = sm_p.N;
94  ip_data.dNdx_p = sm_p.dNdx;
95 
96  ip_data.sigma_eff.setZero(kelvin_vector_size);
97  ip_data.eps.setZero(kelvin_vector_size);
98 
99  ip_data.sigma_eff_prev.resize(kelvin_vector_size);
100  ip_data.eps_prev.resize(kelvin_vector_size);
101  ip_data.C.resize(kelvin_vector_size, kelvin_vector_size);
102 
103  auto const initial_effective_stress =
104  _process_data.initial_effective_stress(0, x_position);
105  for (unsigned i = 0; i < kelvin_vector_size; i++)
106  {
107  ip_data.sigma_eff[i] = initial_effective_stress[i];
108  ip_data.sigma_eff_prev[i] = initial_effective_stress[i];
109  }
110  }
111 }
112 
113 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
114  typename IntegrationMethod, int GlobalDim>
116  ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
117  GlobalDim>::assembleWithJacobianConcrete(double const t,
118  Eigen::VectorXd const& local_x,
119  Eigen::VectorXd const& local_x_dot,
120  Eigen::VectorXd& local_rhs,
121  Eigen::MatrixXd& local_Jac)
122 {
123  auto p = const_cast<Eigen::VectorXd&>(local_x).segment(pressure_index,
124  pressure_size);
125  auto p_dot = const_cast<Eigen::VectorXd&>(local_x_dot)
126  .segment(pressure_index, pressure_size);
127 
128  if (_process_data.deactivate_matrix_in_flow)
129  {
132  }
133 
134  auto u = local_x.segment(displacement_index, displacement_size);
135  auto u_dot = local_x_dot.segment(displacement_index, displacement_size);
136 
137  auto rhs_p = local_rhs.template segment<pressure_size>(pressure_index);
138  auto rhs_u =
139  local_rhs.template segment<displacement_size>(displacement_index);
140 
141  auto J_pp = local_Jac.template block<pressure_size, pressure_size>(
143  auto J_pu = local_Jac.template block<pressure_size, displacement_size>(
145  auto J_uu = local_Jac.template block<displacement_size, displacement_size>(
147  auto J_up = local_Jac.template block<displacement_size, pressure_size>(
149 
150  assembleBlockMatricesWithJacobian(t, p, p_dot, u, u_dot, rhs_p, rhs_u, J_pp,
151  J_pu, J_uu, J_up);
152 }
153 
154 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
155  typename IntegrationMethod, int GlobalDim>
156 void HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
157  ShapeFunctionPressure,
158  IntegrationMethod, GlobalDim>::
160  double const t,
161  Eigen::Ref<const Eigen::VectorXd> const& p,
162  Eigen::Ref<const Eigen::VectorXd> const& p_dot,
163  Eigen::Ref<const Eigen::VectorXd> const& u,
164  Eigen::Ref<const Eigen::VectorXd> const& u_dot,
165  Eigen::Ref<Eigen::VectorXd>
166  rhs_p,
167  Eigen::Ref<Eigen::VectorXd>
168  rhs_u,
169  Eigen::Ref<Eigen::MatrixXd>
170  J_pp,
171  Eigen::Ref<Eigen::MatrixXd>
172  J_pu,
173  Eigen::Ref<Eigen::MatrixXd>
174  J_uu,
175  Eigen::Ref<Eigen::MatrixXd>
176  J_up)
177 {
178  assert(this->_element.getDimension() == GlobalDim);
179 
180  typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
181  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
182  pressure_size);
183 
184  typename ShapeMatricesTypePressure::NodalMatrixType storage_p =
185  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
186  pressure_size);
187 
188  typename ShapeMatricesTypeDisplacement::template MatrixType<
190  Kup = ShapeMatricesTypeDisplacement::template MatrixType<
191  displacement_size, pressure_size>::Zero(displacement_size,
192  pressure_size);
193 
194  double const& dt = _process_data.dt;
195  auto const& gravity_vec = _process_data.specific_body_force;
196 
197  SpatialPosition x_position;
198  x_position.setElementID(_element.getID());
199 
200  unsigned const n_integration_points = _ip_data.size();
201  for (unsigned ip = 0; ip < n_integration_points; ip++)
202  {
203  x_position.setIntegrationPoint(ip);
204 
205  auto& ip_data = _ip_data[ip];
206  auto const& ip_w = ip_data.integration_weight;
207  auto const& N_u = ip_data.N_u;
208  auto const& dNdx_u = ip_data.dNdx_u;
209  auto const& N_p = ip_data.N_p;
210  auto const& dNdx_p = ip_data.dNdx_p;
211  auto const& H_u = ip_data.H_u;
212 
213  auto const x_coord =
214  interpolateXCoordinate<ShapeFunctionDisplacement,
216  N_u);
217  auto const B =
219  ShapeFunctionDisplacement::NPOINTS,
220  typename BMatricesType::BMatrixType>(
221  dNdx_u, N_u, x_coord, _is_axially_symmetric);
222 
223  auto const& eps_prev = ip_data.eps_prev;
224  auto const& sigma_eff_prev = ip_data.sigma_eff_prev;
225  auto& sigma_eff = ip_data.sigma_eff;
226 
227  auto& eps = ip_data.eps;
228  auto& state = ip_data.material_state_variables;
229 
230  auto q = ip_data.darcy_velocity.head(GlobalDim);
231 
232  auto const alpha = _process_data.biot_coefficient(t, x_position)[0];
233  auto const rho_sr = _process_data.solid_density(t, x_position)[0];
234  auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
235  auto const porosity = _process_data.porosity(t, x_position)[0];
236 
237  double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
238  auto const& identity2 =
240 
241  eps.noalias() = B * u;
242 
243  auto&& solution = _ip_data[ip].solid_material.integrateStress(
244  t, x_position, _process_data.dt, eps_prev, eps, sigma_eff_prev,
245  *state, _process_data.reference_temperature);
246 
247  if (!solution)
248  {
249  OGS_FATAL("Computation of local constitutive relation failed.");
250  }
251 
253  std::tie(sigma_eff, state, C) = std::move(*solution);
254 
255  J_uu.noalias() += B.transpose() * C * B * ip_w;
256 
257  rhs_u.noalias() -= B.transpose() * sigma_eff * ip_w;
258  rhs_u.noalias() -= -H_u.transpose() * rho * gravity_vec * ip_w;
259 
260  //
261  // displacement equation, pressure part
262  //
263  Kup.noalias() += B.transpose() * alpha * identity2 * N_p * ip_w;
264 
265  //
266  // pressure equation, pressure part.
267  //
268  if (!_process_data.deactivate_matrix_in_flow) // Only for hydraulically
269  // active matrix
270  {
271  double const k_over_mu =
272  _process_data.intrinsic_permeability(t, x_position)[0] /
273  _process_data.fluid_viscosity(t, x_position)[0];
274  double const S = _process_data.specific_storage(t, x_position)[0];
275 
276  q.noalias() = -k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
277 
278  laplace_p.noalias() +=
279  dNdx_p.transpose() * k_over_mu * dNdx_p * ip_w;
280  storage_p.noalias() += N_p.transpose() * S * N_p * ip_w;
281 
282  rhs_p.noalias() +=
283  dNdx_p.transpose() * rho_fr * k_over_mu * gravity_vec * ip_w;
284  }
285  }
286 
287  // displacement equation, pressure part
288  J_up.noalias() -= Kup;
289 
290  // pressure equation, pressure part.
291  J_pp.noalias() += laplace_p + storage_p / dt;
292 
293  // pressure equation, displacement part.
294  J_pu.noalias() += Kup.transpose() / dt;
295 
296  // pressure equation
297  rhs_p.noalias() -=
298  laplace_p * p + storage_p * p_dot + Kup.transpose() * u_dot;
299 
300  // displacement equation
301  rhs_u.noalias() -= -Kup * p;
302 }
303 
304 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
305  typename IntegrationMethod, int GlobalDim>
306 void HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
307  ShapeFunctionPressure,
308  IntegrationMethod, GlobalDim>::
310  Eigen::VectorXd const& local_x)
311 {
312  auto p = const_cast<Eigen::VectorXd&>(local_x).segment(pressure_index,
313  pressure_size);
314  if (_process_data.deactivate_matrix_in_flow)
315  {
317  }
318  auto u = local_x.segment(displacement_index, displacement_size);
319 
321 }
322 
323 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
324  typename IntegrationMethod, int GlobalDim>
325 void HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
326  ShapeFunctionPressure,
327  IntegrationMethod, GlobalDim>::
329  double const t,
330  Eigen::Ref<const Eigen::VectorXd> const& p,
331  Eigen::Ref<const Eigen::VectorXd> const& u)
332 {
333  SpatialPosition x_position;
334  x_position.setElementID(_element.getID());
335 
336  unsigned const n_integration_points = _ip_data.size();
337  for (unsigned ip = 0; ip < n_integration_points; ip++)
338  {
339  x_position.setIntegrationPoint(ip);
340 
341  auto& ip_data = _ip_data[ip];
342 
343  auto const& eps_prev = ip_data.eps_prev;
344  auto const& sigma_eff_prev = ip_data.sigma_eff_prev;
345 
346  auto& eps = ip_data.eps;
347  auto& sigma_eff = ip_data.sigma_eff;
348  auto& state = ip_data.material_state_variables;
349 
350  auto const& N_u = ip_data.N_u;
351  auto const& dNdx_u = ip_data.dNdx_u;
352 
353  auto const x_coord =
354  interpolateXCoordinate<ShapeFunctionDisplacement,
356  N_u);
357  auto const B =
359  ShapeFunctionDisplacement::NPOINTS,
360  typename BMatricesType::BMatrixType>(
361  dNdx_u, N_u, x_coord, _is_axially_symmetric);
362 
363  eps.noalias() = B * u;
364 
365  auto&& solution = _ip_data[ip].solid_material.integrateStress(
366  t, x_position, _process_data.dt, eps_prev, eps, sigma_eff_prev,
367  *state, _process_data.reference_temperature);
368 
369  if (!solution)
370  {
371  OGS_FATAL("Computation of local constitutive relation failed.");
372  }
373 
375  std::tie(sigma_eff, state, C) = std::move(*solution);
376 
377  if (!_process_data.deactivate_matrix_in_flow) // Only for hydraulically
378  // active matrix
379  {
380  double const k_over_mu =
381  _process_data.intrinsic_permeability(t, x_position)[0] /
382  _process_data.fluid_viscosity(t, x_position)[0];
383  auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
384  auto const& gravity_vec = _process_data.specific_body_force;
385  auto const& dNdx_p = ip_data.dNdx_p;
386 
387  ip_data.darcy_velocity.head(GlobalDim).noalias() =
388  -k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
389  }
390  }
391 
392  int n = GlobalDim == 2 ? 4 : 6;
393  Eigen::VectorXd ele_stress = Eigen::VectorXd::Zero(n);
394  Eigen::VectorXd ele_strain = Eigen::VectorXd::Zero(n);
395  Eigen::Vector3d ele_velocity = Eigen::Vector3d::Zero();
396 
397  for (auto const& ip_data : _ip_data)
398  {
399  ele_stress += ip_data.sigma_eff;
400  ele_strain += ip_data.eps;
401  ele_velocity += ip_data.darcy_velocity;
402  }
403 
404  ele_stress /= static_cast<double>(n_integration_points);
405  ele_strain /= static_cast<double>(n_integration_points);
406  ele_velocity /= static_cast<double>(n_integration_points);
407 
408  auto const element_id = _element.getID();
409  (*_process_data.mesh_prop_stress_xx)[_element.getID()] = ele_stress[0];
410  (*_process_data.mesh_prop_stress_yy)[_element.getID()] = ele_stress[1];
411  (*_process_data.mesh_prop_stress_zz)[_element.getID()] = ele_stress[2];
412  (*_process_data.mesh_prop_stress_xy)[_element.getID()] = ele_stress[3];
413  if (GlobalDim == 3)
414  {
415  (*_process_data.mesh_prop_stress_yz)[_element.getID()] = ele_stress[4];
416  (*_process_data.mesh_prop_stress_xz)[_element.getID()] = ele_stress[5];
417  }
418 
419  (*_process_data.mesh_prop_strain_xx)[_element.getID()] = ele_strain[0];
420  (*_process_data.mesh_prop_strain_yy)[_element.getID()] = ele_strain[1];
421  (*_process_data.mesh_prop_strain_zz)[_element.getID()] = ele_strain[2];
422  (*_process_data.mesh_prop_strain_xy)[_element.getID()] = ele_strain[3];
423  if (GlobalDim == 3)
424  {
425  (*_process_data.mesh_prop_strain_yz)[_element.getID()] = ele_strain[4];
426  (*_process_data.mesh_prop_strain_xz)[_element.getID()] = ele_strain[5];
427  }
428 
429  for (unsigned i = 0; i < 3; i++)
430  {
431  (*_process_data.mesh_prop_velocity)[element_id * 3 + i] =
432  ele_velocity[i];
433  }
434 
436  ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
437  GlobalDim>(_element, _is_axially_symmetric, p,
438  *_process_data.mesh_prop_nodal_p);
439 }
440 
441 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
442  typename IntegrationMethod, int GlobalDim>
444  ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
445  GlobalDim>::setPressureOfInactiveNodes(double const t,
446  Eigen::Ref<Eigen::VectorXd>
447  p)
448 {
449  SpatialPosition x_position;
450  x_position.setElementID(_element.getID());
451  for (unsigned i = 0; i < pressure_size; i++)
452  {
453  // only inactive nodes
454  if (_process_data.p_element_status->isActiveNode(_element.getNode(i)))
455  {
456  continue;
457  }
458  x_position.setNodeID(_element.getNodeIndex(i));
459  auto const p0 = (*_process_data.p0)(t, x_position)[0];
460  p[i] = p0;
461  }
462 }
463 
464 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
465  typename IntegrationMethod, int GlobalDim>
467  ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
468  GlobalDim>::setPressureDotOfInactiveNodes(Eigen::Ref<Eigen::VectorXd> p_dot)
469 {
470  for (unsigned i = 0; i < pressure_size; i++)
471  {
472  // only inactive nodes
473  if (_process_data.p_element_status->isActiveNode(_element.getNode(i)))
474  {
475  continue;
476  }
477  p_dot[i] = 0;
478  }
479 }
480 
481 } // namespace HydroMechanics
482 } // namespace LIE
483 } // namespace ProcessLib
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool is_axially_symmetric, IntegrationMethod const &integration_method)
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)
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, KelvinVectorDimensions< DisplacementDim >::value, Eigen::RowMajor > KelvinMatrixType
Definition: KelvinVector.h:59
void setElementID(std::size_t element_id)
virtual unsigned getDimension() const =0
Get dimension of the mesh element.
void assembleBlockMatricesWithJacobian(double const t, Eigen::Ref< const Eigen::VectorXd > const &p, Eigen::Ref< const Eigen::VectorXd > const &p_dot, Eigen::Ref< const Eigen::VectorXd > const &u, Eigen::Ref< const Eigen::VectorXd > const &u_dot, Eigen::Ref< Eigen::VectorXd > rhs_p, Eigen::Ref< Eigen::VectorXd > rhs_u, Eigen::Ref< Eigen::MatrixXd > J_pp, Eigen::Ref< Eigen::MatrixXd > J_pu, Eigen::Ref< Eigen::MatrixXd > J_uu, Eigen::Ref< Eigen::MatrixXd > J_up)
void computeSecondaryVariableConcreteWithVector(double const t, Eigen::VectorXd const &local_x) override
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
Definition: LinearBMatrix.h:41
std::vector< IntegrationPointDataType, Eigen::aligned_allocator< IntegrationPointDataType > > _ip_data
ShapeMatrixPolicyType< ShapeFunctionPressure, GlobalDim > ShapeMatricesTypePressure
static const double q
void computeSecondaryVariableConcreteWithBlockVectors(double const t, Eigen::Ref< const Eigen::VectorXd > const &p, Eigen::Ref< const Eigen::VectorXd > const &u)
std::size_t getNodeIndex(unsigned i) const
Definition: Element.cpp:164
void assembleWithJacobianConcrete(double const t, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_dot, Eigen::VectorXd &local_rhs, Eigen::MatrixXd &local_Jac) override
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
void setNodeID(std::size_t node_id)
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:90
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
#define OGS_FATAL(fmt,...)
Definition: Error.h:71
void setIntegrationPoint(unsigned integration_point)
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition: BMatrixPolicy.h:54
Definition of the ElementStatus class.
const Node * getNode(unsigned i) const
Definition: Element.cpp:144
ShapeMatrixPolicyType< ShapeFunctionDisplacement, GlobalDim > ShapeMatricesTypeDisplacement