OGS 6.2.1-97-g73d1aeda3
HydroMechanicsFEM.h
Go to the documentation of this file.
1 
10 #pragma once
11 
12 #include <memory>
13 #include <vector>
14 
17 #include "MathLib/KelvinVector.h"
21 #include "ParameterLib/Parameter.h"
26 
29 
30 namespace ProcessLib
31 {
32 namespace HydroMechanics
33 {
34 template <typename BMatricesType, typename ShapeMatrixTypeDisplacement,
35  typename ShapeMatricesTypePressure, int DisplacementDim, int NPoints>
37 {
41  : solid_material(solid_material),
43  solid_material.createMaterialStateVariables())
44  {
45  }
46 
47  typename ShapeMatrixTypeDisplacement::template MatrixType<
48  DisplacementDim, NPoints * DisplacementDim>
52 
53  typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u;
54  typename ShapeMatrixTypeDisplacement::GlobalDimNodalMatrixType dNdx_u;
55 
56  typename ShapeMatricesTypePressure::NodalRowVectorType N_p;
57  typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p;
58 
60  std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
61  DisplacementDim>::MaterialStateVariables>
64 
66  {
67  eps_prev = eps;
68  sigma_eff_prev = sigma_eff;
69  material_state_variables->pushBackState();
70  }
71 
72  template <typename DisplacementVectorType>
74  double const t,
75  ParameterLib::SpatialPosition const& x_position,
76  double const dt,
77  DisplacementVectorType const& /*u*/,
78  double const T)
79  {
80  auto&& solution = solid_material.integrateStress(
81  t, x_position, dt, eps_prev, eps, sigma_eff_prev,
83 
84  if (!solution)
85  {
86  OGS_FATAL("Computation of local constitutive relation failed.");
87  }
88 
90  std::tie(sigma_eff, material_state_variables, C) = std::move(*solution);
91 
92  return C;
93  }
94 
96 };
97 
100 template <typename ShapeMatrixType>
102 {
103  std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N_u;
104 };
105 
106 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
107  typename IntegrationMethod, int DisplacementDim>
109 {
110 public:
113 
114  // Types for pressure.
115  using ShapeMatricesTypePressure =
117 
118  static int const KelvinVectorSize =
121 
124 
126  MeshLib::Element const& e,
127  std::size_t const /*local_matrix_size*/,
128  bool const is_axially_symmetric,
129  unsigned const integration_order,
131 
132  void assemble(double const /*t*/, std::vector<double> const& /*local_x*/,
133  std::vector<double>& /*local_M_data*/,
134  std::vector<double>& /*local_K_data*/,
135  std::vector<double>& /*local_rhs_data*/) override
136  {
137  OGS_FATAL(
138  "HydroMechanicsLocalAssembler: assembly without jacobian is not "
139  "implemented.");
140  }
141 
142  void assembleWithJacobian(double const t,
143  std::vector<double> const& local_x,
144  std::vector<double> const& local_xdot,
145  const double /*dxdot_dx*/, const double /*dx_dx*/,
146  std::vector<double>& /*local_M_data*/,
147  std::vector<double>& /*local_K_data*/,
148  std::vector<double>& local_rhs_data,
149  std::vector<double>& local_Jac_data) override;
150 
151  void assembleWithJacobianForStaggeredScheme(
152  double const t, std::vector<double> const& local_xdot,
153  const double dxdot_dx, const double dx_dx,
154  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
155  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
156  LocalCoupledSolutions const& local_coupled_solutions) override;
157 
158  void initializeConcrete() override
159  {
160  unsigned const n_integration_points =
161  _integration_method.getNumberOfPoints();
162 
163  for (unsigned ip = 0; ip < n_integration_points; ip++)
164  {
165  _ip_data[ip].pushBackState();
166  }
167  }
168 
169  void postTimestepConcrete(std::vector<double> const& /*local_x*/) override
170  {
171  unsigned const n_integration_points =
172  _integration_method.getNumberOfPoints();
173 
174  for (unsigned ip = 0; ip < n_integration_points; ip++)
175  {
176  _ip_data[ip].pushBackState();
177  }
178  }
179 
180  void computeSecondaryVariableConcrete(
181  double const t, std::vector<double> const& local_x) override;
182  void postNonLinearSolverConcrete(std::vector<double> const& local_x,
183  double const t,
184  bool const use_monolithic_scheme) override;
185 
186  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
187  const unsigned integration_point) const override
188  {
189  auto const& N_u = _secondary_data.N_u[integration_point];
190 
191  // assumes N is stored contiguously in memory
192  return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
193  }
194 
195  std::vector<double> const& getIntPtSigmaXX(
196  const double /*t*/,
197  GlobalVector const& /*current_solution*/,
198  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
199  std::vector<double>& cache) const override
200  {
201  return getIntPtSigma(cache, 0);
202  }
203 
204  std::vector<double> const& getIntPtSigmaYY(
205  const double /*t*/,
206  GlobalVector const& /*current_solution*/,
207  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
208  std::vector<double>& cache) const override
209  {
210  return getIntPtSigma(cache, 1);
211  }
212 
213  std::vector<double> const& getIntPtSigmaZZ(
214  const double /*t*/,
215  GlobalVector const& /*current_solution*/,
216  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
217  std::vector<double>& cache) const override
218  {
219  return getIntPtSigma(cache, 2);
220  }
221 
222  std::vector<double> const& getIntPtSigmaXY(
223  const double /*t*/,
224  GlobalVector const& /*current_solution*/,
225  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
226  std::vector<double>& cache) const override
227  {
228  return getIntPtSigma(cache, 3);
229  }
230 
231  std::vector<double> const& getIntPtSigmaXZ(
232  const double /*t*/,
233  GlobalVector const& /*current_solution*/,
234  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
235  std::vector<double>& cache) const override
236  {
237  assert(DisplacementDim == 3);
238  return getIntPtSigma(cache, 4);
239  }
240 
241  std::vector<double> const& getIntPtSigmaYZ(
242  const double /*t*/,
243  GlobalVector const& /*current_solution*/,
244  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
245  std::vector<double>& cache) const override
246  {
247  assert(DisplacementDim == 3);
248  return getIntPtSigma(cache, 5);
249  }
250 
251  std::vector<double> const& getIntPtEpsilonXX(
252  const double /*t*/,
253  GlobalVector const& /*current_solution*/,
254  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
255  std::vector<double>& cache) const override
256  {
257  return getIntPtEpsilon(cache, 0);
258  }
259 
260  std::vector<double> const& getIntPtEpsilonYY(
261  const double /*t*/,
262  GlobalVector const& /*current_solution*/,
263  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
264  std::vector<double>& cache) const override
265  {
266  return getIntPtEpsilon(cache, 1);
267  }
268 
269  std::vector<double> const& getIntPtEpsilonZZ(
270  const double /*t*/,
271  GlobalVector const& /*current_solution*/,
272  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
273  std::vector<double>& cache) const override
274  {
275  return getIntPtEpsilon(cache, 2);
276  }
277 
278  std::vector<double> const& getIntPtEpsilonXY(
279  const double /*t*/,
280  GlobalVector const& /*current_solution*/,
281  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
282  std::vector<double>& cache) const override
283  {
284  return getIntPtEpsilon(cache, 3);
285  }
286 
287  std::vector<double> const& getIntPtEpsilonXZ(
288  const double /*t*/,
289  GlobalVector const& /*current_solution*/,
290  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
291  std::vector<double>& cache) const override
292  {
293  assert(DisplacementDim == 3);
294  return getIntPtEpsilon(cache, 4);
295  }
296 
297  std::vector<double> const& getIntPtEpsilonYZ(
298  const double /*t*/,
299  GlobalVector const& /*current_solution*/,
300  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
301  std::vector<double>& cache) const override
302  {
303  assert(DisplacementDim == 3);
304  return getIntPtEpsilon(cache, 5);
305  }
306 
307  std::vector<double> const& getIntPtDarcyVelocity(
308  const double t,
309  GlobalVector const& current_solution,
310  NumLib::LocalToGlobalIndexMap const& dof_table,
311  std::vector<double>& cache) const override;
312 
313 private:
314  std::vector<double> const& getIntPtSigma(std::vector<double>& cache,
315  std::size_t const component) const
316  {
317  cache.clear();
318  cache.reserve(_ip_data.size());
319 
320  for (auto const& ip_data : _ip_data)
321  {
322  if (component < 3)
323  { // xx, yy, zz components
324  cache.push_back(ip_data.sigma_eff[component]);
325  }
326  else
327  { // mixed xy, yz, xz components
328  cache.push_back(ip_data.sigma_eff[component] / std::sqrt(2));
329  }
330  }
331 
332  return cache;
333  }
334 
335  std::vector<double> const& getIntPtEpsilon(
336  std::vector<double>& cache, std::size_t const component) const
337  {
338  cache.clear();
339  cache.reserve(_ip_data.size());
340 
341  for (auto const& ip_data : _ip_data)
342  {
343  cache.push_back(ip_data.eps[component]);
344  }
345 
346  return cache;
347  }
348 
374  void assembleWithJacobianForDeformationEquations(
375  double const t, std::vector<double> const& local_xdot,
376  const double dxdot_dx, const double dx_dx,
377  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
378  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
379  LocalCoupledSolutions const& local_coupled_solutions);
380 
409  void assembleWithJacobianForPressureEquations(
410  double const t, std::vector<double> const& local_xdot,
411  const double dxdot_dx, const double dx_dx,
412  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
413  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
414  LocalCoupledSolutions const& local_coupled_solutions);
415 
416 private:
418 
419  using BMatricesType =
421  using IpData =
423  ShapeMatricesTypePressure, DisplacementDim,
424  ShapeFunctionDisplacement::NPOINTS>;
425  std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
426 
427  IntegrationMethod _integration_method;
431  typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType>
433 
434  static const int pressure_index = 0;
435  static const int pressure_size = ShapeFunctionPressure::NPOINTS;
436  static const int displacement_index = ShapeFunctionPressure::NPOINTS;
437  static const int displacement_size =
438  ShapeFunctionDisplacement::NPOINTS * DisplacementDim;
439 };
440 
441 } // namespace HydroMechanics
442 } // namespace ProcessLib
443 
444 #include "HydroMechanicsFEM-impl.h"
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:23
ShapeMatrixTypeDisplacement::GlobalDimNodalMatrixType dNdx_u
std::vector< double > const & getIntPtSigmaXZ(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
std::unique_ptr< typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables > material_state_variables
std::vector< double > const & getIntPtSigmaXX(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
void postTimestepConcrete(std::vector< double > const &) override
void assemble(double const, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
ShapeMatrixTypeDisplacement::template MatrixType< DisplacementDim, NPoints *DisplacementDim > N_u_op
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
std::vector< double > const & getIntPtSigmaXY(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
std::vector< double > const & getIntPtSigmaYY(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
std::vector< double > const & getIntPtEpsilonXY(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, KelvinVectorDimensions< DisplacementDim >::value, Eigen::RowMajor > KelvinMatrixType
Definition: KelvinVector.h:59
std::vector< double > const & getIntPtEpsilonXX(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
std::vector< double > const & getIntPtSigmaYZ(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
std::vector< double > const & getIntPtEpsilon(std::vector< double > &cache, std::size_t const component) const
std::vector< double > const & getIntPtEpsilonZZ(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
HydroMechanicsProcessData< DisplacementDim > & _process_data
std::vector< double > const & getIntPtEpsilonYY(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p
ShapeMatricesTypePressure::NodalRowVectorType N_p
std::vector< double > const & getIntPtEpsilonYZ(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u
ShapeMatrixTypeDisplacement::NodalRowVectorType N_u
boost::optional< std::tuple< KelvinVector, std::unique_ptr< MaterialStateVariables >, KelvinMatrix > > integrateStress(double const t, ParameterLib::SpatialPosition const &x, double const dt, Eigen::Matrix< double, Eigen::Dynamic, 1 > const &eps_prev, Eigen::Matrix< double, Eigen::Dynamic, 1 > const &eps, Eigen::Matrix< double, Eigen::Dynamic, 1 > const &sigma_prev, MaterialStateVariables const &material_state_variables, double const T) const
Definition: MechanicsBase.h:82
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
std::vector< double > const & getIntPtSigma(std::vector< double > &cache, std::size_t const component) const
IntegrationPointData(MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material)
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material
std::vector< double > const & getIntPtSigmaZZ(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
BMatricesType::KelvinMatrixType updateConstitutiveRelation(double const t, ParameterLib::SpatialPosition const &x_position, double const dt, DisplacementVectorType const &, double const T)
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:49
std::vector< double > const & getIntPtEpsilonXZ(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
BMatricesType::KelvinVectorType sigma_eff_prev