OGS
PhaseFieldFEM.h
Go to the documentation of this file.
1 
11 #pragma once
12 
13 #include <memory>
14 #include <vector>
15 
26 #include "PhaseFieldProcessData.h"
30 
31 namespace ProcessLib
32 {
33 namespace PhaseField
34 {
35 template <typename BMatricesType, typename ShapeMatrixType, int DisplacementDim>
37 {
43  solid_material.createMaterialStateVariables())
44  {
45  }
46 
47  typename ShapeMatrixType::NodalRowVectorType N;
48  typename ShapeMatrixType::GlobalDimNodalMatrixType dNdx;
49 
51 
55 
57  std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
58  DisplacementDim>::MaterialStateVariables>
60 
64 
66  {
69  eps_prev = eps;
70  material_state_variables->pushBackState();
71  }
72 
73  template <typename DisplacementVectorType>
74  void updateConstitutiveRelation(double const t,
76  double const /*dt*/,
77  DisplacementVectorType const& /*u*/,
78  double const degradation,
79  EnergySplitModel const energy_split_model)
80  {
81  auto linear_elastic_mp =
83  DisplacementDim> const&>(solid_material)
84  .getMaterialProperties();
85 
86  auto const bulk_modulus = linear_elastic_mp.bulk_modulus(t, x);
87  auto const mu = linear_elastic_mp.mu(t, x);
88 
89  switch (energy_split_model)
90  {
92  {
95  MaterialLib::Solids::Phasefield::
96  calculateIsotropicDegradedStress<DisplacementDim>(
97  degradation, bulk_modulus, mu, eps);
98  break;
99  }
101  {
104  MaterialLib::Solids::Phasefield::
105  calculateVolDevDegradedStress<DisplacementDim>(
106  degradation, bulk_modulus, mu, eps);
107  break;
108  }
110  {
113  elastic_energy) = MaterialLib::Solids::Phasefield::
114  calculateIsotropicDegradedStressWithRankineEnergy<
115  DisplacementDim>(degradation, bulk_modulus, mu, eps);
116  break;
117  }
118  }
119 
122  }
123 };
124 
127 template <typename ShapeMatrixType>
129 {
130  std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N;
131 };
132 
133 template <typename ShapeFunction, int DisplacementDim>
135 {
136 private:
137  static constexpr int displacement_index = 0;
138  static constexpr int displacement_size =
139  ShapeFunction::NPOINTS * DisplacementDim;
140  static constexpr int phasefield_index =
142  static constexpr int phasefield_size = ShapeFunction::NPOINTS;
143 
144 public:
147  // Types for displacement.
148  // (Higher order elements = ShapeFunction).
151 
153 
155  typename ShapeMatricesType::template VectorType<displacement_size>;
157  typename ShapeMatricesType::template VectorType<phasefield_size>;
159  typename ShapeMatricesType::template MatrixType<phasefield_size,
163  using IpData =
165 
168 
170  MeshLib::Element const& e,
171  std::size_t const /*local_matrix_size*/,
172  NumLib::GenericIntegrationMethod const& integration_method,
173  bool const is_axially_symmetric,
175  : _process_data(process_data),
176  _integration_method(integration_method),
177  _element(e),
178  _is_axially_symmetric(is_axially_symmetric)
179  {
180  unsigned const n_integration_points =
182 
183  _ip_data.reserve(n_integration_points);
184  _secondary_data.N.resize(n_integration_points);
185 
186  auto& solid_material =
188  _process_data.solid_materials,
189  _process_data.material_ids,
190  e.getID());
192  DisplacementDim> const*>(&solid_material))
193  {
194  OGS_FATAL(
195  "For phasefield process only linear elastic solid material "
196  "support is implemented.");
197  }
198 
199  auto const shape_matrices =
201  DisplacementDim>(e, is_axially_symmetric,
203 
205  x_position.setElementID(_element.getID());
206 
207  for (unsigned ip = 0; ip < n_integration_points; ip++)
208  {
209  _ip_data.emplace_back(solid_material);
210  auto& ip_data = _ip_data[ip];
211  ip_data.integration_weight =
213  shape_matrices[ip].integralMeasure * shape_matrices[ip].detJ;
214 
215  static const int kelvin_vector_size =
217  DisplacementDim);
218  ip_data.eps.setZero(kelvin_vector_size);
219  ip_data.eps_prev.resize(kelvin_vector_size);
220  ip_data.C_tensile.setZero(kelvin_vector_size, kelvin_vector_size);
221  ip_data.C_compressive.setZero(kelvin_vector_size,
222  kelvin_vector_size);
223  ip_data.sigma_tensile.setZero(kelvin_vector_size);
224  ip_data.sigma_compressive.setZero(kelvin_vector_size);
225  ip_data.sigma.setZero(kelvin_vector_size);
226  ip_data.strain_energy_tensile = 0.0;
227  ip_data.elastic_energy = 0.0;
228 
229  ip_data.N = shape_matrices[ip].N;
230  ip_data.dNdx = shape_matrices[ip].dNdx;
231 
232  _secondary_data.N[ip] = shape_matrices[ip].N;
233  }
234  }
235 
236  void assemble(double const /*t*/, double const /*dt*/,
237  std::vector<double> const& /*local_x*/,
238  std::vector<double> const& /*local_xdot*/,
239  std::vector<double>& /*local_M_data*/,
240  std::vector<double>& /*local_K_data*/,
241  std::vector<double>& /*local_rhs_data*/) override
242  {
243  OGS_FATAL(
244  "PhaseFieldLocalAssembler: assembly without jacobian is not "
245  "implemented.");
246  }
247 
249  double const t, double const dt, Eigen::VectorXd const& local_x,
250  Eigen::VectorXd const& local_xdot, int const process_id,
251  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
252  std::vector<double>& local_b_data,
253  std::vector<double>& local_Jac_data) override;
254 
255  void initializeConcrete() override
256  {
257  unsigned const n_integration_points =
259 
260  for (unsigned ip = 0; ip < n_integration_points; ip++)
261  {
262  _ip_data[ip].pushBackState();
263  }
264  }
265 
266  void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
267  double const /*t*/,
268  double const /*dt*/) override
269  {
270  unsigned const n_integration_points =
272 
273  for (unsigned ip = 0; ip < n_integration_points; ip++)
274  {
275  _ip_data[ip].pushBackState();
276  }
277  }
278 
280  std::size_t mesh_item_id,
281  std::vector<
282  std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
283  dof_tables,
284  GlobalVector const& x, double const t, double& crack_volume,
285  CoupledSolutionsForStaggeredScheme const* const cpl_xs) override;
286 
287  void computeEnergy(
288  std::size_t mesh_item_id,
289  std::vector<
290  std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
291  dof_tables,
292  GlobalVector const& x, double const t, double& elastic_energy,
293  double& surface_energy, double& pressure_work,
294  CoupledSolutionsForStaggeredScheme const* const cpl_xs) override;
295 
296  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
297  const unsigned integration_point) const override
298  {
299  auto const& N = _secondary_data.N[integration_point];
300 
301  // assumes N is stored contiguously in memory
302  return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
303  }
304 
305 private:
306  std::vector<double> const& getIntPtSigma(
307  const double /*t*/,
308  std::vector<GlobalVector*> const& /*x*/,
309  std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
310  std::vector<double>& cache) const override
311  {
312  return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
313  _ip_data, &IpData::sigma, cache);
314  }
315 
316  std::vector<double> const& getIntPtEpsilon(
317  const double /*t*/,
318  std::vector<GlobalVector*> const& /*x*/,
319  std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
320  std::vector<double>& cache) const override
321  {
322  return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
323  _ip_data, &IpData::eps, cache);
324  }
325 
327  double const t, double const dt, Eigen::VectorXd const& local_x,
328  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data);
329 
331  double const t, double const dt, Eigen::VectorXd const& local_x,
332  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data);
333 
335 
336  std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
337 
342 
343  static const int phase_process_id = 1;
344  static const int mechanics_process_id = 0;
345 };
346 
347 } // namespace PhaseField
348 } // namespace ProcessLib
349 
350 #include "PhaseFieldFEM-impl.h"
#define OGS_FATAL(...)
Definition: Error.h:26
Global vector based on Eigen vector.
Definition: EigenVector.h:28
double getWeight() const
Definition: WeightedPoint.h:80
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:89
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
void setElementID(std::size_t element_id)
VectorType< _number_of_dof > NodalForceVectorType
Rhs residual.
Definition: BMatrixPolicy.h:44
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
void postTimestepConcrete(Eigen::VectorXd const &, double const, double const) override
typename ShapeMatricesType::template VectorType< phasefield_size > PhaseFieldVector
PhaseFieldProcessData< DisplacementDim > & _process_data
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 assembleWithJacobianForDeformationEquations(double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
void assembleWithJacobianPhaseFieldEquations(double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
SecondaryData< typename ShapeMatrices::ShapeType > _secondary_data
typename BMatricesType::NodalForceVectorType NodalForceVectorType
NumLib::GenericIntegrationMethod const & _integration_method
typename ShapeMatricesType::template VectorType< displacement_size > DeformationVector
ShapeMatrixPolicyType< ShapeFunction, DisplacementDim > ShapeMatricesType
PhaseFieldLocalAssembler(PhaseFieldLocalAssembler const &)=delete
void assembleWithJacobianForStaggeredScheme(double const 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< IpData, Eigen::aligned_allocator< IpData > > _ip_data
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
void computeEnergy(std::size_t mesh_item_id, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap >> const &dof_tables, GlobalVector const &x, double const t, double &elastic_energy, double &surface_energy, double &pressure_work, CoupledSolutionsForStaggeredScheme const *const cpl_xs) override
typename ShapeMatricesType::template MatrixType< phasefield_size, phasefield_size > PhaseFieldMatrix
PhaseFieldLocalAssembler(MeshLib::Element const &e, std::size_t const, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, PhaseFieldProcessData< DisplacementDim > &process_data)
PhaseFieldLocalAssembler(PhaseFieldLocalAssembler &&)=delete
void computeCrackIntegral(std::size_t mesh_item_id, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap >> const &dof_tables, GlobalVector const &x, double const t, double &crack_volume, CoupledSolutionsForStaggeredScheme const *const cpl_xs) override
std::vector< double > const & getIntPtEpsilon(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
typename ShapeMatricesType::NodalMatrixType NodalMatrixType
typename ShapeMatricesType::NodalVectorType NodalVectorType
std::vector< double > const & getIntPtSigma(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
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
static const double t
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)
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
NumLib::ShapeMatrices< NodalRowVectorType, DimNodalMatrixType, DimMatrixType, GlobalDimNodalMatrixType > ShapeMatrices
VectorType< ShapeFunction::NPOINTS > NodalVectorType
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N
Definition: SecondaryData.h:28
BMatricesType::KelvinVectorType sigma_tensile
Definition: PhaseFieldFEM.h:52
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material
Definition: PhaseFieldFEM.h:56
ShapeMatrixType::NodalRowVectorType N
Definition: PhaseFieldFEM.h:47
BMatricesType::KelvinVectorType eps_prev
Definition: PhaseFieldFEM.h:50
BMatricesType::KelvinMatrixType C_compressive
Definition: PhaseFieldFEM.h:61
BMatricesType::KelvinVectorType sigma_compressive
Definition: PhaseFieldFEM.h:52
ShapeMatrixType::GlobalDimNodalMatrixType dNdx
Definition: PhaseFieldFEM.h:48
BMatricesType::KelvinVectorType eps
Definition: PhaseFieldFEM.h:50
BMatricesType::KelvinVectorType sigma
Definition: PhaseFieldFEM.h:53
void updateConstitutiveRelation(double const t, ParameterLib::SpatialPosition const &x, double const, DisplacementVectorType const &, double const degradation, EnergySplitModel const energy_split_model)
Definition: PhaseFieldFEM.h:74
IntegrationPointData(MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material)
Definition: PhaseFieldFEM.h:38
BMatricesType::KelvinMatrixType C_tensile
Definition: PhaseFieldFEM.h:61
std::unique_ptr< typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables > material_state_variables
Definition: PhaseFieldFEM.h:59
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N