117 * #include <deal.II/base/quadrature_lib.h>
118 * #include <deal.II/base/function.h>
119 * #include <deal.II/base/timer.h>
120 * #include <deal.II/lac/generic_linear_algebra.h>
122 * #define FORCE_USE_OF_TRILINOS
125 * #
if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \
126 * !(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
127 *
using namespace dealii::LinearAlgebraPETSc;
128 * # define USE_PETSC_LA
129 * #elif defined(DEAL_II_WITH_TRILINOS)
130 *
using namespace dealii::LinearAlgebraTrilinos;
132 * # error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
136 * #include <deal.II/lac/vector.h>
137 * #include <deal.II/lac/full_matrix.h>
138 * #include <deal.II/lac/precondition.h>
139 * #include <deal.II/lac/solver_cg.h>
140 * #include <deal.II/lac/affine_constraints.h>
141 * #include <deal.II/lac/dynamic_sparsity_pattern.h>
142 * #include <deal.II/grid/grid_generator.h>
143 * #include <deal.II/dofs/dof_handler.h>
144 * #include <deal.II/dofs/dof_tools.h>
145 * #include <deal.II/dofs/dof_accessor.h>
146 * #include <deal.II/grid/tria_accessor.h>
147 * #include <deal.II/fe/fe_values.h>
148 * #include <deal.II/fe/fe_system.h>
149 * #include <deal.II/fe/fe_q.h>
150 * #include <deal.II/grid/grid_refinement.h>
151 * #include <deal.II/numerics/vector_tools.h>
152 * #include <deal.II/numerics/data_out.h>
153 * #include <deal.II/numerics/error_estimator.h>
154 * #include <deal.II/base/utilities.h>
155 * #include <deal.II/base/conditional_ostream.h>
156 * #include <deal.II/base/index_set.h>
157 * #include <deal.II/lac/sparsity_tools.h>
158 * #include <deal.II/distributed/tria.h>
159 * #include <deal.II/distributed/grid_refinement.h>
160 * #include <deal.II/distributed/solution_transfer.h>
161 * #include <deal.II/base/quadrature_point_data.h>
162 * #include <deal.II/base/tensor_function.h>
164 * #include <iostream>
167 *
namespace FracturePropagation
188 * 1. Mesh and boundary conditions setup at the beginning
192 * setup_mesh_and_bcs ();
196 * 2. Elastic subproblem setup and solution
200 * setup_constraints_elastic (
const unsigned int load_step);
202 * setup_system_elastic (
const unsigned int load_step);
204 * assemble_system_elastic ();
206 * solve_linear_system_elastic ();
208 * solve_elastic_subproblem (
const unsigned int load_step);
212 * 3. Damage subproblem setup and solution
216 * setup_boundary_values_damage ();
218 * setup_system_damage ();
220 * assemble_system_damage ();
222 * solve_linear_system_damage ();
226 * solve_damage_subproblem ();
230 * 4. Convergence
check after each iteration
234 * check_convergence ();
238 * 5. Post-processing: output results,
refine grid, and calculate displacement
242 * output_results (
const unsigned int load_step)
const;
244 * load_disp_calculation (
const unsigned int load_step);
246 * update_history_field ();
248 * refine_grid (
const unsigned int load_step);
258 * Objects
for elasticity
263 *
IndexSet locally_owned_dofs_elastic;
264 *
IndexSet locally_relevant_dofs_elastic;
266 * LA::MPI::SparseMatrix system_matrix_elastic;
267 * LA::MPI::Vector locally_relevant_solution_elastic;
268 * LA::MPI::Vector completely_distributed_solution_elastic_old;
269 * LA::MPI::Vector completely_distributed_solution_elastic;
270 * LA::MPI::Vector system_rhs_elastic;
271 *
const QGauss<3> quadrature_formula_elastic;
280 *
IndexSet locally_owned_dofs_damage;
281 *
IndexSet locally_relevant_dofs_damage;
283 * LA::MPI::SparseMatrix system_matrix_damage;
284 * LA::MPI::Vector locally_relevant_solution_damage;
285 * LA::MPI::Vector completely_distributed_solution_damage_old;
286 * LA::MPI::Vector completely_distributed_solution_damage;
287 * LA::MPI::Vector system_rhs_damage;
288 *
const QGauss<3> quadrature_formula_damage;
290 *
const double ux = 1
e-3;
291 *
const double alpha = 1;
292 *
const double uy = alpha * ux;
293 *
const double l = 0.6;
294 *
const unsigned int num_load_steps = 100;
295 *
const double tol = 1
e-2;
296 *
const double GC = 1
e-3;
297 *
const double E = 37.5;
298 *
const double beta = 25;
299 *
const double nu = 0.25;
303 * Objects
for load-displacement calculation
313 * MyQData () =
default;
315 * ~MyQData () =
default;
324 *
pack_values (std::vector<double> &scalars)
const override
326 *
Assert (scalars.size () == 2, ExcInternalError ());
327 * scalars[0] = value_H;
328 * scalars[1] = value_H_new;
332 *
unpack_values (
const std::vector<double> &scalars)
override
334 *
Assert (scalars.size () == 2, ExcInternalError ());
335 * value_H = scalars[0];
336 * value_H_new = scalars[1];
340 *
double value_H_new;
351 *
return (E * nu) / ((1 + nu) * (1 - 2 * nu));
358 *
return E / (2 * (1 + nu));
362 * gc (
const float GC,
368 *
if (((p[2] - z1) > 1e-6) && ((p[2] - z2) < 1e-6))
375 * is_in_middle_layer (
const Point<3> &cell_center,
const float z1,
const float z2)
377 *
return (cell_center[2] >= z1 && cell_center[2] <= z2);
383 *
const double x_min = 0;
384 *
const double y_min = 0;
385 *
const double z_min = 0;
386 *
const double x_max = 30;
387 *
const double y_max = 30;
388 *
const double z_max = 13;
389 *
const double z1 = 5;
390 *
const double z2 = 8;
393 *
namespace RandomMedium
395 *
class EnergyReleaseRate :
public Function<3>
398 * EnergyReleaseRate ()
406 * value_list method calculates the fracture_energy (Gc) at a set of points and stores the result in a vector of
zero order tensors.
411 * std::vector<double> &values,
412 *
const unsigned int = 0) const override
416 *
for (
unsigned int p = 0; p < points.size (); ++p)
420 *
double fracture_energy = 0;
421 *
for (
unsigned int i = 0; i < centers.size (); ++i)
423 * -(points[p] - centers[i]).norm_square () / (1.5 * 1.5));
425 *
const double normalized_fracture_energy =
std::min (
426 *
std::max (fracture_energy, 4e-5), 4e-4);
428 *
values[p] = normalized_fracture_energy;
433 *
static std::vector<Point<3>> centers;
435 *
static std::vector<Point<3>>
438 *
const unsigned int N = 1000;
440 * std::vector<Point<3>> centers_list (N);
441 *
for (
unsigned int i = 0; i <
N; ++i)
442 *
for (
unsigned int d = 0;
d < 3; ++
d)
443 *
if (d == 0 || d == 1)
445 * centers_list[i][
d] =
446 *
static_cast<double> ((rand ()) / RAND_MAX) * Domain::x_max;
449 * x,y will be between 0 to x_max
455 * centers_list[i][
d] =
static_cast<double> (Domain::z1
456 * + ((rand ()) / RAND_MAX) * (Domain::z2 - Domain::z1));
459 *
return centers_list;
463 * std::vector<Point<3>> EnergyReleaseRate::centers =
464 * EnergyReleaseRate::get_centers ();
470 * Conductivity_damage (
const Point<3> &p)
472 *
return p[0] - p[0] + 1;
476 * right_hand_side_elastic (
const std::vector<
Point<3>> &points,
481 *
for (
unsigned int point_n = 0; point_n < points.size (); ++point_n)
490 * Traction_elastic (
const std::vector<
Point<3>> &points,
495 *
for (
unsigned int point_n = 0; point_n < points.size (); ++point_n)
504 * PhaseField::PhaseField ()
506 * mpi_communicator (MPI_COMM_WORLD),
514 * quadrature_formula_elastic (fe_elastic.degree + 1),
517 * quadrature_formula_damage (fe_elastic.degree + 1),
518 * load_values_x(num_load_steps+1),
519 * load_values_y(num_load_steps+1),
520 * displacement_values(num_load_steps+1)
527 *
double Mac_tr_strain, Mac_first_principal_strain,
528 * Mac_second_principal_strain, Mac_third_principal_strain,
529 * tr_sqr_Mac_Principal_strain;
531 *
const double tr_strain =
trace (strain);
533 * Mac_tr_strain = tr_strain >0 ? tr_strain : 0;
534 *
const std::array<double, 3> Principal_strains =
eigenvalues (strain);
536 * Mac_first_principal_strain = (Principal_strains[0] > 0) ? Principal_strains[0] : 0;
537 * Mac_second_principal_strain = (Principal_strains[1] > 0) ? Principal_strains[1] : 0;
538 * Mac_third_principal_strain = (Principal_strains[2] > 0) ? Principal_strains[2] : 0;
540 * tr_sqr_Mac_Principal_strain =
pow (Mac_first_principal_strain, 2)
541 * +
pow (Mac_second_principal_strain, 2)
542 * +
pow (Mac_third_principal_strain, 2);
545 * H_plus_val = 0.5 *
lambda (E, nu) *
pow (Mac_tr_strain, 2)
546 * + mu (E, nu) * tr_sqr_Mac_Principal_strain;
552 * PhaseField::setup_mesh_and_bcs ()
555 *
const unsigned int nx = 20;
556 *
const unsigned int ny = 20;
557 *
const unsigned int nz = 10;
558 *
const std::vector<unsigned int> repetitions = {nx,ny,nz};
560 *
const Point<3> p1(Domain::x_min,Domain::y_min,Domain::z_min), p2(Domain::x_max,Domain::y_max,Domain::z_max);
567 * 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
573 *
for (
const auto &cell :
triangulation.active_cell_iterators ())
575 * for (const auto &face : cell->face_iterators ())
577 * if (face->at_boundary ())
579 * const auto center = face->center ();
580 *
if (std::fabs (center (0) - (Domain::x_min)) < 1e-12)
581 * face->set_boundary_id (0);
583 *
else if (std::fabs (center (0) - Domain::x_max) < 1e-12)
584 * face->set_boundary_id (1);
586 *
else if (std::fabs (center (1) - (Domain::y_min)) < 1e-12)
587 * face->set_boundary_id (2);
589 *
else if (std::fabs (center (1) - Domain::y_max) < 1e-12)
590 * face->set_boundary_id (3);
595 * pcout <<
"No. of levels in triangulation: "
598 * dof_handler_damage.distribute_dofs (fe_damage);
599 * dof_handler_elastic.distribute_dofs (fe_elastic);
601 * pcout <<
" Number of locally owned cells on the process: "
602 * <<
triangulation.n_locally_owned_active_cells () << std::endl;
604 * pcout <<
"Number of global cells:" <<
triangulation.n_global_active_cells ()
607 * pcout <<
" Total Number of globally active cells: "
609 * <<
" Number of degrees of freedom for elasticity: "
610 * << dof_handler_elastic.n_dofs () << std::endl
611 * <<
" Number of degrees of freedom for damage: "
612 * << dof_handler_damage.n_dofs () << std::endl;
616 * Initialising damage vectors
619 * locally_owned_dofs_damage = dof_handler_damage.locally_owned_dofs ();
621 * dof_handler_damage);
623 * completely_distributed_solution_damage_old.reinit (
624 * locally_owned_dofs_damage, mpi_communicator);
625 * locally_relevant_solution_damage.reinit (locally_owned_dofs_damage,
626 * locally_relevant_dofs_damage, mpi_communicator);
628 * locally_owned_dofs_elastic = dof_handler_elastic.locally_owned_dofs ();
630 * dof_handler_elastic);
632 * completely_distributed_solution_elastic_old.reinit (
633 * locally_owned_dofs_elastic, mpi_communicator);
635 *
for (
const auto &cell :
triangulation.active_cell_iterators ())
637 * if (cell->is_locally_owned ())
639 * quadrature_point_history_field.initialize (cell, 8);
643 *
FEValues<3> fe_values_damage (fe_damage, quadrature_formula_damage,
647 *
for (
const auto &cell :
triangulation.active_cell_iterators ())
648 * if (cell->is_locally_owned ())
650 * const
std::vector<
std::shared_ptr<MyQData>> lqph =
651 * quadrature_point_history_field.get_data (cell);
652 *
for (
const unsigned int q_index : fe_values_damage.quadrature_point_indices ())
654 * lqph[q_index]->value_H = 0.0;
655 * lqph[q_index]->value_H_new = 0.0;
661 * PhaseField::setup_constraints_elastic (
const unsigned int load_step)
663 * constraints_elastic.clear ();
664 * constraints_elastic.reinit (locally_relevant_dofs_elastic);
667 * constraints_elastic);
669 *
for (
const auto &cell : dof_handler_elastic.active_cell_iterators ())
670 * if (cell->is_locally_owned ())
672 * for (const auto &face : cell->face_iterators ())
674 * if (face->at_boundary ())
676 * const auto center = face->center ();
677 *
if (std::fabs (center (0) - Domain::x_min) < 1e-12)
682 * const auto vert = cell->vertex (vertex_number);
683 *
const double z_mid = 0.5 * (Domain::z_max + Domain::z_min);
684 *
if (std::fabs (vert (2) - z_mid) < 1
e-12 && std::fabs (
690 *
const unsigned int z_dof =
691 * cell->vertex_dof_index (vertex_number, 2);
692 * constraints_elastic.add_line (z_dof);
693 * constraints_elastic.set_inhomogeneity (z_dof, 0);
694 *
const unsigned int x_dof =
695 * cell->vertex_dof_index (vertex_number, 0);
696 * constraints_elastic.add_line (x_dof);
697 * constraints_elastic.set_inhomogeneity (x_dof, 0);
699 *
else if (std::fabs (vert (0) - Domain::x_min) < 1
e-12)
702 *
const unsigned int x_dof =
703 * cell->vertex_dof_index (vertex_number, 0);
704 * constraints_elastic.add_line (x_dof);
705 * constraints_elastic.set_inhomogeneity (x_dof, 0);
719 *
const double u_x_values_right = ux * load_step;
720 *
const double u_y_values = uy * load_step;
721 *
const double u_fix = 0.0;
725 * constraints_elastic, u_x_mask);
735 * constraints_elastic.close ();
739 * PhaseField::setup_system_elastic (
const unsigned int load_step)
743 * locally_owned_dofs_elastic = dof_handler_elastic.locally_owned_dofs ();
745 * dof_handler_elastic);
747 * locally_relevant_solution_elastic.reinit (locally_owned_dofs_elastic,
748 * locally_relevant_dofs_elastic, mpi_communicator);
750 * system_rhs_elastic.reinit (locally_owned_dofs_elastic, mpi_communicator);
752 * completely_distributed_solution_elastic.reinit (locally_owned_dofs_elastic,
755 * setup_constraints_elastic (load_step);
760 * constraints_elastic,
false);
763 * dof_handler_elastic.locally_owned_dofs (), mpi_communicator,
764 * locally_relevant_dofs_elastic);
766 * system_matrix_elastic.reinit (locally_owned_dofs_elastic,
767 * locally_owned_dofs_elastic, dsp, mpi_communicator);
771 * PhaseField::assemble_system_elastic ()
776 *
FEValues<3> fe_values_elastic (fe_elastic, quadrature_formula_elastic,
779 *
FEValues<3> fe_values_damage (fe_damage, quadrature_formula_elastic,
782 *
const unsigned int dofs_per_cell = fe_elastic.n_dofs_per_cell ();
783 *
const unsigned int n_q_points = quadrature_formula_elastic.size ();
788 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
790 * std::vector<double> damage_values (n_q_points);
791 * std::vector<Tensor<1, 3>> rhs_values_elastic (n_q_points);
793 *
for (
const auto &cell : dof_handler_elastic.active_cell_iterators ())
794 * if (cell->is_locally_owned ())
797 * cell_matrix_elastic = 0;
798 * cell_rhs_elastic = 0;
802 * dof_handler_damage);
804 * fe_values_damage.reinit (damage_cell);
805 * fe_values_elastic.reinit (cell);
807 * fe_values_damage.get_function_values(locally_relevant_solution_damage,
809 * right_hand_side_elastic (fe_values_elastic.get_quadrature_points (),
810 * rhs_values_elastic);
812 *
for (
const unsigned int q_point : fe_values_elastic.quadrature_point_indices ())
814 * const double
d = damage_values[q_point];
816 *
for (
const unsigned int i : fe_values_elastic.dof_indices ())
819 * const unsigned
int component_i =
820 * fe_elastic.system_to_component_index (i).
first;
822 *
for (
const unsigned int j : fe_values_elastic.dof_indices ())
824 * const unsigned
int component_j =
825 * fe_elastic.system_to_component_index (j).
first;
828 * cell_matrix_elastic (i, j) +=
829 *
pow ((1 - d), 2) * (
830 * (fe_values_elastic.shape_grad (i, q_point)[component_i] *
831 * fe_values_elastic.shape_grad (j, q_point)[component_j]
835 * (fe_values_elastic.shape_grad (i, q_point)[component_j] *
836 * fe_values_elastic.shape_grad (j, q_point)[component_i]
840 * ((component_i == component_j) ?
841 * (fe_values_elastic.shape_grad (i, q_point) *
842 * fe_values_elastic.shape_grad (j, q_point)
848 * fe_values_elastic.JxW (q_point);
854 *
for (
const unsigned int i : fe_values_elastic.dof_indices ())
856 * const unsigned
int component_i =
857 * fe_elastic.system_to_component_index (i).
first;
859 *
for (
const unsigned int q_point : fe_values_elastic.quadrature_point_indices ())
860 * cell_rhs_elastic (i) +=
861 * fe_values_elastic.shape_value (i, q_point) * rhs_values_elastic[q_point][component_i]
862 * * fe_values_elastic.JxW (q_point);
867 * traction contribution to rhs
897 * cell->get_dof_indices (local_dof_indices);
899 * constraints_elastic.distribute_local_to_global (cell_matrix_elastic,
900 * cell_rhs_elastic, local_dof_indices, system_matrix_elastic,
901 * system_rhs_elastic);
909 * PhaseField::solve_linear_system_elastic ()
914 * 1e-12* system_rhs_elastic.l2_norm());
918 * LA::MPI::PreconditionAMG::AdditionalData
data;
919 * #ifdef USE_PETSC_LA
920 *
data.symmetric_operator =
true;
924 * Trilinos defaults are good
928 * LA::MPI::PreconditionAMG preconditioner;
929 * preconditioner.initialize (system_matrix_elastic,
data);
931 * solver.solve (system_matrix_elastic,
932 * completely_distributed_solution_elastic, system_rhs_elastic,
935 * pcout <<
" Solved in " << solver_control.last_step () <<
" iterations."
938 * constraints_elastic.distribute (completely_distributed_solution_elastic);
940 * locally_relevant_solution_elastic = completely_distributed_solution_elastic;
944 * PhaseField::setup_boundary_values_damage ()
948 * constraints_damage.clear ();
949 * constraints_damage.reinit (locally_relevant_dofs_damage);
951 * constraints_damage);
991 * constraints_damage.close ();
996 * PhaseField::setup_system_damage ()
1000 * locally_owned_dofs_damage = dof_handler_damage.locally_owned_dofs ();
1002 * dof_handler_damage);
1004 * locally_relevant_solution_damage.reinit (locally_owned_dofs_damage,
1005 * locally_relevant_dofs_damage, mpi_communicator);
1007 * system_rhs_damage.reinit (locally_owned_dofs_damage, mpi_communicator);
1009 * completely_distributed_solution_damage.reinit (locally_owned_dofs_damage,
1010 * mpi_communicator);
1015 * constraints_damage,
false);
1017 * dof_handler_damage.locally_owned_dofs (), mpi_communicator,
1018 * locally_relevant_dofs_damage);
1020 * system_matrix_damage.reinit (locally_owned_dofs_damage,
1021 * locally_owned_dofs_damage, dsp, mpi_communicator);
1025 * PhaseField::assemble_system_damage ()
1030 *
FEValues<3> fe_values_damage (fe_damage, quadrature_formula_damage,
1033 *
FEValues<3> fe_values_elastic (fe_elastic, quadrature_formula_damage,
1037 *
const unsigned int dofs_per_cell = fe_damage.n_dofs_per_cell ();
1038 *
const unsigned int n_q_points = quadrature_formula_damage.size ();
1043 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1045 *
const RandomMedium::EnergyReleaseRate energy_release_rate;
1046 * std::vector<double> energy_release_rate_values (n_q_points);
1050 * Storing strain tensor
for all Gauss points of a cell in vector strain_values.
1053 * std::vector<SymmetricTensor<2, 3>> strain_values (n_q_points);
1055 *
for (
const auto &cell : dof_handler_damage.active_cell_iterators ())
1056 * if (cell->is_locally_owned ())
1058 *
std::vector<
std::shared_ptr<MyQData>> qpdH =
1059 * quadrature_point_history_field.get_data (cell);
1063 * dof_handler_elastic);
1065 * fe_values_damage.reinit (cell);
1066 * fe_values_elastic.reinit (elastic_cell);
1070 * fe_values_elastic[displacements].get_function_symmetric_gradients (
1071 * locally_relevant_solution_elastic, strain_values);
1073 * cell_matrix_damage = 0;
1074 * cell_rhs_damage = 0;
1076 * energy_release_rate.value_list (
1077 * fe_values_damage.get_quadrature_points (),
1078 * energy_release_rate_values);
1080 *
for (
const unsigned int q_index : fe_values_damage.quadrature_point_indices ())
1082 * const auto &x_q = fe_values_damage.quadrature_point (q_index);
1083 *
double H_call = H_plus (strain_values[q_index]);
1085 *
const double H =
std::max(H_call, qpdH[q_index]->value_H);
1086 * qpdH[q_index]->value_H_new = H;
1087 *
Point<3> cell_center = cell->center ();
1090 *
if (is_in_middle_layer (cell_center,Domain::z1,Domain::z2))
1091 * g_c = energy_release_rate_values[q_index];
1093 * g_c = gc (GC, beta, Domain::z1, Domain::z2, x_q);
1095 *
for (
const unsigned int i : fe_values_damage.dof_indices ())
1097 * for (const unsigned
int j : fe_values_damage.dof_indices ())
1100 * cell_matrix_damage (i, j) +=
1103 * contribution to stiffness from -laplace u term
1106 * Conductivity_damage (x_q) * fe_values_damage.shape_grad (
1109 * fe_values_damage.shape_grad (j, q_index)
1111 * fe_values_damage.JxW (q_index)
1115 * Contribution to stiffness from u term
1121 * ((1 + (2 *
l * H) / g_c) * (1 /
pow (
l, 2))
1122 * * fe_values_damage.shape_value (i, q_index) *
1123 * fe_values_damage.shape_value (j, q_index) *
1124 * fe_values_damage.JxW (q_index));
1128 *
const auto &x_q = fe_values_damage.quadrature_point(q_index);
1131 * cell_rhs_damage (i) += (fe_values_damage.shape_value (i,
1135 * * H * fe_values_damage.JxW (q_index));
1139 * cell->get_dof_indices (local_dof_indices);
1140 * constraints_damage.distribute_local_to_global (cell_matrix_damage,
1141 * cell_rhs_damage, local_dof_indices, system_matrix_damage,
1142 * system_rhs_damage);
1150 * PhaseField::solve_linear_system_damage ()
1155 * 1e-12* system_rhs_damage.l2_norm());
1159 * LA::MPI::PreconditionAMG::AdditionalData
data;
1160 * #ifdef USE_PETSC_LA
1161 *
data.symmetric_operator =
true;
1165 * Trilinos defaults are good
1169 * LA::MPI::PreconditionAMG preconditioner;
1170 * preconditioner.initialize (system_matrix_damage,
data);
1172 * solver.solve (system_matrix_damage, completely_distributed_solution_damage,
1173 * system_rhs_damage, preconditioner);
1175 * pcout <<
" Solved in " << solver_control.last_step () <<
" iterations."
1178 * constraints_damage.distribute (completely_distributed_solution_damage);
1179 * locally_relevant_solution_damage = completely_distributed_solution_damage;
1183 * PhaseField::solve_elastic_subproblem (
const unsigned int load_step)
1185 * setup_system_elastic (load_step);
1186 * assemble_system_elastic ();
1187 * solve_linear_system_elastic ();
1191 * PhaseField::output_results (
const unsigned int load_step)
const
1194 * std::vector<std::string> displacement_names (3,
"displacement");
1195 * std::vector<DataComponentInterpretation::DataComponentInterpretation> displacement_component_interpretation (
1200 * locally_relevant_solution_elastic, displacement_names,
1201 * displacement_component_interpretation);
1202 * data_out_phasefield.add_data_vector (dof_handler_damage,
1203 * locally_relevant_solution_damage,
"damage");
1206 *
for (
unsigned int i = 0; i < subdomain.size (); ++i)
1208 * data_out_phasefield.add_data_vector (subdomain,
"subdomain");
1209 * data_out_phasefield.build_patches ();
1210 * data_out_phasefield.write_vtu_with_pvtu_record (
"./",
"solution", load_step,
1211 * mpi_communicator, 2, 0);
1215 * PhaseField::load_disp_calculation (
const unsigned int load_step)
1221 *
const QGauss<2> face_quadrature (fe_elastic.degree);
1224 * std::vector<SymmetricTensor<2, 3>> strain_values (face_quadrature.size ());
1228 * std::vector<double> damage_values (fe_face_values_damage.n_quadrature_points);
1230 *
for (
const auto &cell : dof_handler_elastic.active_cell_iterators ())
1231 * if (cell->is_locally_owned ())
1233 * const
DoFHandler<3>::active_cell_iterator damage_cell =
1234 *
Triangulation<3>::active_cell_iterator (cell)->as_dof_handler_iterator (
1235 * dof_handler_damage);
1236 *
for (
unsigned int f : cell->face_indices ())
1237 * if (cell->face (f)->at_boundary () && (cell->face (f)->
boundary_id ()
1240 * fe_face_values.
reinit (cell, f);
1241 * fe_face_values[displacements].get_function_symmetric_gradients (
1242 * locally_relevant_solution_elastic, strain_values);
1243 * fe_face_values_damage.reinit(damage_cell, f);
1245 * fe_face_values_damage.get_function_values (locally_relevant_solution_damage, damage_values);
1247 *
for (
unsigned int q = 0; q < fe_face_values.n_quadrature_points; ++q)
1251 *
const double tr_strain = strain[0][0] + strain[1][1] + strain[2][2];
1252 *
const double d = damage_values[q];
1255 * stress[0][0] =
pow ((1 - d), 2)
1256 * * (
lambda (E, nu) * tr_strain + 2 * mu (E, nu)
1258 * stress[0][1] =
pow ((1 - d), 2)
1259 * * (2 * mu (E, nu) * strain[0][1]);
1260 * stress[0][2] =
pow ((1 - d), 2)
1261 * * (2 * mu (E, nu) * strain[0][2]);
1262 * stress[1][1] =
pow ((1 - d), 2)
1263 * * (
lambda (E, nu) * tr_strain + 2 * mu (E, nu)
1265 * stress[1][2] =
pow ((1 - d), 2)
1266 * * (2 * mu (E, nu) * strain[1][2]);
1267 * stress[2][2] =
pow ((1 - d), 2)
1268 * * (
lambda (E, nu) * tr_strain + 2 * mu (E, nu)
1272 * * fe_face_values.normal_vector (q);
1273 * x_max_force += force_density * fe_face_values.JxW (q);
1277 *
else if (cell->face (f)->at_boundary ()
1278 * && (cell->face (f)->boundary_id () == 3))
1280 * fe_face_values.reinit (cell, f);
1281 * fe_face_values_damage.reinit (damage_cell, f);
1282 * fe_face_values[displacements].get_function_symmetric_gradients (
1283 * locally_relevant_solution_elastic, strain_values);
1284 * fe_face_values_damage.get_function_values (locally_relevant_solution_damage, damage_values);
1286 *
for (
unsigned int q = 0; q < fe_face_values.n_quadrature_points;
1291 *
const double tr_strain = strain[0][0] + strain[1][1] + strain[2][2];
1292 *
const double d = damage_values[q];
1295 * stress[0][0] =
pow ((1 - d), 2)
1296 * * (
lambda (E, nu) * tr_strain + 2 * mu (E, nu)
1298 * stress[0][1] =
pow ((1 - d), 2)
1299 * * (2 * mu (E, nu) * strain[0][1]);
1300 * stress[0][2] =
pow ((1 - d), 2)
1301 * * (2 * mu (E, nu) * strain[0][2]);
1302 * stress[1][1] =
pow ((1 - d), 2)
1303 * * (
lambda (E, nu) * tr_strain + 2 * mu (E, nu)
1305 * stress[1][2] =
pow ((1 - d), 2)
1306 * * (2 * mu (E, nu) * strain[1][2]);
1307 * stress[2][2] =
pow ((1 - d), 2)
1308 * * (
lambda (E, nu) * tr_strain + 2 * mu (E, nu)
1312 * * fe_face_values.normal_vector (q);
1313 * y_max_force += force_density * fe_face_values.JxW (q);
1317 *
double x_max_force_x;
1318 * x_max_force_x = x_max_force[0];
1320 * pcout <<
"fx: " << x_max_force_x << std::endl;
1322 *
double y_max_force_y;
1323 * y_max_force_y = y_max_force[1];
1325 * pcout <<
"fy: " << y_max_force_y << std::endl;
1331 * Code to be executed only on process 0
1334 * load_values_x[load_step] = x_max_force_x;
1335 * load_values_y[load_step] = y_max_force_y;
1337 *
const double disp = ux * load_step;
1338 * displacement_values[load_step] = disp;
1342 * load displacement plot
1345 * std::ofstream m_file;
1346 * m_file.open (
"load_displacement.m");
1347 * m_file <<
"% Matlab code generated by dealii to plot load displacement curves"
1349 * m_file <<
"clc" << std::endl;
1350 * m_file <<
"clear" << std::endl;
1351 * m_file <<
"close all" << std::endl;
1352 * m_file <<
"load_x=[" << load_values_x <<
"]" << std::endl;
1353 * m_file <<
"load_y=[" << load_values_y <<
"]" << std::endl;
1354 * m_file <<
"displacement=[" << displacement_values <<
"]" << std::endl;
1355 * m_file <<
"plot( displacement,load_x,'linewidth',2)" << std::endl;
1356 * m_file <<
"xlabel(\"Displacement (mm)\",'Interpreter', 'latex')"
1358 * m_file <<
"ylabel(\"Reaction force (kN)\",'Interpreter', 'latex')"
1360 * m_file <<
"hold on" << std::endl;
1361 * m_file <<
"plot( displacement,load_y,'linewidth',2)" << std::endl;
1362 * m_file <<
"set(gca,'fontname', 'Courier','FontSize',15,'FontWeight','bold','linewidth',1); grid on"
1364 * m_file <<
"h=legend(\"fx \",\"fy\")" << std::endl;
1365 * m_file <<
"set(h,'FontSize',14);" << std::endl;
1366 * m_file <<
"set(gcf,'Position', [250 250 700 500])" << std::endl;
1367 * m_file <<
"xlim([0,0.065])" << std::endl;
1368 * m_file <<
"box on" << std::endl;
1374 * PhaseField::solve_damage_subproblem ()
1376 * setup_boundary_values_damage ();
1377 * setup_system_damage ();
1378 * assemble_system_damage ();
1379 * solve_linear_system_damage ();
1383 * PhaseField::refine_grid (
const unsigned int load_step)
1385 *
FEValues<3> fe_values_damage (fe_damage, quadrature_formula_damage,
1390 * fe_damage, quadrature_formula_damage, quadrature_formula_damage);
1394 * 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,
1395 * it is usually not necessary to compute the error indicator corresponding to the solution on a cell particularly accurately.
1396 * We therefore typically use a vector of floats instead of a vector of doubles to represent error indicators.
1404 * { }, locally_relevant_solution_damage, estimated_error_per_cell);
1426 *
for (
const auto &cell :
triangulation.active_cell_iterators_on_level (3))
1427 * if (cell->is_locally_owned ())
1428 * cell->clear_refine_flag ();
1443 * data_transfer.prepare_for_coarsening_and_refinement (
triangulation,
1444 * quadrature_point_history_field);
1449 * and give the solution vector that we intend to
interpolate later,
1452 * soltransDamage.prepare_for_coarsening_and_refinement (
1453 * locally_relevant_solution_damage);
1454 * soltransElastic.prepare_for_coarsening_and_refinement (
1455 * locally_relevant_solution_elastic);
1461 * redistribute dofs,
1464 * dof_handler_damage.distribute_dofs (fe_damage);
1465 * dof_handler_elastic.distribute_dofs (fe_elastic);
1469 * Recreate locally_owned_dofs and locally_relevant_dofs
index sets
1472 * locally_owned_dofs_damage = dof_handler_damage.locally_owned_dofs ();
1474 * dof_handler_damage);
1476 * completely_distributed_solution_damage_old.reinit (
1477 * locally_owned_dofs_damage, mpi_communicator);
1478 * soltransDamage.interpolate (completely_distributed_solution_damage_old);
1482 * Apply constraints on the interpolated solution to make sure it conforms with the
new mesh
1485 * setup_boundary_values_damage ();
1487 * constraints_damage.distribute (completely_distributed_solution_damage_old);
1491 * Copy completely_distributed_solution_damage_old to locally_relevant_solution_damage
1494 * locally_relevant_solution_damage.reinit (locally_owned_dofs_damage,
1495 * locally_relevant_dofs_damage, mpi_communicator);
1496 * locally_relevant_solution_damage =
1497 * completely_distributed_solution_damage_old;
1501 * Interpolating elastic solution similarly
1504 * locally_owned_dofs_elastic = dof_handler_elastic.locally_owned_dofs ();
1506 * dof_handler_elastic);
1508 * completely_distributed_solution_elastic_old.reinit (
1509 * locally_owned_dofs_elastic, mpi_communicator);
1510 * soltransElastic.interpolate (completely_distributed_solution_elastic_old);
1514 * Apply constraints on the interpolated solution to make sure it conforms with the
new mesh
1517 * setup_constraints_elastic (load_step);
1518 * constraints_elastic.distribute (
1519 * completely_distributed_solution_elastic_old);
1523 * Copy completely_distributed_solution_damage_old to locally_relevant_solution_damage
1526 * locally_relevant_solution_elastic.reinit (locally_owned_dofs_elastic,
1527 * locally_relevant_dofs_elastic, mpi_communicator);
1528 * locally_relevant_solution_elastic =
1529 * completely_distributed_solution_elastic_old;
1531 *
for (
const auto &cell :
triangulation.active_cell_iterators ())
1533 * if (cell->is_locally_owned ())
1534 * quadrature_point_history_field.initialize (cell, 8);
1536 * data_transfer.interpolate ();
1541 * PhaseField::check_convergence ()
1543 * LA::MPI::Vector solution_damage_difference (locally_owned_dofs_damage,
1544 * mpi_communicator);
1545 * LA::MPI::Vector solution_elastic_difference (locally_owned_dofs_elastic,
1546 * mpi_communicator);
1547 * LA::MPI::Vector solution_damage_difference_ghost (locally_owned_dofs_damage,
1548 * locally_relevant_dofs_damage, mpi_communicator);
1549 * LA::MPI::Vector solution_elastic_difference_ghost (
1550 * locally_owned_dofs_elastic, locally_relevant_dofs_elastic,
1551 * mpi_communicator);
1553 * solution_damage_difference = locally_relevant_solution_damage;
1555 * solution_damage_difference -= completely_distributed_solution_damage_old;
1557 * solution_elastic_difference = locally_relevant_solution_elastic;
1559 * solution_elastic_difference -= completely_distributed_solution_elastic_old;
1561 * solution_damage_difference_ghost = solution_damage_difference;
1563 * solution_elastic_difference_ghost = solution_elastic_difference;
1565 *
double error_elastic_solution_numerator, error_elastic_solution_denominator,
1566 * error_damage_solution_numerator, error_damage_solution_denominator;
1568 * error_damage_solution_numerator = solution_damage_difference.l2_norm ();
1569 * error_elastic_solution_numerator = solution_elastic_difference.l2_norm ();
1570 * error_damage_solution_denominator =
1571 * completely_distributed_solution_damage.l2_norm ();
1572 * error_elastic_solution_denominator =
1573 * completely_distributed_solution_elastic.l2_norm ();
1575 *
double error_elastic_solution, error_damage_solution;
1576 * error_damage_solution = error_damage_solution_numerator
1577 * / error_damage_solution_denominator;
1578 * error_elastic_solution = error_elastic_solution_numerator
1579 * / error_elastic_solution_denominator;
1581 *
if ((error_elastic_solution < tol) && (error_damage_solution < tol))
1588 * PhaseField::update_history_field ()
1590 *
FEValues<3> fe_values_damage (fe_damage, quadrature_formula_damage,
1594 *
for (
const auto &cell : dof_handler_damage.active_cell_iterators ())
1595 * if (cell->is_locally_owned ())
1597 * const
std::vector<
std::shared_ptr<MyQData>> lqph =
1598 * quadrature_point_history_field.get_data (cell);
1599 *
for (
unsigned int q_index = 0;
1600 * q_index < quadrature_formula_damage.size (); ++q_index)
1601 * lqph[q_index]->value_H = lqph[q_index]->value_H_new;
1606 * PhaseField::run ()
1610 * pcout <<
"Running with "
1611 * #ifdef USE_PETSC_LA
1620 * setup_mesh_and_bcs ();
1624 * Create the pre-crack in the domain
1627 * LA::MPI::Vector initial_soln_damage (locally_owned_dofs_damage,
1628 * mpi_communicator);
1629 *
for (
const auto &cell : dof_handler_damage.active_cell_iterators ())
1630 * if (cell->is_locally_owned ())
1634 * cell->vertex_dof_index (vertex_number, 0);
1635 * initial_soln_damage[vertex_dof_index] = 0;
1638 * locally_relevant_solution_damage = initial_soln_damage;
1642 * Loop over load steps
1645 *
for (
unsigned int load_step = 1; load_step <= num_load_steps; load_step++)
1647 * pcout <<
" \n \n load increment number : " << load_step << std::endl;
1651 * Loop over staggered iterations
1654 *
unsigned int iteration = 0;
1655 *
bool stoppingCriterion =
false;
1656 *
while (stoppingCriterion ==
false)
1658 * pcout <<
" \n iteration number:" << iteration << std::endl;
1659 * solve_elastic_subproblem (load_step);
1660 * solve_damage_subproblem ();
1662 * locally_relevant_solution_damage.update_ghost_values ();
1663 * locally_relevant_solution_elastic.update_ghost_values ();
1665 *
if (iteration > 0)
1666 * stoppingCriterion = check_convergence ();
1668 * completely_distributed_solution_elastic_old =
1669 * locally_relevant_solution_elastic;
1670 * completely_distributed_solution_damage_old =
1671 * locally_relevant_solution_damage;
1673 *
if (stoppingCriterion ==
false)
1675 * refine_grid (load_step);
1677 * iteration = iteration + 1;
1682 * Once converged,
do some clean-up operations
1685 *
if ((load_step == 1) || (load_step >= 1 && load_step <= num_load_steps
1686 * && std::fabs (load_step % 10) < 1e-6))
1689 * output_results (load_step);
1692 * load_disp_calculation (load_step);
1694 * computing_timer.print_summary ();
1695 * computing_timer.reset ();
1696 * pcout << std::endl;
1698 * update_history_field ();
1702 * pcout <<
"Total run time: " << timer.wall_time () <<
" seconds."
1713 *
using namespace dealii;
1714 *
using namespace FracturePropagation;
1718 * PhaseField phasefield;
1719 * phasefield.run ();
1721 *
catch (std::exception &exc)
1723 * std::cerr << std::endl << std::endl
1724 * <<
"----------------------------------------------------" << std::endl;
1725 * std::cerr <<
"Exception on processing: " << std::endl << exc.what ()
1726 * << std::endl <<
"Aborting!" << std::endl
1727 * <<
"----------------------------------------------------" << std::endl;
1733 * std::cerr << std::endl << std::endl
1734 * <<
"----------------------------------------------------" << std::endl;
1735 * std::cerr <<
"Unknown exception!" << std::endl <<
"Aborting!" << std::endl
1736 * <<
"----------------------------------------------------" << std::endl;
std::vector< bool > component_mask
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={})
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< RangeNumberType > &values, const unsigned int component=0) const
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)
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)
#define AssertDimension(dim1, dim2)
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_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
@ component_is_part_of_vector
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)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
constexpr types::blas_int zero
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
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)
::SolutionTransfer< dim, VectorType, spacedim > SolutionTransfer
::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 > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
unsigned int global_dof_index
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)