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:
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 debug
mode. 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:
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: 
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
#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
}
#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
{
class PhaseField
{
Function declarations
public:
PhaseField ();
void
run ();
private:
- Mesh and boundary conditions setup at the beginning
void
setup_mesh_and_bcs ();
- 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);
- Damage subproblem setup and solution
void
setup_boundary_values_damage ();
void
setup_system_damage ();
void
assemble_system_damage ();
void
solve_linear_system_damage ();
double
void
solve_damage_subproblem ();
- 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;
void
load_disp_calculation (const unsigned int load_step);
void
update_history_field ();
void
refine_grid (const unsigned int load_step);
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Objects for elasticity
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;
Objects for 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 double ux = 1e-3;
const double alpha = 1;
const double uy = alpha * ux;
const double l = 0.6;
const unsigned int num_load_steps = 100;
const double tol = 1e-2;
const double GC = 1e-3;
const double E = 37.5;
const double beta = 25;
const double nu = 0.25;
Objects for load-displacement calculation
{
public:
MyQData () = default;
virtual
~MyQData () = default;
unsigned int
{
return 2;
}
virtual void
pack_values (std::vector<double> &scalars)
const override
{
Assert (scalars.size () == 2, ExcInternalError ());
scalars[0] = value_H;
scalars[1] = value_H_new;
}
virtual void
{
Assert (scalars.size () == 2, ExcInternalError ());
value_H = scalars[0];
value_H_new = scalars[1];
}
double value_H;
double value_H_new;
};
};
inline double
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,
{
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 ()
:
{
}
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)
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 = 0) const override
{
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)
-(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)
{
static_cast<double> ((rand ()) / RAND_MAX) * Domain::x_max;
#define AssertDimension(dim1, dim2)
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 ();
}
inline double
{
return p[0] - p[0] + 1;
}
void
right_hand_side_elastic (
const std::vector<
Point<3>> &points,
{
for (unsigned int point_n = 0; point_n < points.size (); ++point_n)
{
}
}
void
Traction_elastic (
const std::vector<
Point<3>> &points,
{
for (unsigned int point_n = 0; point_n < points.size (); ++point_n)
{
}
}
PhaseField::PhaseField ()
:
mpi_communicator (MPI_COMM_WORLD),
(
Utilities::
MPI::this_mpi_process (mpi_communicator) == 0)),
computing_timer (mpi_communicator, pcout,
TimerOutput::never,
fe_elastic (
FE_Q<3> (1), 3),
quadrature_formula_elastic (fe_elastic.degree + 1),
fe_damage (1),
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
{
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);
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)
::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 &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: "
dof_handler_damage.distribute_dofs (fe_damage);
dof_handler_elastic.distribute_dofs (fe_elastic);
pcout << " Number of locally owned cells on the process: "
pcout <<
"Number of global cells:" <<
triangulation.n_global_active_cells ()
<< std::endl;
pcout << " Total Number of globally active cells: "
<< " 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 ();
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 ();
dof_handler_elastic);
completely_distributed_solution_elastic_old.reinit (
locally_owned_dofs_elastic, mpi_communicator);
{
if (cell->is_locally_owned ())
{
quadrature_point_history_field.initialize (cell, 8);
}
}
FEValues<3> fe_values_damage (fe_damage, quadrature_formula_damage,
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)
{
{
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) < 1
e-12 && std::fabs (
vert (
0)
- Domain::x_min)
{
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) < 1
e-12)
{
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 double u_x_values_right = ux * load_step;
const double u_y_values = uy * load_step;
const double u_fix = 0.0;
constraints_elastic, u_x_mask);
u_y_mask);
u_y_mask);
constraints_elastic.close ();
}
void
PhaseField::setup_system_elastic (const unsigned int load_step)
{
locally_owned_dofs_elastic = dof_handler_elastic.locally_owned_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);
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 ()
{
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 ();
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;
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) +=
(fe_values_elastic.shape_grad (i, q_point)[component_i] *
fe_values_elastic.shape_grad (j, q_point)[component_j]
*
+
(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
TriaActiveIterator< CellAccessor< dim, spacedim > > active_cell_iterator
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.
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
traction contribution to rhs
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);
}
}
void
PhaseField::solve_linear_system_elastic ()
{
1e-12* system_rhs_elastic.l2_norm());
LA::MPI::PreconditionAMG::AdditionalData
data;
#ifdef USE_PETSC_LA
data.symmetric_operator =
true;
#else
std::vector< index_type > data
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 ()
{
constraints_damage.clear ();
constraints_damage.reinit (locally_relevant_dofs_damage);
constraints_damage);
Create initial crack, if any
for (const auto &face : cell->face_iterators())
{
const auto vert = cell->vertex(vertex_number);
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))
{
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 ()
{
locally_owned_dofs_damage = dof_handler_damage.locally_owned_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);
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 ()
{
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 ();
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);
dof_handler_elastic);
fe_values_damage.reinit (cell);
fe_values_elastic.reinit (elastic_cell);
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;
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)
*
fe_values_damage.shape_grad (j, q_index)
*
fe_values_damage.JxW (q_index)
+
Contribution to stiffness from u term
((1 + (2 * l * H) / g_c) * (1 / pow (l, 2))
* fe_values_damage.shape_value (i, q_index) *
fe_values_damage.shape_value (j, q_index) *
fe_values_damage.JxW (q_index));
}
const auto &x_q = fe_values_damage.quadrature_point(q_index);
cell_rhs_damage (i) += (fe_values_damage.shape_value (i,
q_index)
*
(2 / (l * g_c))
* H * fe_values_damage.JxW (q_index));
}
}
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);
}
}
void
PhaseField::solve_linear_system_damage ()
{
1e-12* system_rhs_damage.l2_norm());
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 (
locally_relevant_solution_elastic, displacement_names,
displacement_component_interpretation);
data_out_phasefield.add_data_vector (dof_handler_damage,
locally_relevant_solution_damage, "damage");
for (unsigned int i = 0; i < subdomain.size (); ++i)
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)
{
const QGauss<2> face_quadrature (fe_elastic.degree);
std::vector<SymmetricTensor<2, 3>> strain_values (face_quadrature.size ());
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 double tr_strain = strain[0][0] + strain[1][1] + strain[2][2];
const double d = damage_values[q];
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);
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 double tr_strain = strain[0][0] + strain[1][1] + strain[2][2];
const double d = damage_values[q];
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.
@ component_is_part_of_vector
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
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.
{ }, 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
::SolutionTransfer< dim, VectorType, spacedim > SolutionTransfer
Initialize SolutionTransfer object
0.0);
{
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,
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);
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 ();
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 ();
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;
{
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 ()
{
pcout << "Running with "
#ifdef USE_PETSC_LA
<< "PETSc"
#else
<< "Trilinos"
#endif
<< " on "
<< std::endl;
setup_mesh_and_bcs ();
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
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 ())
{
const
types::global_dof_index vertex_dof_index =
cell->vertex_dof_index (vertex_number, 0);
initial_soln_damage[vertex_dof_index] = 0;
}
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))
{
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 FracturePropagation;
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;
}