OGS
ThermoMechanicalPhaseFieldFEM.h
Go to the documentation of this file.
1 
11 #pragma once
12 
13 #include <memory>
14 #include <vector>
15 
29 
30 namespace ProcessLib
31 {
32 namespace ThermoMechanicalPhaseField
33 {
34 template <typename BMatricesType, typename ShapeMatrixType, int DisplacementDim>
36 {
42  solid_material.createMaterialStateVariables())
43  {
44  }
45 
46  typename ShapeMatrixType::NodalRowVectorType N;
47  typename ShapeMatrixType::GlobalDimNodalMatrixType dNdx;
48 
51 
55  typename ShapeMatrixType::GlobalDimVectorType heatflux;
56 
58  std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
59  DisplacementDim>::MaterialStateVariables>
61 
64 
66  {
67  eps_prev = eps;
68  material_state_variables->pushBackState();
69  }
70 
73 
74  template <typename DisplacementVectorType>
75  void updateConstitutiveRelation(double const t,
77  double const /*dt*/,
78  DisplacementVectorType const& /*u*/,
79  double const alpha,
80  double const delta_T,
81  double const degradation)
82  {
83  eps_m.noalias() = eps - alpha * delta_T * Invariants::identity2;
84 
85  auto linear_elastic_mp =
87  DisplacementDim> const&>(solid_material)
88  .getMaterialProperties();
89 
90  auto const bulk_modulus = linear_elastic_mp.bulk_modulus(t, x);
91  auto const mu = linear_elastic_mp.mu(t, x);
92 
96  DisplacementDim>(degradation, bulk_modulus, mu, eps_m);
97  }
99 };
100 
103 template <typename ShapeMatrixType>
105 {
106  std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N;
107 };
108 
109 template <typename ShapeFunction, typename IntegrationMethod,
110  int DisplacementDim>
113 {
114 private:
115  static constexpr int temperature_index = 0;
116  static constexpr int temperature_size = ShapeFunction::NPOINTS;
117  static constexpr int displacement_index =
119  static constexpr int displacement_size =
120  ShapeFunction::NPOINTS * DisplacementDim;
121  static constexpr int phasefield_index =
123  static constexpr int phasefield_size = ShapeFunction::NPOINTS;
124 
125 public:
128  // Types for displacement.
129  // (Higher order elements = ShapeFunction).
132 
134 
136 
138  typename ShapeMatricesType::template VectorType<temperature_size>;
140  typename ShapeMatricesType::template VectorType<displacement_size>;
142  typename ShapeMatricesType::template VectorType<phasefield_size>;
143  using IpData =
145 
150 
152  MeshLib::Element const& e,
153  std::size_t const /*local_matrix_size*/,
154  bool const is_axially_symmetric,
155  unsigned const integration_order,
157  int const mechanics_related_process_id,
158  int const phase_field_process_id,
159  int const heat_conduction_process_id)
160  : _process_data(process_data),
161  _integration_method(integration_order),
162  _element(e),
163  _is_axially_symmetric(is_axially_symmetric),
164  _mechanics_related_process_id(mechanics_related_process_id),
165  _phase_field_process_id(phase_field_process_id),
166  _heat_conduction_process_id(heat_conduction_process_id)
167  {
168  unsigned const n_integration_points =
169  _integration_method.getNumberOfPoints();
170 
171  _ip_data.reserve(n_integration_points);
172  _secondary_data.N.resize(n_integration_points);
173 
174  auto& solid_material =
176  _process_data.solid_materials,
177  _process_data.material_ids,
178  e.getID());
180  DisplacementDim> const*>(&solid_material))
181  {
182  OGS_FATAL(
183  "For phasefield process only linear elastic solid material "
184  "support is implemented.");
185  }
186 
187  auto const shape_matrices =
189  DisplacementDim>(e, is_axially_symmetric,
191 
193  x_position.setElementID(_element.getID());
194 
195  for (unsigned ip = 0; ip < n_integration_points; ip++)
196  {
197  _ip_data.emplace_back(solid_material);
198  auto& ip_data = _ip_data[ip];
199  ip_data.integration_weight =
200  _integration_method.getWeightedPoint(ip).getWeight() *
201  shape_matrices[ip].integralMeasure * shape_matrices[ip].detJ;
202 
203  static const int kelvin_vector_size =
205  DisplacementDim);
206  ip_data.eps.setZero(kelvin_vector_size);
207  ip_data.eps_prev.resize(kelvin_vector_size);
208  ip_data.eps_m.setZero(kelvin_vector_size);
209  ip_data.C_tensile.setZero(kelvin_vector_size, kelvin_vector_size);
210  ip_data.C_compressive.setZero(kelvin_vector_size,
211  kelvin_vector_size);
212  ip_data.sigma_tensile.setZero(kelvin_vector_size);
213  ip_data.sigma_compressive.setZero(kelvin_vector_size);
214  ip_data.heatflux.setZero(DisplacementDim);
215  ip_data.sigma.setZero(kelvin_vector_size);
216  ip_data.strain_energy_tensile = 0.0;
217  ip_data.elastic_energy = 0.0;
218 
219  ip_data.N = shape_matrices[ip].N;
220  ip_data.dNdx = shape_matrices[ip].dNdx;
221 
222  _secondary_data.N[ip] = shape_matrices[ip].N;
223  }
224  }
225 
226  void assemble(double const /*t*/, double const /*dt*/,
227  std::vector<double> const& /*local_x*/,
228  std::vector<double> const& /*local_xdot*/,
229  std::vector<double>& /*local_M_data*/,
230  std::vector<double>& /*local_K_data*/,
231  std::vector<double>& /*local_rhs_data*/) override
232  {
233  OGS_FATAL(
234  "ThermoMechanicalPhaseFieldLocalAssembler: assembly without "
235  "Jacobian is not implemented.");
236  }
237 
239  double const t, double const dt, Eigen::VectorXd const& local_x,
240  Eigen::VectorXd const& local_xdot, const double dxdot_dx,
241  const double dx_dx, int const process_id,
242  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
243  std::vector<double>& local_b_data,
244  std::vector<double>& local_Jac_data) override;
245 
246  void initializeConcrete() override
247  {
248  unsigned const n_integration_points =
249  _integration_method.getNumberOfPoints();
250 
251  for (unsigned ip = 0; ip < n_integration_points; ip++)
252  {
253  _ip_data[ip].pushBackState();
254  }
255  }
256 
257  void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
258  double const /*t*/, double const /*dt*/) override
259  {
260  unsigned const n_integration_points =
261  _integration_method.getNumberOfPoints();
262 
263  for (unsigned ip = 0; ip < n_integration_points; ip++)
264  {
265  _ip_data[ip].pushBackState();
266  }
267  }
268 
269  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
270  const unsigned integration_point) const override
271  {
272  auto const& N = _secondary_data.N[integration_point];
273 
274  // assumes N is stored contiguously in memory
275  return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
276  }
277 
278 private:
279  std::vector<double> const& getIntPtSigma(
280  const double /*t*/,
281  std::vector<GlobalVector*> const& /*x*/,
282  std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
283  std::vector<double>& cache) const override
284  {
285  return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
286  _ip_data, &IpData::sigma, cache);
287  }
288 
289  std::vector<double> const& getIntPtEpsilon(
290  const double /*t*/,
291  std::vector<GlobalVector*> const& /*x*/,
292  std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
293  std::vector<double>& cache) const override
294  {
295  return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
296  _ip_data, &IpData::eps, cache);
297  }
298 
299  std::vector<double> const& getIntPtHeatFlux(
300  const double /*t*/,
301  std::vector<GlobalVector*> const& /*x*/,
302  std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
303  std::vector<double>& cache) const override
304  {
306 
307  auto const n_integration_points = _ip_data.size();
308 
309  cache.clear();
310  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
311  double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
312  cache, DisplacementDim, n_integration_points);
313 
314  for (unsigned ip = 0; ip < n_integration_points; ++ip)
315  {
316  auto const& heatflux = _ip_data[ip].heatflux;
317 
318  for (typename KelvinVectorType::Index component = 0;
319  component < DisplacementDim;
320  ++component)
321  { // x, y, z components
322  cache_mat(component, ip) = heatflux[component];
323  }
324  }
325 
326  return cache;
327  }
328 
330  double const t, double const dt, Eigen::VectorXd const& local_x,
331  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data);
332 
334  double const t, double const dt, Eigen::VectorXd const& local_x,
335  Eigen::VectorXd const& local_xdot, std::vector<double>& local_b_data,
336  std::vector<double>& local_Jac_data);
337 
339  double const t, Eigen::VectorXd const& local_x,
340  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data);
341 
343 
344  std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
345 
346  IntegrationMethod _integration_method;
350 
353 
356 
359 };
360 
361 } // namespace ThermoMechanicalPhaseField
362 } // namespace ProcessLib
363 
#define OGS_FATAL(...)
Definition: Error.h:26
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:82
void setElementID(std::size_t element_id)
VectorType< _number_of_dof > NodalForceVectorType
Rhs residual.
Definition: BMatrixPolicy.h:44
VectorType< _kelvin_vector_size > KelvinVectorType
Definition: BMatrixPolicy.h:48
void assemble(double const, double const, std::vector< double > const &, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
void assembleWithJacobianForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, const double dxdot_dx, const double dx_dx, 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 & getIntPtEpsilon(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
int const _mechanics_related_process_id
ID of the processes that contains mechanical process.
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)
void assembleWithJacobianForPhaseFieldEquations(double const t, Eigen::VectorXd const &local_x, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
typename ShapeMatricesType::template VectorType< temperature_size > TemperatureVector
typename ShapeMatricesType::template VectorType< displacement_size > DeformationVector
typename ShapeMatricesType::template VectorType< phasefield_size > PhaseFieldVector
std::vector< double > const & getIntPtSigma(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
void assembleWithJacobianForHeatConductionEquations(double const 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)
ThermoMechanicalPhaseFieldLocalAssembler(ThermoMechanicalPhaseFieldLocalAssembler const &)=delete
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 & getIntPtHeatFlux(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
void postTimestepConcrete(Eigen::VectorXd const &, double const, double const) override
ThermoMechanicalPhaseFieldLocalAssembler(ThermoMechanicalPhaseFieldLocalAssembler &&)=delete
void assembleWithJacobianForDeformationEquations(double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
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)
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, 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::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:32
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
NumLib::ShapeMatrices< NodalRowVectorType, DimNodalMatrixType, DimMatrixType, GlobalDimNodalMatrixType > ShapeMatrices
VectorType< GlobalDim > GlobalDimVectorType
static Eigen::Matrix< double, KelvinVectorSize, 1 > const identity2
Kelvin mapping of 2nd order identity tensor.
Definition: KelvinVector.h:76
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N
Definition: SecondaryData.h:28
void updateConstitutiveRelation(double const t, ParameterLib::SpatialPosition const &x, double const, DisplacementVectorType const &, double const alpha, double const delta_T, double const degradation)
std::unique_ptr< typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables > material_state_variables
IntegrationPointData(MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material)
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N