Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-3124-g7c413cc06e 2025-04-23 12:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
Phase_field_fracture_model_in_3D.h
Go to the documentation of this file.
1
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>
121 *  
122 *   #define FORCE_USE_OF_TRILINOS
123 *   namespace LA
124 *   {
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;
131 *   #else
132 *   # error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
133 *   #endif
134 *   } // namespace LA
135 *  
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>
163 *   #include <fstream>
164 *   #include <iostream>
165 *   #include <random>
166 *  
167 *   namespace FracturePropagation
168 *  
169 *   {
170 *   using namespace dealii;
171 *  
172 *   class PhaseField
173 *   {
174 * @endcode
175 *
176 * Function declarations
177 *
178 * @code
179 *   public:
180 *   PhaseField ();
181 *   void
182 *   run ();
183 *  
184 *   private:
185 *  
186 * @endcode
187 *
188 * 1. Mesh and boundary conditions setup at the beginning
189 *
190 * @code
191 *   void
192 *   setup_mesh_and_bcs (); // Called at the start to create a mesh
193 *  
194 * @endcode
195 *
196 * 2. Elastic subproblem setup and solution
197 *
198 * @code
199 *   void
200 *   setup_constraints_elastic (const unsigned int load_step);
201 *   void
202 *   setup_system_elastic (const unsigned int load_step);
203 *   void
204 *   assemble_system_elastic ();
205 *   void
206 *   solve_linear_system_elastic ();
207 *   void
208 *   solve_elastic_subproblem (const unsigned int load_step); // Calls the above functions
209 *  
210 * @endcode
211 *
212 * 3. Damage subproblem setup and solution
213 *
214 * @code
215 *   void
216 *   setup_boundary_values_damage ();
217 *   void
218 *   setup_system_damage ();
219 *   void
220 *   assemble_system_damage ();
221 *   void
222 *   solve_linear_system_damage ();
223 *   double
224 *   H_plus (const SymmetricTensor<2, 3> &strain); // Called during system assembly for the damage subproblem
225 *   void
226 *   solve_damage_subproblem (); // Calls the above functions
227 *  
228 * @endcode
229 *
230 * 4. Convergence check after each iteration
231 *
232 * @code
233 *   bool
234 *   check_convergence ();
235 *  
236 * @endcode
237 *
238 * 5. Post-processing: output results, refine grid, and calculate displacement
239 *
240 * @code
241 *   void
242 *   output_results (const unsigned int load_step) const; // Called if convergence is achieved
243 *   void
244 *   load_disp_calculation (const unsigned int load_step); // Called if convergence is achieved
245 *   void
246 *   update_history_field (); // Called if convergence is achieved
247 *   void
248 *   refine_grid (const unsigned int load_step); // Called if convergence is not achieved
249 *  
250 *   MPI_Comm mpi_communicator;
251 *   ConditionalOStream pcout;
252 *   TimerOutput computing_timer;
253 *  
255 *  
256 * @endcode
257 *
258 * Objects for elasticity
259 *
260 * @code
261 *   const FESystem<3> fe_elastic;
262 *   DoFHandler<3> dof_handler_elastic;
263 *   IndexSet locally_owned_dofs_elastic;
264 *   IndexSet locally_relevant_dofs_elastic;
265 *   AffineConstraints<double> constraints_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;
272 *  
273 * @endcode
274 *
275 * Objects for damage
276 *
277 * @code
278 *   const FE_Q<3> fe_damage;
279 *   DoFHandler<3> dof_handler_damage;
280 *   IndexSet locally_owned_dofs_damage;
281 *   IndexSet locally_relevant_dofs_damage;
282 *   AffineConstraints<double> constraints_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;
289 *  
290 *   const double ux = 1e-3; // increment in loading
291 *   const double alpha = 1;
292 *   const double uy = alpha * ux;
293 *   const double l = 0.6; // length scale parameter
294 *   const unsigned int num_load_steps = 100; // number of load steps
295 *   const double tol = 1e-2; // tolerance for error in solution
296 *   const double GC = 1e-3; //energy release rate
297 *   const double E = 37.5;
298 *   const double beta = 25;
299 *   const double nu = 0.25;
300 *  
301 * @endcode
302 *
303 * Objects for load-displacement calculation
304 *
305 * @code
306 *   Vector<double> load_values_x;
307 *   Vector<double> load_values_y;
308 *   Vector<double> displacement_values;
309 *  
310 *   class MyQData : public TransferableQuadraturePointData
311 *   {
312 *   public:
313 *   MyQData () = default;
314 *   virtual
315 *   ~MyQData () = default;
316 *  
317 *   unsigned int
318 *   number_of_values () const override
319 *   {
320 *   return 2; // Indicate that there are two values to handle
321 *   }
322 *  
323 *   virtual void
324 *   pack_values (std::vector<double> &scalars) const override
325 *   {
326 *   Assert (scalars.size () == 2, ExcInternalError ()); // Ensure the vector has exactly two elements
327 *   scalars[0] = value_H; // Pack the first value
328 *   scalars[1] = value_H_new; // Pack the second value
329 *   }
330 *  
331 *   virtual void
332 *   unpack_values (const std::vector<double> &scalars) override
333 *   {
334 *   Assert (scalars.size () == 2, ExcInternalError ()); // Ensure the vector has exactly two elements
335 *   value_H = scalars[0]; // Unpack the first value
336 *   value_H_new = scalars[1]; // Unpack the second value
337 *   }
338 *  
339 *   double value_H; // First value
340 *   double value_H_new; // Second value
341 *   };
342 *  
343 *   CellDataStorage<typename Triangulation<3>::cell_iterator, MyQData> quadrature_point_history_field;
344 *  
345 *   };
346 *  
347 *   inline double
348 *   lambda (const float E,
349 *   const float nu)
350 *   {
351 *   return (E * nu) / ((1 + nu) * (1 - 2 * nu));
352 *   }
353 *  
354 *   inline double
355 *   mu (const float E,
356 *   const float nu)
357 *   {
358 *   return E / (2 * (1 + nu));
359 *   }
360 *  
361 *   inline double
362 *   gc (const float GC,
363 *   const float beta,
364 *   const float z1,
365 *   const float z2,
366 *   const Point<3> &p)
367 *   {
368 *   if (((p[2] - z1) > 1e-6) && ((p[2] - z2) < 1e-6))
369 *   return GC / beta;
370 *   else
371 *   return GC;
372 *   }
373 *  
374 *   bool
375 *   is_in_middle_layer (const Point<3> &cell_center,const float z1,const float z2)
376 *   {
377 *   return (cell_center[2] >= z1 && cell_center[2] <= z2);
378 *   }
379 *  
380 *   namespace Domain
381 *  
382 *   {
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;
391 *   }
392 *  
393 *   namespace RandomMedium
394 *   {
395 *   class EnergyReleaseRate : public Function<3>
396 *   {
397 *   public:
398 *   EnergyReleaseRate ()
399 *   :
400 *   Function<3> ()
401 *   {
402 *   }
403 *  
404 * @endcode
405 *
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.
407 *
408 * @code
409 *   virtual void
410 *   value_list (const std::vector<Point<3>> &points,
411 *   std::vector<double> &values,
412 *   const unsigned int /*component*/= 0) const override
413 *   {
414 *   AssertDimension (points.size (), values.size ());
415 *  
416 *   for (unsigned int p = 0; p < points.size (); ++p)
417 *   {
418 *   values[p] = 0.0;
419 *  
420 *   double fracture_energy = 0;
421 *   for (unsigned int i = 0; i < centers.size (); ++i)
422 *   fracture_energy += std::exp (
423 *   -(points[p] - centers[i]).norm_square () / (1.5 * 1.5));
424 *  
425 *   const double normalized_fracture_energy = std::min (
426 *   std::max (fracture_energy, 4e-5), 4e-4);
427 *  
428 *   values[p] = normalized_fracture_energy;
429 *   }
430 *   }
431 *  
432 *   private:
433 *   static std::vector<Point<3>> centers;
434 *  
435 *   static std::vector<Point<3>>
436 *   get_centers ()
437 *   {
438 *   const unsigned int N = 1000;
439 *  
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)
444 *   {
445 *   centers_list[i][d] =
446 *   static_cast<double> ((rand ()) / RAND_MAX) * Domain::x_max; //domain_size; //generates a number between 0 and domain_size
447 * @endcode
448 *
449 * x,y will be between 0 to x_max
450 *
451 * @code
452 *   }
453 *   else if (d == 2)
454 *   {
455 *   centers_list[i][d] = static_cast<double> (Domain::z1
456 *   + ((rand ()) / RAND_MAX) * (Domain::z2 - Domain::z1));
457 *   }
458 *  
459 *   return centers_list;
460 *   }
461 *   };
462 *  
463 *   std::vector<Point<3>> EnergyReleaseRate::centers =
464 *   EnergyReleaseRate::get_centers ();
465 *  
466 *   } // namespace RandomMedium
467 *  
468 *  
469 *   inline double
470 *   Conductivity_damage (const Point<3> &p) //const
471 *   {
472 *   return p[0] - p[0] + 1;
473 *   }
474 *  
475 *   void
476 *   right_hand_side_elastic (const std::vector<Point<3>> &points,
477 *   std::vector<Tensor<1, 3>> &values)
478 *   {
479 *   AssertDimension (values.size (), points.size ());
480 *  
481 *   for (unsigned int point_n = 0; point_n < points.size (); ++point_n)
482 *   {
483 *   values[point_n][0] = 0; //x component of body force
484 *   values[point_n][1] = 0; //y component of body force
485 *   values[point_n][2] = 0; //z component of body force
486 *   }
487 *  
488 *   }
489 *   void
490 *   Traction_elastic (const std::vector<Point<3>> &points,
491 *   std::vector<Tensor<1, 3>> &values)
492 *  
493 *   {
494 *   AssertDimension (values.size (), points.size ());
495 *   for (unsigned int point_n = 0; point_n < points.size (); ++point_n)
496 *   {
497 *   values[point_n][0] = 0; //x component of traction
498 *   values[point_n][1] = 0; //y component of traction
499 *   values[point_n][2] = 0; //y component of traction
500 *   }
501 *  
502 *   }
503 *  
504 *   PhaseField::PhaseField ()
505 *   :
506 *   mpi_communicator (MPI_COMM_WORLD),
507 *   pcout (std::cout,
508 *   (Utilities::MPI::this_mpi_process (mpi_communicator) == 0)),
509 *   computing_timer (mpi_communicator, pcout, TimerOutput::never,
511 *   triangulation (mpi_communicator),
512 *   fe_elastic (FE_Q<3> (1), 3),
513 *   dof_handler_elastic (triangulation),
514 *   quadrature_formula_elastic (fe_elastic.degree + 1),
515 *   fe_damage (1),
516 *   dof_handler_damage (triangulation),
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)
521 *   {
522 *   }
523 *  
524 *   double
525 *   PhaseField::H_plus (const SymmetricTensor<2, 3> &strain)
526 *   {
527 *   double Mac_tr_strain, Mac_first_principal_strain,
528 *   Mac_second_principal_strain, Mac_third_principal_strain,
529 *   tr_sqr_Mac_Principal_strain;
530 *  
531 *   const double tr_strain = trace (strain);
532 *  
533 *   Mac_tr_strain = tr_strain >0 ? tr_strain : 0;
534 *   const std::array<double, 3> Principal_strains = eigenvalues (strain);
535 *  
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;
539 *  
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);
543 *  
544 *   double H_plus_val;
545 *   H_plus_val = 0.5 * lambda (E, nu) * pow (Mac_tr_strain, 2)
546 *   + mu (E, nu) * tr_sqr_Mac_Principal_strain;
547 *   return H_plus_val;
548 *   }
549 *  
550 *  
551 *   void
552 *   PhaseField::setup_mesh_and_bcs ()
553 *  
554 *   {
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};
559 *  
560 *   const Point<3> p1(Domain::x_min,Domain::y_min,Domain::z_min), p2(Domain::x_max,Domain::y_max,Domain::z_max);
561 *  
563 *   p2); // create coarse mesh
564 *  
565 * @endcode
566 *
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
568 *
569
570 *
571 *
572 * @code
573 *   for (const auto &cell : triangulation.active_cell_iterators ())
574 *   {
575 *   for (const auto &face : cell->face_iterators ())
576 *   {
577 *   if (face->at_boundary ())
578 *   {
579 *   const auto center = face->center ();
580 *   if (std::fabs (center (0) - (Domain::x_min)) < 1e-12)
581 *   face->set_boundary_id (0);
582 *  
583 *   else if (std::fabs (center (0) - Domain::x_max) < 1e-12)
584 *   face->set_boundary_id (1);
585 *  
586 *   else if (std::fabs (center (1) - (Domain::y_min)) < 1e-12)
587 *   face->set_boundary_id (2);
588 *  
589 *   else if (std::fabs (center (1) - Domain::y_max) < 1e-12)
590 *   face->set_boundary_id (3);
591 *   }
592 *   }
593 *   }
594 *  
595 *   pcout << "No. of levels in triangulation: "
596 *   << triangulation.n_global_levels () << std::endl;
597 *  
598 *   dof_handler_damage.distribute_dofs (fe_damage);
599 *   dof_handler_elastic.distribute_dofs (fe_elastic);
600 *  
601 *   pcout << " Number of locally owned cells on the process: "
602 *   << triangulation.n_locally_owned_active_cells () << std::endl;
603 *  
604 *   pcout << "Number of global cells:" << triangulation.n_global_active_cells ()
605 *   << std::endl;
606 *  
607 *   pcout << " Total Number of globally active cells: "
608 *   << triangulation.n_global_active_cells () << std::endl
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;
613 *  
614 * @endcode
615 *
616 * Initialising damage vectors
617 *
618 * @code
619 *   locally_owned_dofs_damage = dof_handler_damage.locally_owned_dofs ();
620 *   locally_relevant_dofs_damage = DoFTools::extract_locally_relevant_dofs (
621 *   dof_handler_damage);
622 *  
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);
627 *  
628 *   locally_owned_dofs_elastic = dof_handler_elastic.locally_owned_dofs ();
629 *   locally_relevant_dofs_elastic = DoFTools::extract_locally_relevant_dofs (
630 *   dof_handler_elastic);
631 *  
632 *   completely_distributed_solution_elastic_old.reinit (
633 *   locally_owned_dofs_elastic, mpi_communicator);
634 *  
635 *   for (const auto &cell : triangulation.active_cell_iterators ())
636 *   {
637 *   if (cell->is_locally_owned ())
638 *   {
639 *   quadrature_point_history_field.initialize (cell, 8);
640 *   }
641 *   }
642 *  
643 *   FEValues<3> fe_values_damage (fe_damage, quadrature_formula_damage,
646 *  
647 *   for (const auto &cell : triangulation.active_cell_iterators ())
648 *   if (cell->is_locally_owned ())
649 *   {
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 ())
653 *   {
654 *   lqph[q_index]->value_H = 0.0;
655 *   lqph[q_index]->value_H_new = 0.0;
656 *   }
657 *   }
658 *   }
659 *  
660 *   void
661 *   PhaseField::setup_constraints_elastic (const unsigned int load_step)
662 *   {
663 *   constraints_elastic.clear ();
664 *   constraints_elastic.reinit (locally_relevant_dofs_elastic);
665 *  
666 *   DoFTools::make_hanging_node_constraints (dof_handler_elastic,
667 *   constraints_elastic);
668 *  
669 *   for (const auto &cell : dof_handler_elastic.active_cell_iterators ())
670 *   if (cell->is_locally_owned ())
671 *   {
672 *   for (const auto &face : cell->face_iterators ())
673 *   {
674 *   if (face->at_boundary ())
675 *   {
676 *   const auto center = face->center ();
677 *   if (std::fabs (center (0) - Domain::x_min) < 1e-12) //face lies at x=x_min
678 *   {
679 *  
680 *   for (const auto vertex_number : cell->vertex_indices ())
681 *   {
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) < 1e-12 && std::fabs (
685 *   vert (
686 *   0)
687 *   - Domain::x_min)
688 *   < 1e-12) // vertex at x=x_min,z=z_mid;
689 *   {
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);
698 *   }
699 *   else if (std::fabs (vert (0) - Domain::x_min) < 1e-12) // vertex at x_min
700 *  
701 *   {
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);
706 *   }
707 *   }
708 *   }
709 *   }
710 *   }
711 *   }
712 *  
713 *   const FEValuesExtractors::Scalar u_x (0);
714 *   const FEValuesExtractors::Scalar u_y (1);
715 *  
716 *   const ComponentMask u_x_mask = fe_elastic.component_mask (u_x);
717 *   const ComponentMask u_y_mask = fe_elastic.component_mask (u_y);
718 *  
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;
722 *  
723 *   VectorTools::interpolate_boundary_values (dof_handler_elastic, 1,
724 *   Functions::ConstantFunction<3> (u_x_values_right, 3),
725 *   constraints_elastic, u_x_mask);
726 *  
727 *   VectorTools::interpolate_boundary_values (dof_handler_elastic, 2,
728 *   Functions::ConstantFunction<3> (u_fix, 3), constraints_elastic,
729 *   u_y_mask);
730 *  
731 *   VectorTools::interpolate_boundary_values (dof_handler_elastic, 3,
732 *   Functions::ConstantFunction<3> (u_y_values, 3), constraints_elastic,
733 *   u_y_mask);
734 *  
735 *   constraints_elastic.close ();
736 *   }
737 *  
738 *   void
739 *   PhaseField::setup_system_elastic (const unsigned int load_step)
740 *   {
741 *   TimerOutput::Scope ts (computing_timer, "setup_system_elastic");
742 *  
743 *   locally_owned_dofs_elastic = dof_handler_elastic.locally_owned_dofs ();
744 *   locally_relevant_dofs_elastic = DoFTools::extract_locally_relevant_dofs (
745 *   dof_handler_elastic);
746 *  
747 *   locally_relevant_solution_elastic.reinit (locally_owned_dofs_elastic,
748 *   locally_relevant_dofs_elastic, mpi_communicator);
749 *  
750 *   system_rhs_elastic.reinit (locally_owned_dofs_elastic, mpi_communicator);
751 *  
752 *   completely_distributed_solution_elastic.reinit (locally_owned_dofs_elastic,
753 *   mpi_communicator);
754 *  
755 *   setup_constraints_elastic (load_step);
756 *  
757 *   DynamicSparsityPattern dsp (locally_relevant_dofs_elastic);
758 *  
759 *   DoFTools::make_sparsity_pattern (dof_handler_elastic, dsp,
760 *   constraints_elastic, false);
761 *  
763 *   dof_handler_elastic.locally_owned_dofs (), mpi_communicator,
764 *   locally_relevant_dofs_elastic);
765 *  
766 *   system_matrix_elastic.reinit (locally_owned_dofs_elastic,
767 *   locally_owned_dofs_elastic, dsp, mpi_communicator);
768 *   }
769 *  
770 *   void
771 *   PhaseField::assemble_system_elastic ()
772 *  
773 *   {
774 *   TimerOutput::Scope ts (computing_timer, "assembly_elastic");
775 *  
776 *   FEValues<3> fe_values_elastic (fe_elastic, quadrature_formula_elastic,
778 *   | update_JxW_values);
779 *   FEValues<3> fe_values_damage (fe_damage, quadrature_formula_elastic,
780 *   update_values);
781 *  
782 *   const unsigned int dofs_per_cell = fe_elastic.n_dofs_per_cell ();
783 *   const unsigned int n_q_points = quadrature_formula_elastic.size ();
784 *  
785 *   FullMatrix<double> cell_matrix_elastic (dofs_per_cell, dofs_per_cell);
786 *   Vector<double> cell_rhs_elastic (dofs_per_cell);
787 *  
788 *   std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
789 *  
790 *   std::vector<double> damage_values (n_q_points);
791 *   std::vector<Tensor<1, 3>> rhs_values_elastic (n_q_points);
792 *  
793 *   for (const auto &cell : dof_handler_elastic.active_cell_iterators ())
794 *   if (cell->is_locally_owned ())
795 *  
796 *   {
797 *   cell_matrix_elastic = 0;
798 *   cell_rhs_elastic = 0;
799 *  
800 *   const DoFHandler<3>::active_cell_iterator damage_cell =
801 *   Triangulation<3>::active_cell_iterator (cell)->as_dof_handler_iterator (
802 *   dof_handler_damage);
803 *  
804 *   fe_values_damage.reinit (damage_cell);
805 *   fe_values_elastic.reinit (cell);
806 *  
807 *   fe_values_damage.get_function_values(locally_relevant_solution_damage,
808 *   damage_values);
809 *   right_hand_side_elastic (fe_values_elastic.get_quadrature_points (),
810 *   rhs_values_elastic);
811 *  
812 *   for (const unsigned int q_point : fe_values_elastic.quadrature_point_indices ())
813 *   {
814 *   const double d = damage_values[q_point];
815 *  
816 *   for (const unsigned int i : fe_values_elastic.dof_indices ())
817 *  
818 *   {
819 *   const unsigned int component_i =
820 *   fe_elastic.system_to_component_index (i).first;
821 *  
822 *   for (const unsigned int j : fe_values_elastic.dof_indices ())
823 *   {
824 *   const unsigned int component_j =
825 *   fe_elastic.system_to_component_index (j).first;
826 *  
827 *   {
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]
832 *   *
833 *   lambda (E, nu))
834 *   +
835 *   (fe_values_elastic.shape_grad (i, q_point)[component_j] *
836 *   fe_values_elastic.shape_grad (j, q_point)[component_i]
837 *   *
838 *   mu (E, nu))
839 *   +
840 *   ((component_i == component_j) ?
841 *   (fe_values_elastic.shape_grad (i, q_point) *
842 *   fe_values_elastic.shape_grad (j, q_point)
843 *   *
844 *   mu (E, nu)) :
845 *   0)
846 *   )
847 *   *
848 *   fe_values_elastic.JxW (q_point);
849 *   }
850 *   }
851 *   }
852 *   }
853 *  
854 *   for (const unsigned int i : fe_values_elastic.dof_indices ())
855 *   {
856 *   const unsigned int component_i =
857 *   fe_elastic.system_to_component_index (i).first;
858 *  
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);
863 *   }
864 *  
865 * @endcode
866 *
867 * traction contribution to rhs
868 *
869 * @code
870 *   /* for (const auto &face : cell->face_iterators())
871 *   {
872 *   if (face->at_boundary() &&
873 *   face->boundary_id() == 10)
874 *  
875 *   {
876 *   fe_face_values_elastic.reinit(cell, face);
877 *   for (const auto i : fe_face_values_elastic.dof_indices())
878 *  
879 *   {const unsigned int component_i =
880 *   fe_elastic.system_to_component_index(i).first;
881 *   for (const auto face_q_point :
882 *   fe_face_values_elastic.quadrature_point_indices())
883 *   {
884 *  
885 *   cell_rhs_elastic(i) +=
886 *   fe_face_values_elastic.shape_value( i, face_q_point)
887 *   *(traction_values_elastic[face_q_point][component_i])
888 *   *fe_face_values_elastic.JxW(face_q_point);
889 *  
890 *  
891 *  
892 *   }
893 *   }
894 *   }
895 *   }*/
896 *  
897 *   cell->get_dof_indices (local_dof_indices);
898 *  
899 *   constraints_elastic.distribute_local_to_global (cell_matrix_elastic,
900 *   cell_rhs_elastic, local_dof_indices, system_matrix_elastic,
901 *   system_rhs_elastic);
902 *   }
903 *  
904 *   system_matrix_elastic.compress (VectorOperation::add);
905 *   system_rhs_elastic.compress (VectorOperation::add);
906 *   }
907 *  
908 *   void
909 *   PhaseField::solve_linear_system_elastic ()
910 *   {
911 *   TimerOutput::Scope ts (computing_timer, "solve_linear_system_elastic");
912 *  
913 *   SolverControl solver_control (10000,
914 *   1e-12* system_rhs_elastic.l2_norm());
915 *  
916 *   SolverCG<LA::MPI::Vector> solver (solver_control);
917 *  
918 *   LA::MPI::PreconditionAMG::AdditionalData data;
919 *   #ifdef USE_PETSC_LA
920 *   data.symmetric_operator = true;
921 *   #else
922 * @endcode
923 *
924 * Trilinos defaults are good
925 *
926 * @code
927 *   #endif
928 *   LA::MPI::PreconditionAMG preconditioner;
929 *   preconditioner.initialize (system_matrix_elastic, data);
930 *  
931 *   solver.solve (system_matrix_elastic,
932 *   completely_distributed_solution_elastic, system_rhs_elastic,
933 *   preconditioner);
934 *  
935 *   pcout << " Solved in " << solver_control.last_step () << " iterations."
936 *   << std::endl;
937 *  
938 *   constraints_elastic.distribute (completely_distributed_solution_elastic);
939 *  
940 *   locally_relevant_solution_elastic = completely_distributed_solution_elastic;
941 *   }
942 *  
943 *   void
944 *   PhaseField::setup_boundary_values_damage ()
945 *   {
946 *   TimerOutput::Scope ts (computing_timer, "setup_bv_damage");
947 *  
948 *   constraints_damage.clear ();
949 *   constraints_damage.reinit (locally_relevant_dofs_damage);
950 *   DoFTools::make_hanging_node_constraints (dof_handler_damage,
951 *   constraints_damage);
952 *  
953 * @endcode
954 *
955 * Create initial crack, if any
956 *
957 * @code
958 *   /*for (const auto &cell : dof_handler_damage.active_cell_iterators())
959 *   if (cell->is_locally_owned())
960 *  
961 *   {
962 * @endcode
963 *
964 * for (const auto &face : cell->face_iterators())
965 *
966
967 *
968 *
969 * @code
970 *   for (const auto vertex_number : cell->vertex_indices())
971 *   {
972 *   const auto vert = cell->vertex(vertex_number);
973 *   const Point<3>& node = vert;
974 *   if (((std::fabs((Domain::z_max*(node(0)-0.5*Domian::x_max + A)) - 2*A*( node(2))) <10*l)
975 *   && node(1)>=0.9*Domain::y_max)
976 *   || ((std::fabs((node(0)-0.5*Domian::x_max + A)*Domain::z_max+2*A*(node(2)-Domain::z_max))<10*l)
977 *   && node(1)<=0.1*Domain::y_max))
978 *   if ((vert(0) - 0.5*(Domain::x_min+Domian::x_max) < 1e-12) &&
979 *   (std::fabs(vert(1) - 0.5*(Domain::y_min + Domain::y_max)) <= bandwidth)) // nodes on initial damage plane
980 *  
981 *  
982 *   {
983 *   const unsigned int dof = cell->vertex_dof_index(vertex_number, 0);
984 *   constraints_damage.add_line(dof);
985 *   constraints_damage.set_inhomogeneity(dof,1);
986 *   }
987 *  
988 *   }
989 *  
990 *   }*/
991 *   constraints_damage.close ();
992 *  
993 *   }
994 *  
995 *   void
996 *   PhaseField::setup_system_damage ()
997 *   {
998 *   TimerOutput::Scope ts (computing_timer, "setup_system_damage");
999 *  
1000 *   locally_owned_dofs_damage = dof_handler_damage.locally_owned_dofs ();
1001 *   locally_relevant_dofs_damage = DoFTools::extract_locally_relevant_dofs (
1002 *   dof_handler_damage);
1003 *  
1004 *   locally_relevant_solution_damage.reinit (locally_owned_dofs_damage,
1005 *   locally_relevant_dofs_damage, mpi_communicator);
1006 *  
1007 *   system_rhs_damage.reinit (locally_owned_dofs_damage, mpi_communicator);
1008 *  
1009 *   completely_distributed_solution_damage.reinit (locally_owned_dofs_damage,
1010 *   mpi_communicator);
1011 *  
1012 *   DynamicSparsityPattern dsp (locally_relevant_dofs_damage);
1013 *  
1014 *   DoFTools::make_sparsity_pattern (dof_handler_damage, dsp,
1015 *   constraints_damage, false);
1017 *   dof_handler_damage.locally_owned_dofs (), mpi_communicator,
1018 *   locally_relevant_dofs_damage);
1019 *  
1020 *   system_matrix_damage.reinit (locally_owned_dofs_damage,
1021 *   locally_owned_dofs_damage, dsp, mpi_communicator);
1022 *   }
1023 *  
1024 *   void
1025 *   PhaseField::assemble_system_damage ()
1026 *   {
1027 *  
1028 *   TimerOutput::Scope ts (computing_timer, "assembly_damage");
1029 *  
1030 *   FEValues<3> fe_values_damage (fe_damage, quadrature_formula_damage,
1033 *   FEValues<3> fe_values_elastic (fe_elastic, quadrature_formula_damage,
1036 *  
1037 *   const unsigned int dofs_per_cell = fe_damage.n_dofs_per_cell ();
1038 *   const unsigned int n_q_points = quadrature_formula_damage.size ();
1039 *  
1040 *   FullMatrix<double> cell_matrix_damage (dofs_per_cell, dofs_per_cell);
1041 *   Vector<double> cell_rhs_damage (dofs_per_cell);
1042 *  
1043 *   std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1044 *  
1045 *   const RandomMedium::EnergyReleaseRate energy_release_rate;
1046 *   std::vector<double> energy_release_rate_values (n_q_points);
1047 *  
1048 * @endcode
1049 *
1050 * Storing strain tensor for all Gauss points of a cell in vector strain_values.
1051 *
1052 * @code
1053 *   std::vector<SymmetricTensor<2, 3>> strain_values (n_q_points);
1054 *  
1055 *   for (const auto &cell : dof_handler_damage.active_cell_iterators ())
1056 *   if (cell->is_locally_owned ())
1057 *   {
1058 *   std::vector<std::shared_ptr<MyQData>> qpdH =
1059 *   quadrature_point_history_field.get_data (cell);
1060 *  
1061 *   const DoFHandler<3>::active_cell_iterator elastic_cell =
1062 *   Triangulation<3>::active_cell_iterator (cell)->as_dof_handler_iterator (
1063 *   dof_handler_elastic);
1064 *  
1065 *   fe_values_damage.reinit (cell);
1066 *   fe_values_elastic.reinit (elastic_cell);
1067 *  
1068 *   const FEValuesExtractors::Vector displacements (
1069 *   /* first vector component = */0);
1070 *   fe_values_elastic[displacements].get_function_symmetric_gradients (
1071 *   locally_relevant_solution_elastic, strain_values);
1072 *  
1073 *   cell_matrix_damage = 0;
1074 *   cell_rhs_damage = 0;
1075 *  
1076 *   energy_release_rate.value_list (
1077 *   fe_values_damage.get_quadrature_points (),
1078 *   energy_release_rate_values);
1079 *  
1080 *   for (const unsigned int q_index : fe_values_damage.quadrature_point_indices ())
1081 *   {
1082 *   const auto &x_q = fe_values_damage.quadrature_point (q_index);
1083 *   double H_call = H_plus (strain_values[q_index]);
1084 *  
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 ();
1088 *  
1089 *   double g_c;
1090 *   if (is_in_middle_layer (cell_center,Domain::z1,Domain::z2))
1091 *   g_c = energy_release_rate_values[q_index];
1092 *   else
1093 *   g_c = gc (GC, beta, Domain::z1, Domain::z2, x_q);
1094 *  
1095 *   for (const unsigned int i : fe_values_damage.dof_indices ())
1096 *   {
1097 *   for (const unsigned int j : fe_values_damage.dof_indices ())
1098 *   {
1099 *  
1100 *   cell_matrix_damage (i, j) +=
1101 * @endcode
1102 *
1103 * contribution to stiffness from -laplace u term
1104 *
1105 * @code
1106 *   Conductivity_damage (x_q) * fe_values_damage.shape_grad (
1107 *   i, q_index)
1108 *   * // kappa*grad phi_i(x_q)
1109 *   fe_values_damage.shape_grad (j, q_index)
1110 *   * // grad phi_j(x_q)
1111 *   fe_values_damage.JxW (q_index) // dx
1112 *   +
1113 * @endcode
1114 *
1115 * Contribution to stiffness from u term
1116 *
1117
1118 *
1119 *
1120 * @code
1121 *   ((1 + (2 * l * H) / g_c) * (1 / pow (l, 2))
1122 *   * fe_values_damage.shape_value (i, q_index) * // phi_i(x_q)
1123 *   fe_values_damage.shape_value (j, q_index) * // phi_j(x_q)
1124 *   fe_values_damage.JxW (q_index)); // dx
1125 *   }
1126 * @endcode
1127 *
1128 * const auto &x_q = fe_values_damage.quadrature_point(q_index);
1129 *
1130 * @code
1131 *   cell_rhs_damage (i) += (fe_values_damage.shape_value (i,
1132 *   q_index)
1133 *   * // phi_i(x_q)
1134 *   (2 / (l * g_c))
1135 *   * H * fe_values_damage.JxW (q_index)); // dx
1136 *   }
1137 *   }
1138 *  
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);
1143 *   }
1144 *  
1145 *   system_matrix_damage.compress (VectorOperation::add);
1146 *   system_rhs_damage.compress (VectorOperation::add);
1147 *   }
1148 *  
1149 *   void
1150 *   PhaseField::solve_linear_system_damage ()
1151 *   {
1152 *   TimerOutput::Scope ts (computing_timer, "solve_linear_system_damage");
1153 *  
1154 *   SolverControl solver_control (10000,
1155 *   1e-12* system_rhs_damage.l2_norm());
1156 *  
1157 *   SolverCG<LA::MPI::Vector> solver (solver_control);
1158 *  
1159 *   LA::MPI::PreconditionAMG::AdditionalData data;
1160 *   #ifdef USE_PETSC_LA
1161 *   data.symmetric_operator = true;
1162 *   #else
1163 * @endcode
1164 *
1165 * Trilinos defaults are good
1166 *
1167 * @code
1168 *   #endif
1169 *   LA::MPI::PreconditionAMG preconditioner;
1170 *   preconditioner.initialize (system_matrix_damage, data);
1171 *  
1172 *   solver.solve (system_matrix_damage, completely_distributed_solution_damage,
1173 *   system_rhs_damage, preconditioner);
1174 *  
1175 *   pcout << " Solved in " << solver_control.last_step () << " iterations."
1176 *   << std::endl;
1177 *  
1178 *   constraints_damage.distribute (completely_distributed_solution_damage);
1179 *   locally_relevant_solution_damage = completely_distributed_solution_damage;
1180 *   }
1181 *  
1182 *   void
1183 *   PhaseField::solve_elastic_subproblem (const unsigned int load_step)
1184 *   {
1185 *   setup_system_elastic (load_step);
1186 *   assemble_system_elastic ();
1187 *   solve_linear_system_elastic ();
1188 *   }
1189 *  
1190 *   void
1191 *   PhaseField::output_results (const unsigned int load_step) const
1192 *  
1193 *   {
1194 *   std::vector<std::string> displacement_names (3, "displacement");
1195 *   std::vector<DataComponentInterpretation::DataComponentInterpretation> displacement_component_interpretation (
1197 *  
1198 *   DataOut<3> data_out_phasefield;
1199 *   data_out_phasefield.add_data_vector (dof_handler_elastic,
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");
1204 *  
1205 *   Vector<double> subdomain (triangulation.n_active_cells ());
1206 *   for (unsigned int i = 0; i < subdomain.size (); ++i)
1207 *   subdomain (i) = triangulation.locally_owned_subdomain ();
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);
1212 *   }
1213 *  
1214 *   void
1215 *   PhaseField::load_disp_calculation (const unsigned int load_step)
1216 *  
1217 *   {
1218 *   Tensor<1, 3> x_max_force; //force vector on the x_max face
1219 *   Tensor<1, 3> y_max_force; //force vector on the y_max face
1220 *  
1221 *   const QGauss<2> face_quadrature (fe_elastic.degree);
1222 *   FEFaceValues<3> fe_face_values (fe_elastic, face_quadrature,
1224 *   std::vector<SymmetricTensor<2, 3>> strain_values (face_quadrature.size ());
1225 *   const FEValuesExtractors::Vector displacements (0);
1226 *  
1227 *   FEFaceValues<3> fe_face_values_damage (fe_damage, face_quadrature, update_values);
1228 *   std::vector<double> damage_values (fe_face_values_damage.n_quadrature_points);
1229 *  
1230 *   for (const auto &cell : dof_handler_elastic.active_cell_iterators ())
1231 *   if (cell->is_locally_owned ())
1232 *   {
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 ()
1238 *   == 1))
1239 *   {
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);
1244 *  
1245 *   fe_face_values_damage.get_function_values (locally_relevant_solution_damage, damage_values);
1246 *  
1247 *   for (unsigned int q = 0; q < fe_face_values.n_quadrature_points; ++q)
1248 *  
1249 *   {
1250 *   const Tensor<2, 3> strain = strain_values[q]; //strain tensor at a gauss point
1251 *   const double tr_strain = strain[0][0] + strain[1][1] + strain[2][2];
1252 *   const double d = damage_values[q];
1253 *  
1254 *   Tensor<2, 3> stress;
1255 *   stress[0][0] = pow ((1 - d), 2)
1256 *   * (lambda (E, nu) * tr_strain + 2 * mu (E, nu)
1257 *   * strain[0][0]);
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)
1264 *   * strain[1][1]);
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)
1269 *   * strain[2][2]);
1270 *  
1271 *   const Tensor<1, 3> force_density = stress
1272 *   * fe_face_values.normal_vector (q);
1273 *   x_max_force += force_density * fe_face_values.JxW (q);
1274 *   }
1275 *   }
1276 *  
1277 *   else if (cell->face (f)->at_boundary ()
1278 *   && (cell->face (f)->boundary_id () == 3))
1279 *   {
1280 *   fe_face_values.reinit (cell, f);
1281 *   fe_face_values_damage.reinit (damage_cell, f); //Convert 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);
1285 *  
1286 *   for (unsigned int q = 0; q < fe_face_values.n_quadrature_points;
1287 *   ++q)
1288 *  
1289 *   {
1290 *   const Tensor<2, 3> strain = strain_values[q]; //strain tensor at a gauss point
1291 *   const double tr_strain = strain[0][0] + strain[1][1] + strain[2][2];
1292 *   const double d = damage_values[q];
1293 *  
1294 *   Tensor<2, 3> stress;
1295 *   stress[0][0] = pow ((1 - d), 2)
1296 *   * (lambda (E, nu) * tr_strain + 2 * mu (E, nu)
1297 *   * strain[0][0]);
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)
1304 *   * strain[1][1]);
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)
1309 *   * strain[2][2]);
1310 *  
1311 *   const Tensor<1, 3> force_density = stress
1312 *   * fe_face_values.normal_vector (q);
1313 *   y_max_force += force_density * fe_face_values.JxW (q);
1314 *   }
1315 *   }
1316 *   }
1317 *   double x_max_force_x;
1318 *   x_max_force_x = x_max_force[0];
1319 *   x_max_force_x = Utilities::MPI::sum (x_max_force_x, mpi_communicator);
1320 *   pcout << "fx: " << x_max_force_x << std::endl;
1321 *  
1322 *   double y_max_force_y;
1323 *   y_max_force_y = y_max_force[1];
1324 *   y_max_force_y = Utilities::MPI::sum (y_max_force_y, mpi_communicator);
1325 *   pcout << "fy: " << y_max_force_y << std::endl;
1326 *  
1327 *   if (Utilities::MPI::this_mpi_process (mpi_communicator) == 0)
1328 *   {
1329 * @endcode
1330 *
1331 * Code to be executed only on process 0
1332 *
1333 * @code
1334 *   load_values_x[load_step] = x_max_force_x;
1335 *   load_values_y[load_step] = y_max_force_y;
1336 *   {
1337 *   const double disp = ux * load_step;
1338 *   displacement_values[load_step] = disp;
1339 *   }
1340 * @endcode
1341 *
1342 * load displacement plot
1343 *
1344 * @code
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"
1348 *   << std::endl;
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')"
1357 *   << std::endl;
1358 *   m_file << "ylabel(\"Reaction force (kN)\",'Interpreter', 'latex')"
1359 *   << std::endl;
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"
1363 *   << std::endl;
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;
1369 *   m_file.close ();
1370 *   }
1371 *   }
1372 *  
1373 *   void
1374 *   PhaseField::solve_damage_subproblem ()
1375 *   {
1376 *   setup_boundary_values_damage ();
1377 *   setup_system_damage ();
1378 *   assemble_system_damage ();
1379 *   solve_linear_system_damage ();
1380 *   }
1381 *  
1382 *   void
1383 *   PhaseField::refine_grid (const unsigned int load_step)
1384 *   {
1385 *   FEValues<3> fe_values_damage (fe_damage, quadrature_formula_damage,
1388 *  
1390 *   fe_damage, quadrature_formula_damage, quadrature_formula_damage);
1391 *  
1392 * @endcode
1393 *
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.
1397 *
1398 * @code
1399 *   Vector<float> estimated_error_per_cell (
1400 *   triangulation.n_locally_owned_active_cells ());
1401 *  
1402 *   KellyErrorEstimator<3>::estimate (dof_handler_damage,
1403 *   QGauss<2> (fe_damage.degree + 1),
1404 *   { }, locally_relevant_solution_damage, estimated_error_per_cell);
1405 *  
1406 * @endcode
1407 *
1408 * Initialize SolutionTransfer object
1409 *
1410 * @code
1411 *   parallel::distributed::SolutionTransfer<3, LA::MPI::Vector> soltransDamage (dof_handler_damage);
1412 *  
1413 * @endcode
1414 *
1415 * Initialize SolutionTransfer object
1416 *
1417 * @code
1418 *   parallel::distributed::SolutionTransfer<3, LA::MPI::Vector> soltransElastic (dof_handler_elastic);
1419 *  
1421 *   triangulation, estimated_error_per_cell, 0.01, // top 1% cells marked for refinement
1422 *   0.0); // bottom 0 % cells marked for coarsening
1423 *  
1424 *   if (triangulation.n_global_levels () >= 4)
1425 *   {
1426 *   for (const auto &cell : triangulation.active_cell_iterators_on_level (3))
1427 *   if (cell->is_locally_owned ())
1428 *   cell->clear_refine_flag ();
1429 *   }
1430 *  
1431 * @endcode
1432 *
1433 * prepare the triangulation,
1434 *
1435 * @code
1436 *   triangulation.prepare_coarsening_and_refinement ();
1437 *  
1438 * @endcode
1439 *
1440 * prepare CellDataStorage object for refinement
1441 *
1442 * @code
1443 *   data_transfer.prepare_for_coarsening_and_refinement (triangulation,
1444 *   quadrature_point_history_field);
1445 *  
1446 * @endcode
1447 *
1448 * prepare the SolutionTransfer object for coarsening and refinement
1449 * and give the solution vector that we intend to interpolate later,
1450 *
1451 * @code
1452 *   soltransDamage.prepare_for_coarsening_and_refinement (
1453 *   locally_relevant_solution_damage);
1454 *   soltransElastic.prepare_for_coarsening_and_refinement (
1455 *   locally_relevant_solution_elastic);
1456 *  
1457 *   triangulation.execute_coarsening_and_refinement ();
1458 *  
1459 * @endcode
1460 *
1461 * redistribute dofs,
1462 *
1463 * @code
1464 *   dof_handler_damage.distribute_dofs (fe_damage);
1465 *   dof_handler_elastic.distribute_dofs (fe_elastic);
1466 *  
1467 * @endcode
1468 *
1469 * Recreate locally_owned_dofs and locally_relevant_dofs index sets
1470 *
1471 * @code
1472 *   locally_owned_dofs_damage = dof_handler_damage.locally_owned_dofs ();
1473 *   locally_relevant_dofs_damage = DoFTools::extract_locally_relevant_dofs (
1474 *   dof_handler_damage);
1475 *  
1476 *   completely_distributed_solution_damage_old.reinit (
1477 *   locally_owned_dofs_damage, mpi_communicator);
1478 *   soltransDamage.interpolate (completely_distributed_solution_damage_old);
1479 *  
1480 * @endcode
1481 *
1482 * Apply constraints on the interpolated solution to make sure it conforms with the new mesh
1483 *
1484 * @code
1485 *   setup_boundary_values_damage ();
1486 *  
1487 *   constraints_damage.distribute (completely_distributed_solution_damage_old);
1488 *  
1489 * @endcode
1490 *
1491 * Copy completely_distributed_solution_damage_old to locally_relevant_solution_damage
1492 *
1493 * @code
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;
1498 *  
1499 * @endcode
1500 *
1501 * Interpolating elastic solution similarly
1502 *
1503 * @code
1504 *   locally_owned_dofs_elastic = dof_handler_elastic.locally_owned_dofs ();
1505 *   locally_relevant_dofs_elastic = DoFTools::extract_locally_relevant_dofs (
1506 *   dof_handler_elastic);
1507 *  
1508 *   completely_distributed_solution_elastic_old.reinit (
1509 *   locally_owned_dofs_elastic, mpi_communicator);
1510 *   soltransElastic.interpolate (completely_distributed_solution_elastic_old);
1511 *  
1512 * @endcode
1513 *
1514 * Apply constraints on the interpolated solution to make sure it conforms with the new mesh
1515 *
1516 * @code
1517 *   setup_constraints_elastic (load_step);
1518 *   constraints_elastic.distribute (
1519 *   completely_distributed_solution_elastic_old);
1520 *  
1521 * @endcode
1522 *
1523 * Copy completely_distributed_solution_damage_old to locally_relevant_solution_damage
1524 *
1525 * @code
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;
1530 *  
1531 *   for (const auto &cell : triangulation.active_cell_iterators ())
1532 *   {
1533 *   if (cell->is_locally_owned ())
1534 *   quadrature_point_history_field.initialize (cell, 8);
1535 *   }
1536 *   data_transfer.interpolate ();
1537 *  
1538 *   }
1539 *  
1540 *   bool
1541 *   PhaseField::check_convergence ()
1542 *   {
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);
1552 *  
1553 *   solution_damage_difference = locally_relevant_solution_damage;
1554 *  
1555 *   solution_damage_difference -= completely_distributed_solution_damage_old;
1556 *  
1557 *   solution_elastic_difference = locally_relevant_solution_elastic;
1558 *  
1559 *   solution_elastic_difference -= completely_distributed_solution_elastic_old;
1560 *  
1561 *   solution_damage_difference_ghost = solution_damage_difference;
1562 *  
1563 *   solution_elastic_difference_ghost = solution_elastic_difference;
1564 *  
1565 *   double error_elastic_solution_numerator, error_elastic_solution_denominator,
1566 *   error_damage_solution_numerator, error_damage_solution_denominator;
1567 *  
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 ();
1574 *  
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;
1580 *  
1581 *   if ((error_elastic_solution < tol) && (error_damage_solution < tol))
1582 *   return true;
1583 *   else
1584 *   return false;
1585 *   }
1586 *  
1587 *   void
1588 *   PhaseField::update_history_field ()
1589 *   {
1590 *   FEValues<3> fe_values_damage (fe_damage, quadrature_formula_damage,
1593 *  
1594 *   for (const auto &cell : dof_handler_damage.active_cell_iterators ())
1595 *   if (cell->is_locally_owned ())
1596 *   {
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;
1602 *   }
1603 *   }
1604 *  
1605 *   void
1606 *   PhaseField::run ()
1607 *   {
1608 *   Timer timer;
1609 *   timer.start ();
1610 *   pcout << "Running with "
1611 *   #ifdef USE_PETSC_LA
1612 *   << "PETSc"
1613 *   #else
1614 *   << "Trilinos"
1615 *   #endif
1616 *   << " on "
1617 *   << Utilities::MPI::n_mpi_processes (mpi_communicator) << " MPI rank(s)..."
1618 *   << std::endl;
1619 *  
1620 *   setup_mesh_and_bcs ();
1621 *  
1622 * @endcode
1623 *
1624 * Create the pre-crack in the domain
1625 *
1626 * @code
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 ())
1631 *   for (const auto vertex_number : cell->vertex_indices ())
1632 *   {
1633 *   const types::global_dof_index vertex_dof_index =
1634 *   cell->vertex_dof_index (vertex_number, 0);
1635 *   initial_soln_damage[vertex_dof_index] = 0;
1636 *   }
1637 *   initial_soln_damage.compress (VectorOperation::insert);
1638 *   locally_relevant_solution_damage = initial_soln_damage;
1639 *  
1640 * @endcode
1641 *
1642 * Loop over load steps
1643 *
1644 * @code
1645 *   for (unsigned int load_step = 1; load_step <= num_load_steps; load_step++)
1646 *   {
1647 *   pcout << " \n \n load increment number : " << load_step << std::endl;
1648 *  
1649 * @endcode
1650 *
1651 * Loop over staggered iterations
1652 *
1653 * @code
1654 *   unsigned int iteration = 0;
1655 *   bool stoppingCriterion = false;
1656 *   while (stoppingCriterion == false)
1657 *   {
1658 *   pcout << " \n iteration number:" << iteration << std::endl;
1659 *   solve_elastic_subproblem (load_step);
1660 *   solve_damage_subproblem ();
1661 *  
1662 *   locally_relevant_solution_damage.update_ghost_values ();
1663 *   locally_relevant_solution_elastic.update_ghost_values ();
1664 *  
1665 *   if (iteration > 0)
1666 *   stoppingCriterion = check_convergence ();
1667 *  
1668 *   completely_distributed_solution_elastic_old =
1669 *   locally_relevant_solution_elastic;
1670 *   completely_distributed_solution_damage_old =
1671 *   locally_relevant_solution_damage;
1672 *  
1673 *   if (stoppingCriterion == false)
1674 *   {
1675 *   refine_grid (load_step);
1676 *   }
1677 *   iteration = iteration + 1;
1678 *   }
1679 *  
1680 * @endcode
1681 *
1682 * Once converged, do some clean-up operations
1683 *
1684 * @code
1685 *   if ((load_step == 1) || (load_step >= 1 && load_step <= num_load_steps
1686 *   && std::fabs (load_step % 10) < 1e-6))
1687 *   {
1688 *   TimerOutput::Scope ts (computing_timer, "output");
1689 *   output_results (load_step);
1690 *   }
1691 *  
1692 *   load_disp_calculation (load_step);
1693 *  
1694 *   computing_timer.print_summary ();
1695 *   computing_timer.reset ();
1696 *   pcout << std::endl;
1697 *  
1698 *   update_history_field ();
1699 *   }
1700 *  
1701 *   timer.stop ();
1702 *   pcout << "Total run time: " << timer.wall_time () << " seconds."
1703 *   << std::endl;
1704 *   }
1705 *   }
1706 *  
1707 *   int
1708 *   main (int argc,
1709 *   char *argv[])
1710 *   {
1711 *   try
1712 *   {
1713 *   using namespace dealii;
1714 *   using namespace FracturePropagation;
1715 *  
1716 *   Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
1717 *  
1718 *   PhaseField phasefield;
1719 *   phasefield.run ();
1720 *   }
1721 *   catch (std::exception &exc)
1722 *   {
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;
1728 *  
1729 *   return 1;
1730 *   }
1731 *   catch (...)
1732 *   {
1733 *   std::cerr << std::endl << std::endl
1734 *   << "----------------------------------------------------" << std::endl;
1735 *   std::cerr << "Unknown exception!" << std::endl << "Aborting!" << std::endl
1736 *   << "----------------------------------------------------" << std::endl;
1737 *   return 1;
1738 *   }
1739 *  
1740 *   return 0;
1741 *  
1742 *   }
1743 *  
1744 * @endcode
1745
1746
1747*/
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={})
Definition fe_q.h:554
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)
Definition point.h:113
@ wall_times
Definition timer.h:651
Definition timer.h:117
void start()
Definition timer.cc:176
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
Point< 2 > first
Definition grid_out.cc:4632
unsigned int vertex_indices[2]
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
TriaActiveIterator< CellAccessor< dim, spacedim > > active_cell_iterator
Definition tria.h:1581
typename ActiveSelector::active_cell_iterator active_cell_iterator
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_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
Definition mpi.cc:746
const Event initial
Definition event.cc:68
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
void interpolate(const DoFHandler< dim, spacedim > &dof1, const InVector &u1, const DoFHandler< dim, spacedim > &dof2, OutVector &u2)
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 char N
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)
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm mpi_comm, const IndexSet &locally_relevant_rows)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:99
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:114
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask={})
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)
constexpr double E
Definition numbers.h:215
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
STL namespace.
::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)
Definition types.h:32
unsigned int global_dof_index
Definition types.h:94
unsigned int boundary_id
Definition types.h:161
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 > &)