Contents
Phase Field for Brittle Fracture
Introduction
TBA
Theory
TBA
Implementation
TBA
Results
TBA
Complete code
#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>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/petsc_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/matrix_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>
const unsigned int damage_type = 2; // 1 for quadratic // 2 for linear
namespace elas
{
using namespace dealii;
template<int dim>
SymmetricTensor<4, dim> get_stress_strain_tensor(const double lambda, const double mu)
{
SymmetricTensor<4, dim> tmp;
for (unsigned int i = 0; i < dim; ++i)
for (unsigned int j = 0; j < dim; ++j)
for (unsigned int k = 0; k < dim; ++k)
for (unsigned int l = 0; l < dim; ++l)
tmp[i][j][k][l] = (((i == k) && (j == l) ? mu : 0.0) +
((i == l) && (j == k) ? mu : 0.0) +
((i == j) && (k == l) ? lambda : 0.0));
return tmp;
}
template<int dim>
inline SymmetricTensor<2, dim> get_strain(const FEValues<dim>& fe_values,
const unsigned int shape_func,
const unsigned int q_point)
{
SymmetricTensor<2,dim> tmp;
for (unsigned int i = 0; i < dim; ++i)
tmp[i][i] = fe_values.shape_grad_component(shape_func, q_point, i)[i];
for (unsigned int i = 0; i < dim; ++i)
for (unsigned int j = i + 1; j < dim; ++j)
tmp[i][j] =
(fe_values.shape_grad_component(shape_func, q_point, i)[j] +
fe_values.shape_grad_component(shape_func, q_point, j)[i]) /
2;
return tmp;
}
template<int dim>
inline SymmetricTensor<2, dim> get_strain(const std::vector<Tensor<1, dim>>& grad)
{
Assert(grad.size() == dim, ExcInternalError());
SymmetricTensor<2, dim> strain;
for (unsigned int i = 0; i < dim; ++i)
strain[i][i] = grad[i][i];
for (unsigned int i = 0; i < dim; ++i)
for (unsigned int j = i + 1; j < dim; ++j)
strain[i][j] = (grad[i][j] + grad[j][i]) / 2;
return strain;
}
template <int dim>
class Elasticity
{
public:
Elasticity();
void run();
private:
// Mesh generation and refinement
void create_mesh();
void refine_mesh(const unsigned int load_step);
// Elasticity problem functions
void setup_constraints_elastic (const unsigned int load_step); // Posso creare funz che applica BC o mettere a fine assembly matrix
void setup_elasticity(const unsigned int load_step);
void assemble_elasticity(const unsigned int load_step, const bool flag_iter);
void solve_linear_problem(const bool flag_iter);
// Damage problem functions
void setup_constraints_alpha();
void setup_system_alpha();
void assemble_system_alpha(PETScWrappers::MPI::Vector &present_solution_alpha, PETScWrappers::MPI::Vector &system_rhs_alpha);
double gc(const float GC, const float beta, const float x1, const float x2, const Point<dim> &p);
static PetscErrorCode FormFunction(SNES, Vec, Vec, void*);
static PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void*);
// Save the results - Disp, stress, reaction, energies
void output_results(const unsigned int load_step) const;
void output_forces(const unsigned int load_step);
void output_energies(const unsigned int load_step);
// Parallel setup and print on screen
MPI_Comm mpi_communicator;
const unsigned int n_mpi_processes;
const unsigned int this_mpi_process;
parallel::distributed::Triangulation<dim> triangulation;
// Objects for elasticity
const FESystem<dim> fe_elastic; // FE System
DoFHandler<dim> dof_handler_elastic; // DoF Handler
IndexSet locally_owned_dofs_elastic; // IndexSet - Locally Owned
IndexSet locally_relevant_dofs_elastic; // IndexSet - Locally relevant
AffineConstraints<double> constraints_elastic; // Affine Constraint
static const SymmetricTensor<4,dim> stress_strain_tensor;
PETScWrappers::MPI::SparseMatrix system_matrix_elastic; // Elastcitiy Matrix
PETScWrappers::MPI::Vector locally_relevant_solution_elastic; // MPI Split solution
PETScWrappers::MPI::Vector completely_distributed_solution_elastic; // Full Solution
PETScWrappers::MPI::Vector newton_update_solution_elastic; // Full Solution
PETScWrappers::MPI::Vector system_rhs_elastic; // RHS
PETScWrappers::MPI::Vector completely_distributed_solution_elastic_old; // This is for the refinement
const QGauss<dim> quadrature_formula_elastic; // Quadrature formula
// Objects for damage
const FESystem<dim> fe_damage; // FE System
DoFHandler<dim> dof_handler_damage; // DoF Handler
IndexSet locally_owned_dofs_damage; // IndexSet - Locally Owned
IndexSet locally_relevant_dofs_damage; // IndexSet - Locally relevant
AffineConstraints<double> constraints_damage; // Affine Constraint
PETScWrappers::MPI::SparseMatrix system_matrix_damage; // Elastcitiy Matrix
PETScWrappers::MPI::Vector system_rhs_damage; // RHS
PETScWrappers::MPI::Vector locally_relevant_solution_damage; // MPI Split solution
PETScWrappers::MPI::Vector completely_distributed_solution_damage; // Full Solution
PETScWrappers::MPI::Vector completely_distributed_solution_damage_old; // This is for the refinement
PETScWrappers::MPI::Vector locally_damageLB; // MPI Split solution -> mi serve per interpolazione durante il raffittimento
PETScWrappers::MPI::Vector completely_distributed_damageLB; // Full Solution LB -> uso nello SNES
PETScWrappers::MPI::Vector completely_distributed_damageLB_old; // This is for the refinement -> storo la soluzione e la uso per interpolazione
PETScWrappers::MPI::Vector completely_distributed_damageUB; // Full Solution UB -> ricreo e setto 1 ogni volta che setuppo
const QGauss<dim> quadrature_formula_damage; // Quadrature formula
// Domain dimensions
double L = 1;
double H = 0.05;
double W = 0.05;
const unsigned int nx = 25;
const unsigned int ny = 2;
const unsigned int nz = 2;
// Material properties
double E = 100;
double nu = 0.3;
// Damage parameters
double Gc = .01;
double ell = 0.02;
double k_res = 1e-6;
double c_w;
double load_increment = 0.01; // Load increment applied
double applied_load = 0; // Current applied_load += load_increment
double num_load_steps = 15; // Total load increment
ConditionalOStream pcout;
TimerOutput computing_timer;
}; // End of Elast class
// Constructor
template <int dim>
Elasticity<dim>::Elasticity()
: mpi_communicator(MPI_COMM_WORLD)
, n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator))
, this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator))
, triangulation (mpi_communicator)
, fe_elastic (FE_Q<dim> (1), dim)
, dof_handler_elastic (triangulation)
, quadrature_formula_elastic (fe_elastic.degree + 1)
, fe_damage (FE_Q<dim> (1), 1)
, dof_handler_damage (triangulation)
, quadrature_formula_damage (fe_elastic.degree + 1)
, c_w(damage_type == 1 ? 2.0 : 8.0 / 3.0)
, pcout (std::cout, this_mpi_process == 0)
, computing_timer (mpi_communicator, pcout, TimerOutput::never, TimerOutput::wall_times)
{}
// Evaluation of the elastic constant
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));
}
// Evaluation of the damage functions
inline double g_alpha (double &alpha)
{
return (1-alpha)*(1-alpha);
}
inline double g_prime_alpha (double &alpha)
{
return -2.+2.*alpha;
}
inline double g_second_alpha ()
{
return 2.;
}
inline double w_alpha (double &alpha)
{
double tmp;
if (damage_type == 1)
tmp = alpha*alpha;
else if (damage_type == 2)
tmp = alpha;
return tmp;
}
inline double w_prime_alpha (double &alpha)
{
double tmp;
if (damage_type == 1)
tmp = 2.*alpha;
else if (damage_type == 2)
tmp = 1.;
return tmp;
}
inline double w_second_alpha ()
{
double tmp;
if (damage_type == 1)
tmp = 2.;
else if (damage_type == 2)
tmp = 0.;
return tmp;
}
template <int dim>
double Elasticity<dim>::gc(const float GC,
const float beta,
const float x1,
const float x2,
const Point<dim> &p)
{
if (((p[0] - x1) > 1e-6) && ((p[0] - x2) < 1e-6))
return GC * beta;
else
return GC;
}
// -----------------------------
// ----------- MESH ------------
// -----------------------------
template <int dim>
void Elasticity<dim>::create_mesh()
{
// Define the initial coarse subdivision
const std::vector<unsigned int> repetitions = {nx,ny,nz};
// Create Mesh
Point<dim> p1, p2;
if constexpr (dim == 2)
{
p1 = Point<dim>(0, 0);
p2 = Point<dim>(L, H);
}
else if constexpr (dim == 3)
{
p1 = Point<dim>(0, 0, 0);
p2 = Point<dim>(L, H, W);
}
GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions, p1, p2); // create coarse mesh
// Set Fce IDs
for (const auto &cell : triangulation.active_cell_iterators ())
{
if (cell->is_locally_owned())
{
if (cell->center()[0] >= (0.5*L-0.1) && cell->center()[0] <= (0.5*L+0.1))
cell->set_material_id(1);
else
cell->set_material_id(0);
}
for (const auto &face : cell->face_iterators ())
{
if (face->at_boundary ())
{
const auto center = face->center ();
if (std::fabs (center(0) - 0) < 1e-12)
face->set_boundary_id (0);
else if (std::fabs (center(0) - L) < 1e-12)
face->set_boundary_id (1);
else if (std::fabs (center(1) - 0) < 1e-12)
face->set_boundary_id (2);
else if (std::fabs (center(1) - H) < 1e-12)
face->set_boundary_id (3);
// Controlla se siamo in 3D prima di accedere a center(2)
if constexpr (dim == 3)
{
if (std::fabs(center(2) - 0) < 1e-12)
face->set_boundary_id(4);
else if (std::fabs(center(2) - W) < 1e-12)
face->set_boundary_id(5);
}
}
}
}
// Initialize dof handler and extact Owned/Relevant DoFs
dof_handler_damage.distribute_dofs (fe_damage);
dof_handler_elastic.distribute_dofs (fe_elastic);
//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);
// Actual solution -> local; Old solution -> global
locally_relevant_solution_damage.reinit(locally_owned_dofs_damage,
locally_relevant_dofs_damage,
mpi_communicator);
locally_damageLB.reinit (locally_owned_dofs_damage,
locally_relevant_dofs_damage,
mpi_communicator);
completely_distributed_solution_damage_old.reinit (locally_owned_dofs_damage, mpi_communicator);
completely_distributed_damageLB_old.reinit (locally_owned_dofs_damage, mpi_communicator);
completely_distributed_damageLB_old = 0;
// Do the same for elasticity
locally_owned_dofs_elastic = dof_handler_elastic.locally_owned_dofs ();
locally_relevant_dofs_elastic = DoFTools::extract_locally_relevant_dofs (dof_handler_elastic);
// Initialize the old solution for refinement purposes
completely_distributed_solution_elastic_old.reinit (locally_owned_dofs_elastic, mpi_communicator);
pcout << "No. of levels in triangulation: " << triangulation.n_global_levels () << std::endl;
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;
}
template<int dim>
void Elasticity<dim>::refine_mesh(const unsigned int load_step)
{
Vector<float> estimated_error_per_cell (triangulation.n_locally_owned_active_cells ());
KellyErrorEstimator<dim>::estimate (dof_handler_damage,
QGauss<dim-1> (fe_damage.degree + 1),
std::map<types::boundary_id, const Function<dim> *>(), locally_relevant_solution_damage, estimated_error_per_cell);
// Initialize SolutionTransfer object
parallel::distributed::SolutionTransfer<dim, PETScWrappers::MPI::Vector> soltransDamage (dof_handler_damage);
parallel::distributed::SolutionTransfer<dim, PETScWrappers::MPI::Vector> soltransDamageLB (dof_handler_damage);
parallel::distributed::SolutionTransfer<dim, PETScWrappers::MPI::Vector> soltransElastic (dof_handler_elastic);
parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction (
triangulation, estimated_error_per_cell, 0.15, 0.05);
if (triangulation.n_global_levels () >= 7)
{
for (const auto &cell : triangulation.active_cell_iterators_on_level (6))
if (cell->is_locally_owned ())
cell->clear_refine_flag ();
}
// Prepare the triangulation
triangulation.prepare_coarsening_and_refinement ();
// 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);
soltransDamageLB.prepare_for_coarsening_and_refinement (locally_damageLB);
soltransElastic.prepare_for_coarsening_and_refinement (locally_relevant_solution_elastic);
// Do the refinement
triangulation.execute_coarsening_and_refinement ();
// redistribute dofs,
dof_handler_damage.distribute_dofs (fe_damage);
dof_handler_elastic.distribute_dofs (fe_elastic);
// --- DAMAGE
// 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);
completely_distributed_damageLB_old.reinit (locally_owned_dofs_damage, mpi_communicator);
soltransDamageLB.interpolate (completely_distributed_damageLB_old);
// Apply constraints on the interpolated solution to make sure it conforms with the new mesh
setup_constraints_alpha ();
constraints_damage.distribute (completely_distributed_solution_damage_old);
constraints_damage.distribute (completely_distributed_damageLB_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_damageLB.reinit (locally_owned_dofs_damage,
locally_relevant_dofs_damage,
mpi_communicator);
locally_relevant_solution_damage = completely_distributed_solution_damage_old;
locally_damageLB = completely_distributed_damageLB_old;
// --- ELASTICITY
// Recreate locally_owned_dofs and locally_relevant_dofs index sets - ELASTICITY
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 ())
{
if (cell->center()[0] >= (0.5*L-0.1) && cell->center()[0] <= (0.5*L+0.1))
cell->set_material_id(1);
else
cell->set_material_id(0);
}
}
}
// -----------------------------
// -------- ELASTICITY ---------
// -----------------------------
template <int dim>
const SymmetricTensor<4,dim>
Elasticity<dim>::stress_strain_tensor = get_stress_strain_tensor<dim>(/*lambda = */ (100*0.3)/((1+0.3)*(1-2*0.3)),
/*mu = */ 100/(2.*(1+0.3)));
template<int dim>
void Elasticity<dim>::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) - 0.5*L) < 1e-12) // Face at L/2
{
for (const auto vertex_number : cell->vertex_indices ())
{
const auto vert = cell->vertex (vertex_number);
const double y_bot = 0;
if (std::fabs (vert (1) - 0) < 1e-12 && std::fabs (vert (0) - L/2) < 1e-12)
{
const unsigned int y_dof = cell->vertex_dof_index (vertex_number, 1);
constraints_elastic.add_line (y_dof);
constraints_elastic.set_inhomogeneity (y_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);
}
}
}
}
}
}
const FEValuesExtractors::Scalar u_x (0);
const ComponentMask u_x_mask = fe_elastic.component_mask (u_x);
applied_load = load_increment*load_step;
const double u_x_values_right = 0.5*applied_load;
const double u_x_values_left = -0.5*applied_load;
VectorTools::interpolate_boundary_values (dof_handler_elastic,
0,
Functions::ConstantFunction<dim> (u_x_values_left, dim),
constraints_elastic,
u_x_mask);
VectorTools::interpolate_boundary_values (dof_handler_elastic,
1,
Functions::ConstantFunction<dim> (u_x_values_right, dim),
constraints_elastic,
u_x_mask);
constraints_elastic.close (); */
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);
applied_load = load_increment*load_step;
const double u_x_values_right = 0.5*applied_load;
const double u_x_values_left = -0.5*applied_load;
/* 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, 0,
Functions::ConstantFunction<dim> (u_x_values_left, dim),
constraints_elastic, u_x_mask);
VectorTools::interpolate_boundary_values (dof_handler_elastic, 0,
Functions::ConstantFunction<dim> (u_fix, dim), constraints_elastic,
u_y_mask);
VectorTools::interpolate_boundary_values (dof_handler_elastic, 1,
Functions::ConstantFunction<dim> (u_x_values_right, dim), constraints_elastic,
u_x_mask);
constraints_elastic.close ();
}
template<int dim>
void Elasticity<dim>::setup_elasticity(const unsigned int load_step)
{
TimerOutput::Scope ts(computing_timer, "setup_system_elastic");
// Initialize the dofs contained locally as well as those of the near cells
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);
// Initialize the contraints objects for BC and Hanging Nodes during refinement
system_rhs_elastic.reinit (locally_owned_dofs_elastic, mpi_communicator);
completely_distributed_solution_elastic.reinit (locally_owned_dofs_elastic,
mpi_communicator);
newton_update_solution_elastic.reinit(locally_owned_dofs_elastic,
mpi_communicator);
setup_constraints_elastic (load_step); // Qui unizializzo gli hanging nodes constrains (e le BC volendo)
// Sparsity pattern initialization (where matrix has nonzero entries)
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);
// Elasticity Matrix
system_matrix_elastic.reinit (locally_owned_dofs_elastic,
locally_owned_dofs_elastic,
dsp,
mpi_communicator);
}
template<int dim>
void Elasticity<dim>::assemble_elasticity(const unsigned int load_step, const bool flag_iter)
{
TimerOutput::Scope ts (computing_timer, "assembly_elastic");
// Clean the current matrices
system_rhs_elastic = 0;
system_matrix_elastic = 0;
// Initialize the FEValues Objects
FEValues<dim> fe_values_elastic (fe_elastic, quadrature_formula_elastic,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
FEValues<dim> fe_values_damage (fe_damage, quadrature_formula_elastic,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
// Number of DoFs and Gauss points per cell
const unsigned int dofs_per_cell = fe_elastic.n_dofs_per_cell ();
const unsigned int n_q_points = quadrature_formula_elastic.size ();
// Local Matrices and Vector
FullMatrix<double> cell_matrix_elastic (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs_elastic (dofs_per_cell);
// Create vector to store the local dof indices
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
std::vector<double> damage_values (n_q_points);
std::vector<SymmetricTensor<2, dim>> strain_values (n_q_points);
// loop over the cells
for (const auto &cell : dof_handler_elastic.active_cell_iterators ())
{
if (cell->is_locally_owned ())
{
// Reset the local matrix/rhs of the cell
cell_matrix_elastic = 0;
cell_rhs_elastic = 0;
const typename DoFHandler<dim>::active_cell_iterator damage_cell = cell->as_dof_handler_iterator(dof_handler_damage);
// const DoFHandler<dim>::active_cell_iterator damage_cell = Triangulation<dim>::active_cell_iterator (cell)->as_dof_handler_iterator (dof_handler_damage);
fe_values_elastic.reinit (cell);
fe_values_damage.reinit (damage_cell);
const FEValuesExtractors::Vector displacements (/* first vector component = */0);
fe_values_elastic[displacements].get_function_symmetric_gradients (locally_relevant_solution_elastic, strain_values);
fe_values_damage.get_function_values(locally_relevant_solution_damage, damage_values);
// loop over the gauss point of the cell
for (const unsigned int q_point : fe_values_elastic.quadrature_point_indices ())
{
// Estraggo i valori relevant sui punti di gauss
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) += (std::pow((1 - d), 2) + k_res)
* ((fe_values_elastic.shape_grad (i, q_point)[component_i] * lambda (E, nu) *
fe_values_elastic.shape_grad (j, q_point)[component_j])
+
(fe_values_elastic.shape_grad (i, q_point)[component_j] * mu (E, nu) *
fe_values_elastic.shape_grad (j, q_point)[component_i])
+
((component_i == component_j) ? (fe_values_elastic.shape_grad (i, q_point) * mu (E, nu) *
fe_values_elastic.shape_grad (j, q_point))
: 0)
) * fe_values_elastic.JxW (q_point); //
}
cell_rhs_elastic (i) += fe_values_elastic.shape_value (i, q_point) * 0 * fe_values_elastic.JxW (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);
/* // Now set BC
const FEValuesExtractors::Scalar u_x (0);
const FEValuesExtractors::Scalar u_y (1);
const FEValuesExtractors::Scalar u_z (2);
const ComponentMask u_x_mask = fe_elastic.component_mask (u_x);
const ComponentMask u_y_mask = fe_elastic.component_mask (u_y);
const ComponentMask u_z_mask = fe_elastic.component_mask (u_z);
applied_load = load_increment*load_step;
const double u_x_values_right = 0.5*applied_load;
const double u_x_values_left = -0.5*applied_load;
std::map<types::global_dof_index, double> boundary_values;
Point<dim> point_one, point_two;
point_one(0) = 0.5*L;
point_one(1) = 0.;
point_one(2) = 0.;
point_two(0) = 0.5*L;
point_two(1) = H;
point_two(2) = 0.;
typename DoFHandler<dim>::active_cell_iterator cell = dof_handler_elastic.begin_active(),
endc = dof_handler_elastic.end(); */
/* // Set up Punctual BC
bool evaluation_point_found = false;
for (; (cell!=endc) && !evaluation_point_found; ++cell)
{
if (cell->is_locally_owned())
{
for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex)
{
if (cell->vertex(vertex).distance (point_one) < cell->diameter() * 1e-12)
{ // Block displacement on the mid point at y = 0
boundary_values[cell->vertex_dof_index(vertex,0)]=0.;
boundary_values[cell->vertex_dof_index(vertex,1)]=0.;
boundary_values[cell->vertex_dof_index(vertex,2)]=0.;
evaluation_point_found = true;
}
if (cell->vertex(vertex).distance (point_two) < cell->diameter() * 1e-12)
{ // Block only ux, uz at midpoint at y = H
boundary_values[cell->vertex_dof_index(vertex,0)]=0.;
boundary_values[cell->vertex_dof_index(vertex,2)]=0.;
evaluation_point_found = true;
}
}
}
} */
// Set the displacement on the left face
/* if (flag_iter == true)
{
VectorTools::interpolate_boundary_values (dof_handler_elastic,
0,
Functions::ConstantFunction<dim> (u_x_values_left, dim),
boundary_values,
u_x_mask);
VectorTools::interpolate_boundary_values (dof_handler_elastic,
1,
Functions::ConstantFunction<dim> (u_x_values_right, dim),
boundary_values,
u_x_mask);
} else
{
VectorTools::interpolate_boundary_values (dof_handler_elastic,
0,
Functions::ConstantFunction<dim> (0., dim),
boundary_values,
u_x_mask);
VectorTools::interpolate_boundary_values (dof_handler_elastic,
1,
Functions::ConstantFunction<dim> (0., dim),
boundary_values,
u_x_mask);
}
PETScWrappers::MPI::Vector tmp(locally_owned_dofs_elastic, mpi_communicator);
MatrixTools::apply_boundary_values(boundary_values, system_matrix_elastic, tmp, system_rhs_elastic, false);
if (flag_iter == true)
completely_distributed_solution_elastic = tmp;
else
newton_update_solution_elastic = tmp; */
}
template<int dim>
void Elasticity<dim>::solve_linear_problem(const bool solve_step)
{
TimerOutput::Scope ts (computing_timer, "solve_linear_system_elastic");
// CG
SolverControl solver_control (10000, 1e-12* system_rhs_elastic.l2_norm());
SolverCG<PETScWrappers::MPI::Vector> solver (solver_control);
/* PETScWrappers::PreconditionBoomerAMG::AdditionalData data;
PETScWrappers::PreconditionBoomerAMG preconditioner;
preconditioner.initialize (system_matrix_elastic, data); */
PETScWrappers::PreconditionBlockJacobi preconditioner(system_matrix_elastic);
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;
}
// -----------------------------
// ---------- DAMAGE -----------
// -----------------------------
template<int dim>
void Elasticity<dim>::setup_constraints_alpha()
{
constraints_damage.clear ();
constraints_damage.reinit (locally_relevant_dofs_damage);
DoFTools::make_hanging_node_constraints (dof_handler_damage, constraints_damage);
constraints_damage.close ();
}
template<int dim>
void Elasticity<dim>::setup_system_alpha()
{
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);
locally_damageLB.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);
completely_distributed_damageLB.reinit (locally_owned_dofs_damage, mpi_communicator);
completely_distributed_damageLB = completely_distributed_damageLB_old;
completely_distributed_damageUB.reinit (locally_owned_dofs_damage, mpi_communicator);
completely_distributed_damageUB = 1;
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);
}
template<int dim>
void Elasticity<dim>::assemble_system_alpha(PETScWrappers::MPI::Vector &present_solution_alpha, PETScWrappers::MPI::Vector &system_rhs_alpha)
{
TimerOutput::Scope ts (computing_timer, "assembly_elastic");
// Clean the current matrices
system_rhs_damage = 0;
system_matrix_damage = 0;
FEValues<dim> fe_values_damage (fe_damage, quadrature_formula_damage,
update_values | update_gradients |
update_JxW_values | update_quadrature_points);
FEValues<dim> 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);
std::vector<SymmetricTensor<2, dim>> strain_values (n_q_points);
std::vector<double> local_solution_alpha(n_q_points);
std::vector<Tensor<1, dim>> local_grad_alpha (n_q_points);
for (const auto &cell : dof_handler_damage.active_cell_iterators ())
{
if (cell->is_locally_owned ())
{
// Clear degli elementi locali
cell_matrix_damage = 0;
cell_rhs_damage = 0;
// Inizializzo i fe_values sulla cella
fe_values_damage.reinit (cell);
const typename DoFHandler<dim>::active_cell_iterator elastic_cell = cell->as_dof_handler_iterator(dof_handler_elastic);
// const DoFHandler<dim>::active_cell_iterator elastic_cell = Triangulation<dim>::active_cell_iterator (cell)->as_dof_handler_iterator (dof_handler_elastic);
fe_values_elastic.reinit (elastic_cell);
// Estraggo valori danno ed elasticità
fe_values_damage.get_function_values(locally_relevant_solution_damage, local_solution_alpha);
fe_values_damage.get_function_gradients(locally_relevant_solution_damage, local_grad_alpha);
const FEValuesExtractors::Vector displacements (/* first vector component = */0);
fe_values_elastic[displacements].get_function_symmetric_gradients (locally_relevant_solution_elastic, strain_values);
// Loop sui punti di gauss
for (const unsigned int q_index : fe_values_damage.quadrature_point_indices ())
{
// Estraggo i valori relevant sui punti di gauss
const SymmetricTensor<2,dim> eps_u = strain_values[q_index];
double elastic_source=0.;
elastic_source = 0.5* eps_u * stress_strain_tensor * eps_u;
const auto &x_q = fe_values_damage.quadrature_point (q_index);
double Gc_dmg = gc(Gc, 0.95, 0.495, 0.505, x_q);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
cell_matrix_damage(i,j)+= (fe_values_damage.shape_grad(i,q_index) *
fe_values_damage.shape_grad(j,q_index) *
(2.*Gc_dmg*ell/c_w)
+
fe_values_damage.shape_value(i,q_index) *
fe_values_damage.shape_value(j,q_index) *
(Gc_dmg/(ell*c_w)*w_second_alpha() + g_second_alpha()*elastic_source))
* fe_values_damage.JxW(q_index);
}
cell_rhs_damage(i) += (g_prime_alpha(local_solution_alpha[q_index])*elastic_source+
w_prime_alpha(local_solution_alpha[q_index])*Gc_dmg/(ell*c_w)) *
fe_values_damage.shape_value(i,q_index) *
fe_values_damage.JxW(q_index);
cell_rhs_damage(i) += (local_grad_alpha[q_index]*2.*Gc_dmg*ell/c_w) *
fe_values_damage.shape_grad(i,q_index) *
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);
}
}
system_matrix_damage.compress (VectorOperation::add);
system_rhs_damage.compress (VectorOperation::add);
// Set BC on damage
std::map<types::global_dof_index, double> boundary_values;
VectorTools::interpolate_boundary_values(dof_handler_damage,
0,
Functions::ZeroFunction<dim>(),
boundary_values);
VectorTools::interpolate_boundary_values(dof_handler_damage,
1,
Functions::ZeroFunction<dim>(),
boundary_values);
// Qui riuso il vettore petsc MPI cosi da applicarci le BC
PETScWrappers::MPI::Vector tmp(locally_owned_dofs_damage, mpi_communicator);
MatrixTools::apply_boundary_values(boundary_values, system_matrix_damage, tmp, system_rhs_damage, false);
completely_distributed_solution_damage = tmp;
}
// Function used to create the FormFunction Input for the SNESSetFunction
template <int dim>
PetscErrorCode Elasticity<dim>::FormFunction(SNES snes, Vec x, Vec f, void* ctx)
{
const auto p_ctx = static_cast<Elasticity*>(ctx);
p_ctx->assemble_system_alpha(p_ctx->completely_distributed_solution_damage, p_ctx->system_rhs_damage);
return 0;
}
// Function used to create the FormFunction Input for the SNESSetJacobian
template <int dim>
PetscErrorCode Elasticity<dim>::FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void* ctx)
{
const auto p_ctx = static_cast<Elasticity*>(ctx);
p_ctx->system_matrix_damage;
return 0;
}
// -----------------------------
// ---------- OUTPUT -----------
// -----------------------------
template <int dim>
void Elasticity<dim>::output_results (const unsigned int load_step) const
{
std::vector<std::string> displacement_names (dim, "displacement");
std::vector<DataComponentInterpretation::DataComponentInterpretation> displacement_component_interpretation (dim, DataComponentInterpretation::component_is_part_of_vector);
DataOut<dim> 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 ();
const std::string pvtu_filename = data_out_phasefield.write_vtu_with_pvtu_record ("./", "solution", load_step, mpi_communicator, 2, 0);
if (this_mpi_process == 0)
{
static std::vector<std::pair<double, std::string>> times_and_names;
times_and_names.push_back(std::pair<double, std::string>(load_step, pvtu_filename));
std::ofstream pvd_output("./solution_u.pvd");
DataOutBase::write_pvd_record(pvd_output, times_and_names);
}
}
// -----------------------------
// ----------- RUN -------------
// -----------------------------
template <int dim>
void Elasticity<dim>::run()
{
Timer timer;
timer.start ();
pcout << "Running on " << Utilities::MPI::n_mpi_processes (mpi_communicator) << " MPI rank(s)..." << std::endl;
create_mesh();
// Loop over load steps
for (unsigned int load_step = 0; load_step <= num_load_steps; load_step++)
{
pcout << " --- Load increment number : " << load_step << std::endl;
double error_alternate = 1.0;
const double error_toll = 1.e-5;
unsigned int iter_counter_am = 0;
const unsigned int max_iteration = 500;
double error_elastic;
const double error_toll_elastic=1.e-6;
unsigned int iter_elastic;
const unsigned int max_iteration_elastic = 10;
bool solve_step;
// Alternate Minimization
while (error_alternate > error_toll && iter_counter_am <max_iteration)
{
pcout << " ------ Alternate minimization: " << iter_counter_am << std::endl;
// --- Elasticity problem
pcout << " --- Solving Elastic problem..." << std::endl;
iter_elastic = 0;
solve_step = true;
setup_elasticity (load_step);
assemble_elasticity (load_step, solve_step);
solve_linear_problem (solve_step);
// Damage problem
pcout << " --- Solving Damage problem..." << std::endl;
setup_constraints_alpha();
setup_system_alpha();
// Solve with SNES
SNES snes;
PetscErrorCode ierr;
ierr = SNESCreate(mpi_communicator, &snes);
ierr = SNESSetType(snes, SNESVINEWTONRSLS);
ierr = SNESSetFunction(snes, system_rhs_damage, FormFunction, this);
ierr = SNESSetJacobian(snes, system_matrix_damage, system_matrix_damage, FormJacobian, this);
ierr = SNESVISetVariableBounds(snes, completely_distributed_damageLB, completely_distributed_damageUB);
ierr = SNESSetFromOptions(snes);
ierr = SNESSolve(snes, nullptr, completely_distributed_solution_damage);
ierr = SNESDestroy(&snes);
// End SNES
constraints_damage.distribute(completely_distributed_solution_damage);
locally_relevant_solution_damage = completely_distributed_solution_damage;
// Update ghost values after the solution
locally_relevant_solution_elastic.update_ghost_values();
locally_relevant_solution_damage.update_ghost_values();
// Check convergence
PETScWrappers::MPI::Vector temp_damage(locally_owned_dofs_damage, mpi_communicator);
temp_damage = locally_relevant_solution_damage;
temp_damage -= completely_distributed_solution_damage_old;
error_alternate = temp_damage.linfty_norm();
pcout << " Number of iteration: " << iter_counter_am << std::endl;
pcout << " Error_on_alpha: " << error_alternate << std::endl;
pcout << " Alpha_max: " << completely_distributed_solution_damage.linfty_norm() << std::endl;
// Prepare for refinement
completely_distributed_solution_elastic_old = locally_relevant_solution_elastic;
completely_distributed_solution_damage_old = locally_relevant_solution_damage;
if (error_alternate > error_toll)
{
refine_mesh(load_step);
pcout << " Total New Number of globally active cells: " << triangulation.n_global_active_cells () << std::endl
<< " New Number of degrees of freedom for elasticity: " << dof_handler_elastic.n_dofs () << std::endl;
}
output_results (load_step*100 + iter_counter_am);
iter_counter_am++;
}
completely_distributed_damageLB = completely_distributed_solution_damage;
completely_distributed_damageLB_old = completely_distributed_solution_damage;
output_results (load_step*100+max_iteration);
computing_timer.print_summary ();
computing_timer.reset ();
pcout << std::endl;
}
timer.stop ();
pcout << "Total run time: " << timer.wall_time () << " seconds." << std::endl;
}
} // End of namespace Elasticity
// -----------------------------
// ---------- MAIIN ------------
// -----------------------------
int main (int argc, char *argv[])
{
try
{
using namespace dealii;
using namespace elas;
Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
Elasticity<2> elasticity;
elasticity.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;
}