OGS 6.2.1-97-g73d1aeda3
PhaseFieldFEM.h
Go to the documentation of this file.
1 
10 #pragma once
11 
12 #include <memory>
13 #include <vector>
14 
25 
27 #include "PhaseFieldProcessData.h"
28 
29 namespace ProcessLib
30 {
31 namespace PhaseField
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 
49 
51  sigma;
53 
55  std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
56  DisplacementDim>::MaterialStateVariables>
58 
62 
64  {
65  history_variable_prev =
66  std::max(history_variable_prev, history_variable);
67  eps_prev = eps;
68  material_state_variables->pushBackState();
69  }
70 
71  template <typename DisplacementVectorType>
72  void updateConstitutiveRelation(double const t,
74  double const /*dt*/,
75  DisplacementVectorType const& /*u*/,
76  double const degradation)
77  {
78  auto linear_elastic_mp =
80  DisplacementDim> const&>(solid_material)
81  .getMaterialProperties();
82 
83  auto const bulk_modulus = linear_elastic_mp.bulk_modulus(t, x);
84  auto const mu = linear_elastic_mp.mu(t, x);
85 
86  std::tie(sigma, sigma_tensile, C_tensile, C_compressive,
87  strain_energy_tensile, elastic_energy) =
89  DisplacementDim>(degradation, bulk_modulus, mu, eps);
90  }
92 
93  static constexpr int kelvin_vector_size =
95 };
96 
99 template <typename ShapeMatrixType>
101 {
102  std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N;
103 };
104 
105 template <typename ShapeFunction, typename IntegrationMethod,
106  int DisplacementDim>
108 {
109 public:
110  using ShapeMatricesType =
112  // Types for displacement.
113  // (Higher order elements = ShapeFunction).
116 
118 
121 
123  MeshLib::Element const& e,
124  std::size_t const /*local_matrix_size*/,
125  bool const is_axially_symmetric,
126  unsigned const integration_order,
128  : _process_data(process_data),
129  _integration_method(integration_order),
130  _element(e),
131  _is_axially_symmetric(is_axially_symmetric)
132  {
133  unsigned const n_integration_points =
134  _integration_method.getNumberOfPoints();
135 
136  _ip_data.reserve(n_integration_points);
137  _secondary_data.N.resize(n_integration_points);
138 
139  auto& solid_material =
141  _process_data.solid_materials,
142  _process_data.material_ids,
143  e.getID());
145  DisplacementDim> const*>(&solid_material))
146  {
147  OGS_FATAL(
148  "For phasefield process only linear elastic solid material "
149  "support is implemented.");
150  }
151 
152  auto const shape_matrices =
153  initShapeMatrices<ShapeFunction, ShapeMatricesType,
154  IntegrationMethod, DisplacementDim>(
155  e, is_axially_symmetric, _integration_method);
156 
158  x_position.setElementID(_element.getID());
159 
160  for (unsigned ip = 0; ip < n_integration_points; ip++)
161  {
162  _ip_data.emplace_back(solid_material);
163  auto& ip_data = _ip_data[ip];
164  ip_data.integration_weight =
165  _integration_method.getWeightedPoint(ip).getWeight() *
166  shape_matrices[ip].integralMeasure * shape_matrices[ip].detJ;
167 
168  static const int kelvin_vector_size =
170  DisplacementDim>::value;
171  ip_data.eps.setZero(kelvin_vector_size);
172  ip_data.eps_prev.resize(kelvin_vector_size);
173  ip_data.C_tensile.setZero(kelvin_vector_size, kelvin_vector_size);
174  ip_data.C_compressive.setZero(kelvin_vector_size,
175  kelvin_vector_size);
176  ip_data.sigma_tensile.setZero(kelvin_vector_size);
177  ip_data.sigma_compressive.setZero(kelvin_vector_size);
178  ip_data.history_variable =
179  _process_data.history_field(0, x_position)[0];
180  ip_data.history_variable_prev =
181  _process_data.history_field(0, x_position)[0];
182  ip_data.sigma.setZero(kelvin_vector_size);
183  ip_data.strain_energy_tensile = 0.0;
184  ip_data.elastic_energy = 0.0;
185 
186  ip_data.N = shape_matrices[ip].N;
187  ip_data.dNdx = shape_matrices[ip].dNdx;
188 
189  _secondary_data.N[ip] = shape_matrices[ip].N;
190  }
191  }
192 
193  void assemble(double const /*t*/, std::vector<double> const& /*local_x*/,
194  std::vector<double>& /*local_M_data*/,
195  std::vector<double>& /*local_K_data*/,
196  std::vector<double>& /*local_rhs_data*/) override
197  {
198  OGS_FATAL(
199  "PhaseFieldLocalAssembler: assembly without jacobian is not "
200  "implemented.");
201  }
202 
203  void assembleWithJacobianForStaggeredScheme(
204  double const t, std::vector<double> const& local_xdot,
205  const double dxdot_dx, const double dx_dx,
206  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
207  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
208  LocalCoupledSolutions const& local_coupled_solutions) override;
209 
210  void initializeConcrete() override
211  {
212  unsigned const n_integration_points =
213  _integration_method.getNumberOfPoints();
214 
215  for (unsigned ip = 0; ip < n_integration_points; ip++)
216  {
217  _ip_data[ip].pushBackState();
218  }
219  }
220 
221  void postTimestepConcrete(std::vector<double> const& /*local_x*/
222  ) override
223  {
224  unsigned const n_integration_points =
225  _integration_method.getNumberOfPoints();
226 
227  for (unsigned ip = 0; ip < n_integration_points; ip++)
228  {
229  _ip_data[ip].pushBackState();
230  }
231  }
232 
233  void computeCrackIntegral(
234  std::size_t mesh_item_id,
235  std::vector<
236  std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
237  dof_tables,
238  GlobalVector const& x, double const t, double& crack_volume,
239  CoupledSolutionsForStaggeredScheme const* const cpl_xs) override;
240 
241  void computeEnergy(
242  std::size_t mesh_item_id,
243  std::vector<
244  std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
245  dof_tables,
246  GlobalVector const& x, double const t, double& elastic_energy,
247  double& surface_energy, double& pressure_work,
248  CoupledSolutionsForStaggeredScheme const* const cpl_xs) override;
249 
250  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
251  const unsigned integration_point) const override
252  {
253  auto const& N = _secondary_data.N[integration_point];
254 
255  // assumes N is stored contiguously in memory
256  return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
257  }
258 
259 private:
260  std::vector<double> const& getIntPtSigma(
261  const double /*t*/,
262  GlobalVector const& /*current_solution*/,
263  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
264  std::vector<double>& cache) const override
265  {
266  static const int kelvin_vector_size =
268  DisplacementDim>::value;
269  auto const num_intpts = _ip_data.size();
270 
271  cache.clear();
272  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
273  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
274  cache, kelvin_vector_size, num_intpts);
275 
276  for (unsigned ip = 0; ip < num_intpts; ++ip)
277  {
278  auto const& sigma = _ip_data[ip].sigma;
279  cache_mat.col(ip) =
281  }
282 
283  return cache;
284  }
285 
286  std::vector<double> const& getIntPtEpsilon(
287  const double /*t*/,
288  GlobalVector const& /*current_solution*/,
289  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
290  std::vector<double>& cache) const override
291  {
292  auto const kelvin_vector_size =
294  DisplacementDim>::value;
295  auto const num_intpts = _ip_data.size();
296 
297  cache.clear();
298  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
299  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
300  cache, kelvin_vector_size, num_intpts);
301 
302  for (unsigned ip = 0; ip < num_intpts; ++ip)
303  {
304  auto const& eps = _ip_data[ip].eps;
305  cache_mat.col(ip) =
307  }
308 
309  return cache;
310  }
311 
312  void assembleWithJacobianPhaseFiledEquations(
313  double const t, std::vector<double> const& local_xdot,
314  const double dxdot_dx, const double dx_dx,
315  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
316  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
317  LocalCoupledSolutions const& local_coupled_solutions);
318 
319  void assembleWithJacobianForDeformationEquations(
320  double const t, std::vector<double> const& local_xdot,
321  const double dxdot_dx, const double dx_dx,
322  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
323  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
324  LocalCoupledSolutions const& local_coupled_solutions);
325 
327 
328  std::vector<
330  Eigen::aligned_allocator<IntegrationPointData<
331  BMatricesType, ShapeMatricesType, DisplacementDim>>>
333 
334  IntegrationMethod _integration_method;
338 
339  static const int phasefield_index = 0;
340  static const int phasefield_size = ShapeFunction::NPOINTS;
341  static const int displacement_index = ShapeFunction::NPOINTS;
342  static const int displacement_size =
343  ShapeFunction::NPOINTS * DisplacementDim;
344 };
345 
346 } // namespace PhaseField
347 } // namespace ProcessLib
348 
349 #include "PhaseFieldFEM-impl.h"
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)
BMatricesType::KelvinMatrixType C_tensile
Definition: PhaseFieldFEM.h:59
BMatricesType::KelvinMatrixType C_compressive
Definition: PhaseFieldFEM.h:59
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:23
PhaseFieldLocalAssembler(MeshLib::Element const &e, std::size_t const, bool const is_axially_symmetric, unsigned const integration_order, PhaseFieldProcessData< DisplacementDim > &process_data)
void setElementID(std::size_t element_id)
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)
ShapeMatrixType::NodalRowVectorType N
Definition: PhaseFieldFEM.h:45
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, KelvinVectorDimensions< DisplacementDim >::value, Eigen::RowMajor > KelvinMatrixType
Definition: KelvinVector.h:59
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:31
IntegrationPointData(MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material)
Definition: PhaseFieldFEM.h:36
SecondaryData< typename ShapeMatrices::ShapeType > _secondary_data
void updateConstitutiveRelation(double const t, ParameterLib::SpatialPosition const &x, double const, DisplacementVectorType const &, double const degradation)
Definition: PhaseFieldFEM.h:72
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)
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N
std::unique_ptr< typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables > material_state_variables
Definition: PhaseFieldFEM.h:57
void assemble(double const, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
typename BMatricesType::NodalForceVectorType NodalForceVectorType
BMatricesType::KelvinVectorType sigma_tensile
Definition: PhaseFieldFEM.h:50
BMatricesType::KelvinVectorType sigma
Definition: PhaseFieldFEM.h:50
std::vector< double > const & getIntPtSigma(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
BMatricesType::KelvinVectorType sigma_compressive
Definition: PhaseFieldFEM.h:50
BMatricesType::KelvinVectorType eps
Definition: PhaseFieldFEM.h:48
ShapeMatrixType::GlobalDimNodalMatrixType dNdx
Definition: PhaseFieldFEM.h:46
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
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material
Definition: PhaseFieldFEM.h:54
std::vector< IntegrationPointData< BMatricesType, ShapeMatricesType, DisplacementDim >, Eigen::aligned_allocator< IntegrationPointData< BMatricesType, ShapeMatricesType, DisplacementDim > > > _ip_data
void postTimestepConcrete(std::vector< double > const &) override
PhaseFieldProcessData< DisplacementDim > & _process_data
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
BMatricesType::KelvinVectorType eps_prev
Definition: PhaseFieldFEM.h:48