OGS 6.1.0-1721-g6382411ad
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"
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  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  OGS_FATAL("Computation of local constitutive relation failed.");
86 
88  std::tie(sigma_eff, material_state_variables, C) = std::move(*solution);
89 
90  return C;
91  }
92 
94 };
95 
98 template <typename ShapeMatrixType>
100 {
101  std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N_u;
102 };
103 
104 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
105  typename IntegrationMethod, int DisplacementDim>
107 {
108 public:
111 
112  // Types for pressure.
113  using ShapeMatricesTypePressure =
115 
116  static int const KelvinVectorSize =
119 
122 
124  MeshLib::Element const& e,
125  std::size_t const /*local_matrix_size*/,
126  bool const is_axially_symmetric,
127  unsigned const integration_order,
129 
130  void assemble(double const /*t*/, std::vector<double> const& /*local_x*/,
131  std::vector<double>& /*local_M_data*/,
132  std::vector<double>& /*local_K_data*/,
133  std::vector<double>& /*local_rhs_data*/) override
134  {
135  OGS_FATAL(
136  "HydroMechanicsLocalAssembler: assembly without jacobian is not "
137  "implemented.");
138  }
139 
140  void assembleWithJacobian(double const t,
141  std::vector<double> const& local_x,
142  std::vector<double> const& local_xdot,
143  const double /*dxdot_dx*/, const double /*dx_dx*/,
144  std::vector<double>& /*local_M_data*/,
145  std::vector<double>& /*local_K_data*/,
146  std::vector<double>& local_rhs_data,
147  std::vector<double>& local_Jac_data) override;
148 
149  void assembleWithJacobianForStaggeredScheme(
150  double const t, std::vector<double> const& local_xdot,
151  const double dxdot_dx, const double dx_dx,
152  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
153  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
154  LocalCoupledSolutions const& local_coupled_solutions) override;
155 
156  void preTimestepConcrete(std::vector<double> const& /*local_x*/,
157  double const /*t*/,
158  double const /*delta_t*/) 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 computeSecondaryVariableConcrete(
170  double const t, std::vector<double> const& local_x) override;
171  void postNonLinearSolverConcrete(std::vector<double> const& local_x,
172  double const t,
173  bool const use_monolithic_scheme) override;
174 
175  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
176  const unsigned integration_point) const override
177  {
178  auto const& N_u = _secondary_data.N_u[integration_point];
179 
180  // assumes N is stored contiguously in memory
181  return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
182  }
183 
184  std::vector<double> const& getIntPtSigmaXX(
185  const double /*t*/,
186  GlobalVector const& /*current_solution*/,
187  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
188  std::vector<double>& cache) const override
189  {
190  return getIntPtSigma(cache, 0);
191  }
192 
193  std::vector<double> const& getIntPtSigmaYY(
194  const double /*t*/,
195  GlobalVector const& /*current_solution*/,
196  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
197  std::vector<double>& cache) const override
198  {
199  return getIntPtSigma(cache, 1);
200  }
201 
202  std::vector<double> const& getIntPtSigmaZZ(
203  const double /*t*/,
204  GlobalVector const& /*current_solution*/,
205  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
206  std::vector<double>& cache) const override
207  {
208  return getIntPtSigma(cache, 2);
209  }
210 
211  std::vector<double> const& getIntPtSigmaXY(
212  const double /*t*/,
213  GlobalVector const& /*current_solution*/,
214  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
215  std::vector<double>& cache) const override
216  {
217  return getIntPtSigma(cache, 3);
218  }
219 
220  std::vector<double> const& getIntPtSigmaXZ(
221  const double /*t*/,
222  GlobalVector const& /*current_solution*/,
223  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
224  std::vector<double>& cache) const override
225  {
226  assert(DisplacementDim == 3);
227  return getIntPtSigma(cache, 4);
228  }
229 
230  std::vector<double> const& getIntPtSigmaYZ(
231  const double /*t*/,
232  GlobalVector const& /*current_solution*/,
233  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
234  std::vector<double>& cache) const override
235  {
236  assert(DisplacementDim == 3);
237  return getIntPtSigma(cache, 5);
238  }
239 
240  std::vector<double> const& getIntPtEpsilonXX(
241  const double /*t*/,
242  GlobalVector const& /*current_solution*/,
243  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
244  std::vector<double>& cache) const override
245  {
246  return getIntPtEpsilon(cache, 0);
247  }
248 
249  std::vector<double> const& getIntPtEpsilonYY(
250  const double /*t*/,
251  GlobalVector const& /*current_solution*/,
252  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
253  std::vector<double>& cache) const override
254  {
255  return getIntPtEpsilon(cache, 1);
256  }
257 
258  std::vector<double> const& getIntPtEpsilonZZ(
259  const double /*t*/,
260  GlobalVector const& /*current_solution*/,
261  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
262  std::vector<double>& cache) const override
263  {
264  return getIntPtEpsilon(cache, 2);
265  }
266 
267  std::vector<double> const& getIntPtEpsilonXY(
268  const double /*t*/,
269  GlobalVector const& /*current_solution*/,
270  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
271  std::vector<double>& cache) const override
272  {
273  return getIntPtEpsilon(cache, 3);
274  }
275 
276  std::vector<double> const& getIntPtEpsilonXZ(
277  const double /*t*/,
278  GlobalVector const& /*current_solution*/,
279  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
280  std::vector<double>& cache) const override
281  {
282  assert(DisplacementDim == 3);
283  return getIntPtEpsilon(cache, 4);
284  }
285 
286  std::vector<double> const& getIntPtEpsilonYZ(
287  const double /*t*/,
288  GlobalVector const& /*current_solution*/,
289  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
290  std::vector<double>& cache) const override
291  {
292  assert(DisplacementDim == 3);
293  return getIntPtEpsilon(cache, 5);
294  }
295 
296  std::vector<double> const& getIntPtDarcyVelocity(
297  const double t,
298  GlobalVector const& current_solution,
299  NumLib::LocalToGlobalIndexMap const& dof_table,
300  std::vector<double>& cache) const override;
301 
302 private:
303  std::vector<double> const& getIntPtSigma(std::vector<double>& cache,
304  std::size_t const component) const
305  {
306  cache.clear();
307  cache.reserve(_ip_data.size());
308 
309  for (auto const& ip_data : _ip_data)
310  {
311  if (component < 3) // xx, yy, zz components
312  cache.push_back(ip_data.sigma_eff[component]);
313  else // mixed xy, yz, xz components
314  cache.push_back(ip_data.sigma_eff[component] / std::sqrt(2));
315  }
316 
317  return cache;
318  }
319 
320  std::vector<double> const& getIntPtEpsilon(
321  std::vector<double>& cache, std::size_t const component) const
322  {
323  cache.clear();
324  cache.reserve(_ip_data.size());
325 
326  for (auto const& ip_data : _ip_data)
327  {
328  cache.push_back(ip_data.eps[component]);
329  }
330 
331  return cache;
332  }
333 
359  void assembleWithJacobianForDeformationEquations(
360  double const t, std::vector<double> const& local_xdot,
361  const double dxdot_dx, const double dx_dx,
362  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
363  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
364  LocalCoupledSolutions const& local_coupled_solutions);
365 
394  void assembleWithJacobianForPressureEquations(
395  double const t, std::vector<double> const& local_xdot,
396  const double dxdot_dx, const double dx_dx,
397  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
398  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
399  LocalCoupledSolutions const& local_coupled_solutions);
400 
401 private:
403 
404  using BMatricesType =
406  using IpData =
408  ShapeMatricesTypePressure, DisplacementDim,
409  ShapeFunctionDisplacement::NPOINTS>;
410  std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
411 
412  IntegrationMethod _integration_method;
416  typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType>
418 
419  static const int pressure_index = 0;
420  static const int pressure_size = ShapeFunctionPressure::NPOINTS;
421  static const int displacement_index = ShapeFunctionPressure::NPOINTS;
422  static const int displacement_size =
423  ShapeFunctionDisplacement::NPOINTS * DisplacementDim;
424 };
425 
426 } // namespace HydroMechanics
427 } // namespace ProcessLib
428 
429 #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
#define OGS_FATAL(fmt,...)
Definition: Error.h:71
std::vector< double > const & getIntPtSigma(std::vector< double > &cache, std::size_t const component) const
BMatricesType::KelvinMatrixType updateConstitutiveRelation(double const t, SpatialPosition const &x_position, double const dt, DisplacementVectorType const &, double const T)
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
boost::optional< std::tuple< KelvinVector, std::unique_ptr< MaterialStateVariables >, KelvinMatrix > > integrateStress(double const t, ProcessLib::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
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