deal.II version GIT relicensing-2838-gd85d4b70e9 2025-03-13 22:40:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
The 'Phase-field fracture model in 3D' code gallery program

This program was contributed by Wasim Niyaz Munshi <munshiw0@gmail.com>.
It comes without any warranty or support by its authors or the authors of deal.II.

This program is part of the deal.II code gallery and consists of the following files (click to inspect):

Pictures from this code gallery program

Annotated version of README.md

Phase-field Fracture in 3D

Motivation

This program implements the hybrid phase-field model of Ambati et al. [1] in deal.II using adaptive mesh refinement and parallel computing capabilities. Phase-field models have been proven to be effective in modeling complicated fractures because of their ability to model crack branching, merging, and fragmentation. Despite this, these models are mostly limited to 2D because of the high computational cost, as these models require very fine meshes to resolve the diffused representation of the crack. This code explores the use of parallel computing and adaptive mesh refinement to model fracture propagation using the hybrid phase-field model.

Note
If you use this program as a basis for your own work, please consider citing it in your list of references. The initial version of this work was contributed to the deal.II project by the authors listed in the following citation: 10.5281/zenodo.14682543

Governing equations

The model this program solves is that of Ambati et al. [1], see there for more information. In short, the model describes the growth of a crack as an object is successively strained. The deformation of the solid is described by the usual force balance appropriate for quasi-static deformation:

\begin{align*} \nabla\cdot\boldsymbol{\sigma} (\boldsymbol{\varepsilon}(\boldsymbol{u}),d) = \mathbf{0} \end{align*}

with Dirichlet boundary conditions

\begin{align*} \boldsymbol{u} = \boldsymbol{u}_D \text{ on } \Gamma_D \end{align*}

The crack is tracked by a "damage field" \(d\) that is a smoothed version of a history field \(\mathcal{H}^{+}\) that corresponds to accumulated strain at a quadrature point:

\begin{align*} -l^2 \nabla^2 d+d=\frac{2 l}{G_c}(1-d) \mathcal{H}^{+} \end{align*}

with the boundary condition

\begin{align*} \left(G_{c} l\right) \nabla d \cdot \boldsymbol{n}=\mathbf{0}. \end{align*}

Here, \(\boldsymbol{u}, d\) represent the displacement and damage(crack) fields, \(l\) is the length scale parameter, \(G_c\) is the critical energy release rate and \(\mathcal{H}^{+}\) is the history field variable.

To run the Code

After running cmake ., run make release or make debug to switch between release and debugmode. Compile using make. Run the executable by using make run on the command line. Run the executable on 'n' processes using 'mpirun -np \(n\) ./phase_field'. For example 'mpirun -np 40 ./phase_field' runs the program on 40 processes.

Numerical example

The program currently models fracture propagation in a 3-layered material subjected to equibiaxial loading along \(+x\) and \(+y\) directions. The problem setup is shown in the following picture: Setup Results are compared for two different initial meshes, the finer \(80\times80\times40\) mesh and the coarser \(20\times20\times10\) mesh. The following picture shows results for the fracture patterns and the load-displacement curves for the 2 meshes: Load displacement curves

An animation of how the crack system in this setup evolves can be found here (left: coarse mesh; right: fine mesh).

References

[1]

@article{Ambati2015,
title={A review on phase-field models of brittle fracture and a new fast hybrid formulation},
author={Ambati, Marreddy and Gerasimov, Tymofiy and De Lorenzis, Laura},
journal={Computational Mechanics},
volume={55},
pages={383--405},
year={2015},
publisher={Springer},
doi={10.1007/s00466-014-1109-y}
}

Annotated version of phase_field.cc

  /* -----------------------------------------------------------------------------
  *
  * SPDX-License-Identifier: LGPL-2.1-or-later
  * Copyright (C) 2024 by Wasim Niyaz Munshi
  *
  * This file is part of the deal.II code gallery.
  *
  * -----------------------------------------------------------------------------
  */
  #include <deal.II/base/quadrature_lib.h>
  #include <deal.II/base/function.h>
  #include <deal.II/base/timer.h>
  #include <deal.II/lac/generic_linear_algebra.h>
 
  #define FORCE_USE_OF_TRILINOS
  namespace LA
  {
  #if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \
  !(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
  using namespace dealii::LinearAlgebraPETSc;
  # define USE_PETSC_LA
  #elif defined(DEAL_II_WITH_TRILINOS)
  using namespace dealii::LinearAlgebraTrilinos;
  #else
  # error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
  #endif
  } // namespace LA
 
  #include <deal.II/lac/vector.h>
  #include <deal.II/lac/full_matrix.h>
  #include <deal.II/lac/precondition.h>
  #include <deal.II/lac/solver_cg.h>
  #include <deal.II/lac/affine_constraints.h>
  #include <deal.II/lac/dynamic_sparsity_pattern.h>
  #include <deal.II/grid/grid_generator.h>
  #include <deal.II/dofs/dof_handler.h>
  #include <deal.II/dofs/dof_tools.h>
  #include <deal.II/dofs/dof_accessor.h>
  #include <deal.II/grid/tria_accessor.h>
  #include <deal.II/fe/fe_values.h>
  #include <deal.II/fe/fe_system.h>
  #include <deal.II/fe/fe_q.h>
  #include <deal.II/grid/grid_refinement.h>
  #include <deal.II/numerics/vector_tools.h>
  #include <deal.II/numerics/data_out.h>
  #include <deal.II/numerics/error_estimator.h>
  #include <deal.II/base/utilities.h>
  #include <deal.II/base/conditional_ostream.h>
  #include <deal.II/base/index_set.h>
  #include <deal.II/lac/sparsity_tools.h>
  #include <deal.II/distributed/tria.h>
  #include <deal.II/distributed/grid_refinement.h>
  #include <deal.II/distributed/solution_transfer.h>
  #include <deal.II/base/quadrature_point_data.h>
  #include <deal.II/base/tensor_function.h>
  #include <fstream>
  #include <iostream>
  #include <random>
 
  namespace FracturePropagation
 
  {
  using namespace dealii;
 
  class PhaseField
  {

Function declarations

  public:
  PhaseField ();
  void
  run ();
 
  private:
 
  1. Mesh and boundary conditions setup at the beginning
  void
  setup_mesh_and_bcs (); // Called at the start to create a mesh
 
  1. Elastic subproblem setup and solution
  void
  setup_constraints_elastic (const unsigned int load_step);
  void
  setup_system_elastic (const unsigned int load_step);
  void
  assemble_system_elastic ();
  void
  solve_linear_system_elastic ();
  void
  solve_elastic_subproblem (const unsigned int load_step); // Calls the above functions
 
  1. Damage subproblem setup and solution
  void
  setup_boundary_values_damage ();
  void
  setup_system_damage ();
  void
  assemble_system_damage ();
  void
  solve_linear_system_damage ();
  double
  H_plus (const SymmetricTensor<2, 3> &strain); // Called during system assembly for the damage subproblem
  void
  solve_damage_subproblem (); // Calls the above functions
 
  1. Convergence check after each iteration
  bool
  check_convergence ();
 
  1. Post-processing: output results, refine grid, and calculate displacement
  void
  output_results (const unsigned int load_step) const; // Called if convergence is achieved
  void
  load_disp_calculation (const unsigned int load_step); // Called if convergence is achieved
  void
  update_history_field (); // Called if convergence is achieved
  void
  refine_grid (const unsigned int load_step); // Called if convergence is not achieved
 
  MPI_Comm mpi_communicator;
  TimerOutput computing_timer;
 
 
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation

Objects for elasticity

  const FESystem<3> fe_elastic;
  DoFHandler<3> dof_handler_elastic;
  IndexSet locally_owned_dofs_elastic;
  IndexSet locally_relevant_dofs_elastic;
  AffineConstraints<double> constraints_elastic;
  LA::MPI::SparseMatrix system_matrix_elastic;
  LA::MPI::Vector locally_relevant_solution_elastic;
  LA::MPI::Vector completely_distributed_solution_elastic_old;
  LA::MPI::Vector completely_distributed_solution_elastic;
  LA::MPI::Vector system_rhs_elastic;
  const QGauss<3> quadrature_formula_elastic;
 

Objects for damage

  const FE_Q<3> fe_damage;
  DoFHandler<3> dof_handler_damage;
  IndexSet locally_owned_dofs_damage;
  IndexSet locally_relevant_dofs_damage;
  AffineConstraints<double> constraints_damage;
  LA::MPI::SparseMatrix system_matrix_damage;
  LA::MPI::Vector locally_relevant_solution_damage;
  LA::MPI::Vector completely_distributed_solution_damage_old;
  LA::MPI::Vector completely_distributed_solution_damage;
  LA::MPI::Vector system_rhs_damage;
  const QGauss<3> quadrature_formula_damage;
 
  const double ux = 1e-3; // increment in loading
  const double alpha = 1;
  const double uy = alpha * ux;
  const double l = 0.6; // length scale parameter
  const unsigned int num_load_steps = 100; // number of load steps
  const double tol = 1e-2; // tolerance for error in solution
  const double GC = 1e-3; //energy release rate
  const double E = 37.5;
  const double beta = 25;
  const double nu = 0.25;
 
Definition fe_q.h:554

Objects for load-displacement calculation

  Vector<double> load_values_x;
  Vector<double> load_values_y;
  Vector<double> displacement_values;
 
  class MyQData : public TransferableQuadraturePointData
  {
  public:
  MyQData () = default;
  virtual
  ~MyQData () = default;
 
  unsigned int
  number_of_values () const override
  {
  return 2; // Indicate that there are two values to handle
  }
 
  virtual void
  pack_values (std::vector<double> &scalars) const override
  {
  Assert (scalars.size () == 2, ExcInternalError ()); // Ensure the vector has exactly two elements
  scalars[0] = value_H; // Pack the first value
  scalars[1] = value_H_new; // Pack the second value
  }
 
  virtual void
  unpack_values (const std::vector<double> &scalars) override
  {
  Assert (scalars.size () == 2, ExcInternalError ()); // Ensure the vector has exactly two elements
  value_H = scalars[0]; // Unpack the first value
  value_H_new = scalars[1]; // Unpack the second value
  }
 
  double value_H; // First value
  double value_H_new; // Second value
  };
 
  CellDataStorage<typename Triangulation<3>::cell_iterator, MyQData> quadrature_point_history_field;
 
  };
 
  inline double
  lambda (const float E,
  const float nu)
  {
  return (E * nu) / ((1 + nu) * (1 - 2 * nu));
  }
 
  inline double
  mu (const float E,
  const float nu)
  {
  return E / (2 * (1 + nu));
  }
 
  inline double
  gc (const float GC,
  const float beta,
  const float z1,
  const float z2,
  const Point<3> &p)
  {
  if (((p[2] - z1) > 1e-6) && ((p[2] - z2) < 1e-6))
  return GC / beta;
  else
  return GC;
  }
 
  bool
  is_in_middle_layer (const Point<3> &cell_center,const float z1,const float z2)
  {
  return (cell_center[2] >= z1 && cell_center[2] <= z2);
  }
 
  namespace Domain
 
  {
  const double x_min = 0;
  const double y_min = 0;
  const double z_min = 0;
  const double x_max = 30;
  const double y_max = 30;
  const double z_max = 13;
  const double z1 = 5;
  const double z2 = 8;
  }
 
  namespace RandomMedium
  {
  class EnergyReleaseRate : public Function<3>
  {
  public:
  EnergyReleaseRate ()
  :
  Function<3> ()
  {
  }
 
Definition point.h:113
virtual unsigned int number_of_values() const =0
virtual void unpack_values(const std::vector< double > &values)=0
virtual void pack_values(std::vector< double > &values) const =0
#define Assert(cond, exc)
constexpr double E
Definition numbers.h:214

value_list method calculates the fracture_energy (Gc) at a set of points and stores the result in a vector of zero order tensors.

  virtual void
  value_list (const std::vector<Point<3>> &points,
  std::vector<double> &values,
  const unsigned int /*component*/= 0) const override
  {
  AssertDimension (points.size (), values.size ());
 
  for (unsigned int p = 0; p < points.size (); ++p)
  {
  values[p] = 0.0;
 
  double fracture_energy = 0;
  for (unsigned int i = 0; i < centers.size (); ++i)
  fracture_energy += std::exp (
  -(points[p] - centers[i]).norm_square () / (1.5 * 1.5));
 
  const double normalized_fracture_energy = std::min (
  std::max (fracture_energy, 4e-5), 4e-4);
 
  values[p] = normalized_fracture_energy;
  }
  }
 
  private:
  static std::vector<Point<3>> centers;
 
  static std::vector<Point<3>>
  get_centers ()
  {
  const unsigned int N = 1000;
 
  std::vector<Point<3>> centers_list (N);
  for (unsigned int i = 0; i < N; ++i)
  for (unsigned int d = 0; d < 3; ++d)
  if (d == 0 || d == 1)
  {
  centers_list[i][d] =
  static_cast<double> ((rand ()) / RAND_MAX) * Domain::x_max; //domain_size; //generates a number between 0 and domain_size
#define AssertDimension(dim1, dim2)
constexpr char N
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)

x,y will be between 0 to x_max

  }
  else if (d == 2)
  {
  centers_list[i][d] = static_cast<double> (Domain::z1
  + ((rand ()) / RAND_MAX) * (Domain::z2 - Domain::z1));
  }
 
  return centers_list;
  }
  };
 
  std::vector<Point<3>> EnergyReleaseRate::centers =
  EnergyReleaseRate::get_centers ();
 
  } // namespace RandomMedium
 
 
  inline double
  Conductivity_damage (const Point<3> &p) //const
  {
  return p[0] - p[0] + 1;
  }
 
  void
  right_hand_side_elastic (const std::vector<Point<3>> &points,
  std::vector<Tensor<1, 3>> &values)
  {
  AssertDimension (values.size (), points.size ());
 
  for (unsigned int point_n = 0; point_n < points.size (); ++point_n)
  {
  values[point_n][0] = 0; //x component of body force
  values[point_n][1] = 0; //y component of body force
  values[point_n][2] = 0; //z component of body force
  }
 
  }
  void
  Traction_elastic (const std::vector<Point<3>> &points,
  std::vector<Tensor<1, 3>> &values)
 
  {
  AssertDimension (values.size (), points.size ());
  for (unsigned int point_n = 0; point_n < points.size (); ++point_n)
  {
  values[point_n][0] = 0; //x component of traction
  values[point_n][1] = 0; //y component of traction
  values[point_n][2] = 0; //y component of traction
  }
 
  }
 
  PhaseField::PhaseField ()
  :
  mpi_communicator (MPI_COMM_WORLD),
  pcout (std::cout,
  (Utilities::MPI::this_mpi_process (mpi_communicator) == 0)),
  computing_timer (mpi_communicator, pcout, TimerOutput::never,
  TimerOutput::wall_times),
  triangulation (mpi_communicator),
  fe_elastic (FE_Q<3> (1), 3),
  dof_handler_elastic (triangulation),
  quadrature_formula_elastic (fe_elastic.degree + 1),
  fe_damage (1),
  dof_handler_damage (triangulation),
  quadrature_formula_damage (fe_elastic.degree + 1),
  load_values_x(num_load_steps+1),
  load_values_y(num_load_steps+1),
  displacement_values(num_load_steps+1)
  {
  }
 
  double
  PhaseField::H_plus (const SymmetricTensor<2, 3> &strain)
  {
  double Mac_tr_strain, Mac_first_principal_strain,
  Mac_second_principal_strain, Mac_third_principal_strain,
  tr_sqr_Mac_Principal_strain;
 
  const double tr_strain = trace (strain);
 
  Mac_tr_strain = tr_strain >0 ? tr_strain : 0;
  const std::array<double, 3> Principal_strains = eigenvalues (strain);
 
  Mac_first_principal_strain = (Principal_strains[0] > 0) ? Principal_strains[0] : 0;
  Mac_second_principal_strain = (Principal_strains[1] > 0) ? Principal_strains[1] : 0;
  Mac_third_principal_strain = (Principal_strains[2] > 0) ? Principal_strains[2] : 0;
 
  tr_sqr_Mac_Principal_strain = pow (Mac_first_principal_strain, 2)
  + pow (Mac_second_principal_strain, 2)
  + pow (Mac_third_principal_strain, 2);
 
  double H_plus_val;
  H_plus_val = 0.5 * lambda (E, nu) * pow (Mac_tr_strain, 2)
  + mu (E, nu) * tr_sqr_Mac_Principal_strain;
  return H_plus_val;
  }
 
 
  void
  PhaseField::setup_mesh_and_bcs ()
 
  {
  const unsigned int nx = 20;
  const unsigned int ny = 20;
  const unsigned int nz = 10;
  const std::vector<unsigned int> repetitions = {nx,ny,nz};
 
  const Point<3> p1(Domain::x_min,Domain::y_min,Domain::z_min), p2(Domain::x_max,Domain::y_max,Domain::z_max);
 
  p2); // create coarse mesh
 
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
STL namespace.
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)

The boundary ids need to be setup right after the mesh is generated (before any refinement) and ids need to be setup using all cells and not just the locally owned cells

  for (const auto &cell : triangulation.active_cell_iterators ())
  {
  for (const auto &face : cell->face_iterators ())
  {
  if (face->at_boundary ())
  {
  const auto center = face->center ();
  if (std::fabs (center (0) - (Domain::x_min)) < 1e-12)
  face->set_boundary_id (0);
 
  else if (std::fabs (center (0) - Domain::x_max) < 1e-12)
  face->set_boundary_id (1);
 
  else if (std::fabs (center (1) - (Domain::y_min)) < 1e-12)
  face->set_boundary_id (2);
 
  else if (std::fabs (center (1) - Domain::y_max) < 1e-12)
  face->set_boundary_id (3);
  }
  }
  }
 
  pcout << "No. of levels in triangulation: "
  << triangulation.n_global_levels () << std::endl;
 
  dof_handler_damage.distribute_dofs (fe_damage);
  dof_handler_elastic.distribute_dofs (fe_elastic);
 
  pcout << " Number of locally owned cells on the process: "
  << triangulation.n_locally_owned_active_cells () << std::endl;
 
  pcout << "Number of global cells:" << triangulation.n_global_active_cells ()
  << std::endl;
 
  pcout << " Total Number of globally active cells: "
  << triangulation.n_global_active_cells () << std::endl
  << " Number of degrees of freedom for elasticity: "
  << dof_handler_elastic.n_dofs () << std::endl
  << " Number of degrees of freedom for damage: "
  << dof_handler_damage.n_dofs () << std::endl;
 

Initialising damage vectors

  locally_owned_dofs_damage = dof_handler_damage.locally_owned_dofs ();
  locally_relevant_dofs_damage = DoFTools::extract_locally_relevant_dofs (
  dof_handler_damage);
 
  completely_distributed_solution_damage_old.reinit (
  locally_owned_dofs_damage, mpi_communicator);
  locally_relevant_solution_damage.reinit (locally_owned_dofs_damage,
  locally_relevant_dofs_damage, mpi_communicator);
 
  locally_owned_dofs_elastic = dof_handler_elastic.locally_owned_dofs ();
  locally_relevant_dofs_elastic = DoFTools::extract_locally_relevant_dofs (
  dof_handler_elastic);
 
  completely_distributed_solution_elastic_old.reinit (
  locally_owned_dofs_elastic, mpi_communicator);
 
  for (const auto &cell : triangulation.active_cell_iterators ())
  {
  if (cell->is_locally_owned ())
  {
  quadrature_point_history_field.initialize (cell, 8);
  }
  }
 
  FEValues<3> fe_values_damage (fe_damage, quadrature_formula_damage,
 
  for (const auto &cell : triangulation.active_cell_iterators ())
  if (cell->is_locally_owned ())
  {
  const std::vector<std::shared_ptr<MyQData>> lqph =
  quadrature_point_history_field.get_data (cell);
  for (const unsigned int q_index : fe_values_damage.quadrature_point_indices ())
  {
  lqph[q_index]->value_H = 0.0;
  lqph[q_index]->value_H_new = 0.0;
  }
  }
  }
 
  void
  PhaseField::setup_constraints_elastic (const unsigned int load_step)
  {
  constraints_elastic.clear ();
  constraints_elastic.reinit (locally_relevant_dofs_elastic);
 
  constraints_elastic);
 
  for (const auto &cell : dof_handler_elastic.active_cell_iterators ())
  if (cell->is_locally_owned ())
  {
  for (const auto &face : cell->face_iterators ())
  {
  if (face->at_boundary ())
  {
  const auto center = face->center ();
  if (std::fabs (center (0) - Domain::x_min) < 1e-12) //face lies at x=x_min
  {
 
  for (const auto vertex_number : cell->vertex_indices ())
  {
  const auto vert = cell->vertex (vertex_number);
  const double z_mid = 0.5 * (Domain::z_max + Domain::z_min);
  if (std::fabs (vert (2) - z_mid) < 1e-12 && std::fabs (
  vert (
  0)
  - Domain::x_min)
  < 1e-12) // vertex at x=x_min,z=z_mid;
  {
  const unsigned int z_dof =
  cell->vertex_dof_index (vertex_number, 2);
  constraints_elastic.add_line (z_dof);
  constraints_elastic.set_inhomogeneity (z_dof, 0);
  const unsigned int x_dof =
  cell->vertex_dof_index (vertex_number, 0);
  constraints_elastic.add_line (x_dof);
  constraints_elastic.set_inhomogeneity (x_dof, 0);
  }
  else if (std::fabs (vert (0) - Domain::x_min) < 1e-12) // vertex at x_min
 
  {
  const unsigned int x_dof =
  cell->vertex_dof_index (vertex_number, 0);
  constraints_elastic.add_line (x_dof);
  constraints_elastic.set_inhomogeneity (x_dof, 0);
  }
  }
  }
  }
  }
  }
 
  const FEValuesExtractors::Scalar u_x (0);
  const FEValuesExtractors::Scalar u_y (1);
 
  const ComponentMask u_x_mask = fe_elastic.component_mask (u_x);
  const ComponentMask u_y_mask = fe_elastic.component_mask (u_y);
 
  const double u_x_values_right = ux * load_step;
  const double u_y_values = uy * load_step;
  const double u_fix = 0.0;
 
  VectorTools::interpolate_boundary_values (dof_handler_elastic, 1,
  Functions::ConstantFunction<3> (u_x_values_right, 3),
  constraints_elastic, u_x_mask);
 
  VectorTools::interpolate_boundary_values (dof_handler_elastic, 2,
  Functions::ConstantFunction<3> (u_fix, 3), constraints_elastic,
  u_y_mask);
 
  VectorTools::interpolate_boundary_values (dof_handler_elastic, 3,
  Functions::ConstantFunction<3> (u_y_values, 3), constraints_elastic,
  u_y_mask);
 
  constraints_elastic.close ();
  }
 
  void
  PhaseField::setup_system_elastic (const unsigned int load_step)
  {
  TimerOutput::Scope ts (computing_timer, "setup_system_elastic");
 
  locally_owned_dofs_elastic = dof_handler_elastic.locally_owned_dofs ();
  locally_relevant_dofs_elastic = DoFTools::extract_locally_relevant_dofs (
  dof_handler_elastic);
 
  locally_relevant_solution_elastic.reinit (locally_owned_dofs_elastic,
  locally_relevant_dofs_elastic, mpi_communicator);
 
  system_rhs_elastic.reinit (locally_owned_dofs_elastic, mpi_communicator);
 
  completely_distributed_solution_elastic.reinit (locally_owned_dofs_elastic,
  mpi_communicator);
 
  setup_constraints_elastic (load_step);
 
  DynamicSparsityPattern dsp (locally_relevant_dofs_elastic);
 
  DoFTools::make_sparsity_pattern (dof_handler_elastic, dsp,
  constraints_elastic, false);
 
  dof_handler_elastic.locally_owned_dofs (), mpi_communicator,
  locally_relevant_dofs_elastic);
 
  system_matrix_elastic.reinit (locally_owned_dofs_elastic,
  locally_owned_dofs_elastic, dsp, mpi_communicator);
  }
 
  void
  PhaseField::assemble_system_elastic ()
 
  {
  TimerOutput::Scope ts (computing_timer, "assembly_elastic");
 
  FEValues<3> fe_values_elastic (fe_elastic, quadrature_formula_elastic,
  FEValues<3> fe_values_damage (fe_damage, quadrature_formula_elastic,
 
  const unsigned int dofs_per_cell = fe_elastic.n_dofs_per_cell ();
  const unsigned int n_q_points = quadrature_formula_elastic.size ();
 
  FullMatrix<double> cell_matrix_elastic (dofs_per_cell, dofs_per_cell);
  Vector<double> cell_rhs_elastic (dofs_per_cell);
 
  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
 
  std::vector<double> damage_values (n_q_points);
  std::vector<Tensor<1, 3>> rhs_values_elastic (n_q_points);
 
  for (const auto &cell : dof_handler_elastic.active_cell_iterators ())
  if (cell->is_locally_owned ())
 
  {
  cell_matrix_elastic = 0;
  cell_rhs_elastic = 0;
 
  Triangulation<3>::active_cell_iterator (cell)->as_dof_handler_iterator (
  dof_handler_damage);
 
  fe_values_damage.reinit (damage_cell);
  fe_values_elastic.reinit (cell);
 
  fe_values_damage.get_function_values(locally_relevant_solution_damage,
  damage_values);
  right_hand_side_elastic (fe_values_elastic.get_quadrature_points (),
  rhs_values_elastic);
 
  for (const unsigned int q_point : fe_values_elastic.quadrature_point_indices ())
  {
  const double d = damage_values[q_point];
 
  for (const unsigned int i : fe_values_elastic.dof_indices ())
 
  {
  const unsigned int component_i =
  fe_elastic.system_to_component_index (i).first;
 
  for (const unsigned int j : fe_values_elastic.dof_indices ())
  {
  const unsigned int component_j =
  fe_elastic.system_to_component_index (j).first;
 
  {
  cell_matrix_elastic (i, j) +=
  pow ((1 - d), 2) * (
  (fe_values_elastic.shape_grad (i, q_point)[component_i] *
  fe_values_elastic.shape_grad (j, q_point)[component_j]
  *
  lambda (E, nu))
  +
  (fe_values_elastic.shape_grad (i, q_point)[component_j] *
  fe_values_elastic.shape_grad (j, q_point)[component_i]
  *
  mu (E, nu))
  +
  ((component_i == component_j) ?
  (fe_values_elastic.shape_grad (i, q_point) *
  fe_values_elastic.shape_grad (j, q_point)
  *
  mu (E, nu)) :
  0)
  )
  *
  fe_values_elastic.JxW (q_point);
  }
  }
  }
  }
 
  for (const unsigned int i : fe_values_elastic.dof_indices ())
  {
  const unsigned int component_i =
  fe_elastic.system_to_component_index (i).first;
 
  for (const unsigned int q_point : fe_values_elastic.quadrature_point_indices ())
  cell_rhs_elastic (i) +=
  fe_values_elastic.shape_value (i, q_point) * rhs_values_elastic[q_point][component_i]
  * fe_values_elastic.JxW (q_point);
  }
 
std::vector< bool > component_mask
Point< 2 > first
Definition grid_out.cc:4632
unsigned int vertex_indices[2]
TriaActiveIterator< CellAccessor< dim, spacedim > > active_cell_iterator
Definition tria.h:1581
typename ActiveSelector::active_cell_iterator active_cell_iterator
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm mpi_comm, const IndexSet &locally_relevant_rows)
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask={})

traction contribution to rhs

  /* for (const auto &face : cell->face_iterators())
  {
  if (face->at_boundary() &&
  face->boundary_id() == 10)
 
  {
  fe_face_values_elastic.reinit(cell, face);
  for (const auto i : fe_face_values_elastic.dof_indices())
 
  {const unsigned int component_i =
  fe_elastic.system_to_component_index(i).first;
  for (const auto face_q_point :
  fe_face_values_elastic.quadrature_point_indices())
  {
 
  cell_rhs_elastic(i) +=
  fe_face_values_elastic.shape_value( i, face_q_point)
  *(traction_values_elastic[face_q_point][component_i])
  *fe_face_values_elastic.JxW(face_q_point);
 
 
 
  }
  }
  }
  }*/
 
  cell->get_dof_indices (local_dof_indices);
 
  constraints_elastic.distribute_local_to_global (cell_matrix_elastic,
  cell_rhs_elastic, local_dof_indices, system_matrix_elastic,
  system_rhs_elastic);
  }
 
  system_matrix_elastic.compress (VectorOperation::add);
  system_rhs_elastic.compress (VectorOperation::add);
  }
 
  void
  PhaseField::solve_linear_system_elastic ()
  {
  TimerOutput::Scope ts (computing_timer, "solve_linear_system_elastic");
 
  SolverControl solver_control (10000,
  1e-12* system_rhs_elastic.l2_norm());
 
  SolverCG<LA::MPI::Vector> solver (solver_control);
 
  LA::MPI::PreconditionAMG::AdditionalData data;
  #ifdef USE_PETSC_LA
  data.symmetric_operator = true;
  #else
std::vector< index_type > data
Definition mpi.cc:746

Trilinos defaults are good

  #endif
  LA::MPI::PreconditionAMG preconditioner;
  preconditioner.initialize (system_matrix_elastic, data);
 
  solver.solve (system_matrix_elastic,
  completely_distributed_solution_elastic, system_rhs_elastic,
  preconditioner);
 
  pcout << " Solved in " << solver_control.last_step () << " iterations."
  << std::endl;
 
  constraints_elastic.distribute (completely_distributed_solution_elastic);
 
  locally_relevant_solution_elastic = completely_distributed_solution_elastic;
  }
 
  void
  PhaseField::setup_boundary_values_damage ()
  {
  TimerOutput::Scope ts (computing_timer, "setup_bv_damage");
 
  constraints_damage.clear ();
  constraints_damage.reinit (locally_relevant_dofs_damage);
  constraints_damage);
 

Create initial crack, if any

  /*for (const auto &cell : dof_handler_damage.active_cell_iterators())
  if (cell->is_locally_owned())
 
  {

for (const auto &face : cell->face_iterators())

  for (const auto vertex_number : cell->vertex_indices())
  {
  const auto vert = cell->vertex(vertex_number);
  const Point<3>& node = vert;
  if (((std::fabs((Domain::z_max*(node(0)-0.5*Domian::x_max + A)) - 2*A*( node(2))) <10*l)
  && node(1)>=0.9*Domain::y_max)
  || ((std::fabs((node(0)-0.5*Domian::x_max + A)*Domain::z_max+2*A*(node(2)-Domain::z_max))<10*l)
  && node(1)<=0.1*Domain::y_max))
  if ((vert(0) - 0.5*(Domain::x_min+Domian::x_max) < 1e-12) &&
  (std::fabs(vert(1) - 0.5*(Domain::y_min + Domain::y_max)) <= bandwidth)) // nodes on initial damage plane
 
 
  {
  const unsigned int dof = cell->vertex_dof_index(vertex_number, 0);
  constraints_damage.add_line(dof);
  constraints_damage.set_inhomogeneity(dof,1);
  }
 
  }
 
  }*/
  constraints_damage.close ();
 
  }
 
  void
  PhaseField::setup_system_damage ()
  {
  TimerOutput::Scope ts (computing_timer, "setup_system_damage");
 
  locally_owned_dofs_damage = dof_handler_damage.locally_owned_dofs ();
  locally_relevant_dofs_damage = DoFTools::extract_locally_relevant_dofs (
  dof_handler_damage);
 
  locally_relevant_solution_damage.reinit (locally_owned_dofs_damage,
  locally_relevant_dofs_damage, mpi_communicator);
 
  system_rhs_damage.reinit (locally_owned_dofs_damage, mpi_communicator);
 
  completely_distributed_solution_damage.reinit (locally_owned_dofs_damage,
  mpi_communicator);
 
  DynamicSparsityPattern dsp (locally_relevant_dofs_damage);
 
  DoFTools::make_sparsity_pattern (dof_handler_damage, dsp,
  constraints_damage, false);
  dof_handler_damage.locally_owned_dofs (), mpi_communicator,
  locally_relevant_dofs_damage);
 
  system_matrix_damage.reinit (locally_owned_dofs_damage,
  locally_owned_dofs_damage, dsp, mpi_communicator);
  }
 
  void
  PhaseField::assemble_system_damage ()
  {
 
  TimerOutput::Scope ts (computing_timer, "assembly_damage");
 
  FEValues<3> fe_values_damage (fe_damage, quadrature_formula_damage,
  FEValues<3> fe_values_elastic (fe_elastic, quadrature_formula_damage,
 
  const unsigned int dofs_per_cell = fe_damage.n_dofs_per_cell ();
  const unsigned int n_q_points = quadrature_formula_damage.size ();
 
  FullMatrix<double> cell_matrix_damage (dofs_per_cell, dofs_per_cell);
  Vector<double> cell_rhs_damage (dofs_per_cell);
 
  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
 
  const RandomMedium::EnergyReleaseRate energy_release_rate;
  std::vector<double> energy_release_rate_values (n_q_points);
 

Storing strain tensor for all Gauss points of a cell in vector strain_values.

  std::vector<SymmetricTensor<2, 3>> strain_values (n_q_points);
 
  for (const auto &cell : dof_handler_damage.active_cell_iterators ())
  if (cell->is_locally_owned ())
  {
  std::vector<std::shared_ptr<MyQData>> qpdH =
  quadrature_point_history_field.get_data (cell);
 
  const DoFHandler<3>::active_cell_iterator elastic_cell =
  Triangulation<3>::active_cell_iterator (cell)->as_dof_handler_iterator (
  dof_handler_elastic);
 
  fe_values_damage.reinit (cell);
  fe_values_elastic.reinit (elastic_cell);
 
  const FEValuesExtractors::Vector displacements (
  /* first vector component = */0);
  fe_values_elastic[displacements].get_function_symmetric_gradients (
  locally_relevant_solution_elastic, strain_values);
 
  cell_matrix_damage = 0;
  cell_rhs_damage = 0;
 
  energy_release_rate.value_list (
  fe_values_damage.get_quadrature_points (),
  energy_release_rate_values);
 
  for (const unsigned int q_index : fe_values_damage.quadrature_point_indices ())
  {
  const auto &x_q = fe_values_damage.quadrature_point (q_index);
  double H_call = H_plus (strain_values[q_index]);
 
  const double H = std::max(H_call, qpdH[q_index]->value_H);
  qpdH[q_index]->value_H_new = H;
  Point<3> cell_center = cell->center ();
 
  double g_c;
  if (is_in_middle_layer (cell_center,Domain::z1,Domain::z2))
  g_c = energy_release_rate_values[q_index];
  else
  g_c = gc (GC, beta, Domain::z1, Domain::z2, x_q);
 
  for (const unsigned int i : fe_values_damage.dof_indices ())
  {
  for (const unsigned int j : fe_values_damage.dof_indices ())
  {
 
  cell_matrix_damage (i, j) +=

contribution to stiffness from -laplace u term

  Conductivity_damage (x_q) * fe_values_damage.shape_grad (
  i, q_index)
  * // kappa*grad phi_i(x_q)
  fe_values_damage.shape_grad (j, q_index)
  * // grad phi_j(x_q)
  fe_values_damage.JxW (q_index) // dx
  +

Contribution to stiffness from u term

  ((1 + (2 * l * H) / g_c) * (1 / pow (l, 2))
  * fe_values_damage.shape_value (i, q_index) * // phi_i(x_q)
  fe_values_damage.shape_value (j, q_index) * // phi_j(x_q)
  fe_values_damage.JxW (q_index)); // dx
  }

const auto &x_q = fe_values_damage.quadrature_point(q_index);

  cell_rhs_damage (i) += (fe_values_damage.shape_value (i,
  q_index)
  * // phi_i(x_q)
  (2 / (l * g_c))
  * H * fe_values_damage.JxW (q_index)); // dx
  }
  }
 
  cell->get_dof_indices (local_dof_indices);
  constraints_damage.distribute_local_to_global (cell_matrix_damage,
  cell_rhs_damage, local_dof_indices, system_matrix_damage,
  system_rhs_damage);
  }
 
  system_matrix_damage.compress (VectorOperation::add);
  system_rhs_damage.compress (VectorOperation::add);
  }
 
  void
  PhaseField::solve_linear_system_damage ()
  {
  TimerOutput::Scope ts (computing_timer, "solve_linear_system_damage");
 
  SolverControl solver_control (10000,
  1e-12* system_rhs_damage.l2_norm());
 
  SolverCG<LA::MPI::Vector> solver (solver_control);
 
  LA::MPI::PreconditionAMG::AdditionalData data;
  #ifdef USE_PETSC_LA
  data.symmetric_operator = true;
  #else

Trilinos defaults are good

  #endif
  LA::MPI::PreconditionAMG preconditioner;
  preconditioner.initialize (system_matrix_damage, data);
 
  solver.solve (system_matrix_damage, completely_distributed_solution_damage,
  system_rhs_damage, preconditioner);
 
  pcout << " Solved in " << solver_control.last_step () << " iterations."
  << std::endl;
 
  constraints_damage.distribute (completely_distributed_solution_damage);
  locally_relevant_solution_damage = completely_distributed_solution_damage;
  }
 
  void
  PhaseField::solve_elastic_subproblem (const unsigned int load_step)
  {
  setup_system_elastic (load_step);
  assemble_system_elastic ();
  solve_linear_system_elastic ();
  }
 
  void
  PhaseField::output_results (const unsigned int load_step) const
 
  {
  std::vector<std::string> displacement_names (3, "displacement");
  std::vector<DataComponentInterpretation::DataComponentInterpretation> displacement_component_interpretation (
 
  DataOut<3> data_out_phasefield;
  data_out_phasefield.add_data_vector (dof_handler_elastic,
  locally_relevant_solution_elastic, displacement_names,
  displacement_component_interpretation);
  data_out_phasefield.add_data_vector (dof_handler_damage,
  locally_relevant_solution_damage, "damage");
 
  Vector<double> subdomain (triangulation.n_active_cells ());
  for (unsigned int i = 0; i < subdomain.size (); ++i)
  subdomain (i) = triangulation.locally_owned_subdomain ();
  data_out_phasefield.add_data_vector (subdomain, "subdomain");
  data_out_phasefield.build_patches ();
  data_out_phasefield.write_vtu_with_pvtu_record ("./", "solution", load_step,
  mpi_communicator, 2, 0);
  }
 
  void
  PhaseField::load_disp_calculation (const unsigned int load_step)
 
  {
  Tensor<1, 3> x_max_force; //force vector on the x_max face
  Tensor<1, 3> y_max_force; //force vector on the y_max face
 
  const QGauss<2> face_quadrature (fe_elastic.degree);
  FEFaceValues<3> fe_face_values (fe_elastic, face_quadrature,
  std::vector<SymmetricTensor<2, 3>> strain_values (face_quadrature.size ());
  const FEValuesExtractors::Vector displacements (0);
 
  FEFaceValues<3> fe_face_values_damage (fe_damage, face_quadrature, update_values);
  std::vector<double> damage_values (fe_face_values_damage.n_quadrature_points);
 
  for (const auto &cell : dof_handler_elastic.active_cell_iterators ())
  if (cell->is_locally_owned ())
  {
  const DoFHandler<3>::active_cell_iterator damage_cell =
  Triangulation<3>::active_cell_iterator (cell)->as_dof_handler_iterator (
  dof_handler_damage);
  for (unsigned int f : cell->face_indices ())
  if (cell->face (f)->at_boundary () && (cell->face (f)->boundary_id ()
  == 1))
  {
  fe_face_values.reinit (cell, f);
  fe_face_values[displacements].get_function_symmetric_gradients (
  locally_relevant_solution_elastic, strain_values);
  fe_face_values_damage.reinit(damage_cell, f);
 
  fe_face_values_damage.get_function_values (locally_relevant_solution_damage, damage_values);
 
  for (unsigned int q = 0; q < fe_face_values.n_quadrature_points; ++q)
 
  {
  const Tensor<2, 3> strain = strain_values[q]; //strain tensor at a gauss point
  const double tr_strain = strain[0][0] + strain[1][1] + strain[2][2];
  const double d = damage_values[q];
 
  Tensor<2, 3> stress;
  stress[0][0] = pow ((1 - d), 2)
  * (lambda (E, nu) * tr_strain + 2 * mu (E, nu)
  * strain[0][0]);
  stress[0][1] = pow ((1 - d), 2)
  * (2 * mu (E, nu) * strain[0][1]);
  stress[0][2] = pow ((1 - d), 2)
  * (2 * mu (E, nu) * strain[0][2]);
  stress[1][1] = pow ((1 - d), 2)
  * (lambda (E, nu) * tr_strain + 2 * mu (E, nu)
  * strain[1][1]);
  stress[1][2] = pow ((1 - d), 2)
  * (2 * mu (E, nu) * strain[1][2]);
  stress[2][2] = pow ((1 - d), 2)
  * (lambda (E, nu) * tr_strain + 2 * mu (E, nu)
  * strain[2][2]);
 
  const Tensor<1, 3> force_density = stress
  * fe_face_values.normal_vector (q);
  x_max_force += force_density * fe_face_values.JxW (q);
  }
  }
 
  else if (cell->face (f)->at_boundary ()
  && (cell->face (f)->boundary_id () == 3))
  {
  fe_face_values.reinit (cell, f);
  fe_face_values_damage.reinit (damage_cell, f); //Convert f?
  fe_face_values[displacements].get_function_symmetric_gradients (
  locally_relevant_solution_elastic, strain_values);
  fe_face_values_damage.get_function_values (locally_relevant_solution_damage, damage_values);
 
  for (unsigned int q = 0; q < fe_face_values.n_quadrature_points;
  ++q)
 
  {
  const Tensor<2, 3> strain = strain_values[q]; //strain tensor at a gauss point
  const double tr_strain = strain[0][0] + strain[1][1] + strain[2][2];
  const double d = damage_values[q];
 
  Tensor<2, 3> stress;
  stress[0][0] = pow ((1 - d), 2)
  * (lambda (E, nu) * tr_strain + 2 * mu (E, nu)
  * strain[0][0]);
  stress[0][1] = pow ((1 - d), 2)
  * (2 * mu (E, nu) * strain[0][1]);
  stress[0][2] = pow ((1 - d), 2)
  * (2 * mu (E, nu) * strain[0][2]);
  stress[1][1] = pow ((1 - d), 2)
  * (lambda (E, nu) * tr_strain + 2 * mu (E, nu)
  * strain[1][1]);
  stress[1][2] = pow ((1 - d), 2)
  * (2 * mu (E, nu) * strain[1][2]);
  stress[2][2] = pow ((1 - d), 2)
  * (lambda (E, nu) * tr_strain + 2 * mu (E, nu)
  * strain[2][2]);
 
  const Tensor<1, 3> force_density = stress
  * fe_face_values.normal_vector (q);
  y_max_force += force_density * fe_face_values.JxW (q);
  }
  }
  }
  double x_max_force_x;
  x_max_force_x = x_max_force[0];
  x_max_force_x = Utilities::MPI::sum (x_max_force_x, mpi_communicator);
  pcout << "fx: " << x_max_force_x << std::endl;
 
  double y_max_force_y;
  y_max_force_y = y_max_force[1];
  y_max_force_y = Utilities::MPI::sum (y_max_force_y, mpi_communicator);
  pcout << "fy: " << y_max_force_y << std::endl;
 
  if (Utilities::MPI::this_mpi_process (mpi_communicator) == 0)
  {
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
@ update_normal_vectors
Normal vectors.
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
unsigned int boundary_id
Definition types.h:153

Code to be executed only on process 0

  load_values_x[load_step] = x_max_force_x;
  load_values_y[load_step] = y_max_force_y;
  {
  const double disp = ux * load_step;
  displacement_values[load_step] = disp;
  }

load displacement plot

  std::ofstream m_file;
  m_file.open ("load_displacement.m");
  m_file << "% Matlab code generated by dealii to plot load displacement curves"
  << std::endl;
  m_file << "clc" << std::endl;
  m_file << "clear" << std::endl;
  m_file << "close all" << std::endl;
  m_file << "load_x=[" << load_values_x << "]" << std::endl;
  m_file << "load_y=[" << load_values_y << "]" << std::endl;
  m_file << "displacement=[" << displacement_values << "]" << std::endl;
  m_file << "plot( displacement,load_x,'linewidth',2)" << std::endl;
  m_file << "xlabel(\"Displacement (mm)\",'Interpreter', 'latex')"
  << std::endl;
  m_file << "ylabel(\"Reaction force (kN)\",'Interpreter', 'latex')"
  << std::endl;
  m_file << "hold on" << std::endl;
  m_file << "plot( displacement,load_y,'linewidth',2)" << std::endl;
  m_file << "set(gca,'fontname', 'Courier','FontSize',15,'FontWeight','bold','linewidth',1); grid on"
  << std::endl;
  m_file << "h=legend(\"fx \",\"fy\")" << std::endl;
  m_file << "set(h,'FontSize',14);" << std::endl;
  m_file << "set(gcf,'Position', [250 250 700 500])" << std::endl;
  m_file << "xlim([0,0.065])" << std::endl;
  m_file << "box on" << std::endl;
  m_file.close ();
  }
  }
 
  void
  PhaseField::solve_damage_subproblem ()
  {
  setup_boundary_values_damage ();
  setup_system_damage ();
  assemble_system_damage ();
  solve_linear_system_damage ();
  }
 
  void
  PhaseField::refine_grid (const unsigned int load_step)
  {
  FEValues<3> fe_values_damage (fe_damage, quadrature_formula_damage,
 
  fe_damage, quadrature_formula_damage, quadrature_formula_damage);
 

The output is a vector of values for all active cells. While it may make sense to compute the value of a solution degree of freedom very accurately, it is usually not necessary to compute the error indicator corresponding to the solution on a cell particularly accurately. We therefore typically use a vector of floats instead of a vector of doubles to represent error indicators.

  Vector<float> estimated_error_per_cell (
  triangulation.n_locally_owned_active_cells ());
 
  KellyErrorEstimator<3>::estimate (dof_handler_damage,
  QGauss<2> (fe_damage.degree + 1),
  { }, locally_relevant_solution_damage, estimated_error_per_cell);
 
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)

Initialize SolutionTransfer object

  parallel::distributed::SolutionTransfer<3, LA::MPI::Vector> soltransDamage (dof_handler_damage);
 
::SolutionTransfer< dim, VectorType, spacedim > SolutionTransfer

Initialize SolutionTransfer object

  parallel::distributed::SolutionTransfer<3, LA::MPI::Vector> soltransElastic (dof_handler_elastic);
 
  triangulation, estimated_error_per_cell, 0.01, // top 1% cells marked for refinement
  0.0); // bottom 0 % cells marked for coarsening
 
  if (triangulation.n_global_levels () >= 4)
  {
  for (const auto &cell : triangulation.active_cell_iterators_on_level (3))
  if (cell->is_locally_owned ())
  cell->clear_refine_flag ();
  }
 
void refine_and_coarsen_fixed_fraction(::Triangulation< dim, spacedim > &tria, const ::Vector< Number > &criteria, const double top_fraction_of_error, const double bottom_fraction_of_error, const VectorTools::NormType norm_type=VectorTools::L1_norm)

prepare the triangulation,

  triangulation.prepare_coarsening_and_refinement ();
 

prepare CellDataStorage object for refinement

  data_transfer.prepare_for_coarsening_and_refinement (triangulation,
  quadrature_point_history_field);
 

prepare the SolutionTransfer object for coarsening and refinement and give the solution vector that we intend to interpolate later,

  soltransDamage.prepare_for_coarsening_and_refinement (
  locally_relevant_solution_damage);
  soltransElastic.prepare_for_coarsening_and_refinement (
  locally_relevant_solution_elastic);
 
  triangulation.execute_coarsening_and_refinement ();
 

redistribute dofs,

  dof_handler_damage.distribute_dofs (fe_damage);
  dof_handler_elastic.distribute_dofs (fe_elastic);
 

Recreate locally_owned_dofs and locally_relevant_dofs index sets

  locally_owned_dofs_damage = dof_handler_damage.locally_owned_dofs ();
  locally_relevant_dofs_damage = DoFTools::extract_locally_relevant_dofs (
  dof_handler_damage);
 
  completely_distributed_solution_damage_old.reinit (
  locally_owned_dofs_damage, mpi_communicator);
  soltransDamage.interpolate (completely_distributed_solution_damage_old);
 

Apply constraints on the interpolated solution to make sure it conforms with the new mesh

  setup_boundary_values_damage ();
 
  constraints_damage.distribute (completely_distributed_solution_damage_old);
 

Copy completely_distributed_solution_damage_old to locally_relevant_solution_damage

  locally_relevant_solution_damage.reinit (locally_owned_dofs_damage,
  locally_relevant_dofs_damage, mpi_communicator);
  locally_relevant_solution_damage =
  completely_distributed_solution_damage_old;
 

Interpolating elastic solution similarly

  locally_owned_dofs_elastic = dof_handler_elastic.locally_owned_dofs ();
  locally_relevant_dofs_elastic = DoFTools::extract_locally_relevant_dofs (
  dof_handler_elastic);
 
  completely_distributed_solution_elastic_old.reinit (
  locally_owned_dofs_elastic, mpi_communicator);
  soltransElastic.interpolate (completely_distributed_solution_elastic_old);
 

Apply constraints on the interpolated solution to make sure it conforms with the new mesh

  setup_constraints_elastic (load_step);
  constraints_elastic.distribute (
  completely_distributed_solution_elastic_old);
 

Copy completely_distributed_solution_damage_old to locally_relevant_solution_damage

  locally_relevant_solution_elastic.reinit (locally_owned_dofs_elastic,
  locally_relevant_dofs_elastic, mpi_communicator);
  locally_relevant_solution_elastic =
  completely_distributed_solution_elastic_old;
 
  for (const auto &cell : triangulation.active_cell_iterators ())
  {
  if (cell->is_locally_owned ())
  quadrature_point_history_field.initialize (cell, 8);
  }
  data_transfer.interpolate ();
 
  }
 
  bool
  PhaseField::check_convergence ()
  {
  LA::MPI::Vector solution_damage_difference (locally_owned_dofs_damage,
  mpi_communicator);
  LA::MPI::Vector solution_elastic_difference (locally_owned_dofs_elastic,
  mpi_communicator);
  LA::MPI::Vector solution_damage_difference_ghost (locally_owned_dofs_damage,
  locally_relevant_dofs_damage, mpi_communicator);
  LA::MPI::Vector solution_elastic_difference_ghost (
  locally_owned_dofs_elastic, locally_relevant_dofs_elastic,
  mpi_communicator);
 
  solution_damage_difference = locally_relevant_solution_damage;
 
  solution_damage_difference -= completely_distributed_solution_damage_old;
 
  solution_elastic_difference = locally_relevant_solution_elastic;
 
  solution_elastic_difference -= completely_distributed_solution_elastic_old;
 
  solution_damage_difference_ghost = solution_damage_difference;
 
  solution_elastic_difference_ghost = solution_elastic_difference;
 
  double error_elastic_solution_numerator, error_elastic_solution_denominator,
  error_damage_solution_numerator, error_damage_solution_denominator;
 
  error_damage_solution_numerator = solution_damage_difference.l2_norm ();
  error_elastic_solution_numerator = solution_elastic_difference.l2_norm ();
  error_damage_solution_denominator =
  completely_distributed_solution_damage.l2_norm ();
  error_elastic_solution_denominator =
  completely_distributed_solution_elastic.l2_norm ();
 
  double error_elastic_solution, error_damage_solution;
  error_damage_solution = error_damage_solution_numerator
  / error_damage_solution_denominator;
  error_elastic_solution = error_elastic_solution_numerator
  / error_elastic_solution_denominator;
 
  if ((error_elastic_solution < tol) && (error_damage_solution < tol))
  return true;
  else
  return false;
  }
 
  void
  PhaseField::update_history_field ()
  {
  FEValues<3> fe_values_damage (fe_damage, quadrature_formula_damage,
 
  for (const auto &cell : dof_handler_damage.active_cell_iterators ())
  if (cell->is_locally_owned ())
  {
  const std::vector<std::shared_ptr<MyQData>> lqph =
  quadrature_point_history_field.get_data (cell);
  for (unsigned int q_index = 0;
  q_index < quadrature_formula_damage.size (); ++q_index)
  lqph[q_index]->value_H = lqph[q_index]->value_H_new;
  }
  }
 
  void
  PhaseField::run ()
  {
  Timer timer;
  timer.start ();
  pcout << "Running with "
  #ifdef USE_PETSC_LA
  << "PETSc"
  #else
  << "Trilinos"
  #endif
  << " on "
  << Utilities::MPI::n_mpi_processes (mpi_communicator) << " MPI rank(s)..."
  << std::endl;
 
  setup_mesh_and_bcs ();
 
Definition timer.h:117
void start()
Definition timer.cc:176
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:99

Create the pre-crack in the domain

  LA::MPI::Vector initial_soln_damage (locally_owned_dofs_damage,
  mpi_communicator);
  for (const auto &cell : dof_handler_damage.active_cell_iterators ())
  if (cell->is_locally_owned ())
  for (const auto vertex_number : cell->vertex_indices ())
  {
  const types::global_dof_index vertex_dof_index =
  cell->vertex_dof_index (vertex_number, 0);
  initial_soln_damage[vertex_dof_index] = 0;
  }
  initial_soln_damage.compress (VectorOperation::insert);
  locally_relevant_solution_damage = initial_soln_damage;
 
Definition types.h:32

Loop over load steps

  for (unsigned int load_step = 1; load_step <= num_load_steps; load_step++)
  {
  pcout << " \n \n load increment number : " << load_step << std::endl;
 

Loop over staggered iterations

  unsigned int iteration = 0;
  bool stoppingCriterion = false;
  while (stoppingCriterion == false)
  {
  pcout << " \n iteration number:" << iteration << std::endl;
  solve_elastic_subproblem (load_step);
  solve_damage_subproblem ();
 
  locally_relevant_solution_damage.update_ghost_values ();
  locally_relevant_solution_elastic.update_ghost_values ();
 
  if (iteration > 0)
  stoppingCriterion = check_convergence ();
 
  completely_distributed_solution_elastic_old =
  locally_relevant_solution_elastic;
  completely_distributed_solution_damage_old =
  locally_relevant_solution_damage;
 
  if (stoppingCriterion == false)
  {
  refine_grid (load_step);
  }
  iteration = iteration + 1;
  }
 

Once converged, do some clean-up operations

  if ((load_step == 1) || (load_step >= 1 && load_step <= num_load_steps
  && std::fabs (load_step % 10) < 1e-6))
  {
  TimerOutput::Scope ts (computing_timer, "output");
  output_results (load_step);
  }
 
  load_disp_calculation (load_step);
 
  computing_timer.print_summary ();
  computing_timer.reset ();
  pcout << std::endl;
 
  update_history_field ();
  }
 
  timer.stop ();
  pcout << "Total run time: " << timer.wall_time () << " seconds."
  << std::endl;
  }
  }
 
  int
  main (int argc,
  char *argv[])
  {
  try
  {
  using namespace dealii;
  using namespace FracturePropagation;
 
  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
 
  PhaseField phasefield;
  phasefield.run ();
  }
  catch (std::exception &exc)
  {
  std::cerr << std::endl << std::endl
  << "----------------------------------------------------" << std::endl;
  std::cerr << "Exception on processing: " << std::endl << exc.what ()
  << std::endl << "Aborting!" << std::endl
  << "----------------------------------------------------" << std::endl;
 
  return 1;
  }
  catch (...)
  {
  std::cerr << std::endl << std::endl
  << "----------------------------------------------------" << std::endl;
  std::cerr << "Unknown exception!" << std::endl << "Aborting!" << std::endl
  << "----------------------------------------------------" << std::endl;
  return 1;
  }
 
  return 0;
 
  }