OGS 6.2.0-97-g4a610c866
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 preTimestepConcrete(std::vector<double> const& /*local_x*/,
159  double const /*t*/,
160  double const /*delta_t*/) override
161  {
162  unsigned const n_integration_points =
163  _integration_method.getNumberOfPoints();
164 
165  for (unsigned ip = 0; ip < n_integration_points; ip++)
166  {
167  _ip_data[ip].pushBackState();
168  }
169  }
170 
171  void computeSecondaryVariableConcrete(
172  double const t, std::vector<double> const& local_x) override;
173  void postNonLinearSolverConcrete(std::vector<double> const& local_x,
174  double const t,
175  bool const use_monolithic_scheme) override;
176 
177  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
178  const unsigned integration_point) const override
179  {
180  auto const& N_u = _secondary_data.N_u[integration_point];
181 
182  // assumes N is stored contiguously in memory
183  return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
184  }
185 
186  std::vector<double> const& getIntPtSigmaXX(
187  const double /*t*/,
188  GlobalVector const& /*current_solution*/,
189  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
190  std::vector<double>& cache) const override
191  {
192  return getIntPtSigma(cache, 0);
193  }
194 
195  std::vector<double> const& getIntPtSigmaYY(
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, 1);
202  }
203 
204  std::vector<double> const& getIntPtSigmaZZ(
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, 2);
211  }
212 
213  std::vector<double> const& getIntPtSigmaXY(
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, 3);
220  }
221 
222  std::vector<double> const& getIntPtSigmaXZ(
223  const double /*t*/,
224  GlobalVector const& /*current_solution*/,
225  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
226  std::vector<double>& cache) const override
227  {
228  assert(DisplacementDim == 3);
229  return getIntPtSigma(cache, 4);
230  }
231 
232  std::vector<double> const& getIntPtSigmaYZ(
233  const double /*t*/,
234  GlobalVector const& /*current_solution*/,
235  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
236  std::vector<double>& cache) const override
237  {
238  assert(DisplacementDim == 3);
239  return getIntPtSigma(cache, 5);
240  }
241 
242  std::vector<double> const& getIntPtEpsilonXX(
243  const double /*t*/,
244  GlobalVector const& /*current_solution*/,
245  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
246  std::vector<double>& cache) const override
247  {
248  return getIntPtEpsilon(cache, 0);
249  }
250 
251  std::vector<double> const& getIntPtEpsilonYY(
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, 1);
258  }
259 
260  std::vector<double> const& getIntPtEpsilonZZ(
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, 2);
267  }
268 
269  std::vector<double> const& getIntPtEpsilonXY(
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, 3);
276  }
277 
278  std::vector<double> const& getIntPtEpsilonXZ(
279  const double /*t*/,
280  GlobalVector const& /*current_solution*/,
281  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
282  std::vector<double>& cache) const override
283  {
284  assert(DisplacementDim == 3);
285  return getIntPtEpsilon(cache, 4);
286  }
287 
288  std::vector<double> const& getIntPtEpsilonYZ(
289  const double /*t*/,
290  GlobalVector const& /*current_solution*/,
291  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
292  std::vector<double>& cache) const override
293  {
294  assert(DisplacementDim == 3);
295  return getIntPtEpsilon(cache, 5);
296  }
297 
298  std::vector<double> const& getIntPtDarcyVelocity(
299  const double t,
300  GlobalVector const& current_solution,
301  NumLib::LocalToGlobalIndexMap const& dof_table,
302  std::vector<double>& cache) const override;
303 
304 private:
305  std::vector<double> const& getIntPtSigma(std::vector<double>& cache,
306  std::size_t const component) const
307  {
308  cache.clear();
309  cache.reserve(_ip_data.size());
310 
311  for (auto const& ip_data : _ip_data)
312  {
313  if (component < 3)
314  { // xx, yy, zz components
315  cache.push_back(ip_data.sigma_eff[component]);
316  }
317  else
318  { // mixed xy, yz, xz components
319  cache.push_back(ip_data.sigma_eff[component] / std::sqrt(2));
320  }
321  }
322 
323  return cache;
324  }
325 
326  std::vector<double> const& getIntPtEpsilon(
327  std::vector<double>& cache, std::size_t const component) const
328  {
329  cache.clear();
330  cache.reserve(_ip_data.size());
331 
332  for (auto const& ip_data : _ip_data)
333  {
334  cache.push_back(ip_data.eps[component]);
335  }
336 
337  return cache;
338  }
339 
365  void assembleWithJacobianForDeformationEquations(
366  double const t, std::vector<double> const& local_xdot,
367  const double dxdot_dx, const double dx_dx,
368  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
369  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
370  LocalCoupledSolutions const& local_coupled_solutions);
371 
400  void assembleWithJacobianForPressureEquations(
401  double const t, std::vector<double> const& local_xdot,
402  const double dxdot_dx, const double dx_dx,
403  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
404  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
405  LocalCoupledSolutions const& local_coupled_solutions);
406 
407 private:
409 
410  using BMatricesType =
412  using IpData =
414  ShapeMatricesTypePressure, DisplacementDim,
415  ShapeFunctionDisplacement::NPOINTS>;
416  std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
417 
418  IntegrationMethod _integration_method;
422  typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType>
424 
425  static const int pressure_index = 0;
426  static const int pressure_size = ShapeFunctionPressure::NPOINTS;
427  static const int displacement_index = ShapeFunctionPressure::NPOINTS;
428  static const int displacement_size =
429  ShapeFunctionDisplacement::NPOINTS * DisplacementDim;
430 };
431 
432 } // namespace HydroMechanics
433 } // namespace ProcessLib
434 
435 #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 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
void preTimestepConcrete(std::vector< double > const &, double const, double const) override
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