OGS
PhaseFieldFEM-impl.h
Go to the documentation of this file.
1 
11 #pragma once
12 
15 
16 namespace ProcessLib
17 {
18 namespace PhaseField
19 {
20 template <typename ShapeFunction, int DisplacementDim>
23  double const t, double const dt, Eigen::VectorXd const& local_x,
24  Eigen::VectorXd const& /*local_xdot*/, int const process_id,
25  std::vector<double>& /*local_M_data*/,
26  std::vector<double>& /*local_K_data*/,
27  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
28 {
29  // For the equations with phase field.
30  if (process_id == phase_process_id)
31  {
32  assembleWithJacobianPhaseFieldEquations(t, dt, local_x, local_b_data,
33  local_Jac_data);
34  return;
35  }
36 
37  // For the equations with deformation
38  assembleWithJacobianForDeformationEquations(t, dt, local_x, local_b_data,
39  local_Jac_data);
40 }
41 
42 template <typename ShapeFunction, int DisplacementDim>
45  double const t, double const dt, Eigen::VectorXd const& local_x,
46  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
47 {
48  using DeformationMatrix =
49  typename ShapeMatricesType::template MatrixType<displacement_size,
50  displacement_size>;
51  auto const d = local_x.template segment<phasefield_size>(phasefield_index);
52  auto const u =
53  local_x.template segment<displacement_size>(displacement_index);
54 
55  auto local_Jac = MathLib::createZeroedMatrix<DeformationMatrix>(
56  local_Jac_data, displacement_size, displacement_size);
57 
58  auto local_rhs = MathLib::createZeroedVector<DeformationVector>(
59  local_b_data, displacement_size);
60 
62  x_position.setElementID(_element.getID());
63 
64  auto local_pressure = 0.0;
65  if (_process_data.crack_pressure)
66  {
67  local_pressure = _process_data.unity_pressure;
68  }
69 
70  int const n_integration_points = _integration_method.getNumberOfPoints();
71  for (int ip = 0; ip < n_integration_points; ip++)
72  {
73  x_position.setIntegrationPoint(ip);
74  auto const& w = _ip_data[ip].integration_weight;
75  auto const& N = _ip_data[ip].N;
76  auto const& dNdx = _ip_data[ip].dNdx;
77 
78  auto const x_coord =
79  NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
80  _element, N);
81 
82  auto const& B =
83  LinearBMatrix::computeBMatrix<DisplacementDim,
84  ShapeFunction::NPOINTS,
86  dNdx, N, x_coord, _is_axially_symmetric);
87 
88  auto& eps = _ip_data[ip].eps;
89  eps.noalias() = B * u;
90  double const k = _process_data.residual_stiffness(t, x_position)[0];
91  double const ls = _process_data.crack_length_scale(t, x_position)[0];
92  double const d_ip = N.dot(d);
93  double const degradation =
94  _process_data.degradation_derivative->degradation(d_ip, k, ls);
95  _ip_data[ip].updateConstitutiveRelation(
96  t, x_position, dt, u, degradation,
97  _process_data.energy_split_model);
98 
99  auto& sigma = _ip_data[ip].sigma;
100  auto const& C_tensile = _ip_data[ip].C_tensile;
101  auto const& C_compressive = _ip_data[ip].C_compressive;
102 
103  typename ShapeMatricesType::template MatrixType<DisplacementDim,
104  displacement_size>
105  N_u = ShapeMatricesType::template MatrixType<
106  DisplacementDim, displacement_size>::Zero(DisplacementDim,
107  displacement_size);
108 
109  for (int i = 0; i < DisplacementDim; ++i)
110  {
111  N_u.template block<1, displacement_size / DisplacementDim>(
112  i, i * displacement_size / DisplacementDim)
113  .noalias() = N;
114  }
115 
116  auto const rho_sr = _process_data.solid_density(t, x_position)[0];
117  auto const& b = _process_data.specific_body_force;
118 
119  local_rhs.noalias() -=
120  (B.transpose() * sigma - N_u.transpose() * rho_sr * b -
121  local_pressure * N_u.transpose() * dNdx * d) *
122  w;
123  local_Jac.noalias() +=
124  B.transpose() * (degradation * C_tensile + C_compressive) * B * w;
125  }
126 }
127 
128 template <typename ShapeFunction, int DisplacementDim>
130  assembleWithJacobianPhaseFieldEquations(double const t, double const dt,
131  Eigen::VectorXd const& local_x,
132  std::vector<double>& local_b_data,
133  std::vector<double>& local_Jac_data)
134 {
135  auto const d = local_x.template segment<phasefield_size>(phasefield_index);
136  auto const u =
137  local_x.template segment<displacement_size>(displacement_index);
138 
139  auto local_Jac = MathLib::createZeroedMatrix<PhaseFieldMatrix>(
140  local_Jac_data, phasefield_size, phasefield_size);
141  auto local_rhs = MathLib::createZeroedVector<PhaseFieldVector>(
142  local_b_data, phasefield_size);
143 
145  x_position.setElementID(_element.getID());
146 
147  auto local_pressure = 0.0;
148  if (_process_data.crack_pressure)
149  {
150  local_pressure = _process_data.unity_pressure;
151  }
152  else if (_process_data.hydro_crack)
153  {
154  local_pressure = _process_data.pressure;
155  }
156 
157  double const k = _process_data.residual_stiffness(t, x_position)[0];
158  double const ls = _process_data.crack_length_scale(t, x_position)[0];
159  double const gc = _process_data.crack_resistance(t, x_position)[0];
160 
161  int const n_integration_points = _integration_method.getNumberOfPoints();
162  for (int ip = 0; ip < n_integration_points; ip++)
163  {
164  x_position.setIntegrationPoint(ip);
165  auto const& w = _ip_data[ip].integration_weight;
166  auto const& N = _ip_data[ip].N;
167  auto const& dNdx = _ip_data[ip].dNdx;
168 
169  double const d_ip = N.dot(d);
170  double const degradation =
171  _process_data.degradation_derivative->degradation(d_ip, k, ls);
172  double const degradation_df1 =
173  _process_data.degradation_derivative->degradation_df1(d_ip, ls);
174  double const degradation_df2 =
175  _process_data.degradation_derivative->degradation_df2(d_ip, ls);
176  auto const x_coord =
177  NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
178  _element, N);
179  auto const& B =
180  LinearBMatrix::computeBMatrix<DisplacementDim,
181  ShapeFunction::NPOINTS,
182  typename BMatricesType::BMatrixType>(
183  dNdx, N, x_coord, _is_axially_symmetric);
184 
185  auto& eps = _ip_data[ip].eps;
186  eps.noalias() = B * u;
187  _ip_data[ip].updateConstitutiveRelation(
188  t, x_position, dt, u, degradation,
189  _process_data.energy_split_model);
190 
191  auto const& strain_energy_tensile = _ip_data[ip].strain_energy_tensile;
192 
193  auto& ip_data = _ip_data[ip];
194  ip_data.strain_energy_tensile = strain_energy_tensile;
195 
196  typename ShapeMatricesType::template MatrixType<DisplacementDim,
197  displacement_size>
198  N_u = ShapeMatricesType::template MatrixType<
199  DisplacementDim, displacement_size>::Zero(DisplacementDim,
200  displacement_size);
201 
202  for (int i = 0; i < DisplacementDim; ++i)
203  {
204  N_u.template block<1, displacement_size / DisplacementDim>(
205  i, i * displacement_size / DisplacementDim)
206  .noalias() = N;
207  }
208 
209  local_Jac.noalias() +=
210  (N.transpose() * N * degradation_df2 * strain_energy_tensile) * w;
211 
212  local_rhs.noalias() -=
213  (N.transpose() * degradation_df1 * strain_energy_tensile -
214  local_pressure * dNdx.transpose() * N_u * u) *
215  w;
216 
217  switch (_process_data.phasefield_model)
218  {
220  {
221  auto const local_Jac_AT1 =
222  (gc * 0.75 * dNdx.transpose() * dNdx * ls * w).eval();
223  local_Jac.noalias() += local_Jac_AT1;
224 
225  local_rhs.noalias() -=
226  gc * (-0.375 * N.transpose() / ls) * w + local_Jac_AT1 * d;
227  break;
228  }
230  {
231  auto const local_Jac_AT2 =
232  (gc *
233  (N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) *
234  w)
235  .eval();
236  local_Jac.noalias() += local_Jac_AT2;
237 
238  local_rhs.noalias() -=
239  local_Jac_AT2 * d - gc * (N.transpose() / ls) * w;
240  break;
241  }
243  {
244  auto const local_Jac_COHESIVE =
245  (2.0 / boost::math::double_constants::pi * gc *
246  (-N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) *
247  w)
248  .eval();
249 
250  local_Jac.noalias() += local_Jac_COHESIVE;
251 
252  local_rhs.noalias() -= local_Jac_COHESIVE * d;
253  break;
254  }
255  }
256  }
257 }
258 
259 template <typename ShapeFunction, int DisplacementDim>
261  computeCrackIntegral(std::size_t mesh_item_id,
262  std::vector<std::reference_wrapper<
263  NumLib::LocalToGlobalIndexMap>> const& dof_tables,
264  GlobalVector const& /*x*/, double const /*t*/,
265  double& crack_volume,
266  CoupledSolutionsForStaggeredScheme const* const cpl_xs)
267 {
268  assert(cpl_xs != nullptr);
269 
270  std::vector<std::vector<GlobalIndexType>> indices_of_processes;
271  indices_of_processes.reserve(dof_tables.size());
272  std::transform(dof_tables.begin(), dof_tables.end(),
273  std::back_inserter(indices_of_processes),
274  [&](NumLib::LocalToGlobalIndexMap const& dof_table)
275  { return NumLib::getIndices(mesh_item_id, dof_table); });
276 
277  auto local_coupled_xs =
278  getCoupledLocalSolutions(cpl_xs->coupled_xs, indices_of_processes);
279  assert(local_coupled_xs.size() == displacement_size + phasefield_size);
280 
281  auto const d = Eigen::Map<PhaseFieldVector const>(
282  &local_coupled_xs[phasefield_index], phasefield_size);
283  auto const u = Eigen::Map<DeformationVector const>(
284  &local_coupled_xs[displacement_index], displacement_size);
285 
287  x_position.setElementID(_element.getID());
288 
289  int const n_integration_points = _integration_method.getNumberOfPoints();
290  for (int ip = 0; ip < n_integration_points; ip++)
291  {
292  x_position.setIntegrationPoint(ip);
293  auto const& w = _ip_data[ip].integration_weight;
294  auto const& N = _ip_data[ip].N;
295  auto const& dNdx = _ip_data[ip].dNdx;
296 
297  typename ShapeMatricesType::template MatrixType<DisplacementDim,
298  displacement_size>
299  N_u = ShapeMatricesType::template MatrixType<
300  DisplacementDim, displacement_size>::Zero(DisplacementDim,
301  displacement_size);
302 
303  for (int i = 0; i < DisplacementDim; ++i)
304  {
305  N_u.template block<1, displacement_size / DisplacementDim>(
306  i, i * displacement_size / DisplacementDim)
307  .noalias() = N;
308  }
309 
310  crack_volume += (N_u * u).dot(dNdx * d) * w;
311  }
312 }
313 
314 template <typename ShapeFunction, int DisplacementDim>
316  std::size_t mesh_item_id,
317  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
318  dof_tables,
319  GlobalVector const& /*x*/, double const t, double& elastic_energy,
320  double& surface_energy, double& pressure_work,
321  CoupledSolutionsForStaggeredScheme const* const cpl_xs)
322 {
323  assert(cpl_xs != nullptr);
324 
325  std::vector<std::vector<GlobalIndexType>> indices_of_processes;
326  indices_of_processes.reserve(dof_tables.size());
327  std::transform(dof_tables.begin(), dof_tables.end(),
328  std::back_inserter(indices_of_processes),
329  [&](NumLib::LocalToGlobalIndexMap const& dof_table)
330  { return NumLib::getIndices(mesh_item_id, dof_table); });
331 
332  auto const local_coupled_xs =
333  getCoupledLocalSolutions(cpl_xs->coupled_xs, indices_of_processes);
334  assert(local_coupled_xs.size() == displacement_size + phasefield_size);
335 
336  auto const d = Eigen::Map<PhaseFieldVector const>(
337  &local_coupled_xs[phasefield_index], phasefield_size);
338  auto const u = Eigen::Map<DeformationVector const>(
339  &local_coupled_xs[displacement_index], displacement_size);
340 
342  x_position.setElementID(_element.getID());
343 
344  double element_elastic_energy = 0.0;
345  double element_surface_energy = 0.0;
346  double element_pressure_work = 0.0;
347 
348  double const gc = _process_data.crack_resistance(t, x_position)[0];
349  double const ls = _process_data.crack_length_scale(t, x_position)[0];
350  int const n_integration_points = _integration_method.getNumberOfPoints();
351  for (int ip = 0; ip < n_integration_points; ip++)
352  {
353  x_position.setIntegrationPoint(ip);
354  auto const& w = _ip_data[ip].integration_weight;
355  auto const& N = _ip_data[ip].N;
356  auto const& dNdx = _ip_data[ip].dNdx;
357  auto pressure_ip = _process_data.pressure;
358 
359  typename ShapeMatricesType::template MatrixType<DisplacementDim,
360  displacement_size>
361  N_u = ShapeMatricesType::template MatrixType<
362  DisplacementDim, displacement_size>::Zero(DisplacementDim,
363  displacement_size);
364 
365  for (int i = 0; i < DisplacementDim; ++i)
366  {
367  N_u.template block<1, displacement_size / DisplacementDim>(
368  i, i * displacement_size / DisplacementDim)
369  .noalias() = N;
370  }
371 
372  element_elastic_energy += _ip_data[ip].elastic_energy * w;
373 
374  double const d_ip = N.dot(d);
375  switch (_process_data.phasefield_model)
376  {
378  {
379  element_surface_energy +=
380  gc * 0.375 *
381  ((1 - d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) * w;
382 
383  break;
384  }
386  {
387  element_surface_energy += 0.5 * gc *
388  ((1 - d_ip) * (1 - d_ip) / ls +
389  (dNdx * d).dot((dNdx * d)) * ls) *
390  w;
391  break;
392  }
394  {
395  element_surface_energy +=
396  gc / boost::math::double_constants::pi *
397  ((1 - d_ip * d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) *
398  w;
399  break;
400  }
401  }
402 
403  if (_process_data.crack_pressure)
404  {
405  element_pressure_work += pressure_ip * (N_u * u).dot(dNdx * d) * w;
406  }
407  }
408 
409 #ifdef USE_PETSC
410  int const n_all_nodes = indices_of_processes[1].size();
411  int const n_regular_nodes = std::count_if(
412  begin(indices_of_processes[1]), end(indices_of_processes[1]),
413  [](GlobalIndexType const& index) { return index >= 0; });
414  if (n_all_nodes != n_regular_nodes)
415  {
416  element_elastic_energy *=
417  static_cast<double>(n_regular_nodes) / n_all_nodes;
418  element_surface_energy *=
419  static_cast<double>(n_regular_nodes) / n_all_nodes;
420  element_pressure_work *=
421  static_cast<double>(n_regular_nodes) / n_all_nodes;
422  }
423 #endif // USE_PETSC
424  elastic_energy += element_elastic_energy;
425  surface_energy += element_surface_energy;
426  pressure_work += element_pressure_work;
427 }
428 } // namespace PhaseField
429 } // namespace ProcessLib
GlobalMatrix::IndexType GlobalIndexType
Global vector based on Eigen vector.
Definition: EigenVector.h:28
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition: BMatrixPolicy.h:55
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)
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
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
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
static const double u
static const double t
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
Definition: LinearBMatrix.h:42
std::vector< double > getCoupledLocalSolutions(std::vector< GlobalVector * > const &global_solutions, std::vector< std::vector< GlobalIndexType >> const &indices)
std::vector< GlobalVector * > const & coupled_xs
References to the current solutions of the coupled processes.