OGS 6.2.0-97-g4a610c866
ThermoMechanicalPhaseFieldFEM.h
Go to the documentation of this file.
1 
10 #pragma once
11 
12 #include <memory>
13 #include <vector>
14 
25 
28 
29 namespace ProcessLib
30 {
31 namespace ThermoMechanicalPhaseField
32 {
33 template <typename BMatricesType, typename ShapeMatrixType, int DisplacementDim>
35 {
39  : solid_material(solid_material),
41  solid_material.createMaterialStateVariables())
42  {
43  }
44 
45  typename ShapeMatrixType::NodalRowVectorType N;
46  typename ShapeMatrixType::GlobalDimNodalMatrixType dNdx;
47 
50 
52  sigma;
54  typename ShapeMatrixType::GlobalDimVectorType heatflux;
55 
57  std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
58  DisplacementDim>::MaterialStateVariables>
60 
63 
65  {
66  eps_prev = eps;
67  material_state_variables->pushBackState();
68  }
69 
72 
73  template <typename DisplacementVectorType>
74  void updateConstitutiveRelation(double const t,
76  double const /*dt*/,
77  DisplacementVectorType const& /*u*/,
78  double const alpha,
79  double const delta_T,
80  double const degradation)
81  {
82  eps_m.noalias() = eps - alpha * delta_T * Invariants::identity2;
83 
84  auto linear_elastic_mp =
86  DisplacementDim> const&>(solid_material)
87  .getMaterialProperties();
88 
89  auto const bulk_modulus = linear_elastic_mp.bulk_modulus(t, x);
90  auto const mu = linear_elastic_mp.mu(t, x);
91 
92  std::tie(sigma, sigma_tensile, C_tensile, C_compressive,
93  strain_energy_tensile, elastic_energy) =
95  DisplacementDim>(degradation, bulk_modulus, mu, eps_m);
96  }
98 };
99 
102 template <typename ShapeMatrixType>
104 {
105  std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N;
106 };
107 
108 template <typename ShapeFunction, typename IntegrationMethod,
109  int DisplacementDim>
112 {
113 public:
114  using ShapeMatricesType =
116  // Types for displacement.
117  // (Higher order elements = ShapeFunction).
120 
122 
124 
129 
131  MeshLib::Element const& e,
132  std::size_t const /*local_matrix_size*/,
133  bool const is_axially_symmetric,
134  unsigned const integration_order,
136  int const mechanics_related_process_id,
137  int const phase_field_process_id,
138  int const heat_conduction_process_id)
139  : _process_data(process_data),
140  _integration_method(integration_order),
141  _element(e),
142  _is_axially_symmetric(is_axially_symmetric),
143  _mechanics_related_process_id(mechanics_related_process_id),
144  _phase_field_process_id(phase_field_process_id),
145  _heat_conduction_process_id(heat_conduction_process_id)
146  {
147  unsigned const n_integration_points =
148  _integration_method.getNumberOfPoints();
149 
150  _ip_data.reserve(n_integration_points);
151  _secondary_data.N.resize(n_integration_points);
152 
153  auto& solid_material =
155  _process_data.solid_materials,
156  _process_data.material_ids,
157  e.getID());
159  DisplacementDim> const*>(&solid_material))
160  {
161  OGS_FATAL(
162  "For phasefield process only linear elastic solid material "
163  "support is implemented.");
164  }
165 
166 
167  auto const shape_matrices =
168  initShapeMatrices<ShapeFunction, ShapeMatricesType,
169  IntegrationMethod, DisplacementDim>(
170  e, is_axially_symmetric, _integration_method);
171 
173  x_position.setElementID(_element.getID());
174 
175  for (unsigned ip = 0; ip < n_integration_points; ip++)
176  {
177  _ip_data.emplace_back(solid_material);
178  auto& ip_data = _ip_data[ip];
179  ip_data.integration_weight =
180  _integration_method.getWeightedPoint(ip).getWeight() *
181  shape_matrices[ip].integralMeasure * shape_matrices[ip].detJ;
182 
183  static const int kelvin_vector_size =
185  DisplacementDim>::value;
186  ip_data.eps.setZero(kelvin_vector_size);
187  ip_data.eps_prev.resize(kelvin_vector_size);
188  ip_data.eps_m.setZero(kelvin_vector_size);
189  ip_data.C_tensile.setZero(kelvin_vector_size, kelvin_vector_size);
190  ip_data.C_compressive.setZero(kelvin_vector_size,
191  kelvin_vector_size);
192  ip_data.sigma_tensile.setZero(kelvin_vector_size);
193  ip_data.sigma_compressive.setZero(kelvin_vector_size);
194  ip_data.heatflux.setZero(DisplacementDim);
195  ip_data.sigma.setZero(kelvin_vector_size);
196  ip_data.strain_energy_tensile = 0.0;
197  ip_data.elastic_energy = 0.0;
198 
199  ip_data.N = shape_matrices[ip].N;
200  ip_data.dNdx = shape_matrices[ip].dNdx;
201 
202  _secondary_data.N[ip] = shape_matrices[ip].N;
203  }
204  }
205 
206  void assemble(double const /*t*/, std::vector<double> const& /*local_x*/,
207  std::vector<double>& /*local_M_data*/,
208  std::vector<double>& /*local_K_data*/,
209  std::vector<double>& /*local_rhs_data*/) override
210  {
211  OGS_FATAL(
212  "ThermoMechanicalPhaseFieldLocalAssembler: assembly without "
213  "Jacobian is not implemented.");
214  }
215 
216  void assembleWithJacobianForStaggeredScheme(
217  double const t, std::vector<double> const& local_xdot,
218  const double dxdot_dx, const double dx_dx,
219  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
220  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
221  LocalCoupledSolutions const& local_coupled_solutions) override;
222 
223  void preTimestepConcrete(std::vector<double> const& /*local_x*/,
224  double const /*t*/,
225  double const /*delta_t*/) override
226  {
227  unsigned const n_integration_points =
228  _integration_method.getNumberOfPoints();
229 
230  for (unsigned ip = 0; ip < n_integration_points; ip++)
231  {
232  _ip_data[ip].pushBackState();
233  }
234  }
235 
236  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
237  const unsigned integration_point) const override
238  {
239  auto const& N = _secondary_data.N[integration_point];
240 
241  // assumes N is stored contiguously in memory
242  return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
243  }
244 
245 private:
246  std::vector<double> const& getIntPtSigma(
247  const double /*t*/,
248  GlobalVector const& /*current_solution*/,
249  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
250  std::vector<double>& cache) const override
251  {
252  static const int kelvin_vector_size =
254  DisplacementDim>::value;
255  auto const num_intpts = _ip_data.size();
256 
257  cache.clear();
258  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
259  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
260  cache, kelvin_vector_size, num_intpts);
261 
262  for (unsigned ip = 0; ip < num_intpts; ++ip)
263  {
264  auto const& sigma = _ip_data[ip].sigma;
265  cache_mat.col(ip) =
267  }
268 
269  return cache;
270  }
271 
272  std::vector<double> const& getIntPtEpsilon(
273  const double /*t*/,
274  GlobalVector const& /*current_solution*/,
275  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
276  std::vector<double>& cache) const override
277  {
278  auto const kelvin_vector_size =
280  DisplacementDim>::value;
281  auto const num_intpts = _ip_data.size();
282 
283  cache.clear();
284  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
285  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
286  cache, kelvin_vector_size, num_intpts);
287 
288  for (unsigned ip = 0; ip < num_intpts; ++ip)
289  {
290  auto const& eps = _ip_data[ip].eps;
291  cache_mat.col(ip) =
293  }
294 
295  return cache;
296  }
297 
298  std::vector<double> const& getIntPtHeatFlux(
299  const double /*t*/,
300  GlobalVector const& /*current_solution*/,
301  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
302  std::vector<double>& cache) const override
303  {
305 
306  auto const num_intpts = _ip_data.size();
307 
308  cache.clear();
309  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
310  double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
311  cache, DisplacementDim, num_intpts);
312 
313  for (unsigned ip = 0; ip < num_intpts; ++ip)
314  {
315  auto const& heatflux = _ip_data[ip].heatflux;
316 
317  for (typename KelvinVectorType::Index component = 0;
318  component < DisplacementDim;
319  ++component)
320  { // x, y, z components
321  cache_mat(component, ip) = heatflux[component];
322  }
323  }
324 
325  return cache;
326  }
327 
328  void assembleWithJacobianForDeformationEquations(
329  double const t, std::vector<double> const& local_xdot,
330  const double dxdot_dx, const double dx_dx,
331  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
332  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
333  LocalCoupledSolutions const& local_coupled_solutions);
334 
335  void assembleWithJacobianForHeatConductionEquations(
336  double const t, std::vector<double> const& local_xdot,
337  const double dxdot_dx, const double dx_dx,
338  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
339  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
340  LocalCoupledSolutions const& local_coupled_solutions);
341 
342  void assembleWithJacobianForPhaseFieldEquations(
343  double const t, std::vector<double> const& local_xdot,
344  const double dxdot_dx, const double dx_dx,
345  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
346  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
347  LocalCoupledSolutions const& local_coupled_solutions);
348 
350 
351  std::vector<
353  Eigen::aligned_allocator<IntegrationPointData<
354  BMatricesType, ShapeMatricesType, DisplacementDim>>>
356 
357  IntegrationMethod _integration_method;
361 
362  static const int temperature_index = 0;
363  static const int temperature_size = ShapeFunction::NPOINTS;
364  static const int phasefield_index = ShapeFunction::NPOINTS;
365  static const int phasefield_size = ShapeFunction::NPOINTS;
366  static const int displacement_index = 2 * ShapeFunction::NPOINTS;
367  static const int displacement_size =
368  ShapeFunction::NPOINTS * DisplacementDim;
369 
372 
375 
378 };
379 
380 } // namespace ThermoMechanicalPhaseField
381 } // namespace ProcessLib
382 
void assemble(double const, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
static Eigen::Matrix< double, KelvinVectorSize, 1 > const identity2
Kelvin mapping of 2nd order identity tensor.
Definition: KelvinVector.h:79
std::tuple< MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >, double, double > calculateDegradedStress(double const degradation, double const bulk_modulus, double const mu, MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &eps)
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:23
void setElementID(std::size_t element_id)
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
void preTimestepConcrete(std::vector< double > const &, double const, double const) override
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material
IntegrationPointData(MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material)
VectorType< _number_of_dof > NodalForceVectorType
Rhs residual.
Definition: BMatrixPolicy.h:43
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool is_axially_symmetric, IntegrationMethod const &integration_method)
MechanicsBase< DisplacementDim > & selectSolidConstitutiveRelation(std::map< int, std::unique_ptr< MechanicsBase< DisplacementDim >>> const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, KelvinVectorDimensions< DisplacementDim >::value, Eigen::RowMajor > KelvinMatrixType
Definition: KelvinVector.h:59
ThermoMechanicalPhaseFieldLocalAssembler(MeshLib::Element const &e, std::size_t const, bool const is_axially_symmetric, unsigned const integration_order, ThermoMechanicalPhaseFieldProcessData< DisplacementDim > &process_data, int const mechanics_related_process_id, int const phase_field_process_id, int const heat_conduction_process_id)
std::unique_ptr< typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables > material_state_variables
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:31
int const _mechanics_related_process_id
ID of the processes that contains mechanical process.
Coordinates mapping matrices at particular location.
Definition: ShapeMatrices.h:45
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
VectorType< GlobalDim > GlobalDimVectorType
std::vector< IntegrationPointData< BMatricesType, ShapeMatricesType, DisplacementDim >, Eigen::aligned_allocator< IntegrationPointData< BMatricesType, ShapeMatricesType, DisplacementDim > > > _ip_data
std::vector< double > const & getIntPtHeatFlux(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:90
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
std::vector< double > const & getIntPtSigma(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
void updateConstitutiveRelation(double const t, ParameterLib::SpatialPosition const &x, double const, DisplacementVectorType const &, double const alpha, double const delta_T, double const degradation)
std::vector< double > const & getIntPtEpsilon(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:49