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