OGS
HydroMechanicsFEM.h
Go to the documentation of this file.
1 
11 #pragma once
12 
13 #include <memory>
14 #include <vector>
15 
20 #include "MathLib/KelvinVector.h"
26 #include "ParameterLib/Parameter.h"
31 
32 namespace ProcessLib
33 {
34 namespace HydroMechanics
35 {
36 namespace MPL = MaterialPropertyLib;
37 
38 template <typename BMatricesType, typename ShapeMatricesTypeDisplacement,
39  typename ShapeMatricesTypePressure, int DisplacementDim, int NPoints>
41 {
47  solid_material.createMaterialStateVariables())
48  {
49  }
50 
51  typename ShapeMatricesTypeDisplacement::template MatrixType<
52  DisplacementDim, NPoints * DisplacementDim>
57 
58  typename ShapeMatricesTypeDisplacement::NodalRowVectorType N_u;
59  typename ShapeMatricesTypeDisplacement::GlobalDimNodalMatrixType dNdx_u;
60 
61  typename ShapeMatricesTypePressure::NodalRowVectorType N_p;
62  typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p;
63 
65  std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
66  DisplacementDim>::MaterialStateVariables>
69 
70  // TODO disable in monolithic scheme to save memory
71  double coupling_pressure = std::numeric_limits<double>::quiet_NaN();
75  {
76  eps_prev = eps;
78  material_state_variables->pushBackState();
79  }
80 
83  double const t,
84  ParameterLib::SpatialPosition const& x_position,
85  double const dt,
86  double const temperature)
87  {
88  namespace MPL = MaterialPropertyLib;
89 
90  MPL::VariableArray variable_array;
91  MPL::VariableArray variable_array_prev;
92 
93  auto const null_state = solid_material.createMaterialStateVariables();
94 
96 
97  variable_array[static_cast<int>(MPL::Variable::stress)].emplace<KV>(
98  KV::Zero());
99  variable_array[static_cast<int>(MPL::Variable::mechanical_strain)]
100  .emplace<KV>(KV::Zero());
101  variable_array[static_cast<int>(MPL::Variable::temperature)]
102  .emplace<double>(temperature);
103 
104  variable_array_prev[static_cast<int>(MPL::Variable::stress)]
105  .emplace<KV>(KV::Zero());
106  variable_array_prev[static_cast<int>(MPL::Variable::mechanical_strain)]
107  .emplace<KV>(KV::Zero());
108  variable_array_prev[static_cast<int>(MPL::Variable::temperature)]
109  .emplace<double>(temperature);
110 
111  auto&& solution =
112  solid_material.integrateStress(variable_array_prev, variable_array,
113  t, x_position, dt, *null_state);
114 
115  if (!solution)
116  {
117  OGS_FATAL("Computation of elastic tangent stiffness failed.");
118  }
119 
121  std::move(std::get<2>(*solution));
122 
123  return C;
124  }
125 
126  template <typename DisplacementVectorType>
128  MPL::VariableArray const& variable_array,
129  double const t,
130  ParameterLib::SpatialPosition const& x_position,
131  double const dt,
132  DisplacementVectorType const& /*u*/,
133  double const T)
134  {
135  MaterialPropertyLib::VariableArray variable_array_prev;
136  variable_array_prev[static_cast<int>(
140  variable_array_prev[static_cast<int>(MaterialPropertyLib::Variable::
143  eps_prev);
144  variable_array_prev[static_cast<int>(
146  .emplace<double>(T);
147 
148  auto&& solution = solid_material.integrateStress(
149  variable_array_prev, variable_array, t, x_position, dt,
151 
152  if (!solution)
153  {
154  OGS_FATAL("Computation of local constitutive relation failed.");
155  }
156 
158  std::tie(sigma_eff, material_state_variables, C) = std::move(*solution);
159 
160  return C;
161  }
162 
164 };
165 
168 template <typename ShapeMatrixType>
170 {
171  std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N_u;
172 };
173 
174 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
175  int DisplacementDim>
177  : public LocalAssemblerInterface<DisplacementDim>
178 {
179 public:
182 
183  // Types for pressure.
186 
187  static int const KelvinVectorSize =
190 
191  using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>;
192 
195 
197  MeshLib::Element const& e,
198  std::size_t const /*local_matrix_size*/,
199  NumLib::GenericIntegrationMethod const& integration_method,
200  bool const is_axially_symmetric,
202 
204  std::size_t setIPDataInitialConditions(
205  std::string const& name,
206  double const* values,
207  int const integration_order) override;
208 
209  void assemble(double const /*t*/, double const /*dt*/,
210  std::vector<double> const& /*local_x*/,
211  std::vector<double> const& /*local_xdot*/,
212  std::vector<double>& /*local_M_data*/,
213  std::vector<double>& /*local_K_data*/,
214  std::vector<double>& /*local_rhs_data*/) override
215  {
216  OGS_FATAL(
217  "HydroMechanicsLocalAssembler: assembly without jacobian is not "
218  "implemented.");
219  }
220 
221  void assembleWithJacobian(double const t, double const dt,
222  std::vector<double> const& local_x,
223  std::vector<double> const& local_xdot,
224  std::vector<double>& /*local_M_data*/,
225  std::vector<double>& /*local_K_data*/,
226  std::vector<double>& local_rhs_data,
227  std::vector<double>& local_Jac_data) override;
228 
230  const double t, double const dt, Eigen::VectorXd const& local_x,
231  Eigen::VectorXd const& local_xdot, int const process_id,
232  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
233  std::vector<double>& local_b_data,
234  std::vector<double>& local_Jac_data) override;
235 
236  void initializeConcrete() override
237  {
238  unsigned const n_integration_points =
240 
241  for (unsigned ip = 0; ip < n_integration_points; ip++)
242  {
243  auto& ip_data = _ip_data[ip];
244 
246  if (_process_data.initial_stress != nullptr)
247  {
248  ParameterLib::SpatialPosition const x_position{
249  std::nullopt, _element.getID(), ip,
251  ShapeFunctionDisplacement,
253  _element, ip_data.N_u))};
254 
255  ip_data.sigma_eff =
257  DisplacementDim>((*_process_data.initial_stress)(
258  std::numeric_limits<
259  double>::quiet_NaN() /* time independent */,
260  x_position));
261  }
262 
263  ip_data.pushBackState();
264  }
265  }
266 
267  void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
268  double const /*t*/,
269  double const /*dt*/) override
270  {
271  unsigned const n_integration_points =
273 
274  for (unsigned ip = 0; ip < n_integration_points; ip++)
275  {
276  _ip_data[ip].pushBackState();
277  }
278  }
279 
281  double const t, double const dt, Eigen::VectorXd const& local_xs,
282  Eigen::VectorXd const& local_x_dot) override;
283 
284  void postNonLinearSolverConcrete(std::vector<double> const& local_x,
285  std::vector<double> const& local_xdot,
286  double const t, double const dt,
287  bool const use_monolithic_scheme,
288  int const process_id) override;
289 
290  void setInitialConditionsConcrete(std::vector<double> const& local_x,
291  double const t,
292  bool const use_monolithic_scheme,
293  int const process_id) override;
294 
295  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
296  const unsigned integration_point) const override
297  {
298  auto const& N_u = _secondary_data.N_u[integration_point];
299 
300  // assumes N is stored contiguously in memory
301  return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
302  }
303 
304  // TODO (naumov) This method is same as getIntPtSigma but for arguments and
305  // the ordering of the cache_mat.
306  // There should be only one.
307  std::vector<double> getSigma() const override;
308 
309  std::vector<double> getEpsilon() const override;
310 
311  std::vector<double> const& getIntPtDarcyVelocity(
312  const double t,
313  std::vector<GlobalVector*> const& x,
314  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
315  std::vector<double>& cache) const override;
316 
317  std::vector<double> const& getIntPtSigma(
318  const double /*t*/,
319  std::vector<GlobalVector*> const& /*x*/,
320  std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
321  std::vector<double>& cache) const override
322  {
323  return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
324  _ip_data, &IpData::sigma_eff, cache);
325  }
326 
327  std::vector<double> const& getIntPtEpsilon(
328  const double /*t*/,
329  std::vector<GlobalVector*> const& /*x*/,
330  std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
331  std::vector<double>& cache) const override
332  {
333  return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
334  _ip_data, &IpData::eps, cache);
335  }
336 
337 private:
357  const double t, double const dt, Eigen::VectorXd const& local_x,
358  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data);
359 
384  const double t, double const dt, Eigen::VectorXd const& local_x,
385  Eigen::VectorXd const& local_xdot, std::vector<double>& local_b_data,
386  std::vector<double>& local_Jac_data);
387 
388  unsigned getNumberOfIntegrationPoints() const override;
389 
391  DisplacementDim>::MaterialStateVariables const&
392  getMaterialStateVariablesAt(unsigned integration_point) const override;
393 
394 private:
396 
399  using IpData =
401  ShapeMatricesTypePressure, DisplacementDim,
402  ShapeFunctionDisplacement::NPOINTS>;
403  std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
404 
411 
412  static const int pressure_index = 0;
413  static const int pressure_size = ShapeFunctionPressure::NPOINTS;
414  static const int displacement_index = ShapeFunctionPressure::NPOINTS;
415  static const int displacement_size =
416  ShapeFunctionDisplacement::NPOINTS * DisplacementDim;
417 };
418 
419 } // namespace HydroMechanics
420 } // namespace ProcessLib
421 
422 #include "HydroMechanicsFEM-impl.h"
#define OGS_FATAL(...)
Definition: Error.h:26
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:89
void computeSecondaryVariableConcrete(double const t, double const dt, Eigen::VectorXd const &local_xs, Eigen::VectorXd const &local_x_dot) override
NumLib::GenericIntegrationMethod const & _integration_method
void assemble(double const, double const, std::vector< double > const &, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
HydroMechanicsLocalAssembler(HydroMechanicsLocalAssembler &&)=delete
std::size_t setIPDataInitialConditions(std::string const &name, double const *values, int const integration_order) override
Returns number of read integration points.
Eigen::Matrix< double, KelvinVectorSize, 1 > SymmetricTensor
void postTimestepConcrete(Eigen::VectorXd const &, double const, double const) override
void assembleWithJacobianForPressureEquations(const double t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt(unsigned integration_point) const override
BMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > BMatricesType
void assembleWithJacobianForDeformationEquations(const double t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
std::vector< double > const & getIntPtDarcyVelocity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
void postNonLinearSolverConcrete(std::vector< double > const &local_x, std::vector< double > const &local_xdot, double const t, double const dt, bool const use_monolithic_scheme, int const process_id) override
std::vector< double > const & getIntPtEpsilon(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
HydroMechanicsLocalAssembler(HydroMechanicsLocalAssembler const &)=delete
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
void assembleWithJacobian(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_xdot, std::vector< double > &, std::vector< double > &, std::vector< double > &local_rhs_data, std::vector< double > &local_Jac_data) override
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
HydroMechanicsProcessData< DisplacementDim > & _process_data
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
void assembleWithJacobianForStaggeredScheme(const double t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
std::vector< double > const & getIntPtSigma(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
void setInitialConditionsConcrete(std::vector< double > const &local_x, double const t, bool const use_monolithic_scheme, int const process_id) override
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
Definition: VariableType.h:110
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:48
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:23
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType
Definition: KelvinVector.h:56
Eigen::Matrix< double, Eigen::MatrixBase< Derived >::RowsAtCompileTime, 1 > symmetricTensorToKelvinVector(Eigen::MatrixBase< Derived > const &v)
Definition: KelvinVector.h:170
static const double t
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
BMatricesType::KelvinMatrixType updateConstitutiveRelation(MPL::VariableArray const &variable_array, double const t, ParameterLib::SpatialPosition const &x_position, double const dt, DisplacementVectorType const &, double const T)
ShapeMatricesTypeDisplacement::GlobalDimNodalMatrixType dNdx_u
std::unique_ptr< typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables > material_state_variables
ShapeMatricesTypeDisplacement::NodalRowVectorType N_u
ShapeMatricesTypeDisplacement::template MatrixType< DisplacementDim, NPoints *DisplacementDim > N_u_op
ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > computeElasticTangentStiffness(double const t, ParameterLib::SpatialPosition const &x_position, double const dt, double const temperature)
ShapeMatricesTypePressure::NodalRowVectorType N_p
BMatricesType::KelvinVectorType sigma_eff_prev
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material
IntegrationPointData(MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material)
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u