OGS
ThermoRichardsMechanicsFEM.h
Go to the documentation of this file.
1 
11 #pragma once
12 
13 #include <memory>
14 #include <vector>
15 
16 #include "ConstitutiveSetting.h"
17 #include "IntegrationPointData.h"
20 #include "MathLib/KelvinVector.h"
26 #include "ParameterLib/Parameter.h"
31 
32 namespace ProcessLib
33 {
34 namespace ThermoRichardsMechanics
35 {
36 namespace MPL = MaterialPropertyLib;
37 
40 template <typename ShapeMatrixType>
42 {
43  std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N_u;
44 };
45 
46 template <typename ShapeFunctionDisplacement, typename ShapeFunction,
47  int DisplacementDim>
49  : public LocalAssemblerInterface<DisplacementDim>
50 {
51  static constexpr int temperature_index = 0;
52  static constexpr int temperature_size = ShapeFunction::NPOINTS;
53  static constexpr int pressure_index = temperature_size;
54  static constexpr int pressure_size = ShapeFunction::NPOINTS;
55  static constexpr int displacement_index = 2 * ShapeFunction::NPOINTS;
56  static constexpr int displacement_size =
57  ShapeFunctionDisplacement::NPOINTS * DisplacementDim;
58 
59 public:
62  // Note: temperature variable uses the same shape functions as that are used
63  // by pressure variable.
66 
69 
70  using BMatricesType =
73 
75  ShapeMatricesType, DisplacementDim,
76  ShapeFunctionDisplacement::NPOINTS>;
77 
78  static int const KelvinVectorSize =
81 
82  using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>;
83 
88 
90  MeshLib::Element const& e,
91  std::size_t const /*local_matrix_size*/,
92  NumLib::GenericIntegrationMethod const& integration_method,
93  bool const is_axially_symmetric,
95 
97  std::size_t setIPDataInitialConditions(
98  std::string const& name,
99  double const* values,
100  int const integration_order) override;
101 
102  void setInitialConditionsConcrete(std::vector<double> const& local_x,
103  double const t,
104  bool const use_monolithic_scheme,
105  int const process_id) override;
106 
108  {
110 
111  static auto constexpr local_matrix_dim =
113 
114  template <Eigen::Index rows, Eigen::Index cols>
115  using Mat =
116  typename ShapeMatricesTypeDisplacement::template MatrixType<rows,
117  cols>;
118  using Vec = typename ShapeMatricesTypeDisplacement::template VectorType<
120 
121  public:
122  void setZero()
123  {
124  M_TT = NodalMatrix::Zero(temperature_size, temperature_size);
125  M_Tp = NodalMatrix::Zero(temperature_size, pressure_size);
126  K_TT = NodalMatrix::Zero(temperature_size, temperature_size);
127  K_Tp = NodalMatrix::Zero(temperature_size, pressure_size);
128  dK_TT_dp = NodalMatrix::Zero(temperature_size, pressure_size);
129 
132 
133  M_pT = NodalMatrix::Zero(pressure_size, temperature_size);
134 
135  K_pp = NodalMatrix::Zero(pressure_size, pressure_size);
136  K_pT = NodalMatrix::Zero(pressure_size, temperature_size);
137 
138  storage_p_a_p = NodalMatrix::Zero(pressure_size, pressure_size);
139  storage_p_a_S_Jpp = NodalMatrix::Zero(pressure_size, pressure_size);
140  storage_p_a_S = NodalMatrix::Zero(pressure_size, pressure_size);
141 
144  res = Vec::Zero(local_matrix_dim);
145  }
146 
148  {
149  M_TT += other.M_TT;
150  M_Tp += other.M_Tp;
151  K_TT += other.K_TT;
152  K_Tp += other.K_Tp;
153  dK_TT_dp += other.dK_TT_dp;
154 
155  M_pu += other.M_pu;
156 
157  M_pT += other.M_pT;
158 
159  K_pp += other.K_pp;
160  K_pT += other.K_pT;
161 
162  storage_p_a_p += other.storage_p_a_p;
164  storage_p_a_S += other.storage_p_a_S;
165 
166  Jac += other.Jac;
167  res += other.res;
168 
169  return *this;
170  }
171 
172  LocalMatrices& operator*=(double const a)
173  {
174  M_TT *= a;
175  M_Tp *= a;
176  K_TT *= a;
177  K_Tp *= a;
178  dK_TT_dp *= a;
179 
180  M_pu *= a;
181 
182  M_pT *= a;
183 
184  K_pp *= a;
185  K_pT *= a;
186 
187  storage_p_a_p *= a;
188  storage_p_a_S_Jpp *= a;
189  storage_p_a_S *= a;
190 
191  Jac *= a;
192  res *= a;
193 
194  return *this;
195  }
196 
202 
204 
206 
209 
213 
217 
221  };
222 
223  void assembleWithJacobian(double const t, double const dt,
224  std::vector<double> const& local_x,
225  std::vector<double> const& local_xdot,
226  std::vector<double>& /*local_M_data*/,
227  std::vector<double>& /*local_K_data*/,
228  std::vector<double>& local_rhs_data,
229  std::vector<double>& local_Jac_data) override;
230 
231 private:
233  double const t, double const dt,
234  ParameterLib::SpatialPosition const& x_position,
235  std::vector<double> const& local_x,
236  std::vector<double> const& local_xdot, IpData const& ip_data,
239  StatefulData<DisplacementDim>& current_state,
240  StatefulData<DisplacementDim> const& prev_state,
242  OutputData<DisplacementDim>& output_data) const;
243 
244  void addToLocalMatrixData(double const dt,
245  std::vector<double> const& local_x,
246  std::vector<double> const& local_xdot,
247  LocalMatrices const& loc_mat,
248  std::vector<double>& local_rhs_data,
249  std::vector<double>& local_Jac_data) const;
250 
251  void massLumping(LocalMatrices& loc_mat) const;
252 
255  auto localDOF(std::vector<double> const& local_dof_data) const
256  {
257  static_assert(temperature_size == pressure_size);
258 
259  using NodalTOrPVec =
260  typename ShapeMatricesType::template VectorType<temperature_size>;
261  using NodalDispVec =
262  typename ShapeMatricesTypeDisplacement::template VectorType<
264 
265  return std::tuple<Eigen::Map<NodalTOrPVec const>,
266  Eigen::Map<NodalTOrPVec const>,
267  Eigen::Map<NodalDispVec const>>(
268  {local_dof_data.data() + temperature_index, temperature_size},
269  {local_dof_data.data() + pressure_index, pressure_size},
270  {local_dof_data.data() + displacement_index, displacement_size});
271  };
272 
273 public:
274  void initializeConcrete() override
275  {
276  unsigned const n_integration_points =
278 
279  for (unsigned ip = 0; ip < n_integration_points; ip++)
280  {
281  // Set initial stress from parameter.
282  if (process_data_.initial_stress != nullptr)
283  {
284  ParameterLib::SpatialPosition const x_position{
285  std::nullopt, element_.getID(), ip,
287  ShapeFunctionDisplacement,
289  element_, ip_data_[ip].N_u))};
290 
291  current_states_[ip].s_mech_data.sigma_eff =
293  DisplacementDim>((*process_data_.initial_stress)(
294  std::numeric_limits<
295  double>::quiet_NaN() /* time independent */,
296  x_position));
297  }
298 
299  material_states_[ip].pushBackState();
300  }
301 
303  }
304 
305  void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
306  double const /*t*/,
307  double const /*dt*/) override
308  {
309  unsigned const n_integration_points =
311 
312  for (unsigned ip = 0; ip < n_integration_points; ip++)
313  {
314  // TODO re-evaluate part of the assembly in order to be consistent?
315  material_states_[ip].pushBackState();
316  }
317 
319  }
320 
322  double const t, double const dt, Eigen::VectorXd const& local_x,
323  Eigen::VectorXd const& local_x_dot) override;
324 
325  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
326  const unsigned integration_point) const override
327  {
328  auto const& N_u = secondary_data_.N_u[integration_point];
329 
330  // assumes N is stored contiguously in memory
331  return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
332  }
333 
334  std::vector<double> getSigma() const override;
335 
336  std::vector<double> const& getIntPtDarcyVelocity(
337  const double t,
338  std::vector<GlobalVector*> const& x,
339  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
340  std::vector<double>& cache) const override;
341 
342  std::vector<double> getSaturation() const override;
343  std::vector<double> const& getIntPtSaturation(
344  const double t,
345  std::vector<GlobalVector*> const& x,
346  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
347  std::vector<double>& cache) const override;
348 
349  std::vector<double> getPorosity() const override;
350  std::vector<double> const& getIntPtPorosity(
351  const double t,
352  std::vector<GlobalVector*> const& x,
353  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
354  std::vector<double>& cache) const override;
355 
356  std::vector<double> getTransportPorosity() const override;
357  std::vector<double> const& getIntPtTransportPorosity(
358  const double t,
359  std::vector<GlobalVector*> const& x,
360  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
361  std::vector<double>& cache) const override;
362 
363  std::vector<double> const& getIntPtSigma(
364  const double t,
365  std::vector<GlobalVector*> const& x,
366  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
367  std::vector<double>& cache) const override;
368 
369  std::vector<double> getSwellingStress() const override;
370  std::vector<double> const& getIntPtSwellingStress(
371  const double t,
372  std::vector<GlobalVector*> const& x,
373  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
374  std::vector<double>& cache) const override;
375 
376  std::vector<double> getEpsilon() const override;
377  std::vector<double> const& getIntPtEpsilon(
378  const double t,
379  std::vector<GlobalVector*> const& x,
380  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
381  std::vector<double>& cache) const override;
382 
383  std::vector<double> const& getIntPtDryDensitySolid(
384  const double t,
385  std::vector<GlobalVector*> const& x,
386  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
387  std::vector<double>& cache) const override;
388 
389  std::vector<double> const& getIntPtLiquidDensity(
390  const double t,
391  std::vector<GlobalVector*> const& x,
392  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
393  std::vector<double>& cache) const override;
394 
395  std::vector<double> const& getIntPtViscosity(
396  const double t,
397  std::vector<GlobalVector*> const& x,
398  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
399  std::vector<double>& cache) const override;
400 
401 private:
402  unsigned getNumberOfIntegrationPoints() const override;
403 
405  DisplacementDim>::MaterialStateVariables const&
406  getMaterialStateVariablesAt(unsigned integration_point) const override;
407 
408 private:
410 
411  std::vector<StatefulData<DisplacementDim>>
412  current_states_; // TODO maybe do not store but rather re-evaluate for
413  // state update
414  std::vector<StatefulData<DisplacementDim>> prev_states_;
415 
416  // Material state is special, because it contains both the current and the
417  // old state.
418  std::vector<MaterialStateData<DisplacementDim>> material_states_;
419 
420  std::vector<IpData> ip_data_;
421 
428 
430 
431  std::vector<OutputData<DisplacementDim>> output_data_;
432 
433  static auto block_uu(auto& mat)
434  {
435  return mat.template block<displacement_size, displacement_size>(
437  }
438  static auto block_up(auto& mat)
439  {
440  return mat.template block<displacement_size, pressure_size>(
442  }
443  static auto block_uT(auto& mat)
444  {
445  return mat.template block<displacement_size, temperature_size>(
447  }
448  static auto block_pu(auto& mat)
449  {
450  return mat.template block<pressure_size, displacement_size>(
452  }
453  static auto block_pp(auto& mat)
454  {
455  return mat.template block<pressure_size, pressure_size>(pressure_index,
457  }
458  static auto block_pT(auto& mat)
459  {
460  return mat.template block<pressure_size, temperature_size>(
462  }
463  static auto block_Tp(auto& mat)
464  {
465  return mat.template block<temperature_size, pressure_size>(
467  }
468  static auto block_TT(auto& mat)
469  {
470  return mat.template block<temperature_size, temperature_size>(
472  }
473 
474  static auto block_u(auto& vec)
475  {
476  return vec.template segment<displacement_size>(displacement_index);
477  }
478  static auto block_p(auto& vec)
479  {
480  return vec.template segment<pressure_size>(pressure_index);
481  }
482  static auto block_T(auto& vec)
483  {
484  return vec.template segment<temperature_size>(temperature_index);
485  }
486 };
487 
488 } // namespace ThermoRichardsMechanics
489 } // namespace ProcessLib
490 
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:89
VectorType< _kelvin_vector_size > KelvinVectorType
Definition: BMatrixPolicy.h:48
typename ShapeMatricesTypeDisplacement::template MatrixType< rows, cols > Mat
typename ShapeMatricesTypeDisplacement::template VectorType< local_matrix_dim > Vec
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt(unsigned integration_point) const override
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > secondary_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
std::vector< double > const & getIntPtTransportPorosity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
void computeSecondaryVariableConcrete(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_dot) override
ThermoRichardsMechanicsProcessData< DisplacementDim > & process_data_
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
std::vector< double > const & getIntPtPorosity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::vector< double > const & getIntPtDryDensitySolid(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::vector< double > const & getIntPtViscosity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
ThermoRichardsMechanicsLocalAssembler(ThermoRichardsMechanicsLocalAssembler &&)=delete
ShapeMatrixPolicyType< ShapeFunction, DisplacementDim > ShapeMatricesType
std::vector< double > const & getIntPtSaturation(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
void setInitialConditionsConcrete(std::vector< double > const &local_x, double const t, bool const use_monolithic_scheme, int const process_id) override
void addToLocalMatrixData(double const dt, std::vector< double > const &local_x, std::vector< double > const &local_xdot, LocalMatrices const &loc_mat, std::vector< double > &local_rhs_data, std::vector< double > &local_Jac_data) const
void postTimestepConcrete(Eigen::VectorXd const &, double const, double const) override
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material_
std::size_t setIPDataInitialConditions(std::string const &name, double const *values, int const integration_order) override
ThermoRichardsMechanicsLocalAssembler(ThermoRichardsMechanicsLocalAssembler const &)=delete
std::vector< double > const & getIntPtSwellingStress(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
std::vector< double > const & getIntPtSigma(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
void assembleWithJacobianSingleIP(double const t, double const dt, ParameterLib::SpatialPosition const &x_position, std::vector< double > const &local_x, std::vector< double > const &local_xdot, IpData const &ip_data, ConstitutiveSetting< DisplacementDim > &CS, MaterialPropertyLib::Medium &medium, LocalMatrices &out, StatefulData< DisplacementDim > &current_state, StatefulData< DisplacementDim > const &prev_state, MaterialStateData< DisplacementDim > &mat_state, OutputData< DisplacementDim > &output_data) const
std::vector< double > const & getIntPtLiquidDensity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::vector< double > const & getIntPtEpsilon(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:23
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)
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
MatrixType< GlobalDim, GlobalDim > GlobalDimMatrixType
VectorType< GlobalDim > GlobalDimVectorType
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u
Data that is needed for output purposes, but not directly for the assembly.
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u
Data whose state must be tracked by the TRM process.