This program was contributed by Wasim Niyaz Munshi <[email protected]>.
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:

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/numerics/solution_transfer.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/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:
  

Mesh and boundary conditions setup at the beginning

        void
        setup_mesh_and_bcs (); // Called at the start to create a mesh
  

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
  

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
  

Convergence check after each iteration

        bool
        check_convergence ();
  

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;
        ConditionalOStream pcout;
        TimerOutput computing_timer;
  
        parallel::distributed::Triangulation<3> 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;
  

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> ()
          {
          }
  

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

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);
  
      GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions, p1,
          p2); // create coarse mesh
  

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,
          update_values | update_gradients | update_JxW_values
          | update_quadrature_points);
  
      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);
  
      DoFTools::make_hanging_node_constraints (dof_handler_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);
  
      SparsityTools::distribute_sparsity_pattern (dsp,
          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,
          update_values | update_gradients | update_quadrature_points
          | update_JxW_values);
      FEValues<3> fe_values_damage (fe_damage, quadrature_formula_elastic,
          update_values);
  
      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;
  
            const DoFHandler<3>::active_cell_iterator damage_cell =
                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);
              }
  

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

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);
      DoFTools::make_hanging_node_constraints (dof_handler_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);
      SparsityTools::distribute_sparsity_pattern (dsp,
          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,
          update_values | update_gradients | update_JxW_values
          | update_quadrature_points);
      FEValues<3> fe_values_elastic (fe_elastic, quadrature_formula_damage,
          update_values | update_gradients | update_JxW_values
          | update_quadrature_points);
  
      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 (
          3, DataComponentInterpretation::component_is_part_of_vector);
  
      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,
          update_gradients | update_normal_vectors);
      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)
        {

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,
          update_values | update_gradients | update_JxW_values
          | update_quadrature_points);
  
      parallel::distributed::ContinuousQuadratureDataTransfer<3, MyQData> data_transfer (
          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);
  

Initialize SolutionTransfer object

SolutionTransfer<3, LA::MPI::Vector> soltransDamage (dof_handler_damage);

SolutionTransfer

Definition solution_transfer.h:304

Initialize SolutionTransfer object

      SolutionTransfer<3, LA::MPI::Vector> soltransElastic (dof_handler_elastic);
  
      parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction (
          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 ();
        }
  

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,
          update_values | update_gradients | update_JxW_values
          | update_quadrature_points);
  
      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 ();
  

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;
  

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;
  
  }
  