OGS 6.1.0-1721-g6382411ad
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,
73  SpatialPosition const& x,
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 
157  SpatialPosition x_position;
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 preTimestepConcrete(std::vector<double> const& /*local_x*/,
211  double const /*t*/,
212  double const /*delta_t*/) override
213  {
214  unsigned const n_integration_points =
215  _integration_method.getNumberOfPoints();
216 
217  for (unsigned ip = 0; ip < n_integration_points; ip++)
218  {
219  _ip_data[ip].pushBackState();
220  }
221  }
222 
223  void computeCrackIntegral(
224  std::size_t mesh_item_id,
225  std::vector<
226  std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
227  dof_tables,
228  GlobalVector const& x, double const t, double& crack_volume,
229  CoupledSolutionsForStaggeredScheme const* const cpl_xs) override;
230 
231  void computeEnergy(
232  std::size_t mesh_item_id,
233  std::vector<
234  std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
235  dof_tables,
236  GlobalVector const& x, double const t, double& elastic_energy,
237  double& surface_energy, double& pressure_work,
238  CoupledSolutionsForStaggeredScheme const* const cpl_xs) override;
239 
240  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
241  const unsigned integration_point) const override
242  {
243  auto const& N = _secondary_data.N[integration_point];
244 
245  // assumes N is stored contiguously in memory
246  return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
247  }
248 
249 private:
250  std::vector<double> const& getIntPtSigma(
251  const double /*t*/,
252  GlobalVector const& /*current_solution*/,
253  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
254  std::vector<double>& cache) const override
255  {
256  static const int kelvin_vector_size =
258  DisplacementDim>::value;
259  auto const num_intpts = _ip_data.size();
260 
261  cache.clear();
262  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
263  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
264  cache, kelvin_vector_size, num_intpts);
265 
266  for (unsigned ip = 0; ip < num_intpts; ++ip)
267  {
268  auto const& sigma = _ip_data[ip].sigma;
269  cache_mat.col(ip) =
271  }
272 
273  return cache;
274  }
275 
276  virtual std::vector<double> const& getIntPtEpsilon(
277  const double /*t*/,
278  GlobalVector const& /*current_solution*/,
279  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
280  std::vector<double>& cache) const override
281  {
282  auto const kelvin_vector_size =
284  DisplacementDim>::value;
285  auto const num_intpts = _ip_data.size();
286 
287  cache.clear();
288  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
289  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
290  cache, kelvin_vector_size, num_intpts);
291 
292  for (unsigned ip = 0; ip < num_intpts; ++ip)
293  {
294  auto const& eps = _ip_data[ip].eps;
295  cache_mat.col(ip) =
297  }
298 
299  return cache;
300  }
301 
302  void assembleWithJacobianPhaseFiledEquations(
303  double const t, std::vector<double> const& local_xdot,
304  const double dxdot_dx, const double dx_dx,
305  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
306  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
307  LocalCoupledSolutions const& local_coupled_solutions);
308 
309  void assembleWithJacobianForDeformationEquations(
310  double const t, std::vector<double> const& local_xdot,
311  const double dxdot_dx, const double dx_dx,
312  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
313  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
314  LocalCoupledSolutions const& local_coupled_solutions);
315 
317 
318  std::vector<
320  Eigen::aligned_allocator<IntegrationPointData<
321  BMatricesType, ShapeMatricesType, DisplacementDim>>>
323 
324  IntegrationMethod _integration_method;
328 
329  static const int phasefield_index = 0;
330  static const int phasefield_size = ShapeFunction::NPOINTS;
331  static const int displacement_index = ShapeFunction::NPOINTS;
332  static const int displacement_size =
333  ShapeFunction::NPOINTS * DisplacementDim;
334 };
335 
336 } // namespace PhaseField
337 } // namespace ProcessLib
338 
339 #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
virtual std::vector< double > const & getIntPtEpsilon(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
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)
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
void preTimestepConcrete(std::vector< double > const &, double const, double const) override
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, KelvinVectorDimensions< DisplacementDim >::value, Eigen::RowMajor > KelvinMatrixType
Definition: KelvinVector.h:59
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.
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
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
void updateConstitutiveRelation(double const t, SpatialPosition const &x, double const, DisplacementVectorType const &, double const degradation)
Definition: PhaseFieldFEM.h:72
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:71
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
PhaseFieldProcessData< DisplacementDim > & _process_data
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:49
BMatricesType::KelvinVectorType eps_prev
Definition: PhaseFieldFEM.h:48