deal.II version GIT relicensing-2173-gae8fc9d14b 2024-11-24 06:40:00+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\}}\)
Loading...
Searching...
No Matches
CeresFE.h
Go to the documentation of this file.
1
129 *  
130 *   /*
131 *   Summary of output files:
132 *  
133 *   One per run:
134 *  
135 *   - initial_mesh.eps : Visualization of initially imported mesh
136 *   - physical_times.txt : Columns are (1) step number corresponding to other files, (2) physical times at the time when each calculation is run in sec, (3) number of the final plasticity iteration in each timestep. Written in do_elastic_steps() for elastic steps and do_flow_step() for viscous steps
137 *  
138 *   One per time step:
139 *  
140 *   - timeXX_elastic_displacements.txt : Vtk-readable file with columns (1) x, (2) y, (3) u_x, (4) u_y, (5) P. Written in output_results() function, which is run immediately after solve().
141 *   - timeXX_baseviscosities.txt : Columns (1) cell x, (2) cell y, (3) base viscosity in Pa s. Written in solution_stresses().
142 *   - timeXX_surface.txt : Surface (defined as where P=0 boundary condition is applied) vertices at the beginning of timestep, except for the final timestep. Written in write_vertices() function, which is called immediately after setup_dofs() except for the final iteration, when it is called after move_mesh()
143 *  
144 *   One per plasticity step:
145 *  
146 *   - timeXX_flowYY.txt : Same as timeXX_elastic_displacements.txt above
147 *   - timeXX_principalstressesYY.txt : Columns with sigma1 and sigma3 at each cell. Same order as timeXX_baseviscosities.txt. Written in solution_stresses().
148 *   - timeXX_stresstensorYY.txt : Columns with components 11, 22, 33, and 13 of stress tensor at each cell. Written in solution_stresses().
149 *   - timeXX_failurelocationsYY.txt : Gives x,y coordinates of all cells where failure occurred. Written in solution_stresses().
150 *   - timeXX_viscositiesregYY.txt : Gives smoothed and regularized (i.e., floor and ceiling-filtered) effective viscosities. Written at end of solution_stresses().
151 *  
152 *   */
153 *  
154 *   #include <deal.II/base/quadrature_lib.h>
155 *   #include <deal.II/base/logstream.h>
156 *   #include <deal.II/base/function.h>
157 *   #include <deal.II/base/utilities.h>
158 *  
159 *   #include <deal.II/lac/block_vector.h>
160 *   #include <deal.II/lac/full_matrix.h>
161 *   #include <deal.II/lac/block_sparse_matrix.h>
162 *   #include <deal.II/lac/solver_cg.h>
163 *   #include <deal.II/lac/precondition.h>
164 *   #include <deal.II/lac/affine_constraints.h>
165 *  
166 *   #include <deal.II/grid/tria.h>
167 *   #include <deal.II/grid/grid_generator.h>
168 *   #include <deal.II/grid/tria_accessor.h>
169 *   #include <deal.II/grid/tria_iterator.h>
170 *   #include <deal.II/grid/grid_tools.h>
171 *   #include <deal.II/grid/manifold_lib.h>
172 *   #include <deal.II/grid/grid_refinement.h>
173 *   #include <deal.II/grid/grid_in.h>
174 *  
175 *   #include <deal.II/dofs/dof_handler.h>
176 *   #include <deal.II/dofs/dof_renumbering.h>
177 *   #include <deal.II/dofs/dof_accessor.h>
178 *   #include <deal.II/dofs/dof_tools.h>
179 *  
180 *   #include <deal.II/fe/fe_q.h>
181 *   #include <deal.II/fe/fe_system.h>
182 *   #include <deal.II/fe/fe_values.h>
183 *  
184 *   #include <deal.II/numerics/vector_tools.h>
185 *   #include <deal.II/numerics/matrix_tools.h>
186 *   #include <deal.II/numerics/data_out.h>
187 *   #include <deal.II/numerics/error_estimator.h>
188 *   #include <deal.II/numerics/derivative_approximation.h>
189 *   #include <deal.II/numerics/fe_field_function.h>
190 *  
191 *   #include <deal.II/lac/sparse_direct.h>
192 *   #include <deal.II/lac/sparse_ilu.h>
193 *  
194 *   #include <iostream>
195 *   #include <fstream>
196 *   #include <sstream>
197 *   #include <string>
198 *   #include <time.h>
199 *   #include <math.h>
200 *  
201 *   DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
202 *   #include <armadillo>
203 *   DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
204 *  
205 *   #include "../support_code/ellipsoid_grav.h"
206 *   #include "../support_code/ellipsoid_fit.h"
207 *   #include "../support_code/config_in.h"
208 *  
209 * @endcode
210 *
211 * As in all programs, the namespace dealii
212 * is included:
213 *
214 * @code
215 *   namespace Step22
216 *   {
217 *  
218 *   using namespace dealii;
219 *   using namespace arma;
220 *  
221 *   template<int dim>
222 *   struct InnerPreconditioner;
223 *  
224 *   template<>
225 *   struct InnerPreconditioner<2>
226 *   {
227 *   typedef SparseDirectUMFPACK type;
228 *   };
229 *  
230 *   template<>
231 *   struct InnerPreconditioner<3>
232 *   {
233 *   typedef SparseILU<double> type;
234 *   };
235 *  
236 * @endcode
237 *
238 * Auxiliary functions
239 *
240
241 *
242 *
243 * @code
244 *   template<int dim>
245 *   class AuxFunctions
246 *   {
247 *   public:
248 *   Tensor<2, 2> get_rotation_matrix(const std::vector<Tensor<1, 2> > &grad_u);
249 *   };
250 *  
251 *   template<int dim>
252 *   Tensor<2, 2> AuxFunctions<dim>::get_rotation_matrix(
253 *   const std::vector<Tensor<1, 2> > &grad_u)
254 *   {
255 *   const double curl = (grad_u[1][0] - grad_u[0][1]);
256 *  
257 *   const double angle = std::atan(curl);
258 *  
259 *   const double t[2][2] = { { cos(angle), sin(angle) }, {
260 *   -sin(angle), cos(
261 *   angle)
262 *   }
263 *   };
264 *   return Tensor<2, 2>(t);
265 *   }
266 *  
267 * @endcode
268 *
269 * Class for remembering material state/properties at each quadrature point
270 *
271
272 *
273 *
274 * @code
275 *   template<int dim>
276 *   struct PointHistory
277 *   {
278 *   SymmetricTensor<2, dim> old_stress;
279 *   double old_phiphi_stress;
280 *   double first_eta;
281 *   double new_eta;
282 *   double G;
283 *   };
284 *  
285 * @endcode
286 *
287 * Primary class of this problem
288 *
289
290 *
291 *
292 * @code
293 *   template<int dim>
294 *   class StokesProblem
295 *   {
296 *   public:
297 *   StokesProblem(const unsigned int degree);
298 *   void run();
299 *  
300 *   private:
301 *   void setup_dofs();
302 *   void assemble_system();
303 *   void solve();
304 *   void output_results() const;
305 *   void refine_mesh();
306 *   void solution_stesses();
307 *   void smooth_eta_field(std::vector<bool> failing_cells);
308 *  
309 *   void setup_initial_mesh();
310 *   void do_elastic_steps();
311 *   void do_flow_step();
312 *   void update_time_interval();
313 *   void initialize_eta_and_G();
314 *   void move_mesh();
315 *   void do_ellipse_fits();
316 *   void append_physical_times(int max_plastic);
317 *   void write_vertices(unsigned char);
318 *   void write_mesh();
319 *   void setup_quadrature_point_history();
320 *   void update_quadrature_point_history();
321 *  
322 *   const unsigned int degree;
323 *  
324 *   Triangulation<dim> triangulation;
325 *   const MappingQ1<dim> mapping;
326 *   FESystem<dim> fe;
327 *   DoFHandler<dim> dof_handler;
328 *   unsigned int n_u = 0, n_p = 0;
329 *   unsigned int plastic_iteration = 0;
330 *   unsigned int last_max_plasticity = 0;
331 *  
332 *   QGauss<dim> quadrature_formula;
333 *   std::vector< std::vector <Vector<double> > > quad_viscosities; // Indices for this object are [cell][q][q coords, eta]
334 *   std::vector<double> cell_viscosities; // This vector is only used for output, not FE computations
335 *   std::vector<PointHistory<dim> > quadrature_point_history;
336 *  
337 *   AffineConstraints<double> constraints;
338 *  
339 *   BlockSparsityPattern sparsity_pattern;
340 *   BlockSparseMatrix<double> system_matrix;
341 *  
342 *   BlockVector<double> solution;
343 *   BlockVector<double> system_rhs;
344 *  
345 *   std::shared_ptr<typename InnerPreconditioner<dim>::type> A_preconditioner;
346 *  
347 *   ellipsoid_fit<dim> ellipsoid;
348 *   };
349 *  
350 * @endcode
351 *
352 * Class for boundary conditions and rhs
353 *
354
355 *
356 *
357 * @code
358 *   template<int dim>
359 *   class BoundaryValuesP: public Function<dim>
360 *   {
361 *   public:
362 *   BoundaryValuesP() :
363 *   Function<dim>(dim + 1)
364 *   {
365 *   }
366 *  
367 *   virtual double value(const Point<dim> &p,
368 *   const unsigned int component = 0) const override;
369 *  
370 *   virtual void vector_value(const Point<dim> &p, Vector<double> &value) const override;
371 *   };
372 *  
373 *   template<int dim>
374 *   double BoundaryValuesP<dim>::value(const Point<dim> &p,
375 *   const unsigned int component) const
376 *   {
377 *   Assert(component < this->n_components,
378 *   ExcIndexRange (component, 0, this->n_components));
379 *   (void)p;
380 *   (void)component;
381 *  
382 *   Assert(p[0] >= -10.0, ExcLowerRangeType<double>(p[0], -10.0)); //value of -10 is to permit some small numerical error moving nodes left of x=0; a << value is in fact sufficient
383 *  
384 *   return 0;
385 *   }
386 *  
387 *   template<int dim>
388 *   void BoundaryValuesP<dim>::vector_value(const Point<dim> &p,
389 *   Vector<double> &values) const
390 *   {
391 *   for (unsigned int c = 0; c < this->n_components; ++c)
392 *   values(c) = BoundaryValuesP<dim>::value(p, c);
393 *   }
394 *  
395 *   template<int dim>
396 *   class RightHandSide: public Function<dim>
397 *   {
398 *   public:
399 *   RightHandSide () : Function<dim>(dim+1) {}
400 *  
401 *   using Function<dim>::value;
402 *   using Function<dim>::vector_value;
403 *   using Function<dim>::vector_value_list;
404 *  
405 *   virtual double value(const Point<dim> &p, const unsigned int component,
406 *   A_Grav_namespace::AnalyticGravity<dim> *aGrav) const;
407 *  
408 *   virtual void vector_value(const Point<dim> &p, Vector<double> &value,
409 *   A_Grav_namespace::AnalyticGravity<dim> *aGrav) const;
410 *  
411 *   virtual void vector_value_list(const std::vector<Point<dim> > &points,
412 *   std::vector<Vector<double> > &values,
413 *   A_Grav_namespace::AnalyticGravity<dim> *aGrav) const;
414 *  
415 *   };
416 *  
417 *   template<int dim>
418 *   double RightHandSide<dim>::value(const Point<dim> &p,
419 *   const unsigned int component,
420 *   A_Grav_namespace::AnalyticGravity<dim> *aGrav) const
421 *   {
422 *  
423 *   std::vector<double> temp_vector(2);
424 *   aGrav->get_gravity(p, temp_vector);
425 *  
426 *   if (component == 0)
427 *   {
428 *   return temp_vector[0] + system_parameters::omegasquared * p[0]; // * 1.2805;
429 *   }
430 *   else
431 *   {
432 *   if (component == 1)
433 *   return temp_vector[1];
434 *   else
435 *   return 0;
436 *   }
437 *   }
438 *  
439 *   template<int dim>
440 *   void RightHandSide<dim>::vector_value(const Point<dim> &p,
441 *   Vector<double> &values,
442 *   A_Grav_namespace::AnalyticGravity<dim> *aGrav) const
443 *   {
444 *   for (unsigned int c = 0; c < this->n_components; ++c)
445 *   values(c) = RightHandSide<dim>::value(p, c, aGrav);
446 *   }
447 *  
448 *   template<int dim>
449 *   void RightHandSide<dim>::vector_value_list(
450 *   const std::vector<Point<dim> > &points,
451 *   std::vector<Vector<double> > &values,
452 *   A_Grav_namespace::AnalyticGravity<dim> *aGrav) const
453 *   {
454 * @endcode
455 *
456 * check whether component is in
457 * the valid range is up to the
458 * derived class
459 *
460 * @code
461 *   Assert(values.size() == points.size(),
462 *   ExcDimensionMismatch(values.size(), points.size()));
463 *  
464 *   for (unsigned int i = 0; i < points.size(); ++i)
465 *   this->vector_value(points[i], values[i], aGrav);
466 *   }
467 *  
468 * @endcode
469 *
470 * Class for linear solvers and preconditioners
471 *
472
473 *
474 *
475 * @code
476 *   template<class Matrix, class Preconditioner>
477 *   class InverseMatrix: public Subscriptor
478 *   {
479 *   public:
480 *   InverseMatrix(const Matrix &m, const Preconditioner &preconditioner);
481 *  
482 *   void vmult(Vector<double> &dst, const Vector<double> &src) const;
483 *  
484 *   private:
485 *   const SmartPointer<const Matrix> matrix;
486 *   const SmartPointer<const Preconditioner> preconditioner;
487 *   };
488 *  
489 *   template<class Matrix, class Preconditioner>
490 *   InverseMatrix<Matrix, Preconditioner>::InverseMatrix(const Matrix &m,
491 *   const Preconditioner &preconditioner) :
492 *   matrix(&m), preconditioner(&preconditioner)
493 *   {
494 *   }
495 *  
496 *   template<class Matrix, class Preconditioner>
497 *   void InverseMatrix<Matrix, Preconditioner>::vmult(Vector<double> &dst,
498 *   const Vector<double> &src) const
499 *   {
500 *   SolverControl solver_control(1000 * src.size(), 1e-9 * src.l2_norm());
501 *  
502 *   SolverCG<> cg(solver_control);
503 *  
504 *   dst = 0;
505 *  
506 *   cg.solve(*matrix, dst, src, *preconditioner);
507 *   }
508 *  
509 * @endcode
510 *
511 * Class for the SchurComplement
512 *
513
514 *
515 *
516 * @code
517 *   template<class Preconditioner>
518 *   class SchurComplement: public Subscriptor
519 *   {
520 *   public:
521 *   SchurComplement(const BlockSparseMatrix<double> &system_matrix,
522 *   const InverseMatrix<SparseMatrix<double>, Preconditioner> &A_inverse);
523 *  
524 *   void vmult(Vector<double> &dst, const Vector<double> &src) const;
525 *  
526 *   private:
527 *   const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
528 *   const SmartPointer<const InverseMatrix<SparseMatrix<double>, Preconditioner> > A_inverse;
529 *  
530 *   mutable Vector<double> tmp1, tmp2;
531 *   };
532 *  
533 *   template<class Preconditioner>
534 *   SchurComplement<Preconditioner>::SchurComplement(
535 *   const BlockSparseMatrix<double> &system_matrix,
536 *   const InverseMatrix<SparseMatrix<double>, Preconditioner> &A_inverse) :
537 *   system_matrix(&system_matrix), A_inverse(&A_inverse), tmp1(
538 *   system_matrix.block(0, 0).m()), tmp2(
539 *   system_matrix.block(0, 0).m())
540 *   {
541 *   }
542 *  
543 *   template<class Preconditioner>
544 *   void SchurComplement<Preconditioner>::vmult(Vector<double> &dst,
545 *   const Vector<double> &src) const
546 *   {
547 *   system_matrix->block(0, 1).vmult(tmp1, src);
548 *   A_inverse->vmult(tmp2, tmp1);
549 *   system_matrix->block(1, 0).vmult(dst, tmp2);
550 *   }
551 *  
552 * @endcode
553 *
554 * StokesProblem::StokesProblem
555 *
556
557 *
558 *
559 * @code
560 *   template<int dim>
561 *   StokesProblem<dim>::StokesProblem(const unsigned int degree) :
562 *   degree(degree),
563 *   triangulation(Triangulation<dim>::maximum_smoothing),
564 *   mapping(),
565 *   fe(FE_Q<dim>(degree + 1), dim, FE_Q<dim>(degree), 1),
566 *   dof_handler(triangulation),
567 *   quadrature_formula(degree + 2),
568 *   ellipsoid(&triangulation)
569 *   {}
570 *  
571 * @endcode
572 *
573 * Set up dofs
574 *
575
576 *
577 *
578 * @code
579 *   template<int dim>
580 *   void StokesProblem<dim>::setup_dofs()
581 *   {
582 *   A_preconditioner.reset();
583 *   system_matrix.clear();
584 *  
585 *   dof_handler.distribute_dofs(fe);
586 *   DoFRenumbering::Cuthill_McKee(dof_handler);
587 *  
588 *   std::vector<unsigned int> block_component(dim + 1, 0);
589 *   block_component[dim] = 1;
590 *   DoFRenumbering::component_wise(dof_handler, block_component);
591 *  
592 * @endcode
593 *
594 * ========================================Apply Boundary Conditions=====================================
595 *
596 * @code
597 *   {
598 *   constraints.clear();
599 *   std::vector<bool> component_maskP(dim + 1, false);
600 *   component_maskP[dim] = true;
601 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
602 *   VectorTools::interpolate_boundary_values(dof_handler, 1,
603 *   BoundaryValuesP<dim>(), constraints, component_maskP);
604 *   }
605 *   {
606 *   std::set<types::boundary_id> no_normal_flux_boundaries;
607 *   no_normal_flux_boundaries.insert(99);
608 *   VectorTools::compute_no_normal_flux_constraints(dof_handler, 0,
609 *   no_normal_flux_boundaries, constraints);
610 *   }
611 *  
612 *   constraints.close();
613 *  
614 *   std::vector<types::global_dof_index> dofs_per_block = DoFTools::count_dofs_per_fe_block(dof_handler, block_component);
615 *   n_u = dofs_per_block[0];
616 *   n_p = dofs_per_block[1];
617 *  
618 *   std::cout << " Number of active cells: " << triangulation.n_active_cells()
619 *   << std::endl << " Number of degrees of freedom: "
620 *   << dof_handler.n_dofs() << " (" << n_u << '+' << n_p << ')'
621 *   << std::endl;
622 *  
623 *   {
624 *   BlockDynamicSparsityPattern csp(2, 2);
625 *  
626 *   csp.block(0, 0).reinit(n_u, n_u);
627 *   csp.block(1, 0).reinit(n_p, n_u);
628 *   csp.block(0, 1).reinit(n_u, n_p);
629 *   csp.block(1, 1).reinit(n_p, n_p);
630 *  
631 *   csp.collect_sizes();
632 *  
633 *   DoFTools::make_sparsity_pattern(dof_handler, csp, constraints, false);
634 *   sparsity_pattern.copy_from(csp);
635 *   }
636 *  
637 *   system_matrix.reinit(sparsity_pattern);
638 *  
639 *   solution.reinit(2);
640 *   solution.block(0).reinit(n_u);
641 *   solution.block(1).reinit(n_p);
642 *   solution.collect_sizes();
643 *  
644 *   system_rhs.reinit(2);
645 *   system_rhs.block(0).reinit(n_u);
646 *   system_rhs.block(1).reinit(n_p);
647 *   system_rhs.collect_sizes();
648 *  
649 *   }
650 *  
651 * @endcode
652 *
653 * Viscosity and Shear modulus functions
654 *
655
656 *
657 *
658 * @code
659 *   template<int dim>
660 *   class Rheology
661 *   {
662 *   public:
663 *   double get_eta(double &r, double &z);
664 *   double get_G(unsigned int mat_id);
665 *  
666 *   private:
667 *   std::vector<double> get_manual_eta_profile();
668 *  
669 *   };
670 *  
671 *   template<int dim>
672 *   std::vector<double> Rheology<dim>::get_manual_eta_profile()
673 *   {
674 *   vector<double> etas;
675 *  
676 *   for (unsigned int i=0; i < system_parameters::sizeof_depths_eta; ++i)
677 *   {
678 *   etas.push_back(system_parameters::depths_eta[i]);
679 *   etas.push_back(system_parameters::eta_kinks[i]);
680 *   }
681 *   return etas;
682 *   }
683 *  
684 *   template<int dim>
685 *   double Rheology<dim>::get_eta(double &r, double &z)
686 *   {
687 * @endcode
688 *
689 * compute local depth
690 *
691 * @code
692 *   double ecc = system_parameters::q_axes[0] / system_parameters::p_axes[0];
693 *   double Rminusr = system_parameters::q_axes[0] - system_parameters::p_axes[0];
694 *   double approx_a = std::sqrt(r * r + z * z * ecc * ecc);
695 *   double group1 = r * r + z * z - Rminusr * Rminusr;
696 *  
697 *   double a0 = approx_a;
698 *   double error = 10000;
699 * @endcode
700 *
701 * While loop finds the a axis of the "isodepth" ellipse for which the input point is on the surface.
702 * An "isodepth" ellipse is defined as one whose axes a,b are related to the global axes A, B by: A-h = B-h
703 *
704
705 *
706 *
707 * @code
708 *   if ((r > system_parameters::q_axes[0] - system_parameters::depths_eta.back()) ||
709 *   (z > system_parameters::p_axes[0] - system_parameters::depths_eta.back()))
710 *   {
711 *  
712 *   double eps = 10.0;
713 *   while (error >= eps)
714 *   {
715 *   double a02 = a0 * a0;
716 *   double a03 = a0 * a02;
717 *   double a04 = a0 * a03;
718 *   double fofa = a04 - (2 * Rminusr * a03) - (group1 * a02)
719 *   + (2 * r * r * Rminusr * a0) - (r * r * Rminusr * Rminusr);
720 *   double fprimeofa = 4 * a03 - (6 * Rminusr * a02) - (2 * group1 * a0)
721 *   + (2 * r * r * Rminusr);
722 *   double deltaa = -fofa / fprimeofa;
723 *   a0 += deltaa;
724 *   error = std::abs(deltaa);
725 * @endcode
726 *
727 * cout << "error = " << error << endl;
728 *
729 * @code
730 *   }
731 *   }
732 *   else
733 *   {
734 *   a0 = 0.0;
735 *   }
736 *  
737 *   double local_depth = system_parameters::q_axes[0] - a0;
738 *   if (local_depth < 0)
739 *   local_depth = 0;
740 *  
741 *   if (local_depth > system_parameters::depths_eta.back())
742 *   {
743 *   if (system_parameters::eta_kinks.back() < system_parameters::eta_floor)
744 *   return system_parameters::eta_floor;
745 *   else if (system_parameters::eta_kinks.back() > system_parameters::eta_ceiling)
746 *   return system_parameters::eta_ceiling;
747 *   else
748 *   return system_parameters::eta_kinks.back();
749 *   }
750 *  
751 *   std::vector<double> viscosity_function = get_manual_eta_profile();
752 *  
753 *   unsigned int n_visc_kinks = viscosity_function.size() / 2;
754 *  
755 * @endcode
756 *
757 * find the correct interval to do the interpolation in
758 *
759 * @code
760 *   int n_minus_one = -1;
761 *   for (unsigned int n = 1; n <= n_visc_kinks; ++n)
762 *   {
763 *   unsigned int ndeep = 2 * n - 2;
764 *   unsigned int nshallow = 2 * n;
765 *   if (local_depth >= viscosity_function[ndeep] && local_depth <= viscosity_function[nshallow])
766 *   n_minus_one = ndeep;
767 *   }
768 *  
769 * @endcode
770 *
771 * find the viscosity interpolation
772 *
773 * @code
774 *   if (n_minus_one == -1)
775 *   return system_parameters::eta_ceiling;
776 *   else
777 *   {
778 *   double visc_exponent =
779 *   (viscosity_function[n_minus_one]
780 *   - local_depth)
781 *   / (viscosity_function[n_minus_one]
782 *   - viscosity_function[n_minus_one + 2]);
783 *   double visc_base = viscosity_function[n_minus_one + 3]
784 *   / viscosity_function[n_minus_one + 1];
785 * @endcode
786 *
787 * This is the true viscosity given the thermal profile
788 *
789 * @code
790 *   double true_eta = viscosity_function[n_minus_one + 1] * std::pow(visc_base, visc_exponent);
791 *  
792 * @endcode
793 *
794 * Implement latitude-dependence viscosity
795 *
796 * @code
797 *   if (system_parameters::lat_dependence)
798 *   {
799 *   double lat = 180 / PI * std::atan(z / r);
800 *   if (lat > 80)
801 *   lat = 80;
802 *   double T_eq = 155;
803 *   double T_surf = T_eq * std::sqrt( std::sqrt( std::cos( PI / 180 * lat ) ) );
804 *   double taper_depth = 40000;
805 *   double surface_taper = (taper_depth - local_depth) / taper_depth;
806 *   if (surface_taper < 0)
807 *   surface_taper = 0;
808 *   double log_eta_contrast = surface_taper * system_parameters::eta_Ea * 52.5365 * (T_eq - T_surf) / T_eq / T_surf;
809 *   true_eta *= std::pow(10, log_eta_contrast);
810 *   }
811 *  
812 *   if (true_eta > system_parameters::eta_ceiling)
813 *   return system_parameters::eta_ceiling;
814 *   else if (true_eta < system_parameters::eta_floor)
815 *   return system_parameters::eta_floor;
816 *   else
817 *   return true_eta;
818 *   }
819 *   }
820 *  
821 *  
822 *   template<int dim>
823 *   double Rheology<dim>::get_G(unsigned int mat_id)
824 *   {
825 *   return system_parameters::G[mat_id];
826 *   }
827 *  
828 *  
829 * @endcode
830 *
831 * Initialize the eta and G parts of the quadrature_point_history object
832 *
833
834 *
835 *
836 * @code
837 *   template<int dim>
838 *   void StokesProblem<dim>::initialize_eta_and_G()
839 *   {
840 *   FEValues<dim> fe_values(fe, quadrature_formula, update_quadrature_points);
841 *  
842 *   const unsigned int n_q_points = quadrature_formula.size();
843 *   Rheology<dim> rheology;
844 *  
845 *   for (typename DoFHandler<dim>::active_cell_iterator cell =
846 *   dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
847 *   {
848 *   PointHistory<dim> *local_quadrature_points_history =
849 *   reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
850 *   Assert(
851 *   local_quadrature_points_history >= &quadrature_point_history.front(),
852 *   ExcInternalError());
853 *   Assert(
854 *   local_quadrature_points_history < &quadrature_point_history.back(),
855 *   ExcInternalError());
856 *   fe_values.reinit(cell);
857 *  
858 *   for (unsigned int q = 0; q < n_q_points; ++q)
859 *   {
860 *  
861 *   double r_value = fe_values.quadrature_point(q)[0];
862 *   double z_value = fe_values.quadrature_point(q)[1];
863 *  
864 * @endcode
865 *
866 * defines local viscosity
867 *
868 * @code
869 *   double local_viscosity = 0;
870 *  
871 *   local_viscosity = rheology.get_eta(r_value, z_value);
872 *  
873 *   local_quadrature_points_history[q].first_eta = local_viscosity;
874 *   local_quadrature_points_history[q].new_eta = local_viscosity;
875 *  
876 * @endcode
877 *
878 * defines local shear modulus
879 *
880 * @code
881 *   double local_G = 0;
882 *  
883 *   unsigned int mat_id = cell->material_id();
884 *  
885 *   local_G = rheology.get_G(mat_id);
886 *   local_quadrature_points_history[q].G = local_G;
887 *  
888 * @endcode
889 *
890 * initializes the phi-phi stress
891 *
892 * @code
893 *   local_quadrature_points_history[q].old_phiphi_stress = 0;
894 *   }
895 *   }
896 *   }
897 *  
898 * @endcode
899 *
900 * ====================== ASSEMBLE THE SYSTEM ======================
901 *
902
903 *
904 *
905 * @code
906 *   template<int dim>
907 *   void StokesProblem<dim>::assemble_system()
908 *   {
909 *   system_matrix = 0;
910 *   system_rhs = 0;
911 *  
912 *   FEValues<dim> fe_values(fe, quadrature_formula,
913 *   update_values | update_quadrature_points | update_JxW_values
914 *   | update_gradients);
915 *  
916 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
917 *  
918 *   const unsigned int n_q_points = quadrature_formula.size();
919 *  
920 *   FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
921 *   Vector<double> local_rhs(dofs_per_cell);
922 *  
923 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
924 *  
925 * @endcode
926 *
927 * runs the gravity script function
928 *
929 * @code
930 *   const RightHandSide<dim> right_hand_side;
931 *  
932 *   A_Grav_namespace::AnalyticGravity<dim> *aGrav =
933 *   new A_Grav_namespace::AnalyticGravity<dim>;
934 *   std::vector<double> grav_parameters;
935 *   grav_parameters.push_back(system_parameters::q_axes[system_parameters::present_timestep * 2 + 0]);
936 *   grav_parameters.push_back(system_parameters::p_axes[system_parameters::present_timestep * 2 + 0]);
937 *   grav_parameters.push_back(system_parameters::q_axes[system_parameters::present_timestep * 2 + 1]);
938 *   grav_parameters.push_back(system_parameters::p_axes[system_parameters::present_timestep * 2 + 1]);
939 *   grav_parameters.push_back(system_parameters::rho[0]);
940 *   grav_parameters.push_back(system_parameters::rho[1]);
941 *  
942 *   std::cout << "Body parameters are: " ;
943 *   for (int i=0; i<6; ++i)
944 *   std::cout << grav_parameters[i] << " ";
945 *   std::cout << endl;
946 *  
947 *   aGrav->setup_vars(grav_parameters);
948 *  
949 *   std::vector<Vector<double> > rhs_values(n_q_points,
950 *   Vector<double>(dim + 1));
951 *  
952 *   const FEValuesExtractors::Vector velocities(0);
953 *   const FEValuesExtractors::Scalar pressure(dim);
954 *  
955 *   std::vector<SymmetricTensor<2, dim> > phi_grads_u(dofs_per_cell);
956 *   std::vector<double> div_phi_u(dofs_per_cell);
957 *   std::vector<Tensor<1, dim> > phi_u(dofs_per_cell);
958 *   std::vector<double> phi_p(dofs_per_cell);
959 *  
960 *   typename DoFHandler<dim>::active_cell_iterator cell =
961 *   dof_handler.begin_active(), first_cell = dof_handler.begin_active(),
962 *   endc = dof_handler.end();
963 *  
964 *   for (; cell != endc; ++cell)
965 *   {
966 *   PointHistory<dim> *local_quadrature_points_history =
967 *   reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
968 *   Assert(
969 *   local_quadrature_points_history >= &quadrature_point_history.front(),
970 *   ExcInternalError());
971 *   Assert(
972 *   local_quadrature_points_history < &quadrature_point_history.back(),
973 *   ExcInternalError());
974 *  
975 *   double cell_area = cell->measure();
976 *  
977 *   if (cell_area<0)
978 *   append_physical_times(-1); // This writes final line to physical_times.txt if step terminates prematurely
979 *   AssertThrow(cell_area > 0
980 *   ,
981 *   ExcInternalError());
982 *  
983 *  
984 *   unsigned int m_id = cell->material_id();
985 *  
986 * @endcode
987 *
988 * initializes the rhs vector to the correct g values
989 *
990 * @code
991 *   fe_values.reinit(cell);
992 *   right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
993 *   rhs_values, aGrav);
994 *  
995 *   std::vector<Vector<double> > new_viscosities(quadrature_formula.size(), Vector<double>(dim + 1));
996 *  
997 * @endcode
998 *
999 * Finds vertices where the radius is zero DIM
1000 *
1001 * @code
1002 *   bool is_singular = false;
1003 *   for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
1004 *   if (cell->face(f)->center()[0] == 0)
1005 *   is_singular = true;
1006 *  
1007 *   if (is_singular == false || system_parameters::cylindrical == false)
1008 *   {
1009 *   local_matrix = 0;
1010 *   local_rhs = 0;
1011 *  
1012 * @endcode
1013 *
1014 * ===== outputs the local gravity
1015 *
1016 * @code
1017 *   std::vector<Point<dim> > quad_points_list(n_q_points);
1018 *   quad_points_list = fe_values.get_quadrature_points();
1019 *  
1020 *   if (plastic_iteration
1021 *   == (system_parameters::max_plastic_iterations - 1))
1022 *   {
1023 *   if (cell != first_cell)
1024 *   {
1025 *   std::ofstream fout("gravity_field.txt", std::ios::app);
1026 *   fout << quad_points_list[0] << " " << rhs_values[0];
1027 *   fout.close();
1028 *   }
1029 *   else
1030 *   {
1031 *   std::ofstream fout("gravity_field.txt");
1032 *   fout << quad_points_list[0] << " " << rhs_values[0];
1033 *   fout.close();
1034 *   }
1035 *   }
1036 *  
1037 *   for (unsigned int q = 0; q < n_q_points; ++q)
1038 *   {
1039 *   SymmetricTensor<2, dim> &old_stress =
1040 *   local_quadrature_points_history[q].old_stress;
1041 *   double &local_old_phiphi_stress =
1042 *   local_quadrature_points_history[q].old_phiphi_stress;
1043 *   double r_value = fe_values.quadrature_point(q)[0];
1044 *  
1045 * @endcode
1046 *
1047 * get local density based on mat id
1048 *
1049 * @code
1050 *   double local_density = system_parameters::rho[m_id];
1051 *  
1052 * @endcode
1053 *
1054 * defines local viscosities
1055 *
1056 * @code
1057 *   double local_viscosity = 0;
1058 *   if (plastic_iteration == 0)
1059 *   local_viscosity = local_quadrature_points_history[q].first_eta;
1060 *   else
1061 *   local_viscosity = local_quadrature_points_history[q].new_eta;
1062 *  
1063 * @endcode
1064 *
1065 * Define the local viscoelastic constants
1066 *
1067 * @code
1068 *   double local_eta_ve = 2
1069 *   / ((1 / local_viscosity)
1070 *   + (1 / local_quadrature_points_history[q].G
1071 *   / system_parameters::current_time_interval));
1072 *   double local_chi_ve = 1
1073 *   / (1
1074 *   + (local_quadrature_points_history[q].G
1075 *   * system_parameters::current_time_interval
1076 *   / local_viscosity));
1077 *  
1078 *   for (unsigned int k = 0; k < dofs_per_cell; ++k)
1079 *   {
1080 *   phi_grads_u[k] = fe_values[velocities].symmetric_gradient(k,
1081 *   q);
1082 *   div_phi_u[k] = (fe_values[velocities].divergence(k, q));
1083 *   phi_u[k] = (fe_values[velocities].value(k, q));
1084 *   if (system_parameters::cylindrical == true)
1085 *   {
1086 *   div_phi_u[k] *= (r_value);
1087 *   div_phi_u[k] += (phi_u[k][0]);
1088 *   }
1089 *   phi_p[k] = fe_values[pressure].value(k, q);
1090 *   }
1091 *  
1092 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1093 *   {
1094 *   for (unsigned int j = 0; j <= i; ++j)
1095 *   {
1096 *   if (system_parameters::cylindrical == true)
1097 *   {
1098 *   local_matrix(i, j) += (phi_grads_u[i]
1099 *   * phi_grads_u[j] * 2 * local_eta_ve
1100 *   * r_value
1101 *   + 2 * phi_u[i][0] * phi_u[j][0]
1102 *   * local_eta_ve / r_value
1103 *   - div_phi_u[i] * phi_p[j]
1104 *   * system_parameters::pressure_scale
1105 *   - phi_p[i] * div_phi_u[j]
1106 *   * system_parameters::pressure_scale
1107 *   + phi_p[i] * phi_p[j] * r_value
1108 *   * system_parameters::pressure_scale)
1109 *   * fe_values.JxW(q);
1110 *   }
1111 *   else
1112 *   {
1113 *   local_matrix(i, j) += (phi_grads_u[i]
1114 *   * phi_grads_u[j] * 2 * local_eta_ve
1115 *   - div_phi_u[i] * phi_p[j]
1116 *   * system_parameters::pressure_scale
1117 *   - phi_p[i] * div_phi_u[j]
1118 *   * system_parameters::pressure_scale
1119 *   + phi_p[i] * phi_p[j]) * fe_values.JxW(q);
1120 *   }
1121 *   }
1122 *   if (system_parameters::cylindrical == true)
1123 *   {
1124 *   const unsigned int component_i =
1125 *   fe.system_to_component_index(i).first;
1126 *   local_rhs(i) += (fe_values.shape_value(i, q)
1127 *   * rhs_values[q](component_i) * r_value
1128 *   * local_density
1129 *   - local_chi_ve * phi_grads_u[i] * old_stress
1130 *   * r_value
1131 *   - local_chi_ve * phi_u[i][0]
1132 *   * local_old_phiphi_stress)
1133 *   * fe_values.JxW(q);
1134 *   }
1135 *   else
1136 *   {
1137 *   const unsigned int component_i =
1138 *   fe.system_to_component_index(i).first;
1139 *   local_rhs(i) += fe_values.shape_value(i, q)
1140 *   * rhs_values[q](component_i) * fe_values.JxW(q)
1141 *   * local_density;
1142 *   }
1143 *   }
1144 *   }
1145 *   } // end of non-singular
1146 *   else
1147 *   {
1148 *   local_matrix = 0;
1149 *   local_rhs = 0;
1150 *  
1151 * @endcode
1152 *
1153 * ===== outputs the local gravity
1154 *
1155 * @code
1156 *   std::vector<Point<dim> > quad_points_list(n_q_points);
1157 *   quad_points_list = fe_values.get_quadrature_points();
1158 *  
1159 *   for (unsigned int q = 0; q < n_q_points; ++q)
1160 *   {
1161 *   const SymmetricTensor<2, dim> &old_stress =
1162 *   local_quadrature_points_history[q].old_stress;
1163 *   double &local_old_phiphi_stress =
1164 *   local_quadrature_points_history[q].old_phiphi_stress;
1165 *   double r_value = fe_values.quadrature_point(q)[0];
1166 *  
1167 *   double local_density = system_parameters::rho[m_id];
1168 *  
1169 * @endcode
1170 *
1171 * defines local viscosities
1172 *
1173 * @code
1174 *   double local_viscosity = 0;
1175 *   if (plastic_iteration == 0)
1176 *   {
1177 *   local_viscosity = local_quadrature_points_history[q].first_eta;
1178 *   }
1179 *   else
1180 *   local_viscosity = local_quadrature_points_history[q].new_eta;
1181 *  
1182 * @endcode
1183 *
1184 * Define the local viscoelastic constants
1185 *
1186 * @code
1187 *   double local_eta_ve = 2
1188 *   / ((1 / local_viscosity)
1189 *   + (1 / local_quadrature_points_history[q].G
1190 *   / system_parameters::current_time_interval));
1191 *   double local_chi_ve = 1
1192 *   / (1
1193 *   + (local_quadrature_points_history[q].G
1194 *   * system_parameters::current_time_interval
1195 *   / local_viscosity));
1196 *  
1197 *   for (unsigned int k = 0; k < dofs_per_cell; ++k)
1198 *   {
1199 *   phi_grads_u[k] = fe_values[velocities].symmetric_gradient(k,
1200 *   q);
1201 *   div_phi_u[k] = (fe_values[velocities].divergence(k, q));
1202 *   phi_u[k] = (fe_values[velocities].value(k, q));
1203 *   if (system_parameters::cylindrical == true)
1204 *   {
1205 *   div_phi_u[k] *= (r_value);
1206 *   div_phi_u[k] += (phi_u[k][0]);
1207 *   }
1208 *   phi_p[k] = fe_values[pressure].value(k, q);
1209 *   }
1210 *  
1211 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1212 *   {
1213 *   for (unsigned int j = 0; j <= i; ++j)
1214 *   {
1215 *   if (system_parameters::cylindrical == true)
1216 *   {
1217 *   local_matrix(i, j) += (phi_grads_u[i]
1218 *   * phi_grads_u[j] * 2 * local_eta_ve
1219 *   * r_value
1220 *   + 2 * phi_u[i][0] * phi_u[j][0]
1221 *   * local_eta_ve / r_value
1222 *   - div_phi_u[i] * phi_p[j]
1223 *   * system_parameters::pressure_scale
1224 *   - phi_p[i] * div_phi_u[j]
1225 *   * system_parameters::pressure_scale
1226 *   + phi_p[i] * phi_p[j] * r_value
1227 *   * system_parameters::pressure_scale)
1228 *   * fe_values.JxW(q);
1229 *   }
1230 *   else
1231 *   {
1232 *   local_matrix(i, j) += (phi_grads_u[i]
1233 *   * phi_grads_u[j] * 2 * local_eta_ve
1234 *   - div_phi_u[i] * phi_p[j]
1235 *   * system_parameters::pressure_scale
1236 *   - phi_p[i] * div_phi_u[j]
1237 *   * system_parameters::pressure_scale
1238 *   + phi_p[i] * phi_p[j]) * fe_values.JxW(q);
1239 *   }
1240 *   }
1241 *   if (system_parameters::cylindrical == true)
1242 *   {
1243 *   const unsigned int component_i =
1244 *   fe.system_to_component_index(i).first;
1245 *   local_rhs(i) += (fe_values.shape_value(i, q)
1246 *   * rhs_values[q](component_i) * r_value
1247 *   * local_density
1248 *   - local_chi_ve * phi_grads_u[i] * old_stress
1249 *   * r_value
1250 *   - local_chi_ve * phi_u[i][0]
1251 *   * local_old_phiphi_stress)
1252 *   * fe_values.JxW(q);
1253 *   }
1254 *   else
1255 *   {
1256 *   const unsigned int component_i =
1257 *   fe.system_to_component_index(i).first;
1258 *   local_rhs(i) += fe_values.shape_value(i, q)
1259 *   * rhs_values[q](component_i) * fe_values.JxW(q)
1260 *   * local_density;
1261 *   }
1262 *   }
1263 *   }
1264 *   } // end of singular
1265 *  
1266 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1267 *   for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
1268 *   local_matrix(i, j) = local_matrix(j, i);
1269 *  
1270 *   cell->get_dof_indices(local_dof_indices);
1271 *   constraints.distribute_local_to_global(local_matrix, local_rhs,
1272 *   local_dof_indices, system_matrix, system_rhs);
1273 *   }
1274 *  
1275 *   std::cout << " Computing preconditioner..." << std::endl << std::flush;
1276 *  
1277 *   A_preconditioner = std::shared_ptr<
1278 *   typename InnerPreconditioner<dim>::type>(
1279 *   new typename InnerPreconditioner<dim>::type());
1280 *   A_preconditioner->initialize(system_matrix.block(0, 0),
1281 *   typename InnerPreconditioner<dim>::type::AdditionalData());
1282 *  
1283 *   delete aGrav;
1284 *   }
1285 *  
1286 * @endcode
1287 *
1288 * ====================== SOLVER ======================
1289 *
1290
1291 *
1292 *
1293 * @code
1294 *   template<int dim>
1295 *   void StokesProblem<dim>::solve()
1296 *   {
1297 *   const InverseMatrix<SparseMatrix<double>,
1298 *   typename InnerPreconditioner<dim>::type> A_inverse(
1299 *   system_matrix.block(0, 0), *A_preconditioner);
1300 *   Vector<double> tmp(solution.block(0).size());
1301 *  
1302 *   {
1303 *   Vector<double> schur_rhs(solution.block(1).size());
1304 *   A_inverse.vmult(tmp, system_rhs.block(0));
1305 *   system_matrix.block(1, 0).vmult(schur_rhs, tmp);
1306 *   schur_rhs -= system_rhs.block(1);
1307 *  
1308 *   SchurComplement<typename InnerPreconditioner<dim>::type> schur_complement(
1309 *   system_matrix, A_inverse);
1310 *  
1311 *   int n_iterations = system_parameters::iteration_coefficient
1312 *   * solution.block(1).size();
1313 *   double tolerance_goal = system_parameters::tolerance_coefficient
1314 *   * schur_rhs.l2_norm();
1315 *  
1316 *   SolverControl solver_control(n_iterations, tolerance_goal);
1317 *   SolverCG<> cg(solver_control);
1318 *  
1319 *   std::cout << "\nMax iterations and tolerance are: " << n_iterations
1320 *   << " and " << tolerance_goal << std::endl;
1321 *  
1322 *   SparseILU<double> preconditioner;
1323 *   preconditioner.initialize(system_matrix.block(1, 1),
1324 *   SparseILU<double>::AdditionalData());
1325 *  
1326 *   InverseMatrix<SparseMatrix<double>, SparseILU<double> > m_inverse(
1327 *   system_matrix.block(1, 1), preconditioner);
1328 *  
1329 *   cg.solve(schur_complement, solution.block(1), schur_rhs, m_inverse);
1330 *  
1331 *   constraints.distribute(solution);
1332 *  
1333 *  
1334 *   std::cout << " " << solver_control.last_step()
1335 *   << " outer CG Schur complement iterations for pressure"
1336 *   << std::endl;
1337 *   }
1338 *  
1339 *   {
1340 *   system_matrix.block(0, 1).vmult(tmp, solution.block(1));
1341 *   tmp *= -1;
1342 *   tmp += system_rhs.block(0);
1343 *  
1344 *   A_inverse.vmult(solution.block(0), tmp);
1345 *   constraints.distribute(solution);
1346 *   solution.block(1) *= (system_parameters::pressure_scale);
1347 *   }
1348 *   }
1349 *  
1350 * @endcode
1351 *
1352 * ====================== OUTPUT RESULTS ======================
1353 *
1354 * @code
1355 *   template<int dim>
1356 *   void StokesProblem<dim>::output_results() const
1357 *   {
1358 *   std::vector < std::string > solution_names(dim, "velocity");
1359 *   solution_names.push_back("pressure");
1360 *  
1361 *   std::vector<DataComponentInterpretation::DataComponentInterpretation> data_component_interpretation(
1362 *   dim, DataComponentInterpretation::component_is_part_of_vector);
1363 *   data_component_interpretation.push_back(
1364 *   DataComponentInterpretation::component_is_scalar);
1365 *  
1366 *   DataOut<dim> data_out;
1367 *   data_out.attach_dof_handler(dof_handler);
1368 *   data_out.add_data_vector(solution, solution_names,
1369 *   DataOut<dim>::type_dof_data, data_component_interpretation);
1370 *   data_out.build_patches();
1371 *  
1372 *   std::ostringstream filename;
1373 *   if (system_parameters::present_timestep < system_parameters::initial_elastic_iterations)
1374 *   {
1375 *   filename << system_parameters::output_folder << "/time"
1376 *   << Utilities::int_to_string(system_parameters::present_timestep, 2)
1377 *   << "_elastic_displacements" << ".txt";
1378 *   }
1379 *   else
1380 *   {
1381 *   filename << system_parameters::output_folder << "/time"
1382 *   << Utilities::int_to_string(system_parameters::present_timestep, 2)
1383 *   << "_flow" << Utilities::int_to_string(plastic_iteration, 2) << ".txt";
1384 *   }
1385 *  
1386 *   std::ofstream output(filename.str().c_str());
1387 *   data_out.write_gnuplot(output);
1388 *   }
1389 *  
1390 * @endcode
1391 *
1392 * ====================== FIND AND WRITE TO FILE THE STRESS TENSOR; IMPLEMENT PLASTICITY ======================
1393 *
1394
1395 *
1396 *
1397 * @code
1398 *   template<int dim>
1399 *   void StokesProblem<dim>::solution_stesses()
1400 *   {
1401 * @endcode
1402 *
1403 * note most of this section only works with dim=2
1404 *
1405
1406 *
1407 * name the output text files
1408 *
1409 * @code
1410 *   std::ostringstream stress_output;
1411 *   stress_output << system_parameters::output_folder << "/time"
1412 *   << Utilities::int_to_string(system_parameters::present_timestep, 2)
1413 *   << "_principalstresses" << Utilities::int_to_string(plastic_iteration, 2)
1414 *   << ".txt";
1415 *   std::ofstream fout_snew(stress_output.str().c_str());
1416 *   fout_snew.close();
1417 *  
1418 *   std::ostringstream stresstensor_output;
1419 *   stresstensor_output << system_parameters::output_folder << "/time"
1420 *   << Utilities::int_to_string(system_parameters::present_timestep, 2)
1421 *   << "_stresstensor" << Utilities::int_to_string(plastic_iteration, 2)
1422 *   << ".txt";
1423 *   std::ofstream fout_sfull(stresstensor_output.str().c_str());
1424 *   fout_sfull.close();
1425 *  
1426 *   std::ostringstream failed_cells_output;
1427 *   failed_cells_output << system_parameters::output_folder << "/time"
1428 *   << Utilities::int_to_string(system_parameters::present_timestep, 2)
1429 *   << "_failurelocations" << Utilities::int_to_string(plastic_iteration, 2)
1430 *   << ".txt";
1431 *   std::ofstream fout_failed_cells(failed_cells_output.str().c_str());
1432 *   fout_failed_cells.close();
1433 *  
1434 *   std::ostringstream plastic_eta_output;
1435 *   plastic_eta_output << system_parameters::output_folder << "/time"
1436 *   << Utilities::int_to_string(system_parameters::present_timestep, 2)
1437 *   << "_viscositiesreg" << Utilities::int_to_string(plastic_iteration, 2)
1438 *   << ".txt";
1439 *   std::ofstream fout_vrnew(plastic_eta_output.str().c_str());
1440 *   fout_vrnew.close();
1441 *  
1442 *   std::ostringstream initial_eta_output;
1443 *   if (plastic_iteration == 0)
1444 *   {
1445 *   initial_eta_output << system_parameters::output_folder << "/time"
1446 *   << Utilities::int_to_string(system_parameters::present_timestep, 2)
1447 *   << "_baseviscosities.txt";
1448 *   std::ofstream fout_baseeta(initial_eta_output.str().c_str());
1449 *   fout_baseeta.close();
1450 *   }
1451 *  
1452 *   std::cout << "Running stress calculations for plasticity iteration "
1453 *   << plastic_iteration << "...\n";
1454 *  
1455 * @endcode
1456 *
1457 * This makes the set of points at which the stress tensor is calculated
1458 *
1459 * @code
1460 *   std::vector<Point<dim> > points_list(0);
1461 *   std::vector<unsigned int> material_list(0);
1462 *   typename DoFHandler<dim>::active_cell_iterator cell =
1463 *   dof_handler.begin_active(), endc = dof_handler.end();
1464 * @endcode
1465 *
1466 * This loop gets the gradients of the velocity field and saves it in the tensor_gradient_? objects DIM
1467 *
1468 * @code
1469 *   for (; cell != endc; ++cell)
1470 *   {
1471 *   points_list.push_back(cell->center());
1472 *   material_list.push_back(cell->material_id());
1473 *   }
1474 * @endcode
1475 *
1476 * Make the FEValues object to evaluate values and derivatives at quadrature points
1477 *
1478 * @code
1479 *   FEValues<dim> fe_values(fe, quadrature_formula,
1480 *   update_values | update_gradients | update_quadrature_points | update_JxW_values);
1481 *  
1482 * @endcode
1483 *
1484 * Make the object that will hold the velocities and velocity gradients at the quadrature points
1485 *
1486 * @code
1487 *   std::vector < std::vector<Tensor<1, dim> >> velocity_grads(quadrature_formula.size(),
1488 *   std::vector < Tensor<1, dim> > (dim + 1));
1489 *   std::vector<Vector<double> > velocities(quadrature_formula.size(),
1490 *   Vector<double>(dim + 1));
1491 * @endcode
1492 *
1493 * Make the object to find rheology
1494 *
1495 * @code
1496 *   Rheology<dim> rheology;
1497 *  
1498 * @endcode
1499 *
1500 * Write the solution flow velocity and derivative for each cell
1501 *
1502 * @code
1503 *   std::vector<Vector<double> > vector_values(0);
1504 *   std::vector < std::vector<Tensor<1, dim> > > gradient_values(0);
1505 *   std::vector<bool> failing_cells;
1506 * @endcode
1507 *
1508 * Write the stresses from the previous step into vectors
1509 *
1510 * @code
1511 *   std::vector<SymmetricTensor<2, dim>> old_stress;
1512 *   std::vector<double> old_phiphi_stress;
1513 *   std::vector<double> cell_Gs;
1514 *   for (typename DoFHandler<dim>::active_cell_iterator cell =
1515 *   dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
1516 *   {
1517 * @endcode
1518 *
1519 * Makes pointer to data in quadrature_point_history
1520 *
1521 * @code
1522 *   PointHistory<dim> *local_quadrature_points_history =
1523 *   reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
1524 *  
1525 *   fe_values.reinit(cell);
1526 *   fe_values.get_function_gradients(solution, velocity_grads);
1527 *   fe_values.get_function_values(solution, velocities);
1528 *   Vector<double> current_cell_velocity(dim+1);
1529 *   std::vector<Tensor<1, dim>> current_cell_grads(dim+1);
1530 *   SymmetricTensor<2, dim> current_cell_old_stress;
1531 *   current_cell_old_stress = 0;
1532 *   double current_cell_old_phiphi_stress = 0;
1533 *   double cell_area = 0;
1534 *  
1535 * @endcode
1536 *
1537 * Averages across each cell to find mean velocities, gradients, and old stresses
1538 *
1539 * @code
1540 *   for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
1541 *   {
1542 *   cell_area += fe_values.JxW(q);
1543 *   velocities[q] *= fe_values.JxW(q);
1544 *   current_cell_velocity += velocities[q];
1545 *   for (unsigned int i = 0; i < (dim+1); ++i)
1546 *   {
1547 *   velocity_grads[q][i] *= fe_values.JxW(q);
1548 *   current_cell_grads[i] += velocity_grads[q][i];
1549 *   }
1550 *   current_cell_old_stress += local_quadrature_points_history[q].old_stress * fe_values.JxW(q);
1551 *   current_cell_old_phiphi_stress += local_quadrature_points_history[q].old_phiphi_stress * fe_values.JxW(q);
1552 *   }
1553 *   current_cell_velocity /= cell_area;
1554 *   for (unsigned int i = 0; i < (dim+1); ++i)
1555 *   current_cell_grads[i] /= cell_area;
1556 *   current_cell_old_stress /= cell_area;
1557 *   current_cell_old_phiphi_stress /= cell_area;
1558 *  
1559 *   vector_values.push_back(current_cell_velocity);
1560 *   gradient_values.push_back(current_cell_grads);
1561 *   old_stress.push_back(current_cell_old_stress);
1562 *   old_phiphi_stress.push_back(current_cell_old_phiphi_stress);
1563 *  
1564 * @endcode
1565 *
1566 * Get cell shear modulus: assumes it's constant for the cell
1567 *
1568 * @code
1569 *   unsigned int mat_id = cell->material_id();
1570 *   double local_G = rheology.get_G(mat_id);
1571 *   cell_Gs.push_back(local_G);
1572 *   }
1573 *  
1574 * @endcode
1575 *
1576 * tracks where failure occurred
1577 *
1578 * @code
1579 *   std::vector<double> reduction_factor;
1580 *   unsigned int total_fails = 0;
1581 *   if (plastic_iteration == 0)
1582 *   cell_viscosities.resize(0);
1583 * @endcode
1584 *
1585 * loop across all the cells to find and adjust eta of failing cells
1586 *
1587 * @code
1588 *   for (unsigned int i = 0; i < triangulation.n_active_cells(); ++i)
1589 *   {
1590 *   double current_cell_viscosity = 0;
1591 *  
1592 * @endcode
1593 *
1594 * Fill viscosities vector, analytically if plastic_iteration == 0 and from previous viscosities for later iteration
1595 *
1596 * @code
1597 *   if (plastic_iteration == 0)
1598 *   {
1599 *   double local_viscosity;
1600 *   local_viscosity = rheology.get_eta(points_list[i][0], points_list[i][1]);
1601 *   current_cell_viscosity = local_viscosity;
1602 *   cell_viscosities.push_back(current_cell_viscosity);
1603 *   }
1604 *   else
1605 *   {
1606 *   current_cell_viscosity = cell_viscosities[i];
1607 *   }
1608 *  
1609 *  
1610 *   double cell_eta_ve = 2
1611 *   / ((1 / current_cell_viscosity)
1612 *   + (1 / cell_Gs[i]
1613 *   / system_parameters::current_time_interval));
1614 *   double cell_chi_ve = 1
1615 *   / (1
1616 *   + (cell_Gs[i]
1617 *   * system_parameters::current_time_interval
1618 *   / current_cell_viscosity));
1619 *  
1620 * @endcode
1621 *
1622 * find local pressure
1623 *
1624 * @code
1625 *   double cell_p = vector_values[i].operator()(2);
1626 * @endcode
1627 *
1628 * find stresses tensor
1629 * makes non-diagonalized local matrix A
1630 *
1631 * @code
1632 *   double sigma13 = 0.5
1633 *   * (gradient_values[i][0][1] + gradient_values[i][1][0]);
1634 *   mat A;
1635 *   A << gradient_values[i][0][0] << 0 << sigma13 << endr
1636 *   << 0 << vector_values[i].operator()(0) / points_list[i].operator()(0)<< 0 << endr
1637 *   << sigma13 << 0 << gradient_values[i][1][1] << endr;
1638 *   mat olddevstress;
1639 *   olddevstress << old_stress[i][0][0] << 0 << old_stress[i][0][1] << endr
1640 *   << 0 << old_phiphi_stress[i] << 0 << endr
1641 *   << old_stress[i][0][1] << 0 << old_stress[i][1][1] << endr;
1642 *   vec P;
1643 *   P << cell_p << cell_p << cell_p;
1644 *   mat Pmat = diagmat(P);
1645 *   mat B;
1646 *   B = (cell_eta_ve * A + cell_chi_ve * olddevstress) - Pmat;
1647 *  
1648 * @endcode
1649 *
1650 * finds principal stresses
1651 *
1652 * @code
1653 *   vec eigval;
1654 *   mat eigvec;
1655 *   eig_sym(eigval, eigvec, B);
1656 *   double sigma1 = -min(eigval);
1657 *   double sigma3 = -max(eigval);
1658 *  
1659 * @endcode
1660 *
1661 * Writes text files for principal stresses, full stress tensor, base viscosities
1662 *
1663 * @code
1664 *   std::ofstream fout_snew(stress_output.str().c_str(), std::ios::app);
1665 *   fout_snew << " " << sigma1 << " " << sigma3 << "\n";
1666 *   fout_snew.close();
1667 *  
1668 *   std::ofstream fout_sfull(stresstensor_output.str().c_str(), std::ios::app);
1669 *   fout_sfull << A(0,0) << " " << A(1,1) << " " << A(2,2) << " " << A(0,2) << "\n";
1670 *   fout_sfull.close();
1671 *  
1672 *   if (plastic_iteration == 0)
1673 *   {
1674 *   std::ofstream fout_baseeta(initial_eta_output.str().c_str(), std::ios::app);
1675 *   fout_baseeta << points_list[i]<< " " << current_cell_viscosity << "\n";
1676 *   fout_baseeta.close();
1677 *   }
1678 *  
1679 * @endcode
1680 *
1681 * Finds adjusted effective viscosity
1682 *
1683 * @code
1684 *   if (system_parameters::plasticity_on)
1685 *   {
1686 *   if (system_parameters::failure_criterion == 0) //Apply Byerlee's rule
1687 *   {
1688 *   if (sigma1 >= 5 * sigma3) // this guarantees that viscosities only go down, never up
1689 *   {
1690 *   failing_cells.push_back(true);
1691 *   double temp_reductionfactor = 1;
1692 *   if (sigma3 < 0)
1693 *   temp_reductionfactor = 100;
1694 *   else
1695 *   temp_reductionfactor = 1.9 * sigma1 / 5 / sigma3;
1696 *  
1697 *   reduction_factor.push_back(temp_reductionfactor);
1698 *   total_fails++;
1699 *  
1700 * @endcode
1701 *
1702 * Text file of all failure locations
1703 *
1704 * @code
1705 *   std::ofstream fout_failed_cells(failed_cells_output.str().c_str(), std::ios::app);
1706 *   fout_failed_cells << points_list[i] << "\n";
1707 *   fout_failed_cells.close();
1708 *   }
1709 *   else
1710 *   {
1711 *   reduction_factor.push_back(1);
1712 *   failing_cells.push_back(false);
1713 *   }
1714 *   }
1715 *   else
1716 *   {
1717 *   if (system_parameters::failure_criterion == 1) //Apply Schultz criterion for frozen sand, RMR=45
1718 *   {
1719 *   double temp_reductionfactor = 1;
1720 *   if (sigma3 < -114037)
1721 *   {
1722 * @endcode
1723 *
1724 * std::cout << " ext ";
1725 *
1726 * @code
1727 *   failing_cells.push_back(true);
1728 *   temp_reductionfactor = 10;
1729 *   reduction_factor.push_back(temp_reductionfactor);
1730 *   total_fails++;
1731 *  
1732 * @endcode
1733 *
1734 * Text file of all failure locations
1735 *
1736 * @code
1737 *   std::ofstream fout_failed_cells(failed_cells_output.str().c_str(), std::ios::app);
1738 *   fout_failed_cells << points_list[i] << "\n";
1739 *   fout_failed_cells.close();
1740 *   }
1741 *   else
1742 *   {
1743 *   double sigma_c = 160e6; //Unconfined compressive strength
1744 *   double yield_sigma1 = sigma3 + std::sqrt( (3.086 * sigma_c * sigma3) + (0.002 * sigma3 * sigma3) );
1745 *   if (sigma1 >= yield_sigma1)
1746 *   {
1747 * @endcode
1748 *
1749 * std::cout << " comp ";
1750 *
1751 * @code
1752 *   failing_cells.push_back(true);
1753 *   temp_reductionfactor = 1.0 * sigma1 / 5 / sigma3;
1754 *  
1755 *   reduction_factor.push_back(temp_reductionfactor);
1756 *   total_fails++;
1757 *  
1758 * @endcode
1759 *
1760 * Text file of all failure locations
1761 *
1762 * @code
1763 *   std::ofstream fout_failed_cells(failed_cells_output.str().c_str(), std::ios::app);
1764 *   fout_failed_cells << points_list[i] << "\n";
1765 *   fout_failed_cells.close();
1766 *   }
1767 *   else
1768 *   {
1769 *   reduction_factor.push_back(1);
1770 *   failing_cells.push_back(false);
1771 *   }
1772 *   }
1773 *   }
1774 *   else
1775 *   {
1776 *   std::cout << "Specified failure criterion not found\n";
1777 *   }
1778 *   }
1779 *   }
1780 *   else
1781 *   reduction_factor.push_back(1);
1782 *   }
1783 *  
1784 * @endcode
1785 *
1786 * If there are enough failed cells, update eta at all quadrature points and perform smoothing
1787 *
1788 * @code
1789 *   std::cout << " Number of failing cells: " << total_fails << "\n";
1790 *   double last_max_plasticity_double = last_max_plasticity;
1791 *   double total_fails_double = total_fails;
1792 *   double decrease_in_plasticity = ((last_max_plasticity_double - total_fails_double) / last_max_plasticity_double);
1793 *   if (plastic_iteration == 0)
1794 *   decrease_in_plasticity = 1;
1795 *   last_max_plasticity = total_fails;
1796 *   if (total_fails <= 100 || decrease_in_plasticity <= 0.2)
1797 *   {
1798 *   system_parameters::continue_plastic_iterations = false;
1799 *   for (unsigned int j=0; j < triangulation.n_active_cells(); ++j)
1800 *   {
1801 * @endcode
1802 *
1803 * Writes to file the undisturbed cell viscosities
1804 *
1805 * @code
1806 *   std::ofstream fout_vrnew(plastic_eta_output.str().c_str(), std::ios::app);
1807 *   fout_vrnew << " " << cell_viscosities[j] << "\n";
1808 *   fout_vrnew.close();
1809 *   }
1810 *   }
1811 *   else
1812 *   {
1813 *   quad_viscosities.resize(triangulation.n_active_cells());
1814 * @endcode
1815 *
1816 * Decrease the eta at quadrature points in failing cells
1817 *
1818 * @code
1819 *   unsigned int cell_no = 0;
1820 *   for (typename DoFHandler<dim>::active_cell_iterator cell =
1821 *   dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
1822 *   {
1823 *   fe_values.reinit(cell);
1824 * @endcode
1825 *
1826 * Make local_quadrature_points_history pointer to the cell data
1827 *
1828 * @code
1829 *   PointHistory<dim> *local_quadrature_points_history =
1830 *   reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
1831 *   Assert(
1832 *   local_quadrature_points_history >= &quadrature_point_history.front(),
1833 *   ExcInternalError());
1834 *   Assert(
1835 *   local_quadrature_points_history < &quadrature_point_history.back(),
1836 *   ExcInternalError());
1837 *  
1838 *   quad_viscosities[cell_no].resize(quadrature_formula.size());
1839 *  
1840 *   for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
1841 *   {
1842 *   if (plastic_iteration == 0)
1843 *   local_quadrature_points_history[q].new_eta = local_quadrature_points_history[q].first_eta;
1844 *   local_quadrature_points_history[q].new_eta /= reduction_factor[cell_no];
1845 * @endcode
1846 *
1847 * Prevents viscosities from dropping below the floor necessary for numerical stability
1848 *
1849 * @code
1850 *   if (local_quadrature_points_history[q].new_eta < system_parameters::eta_floor)
1851 *   local_quadrature_points_history[q].new_eta = system_parameters::eta_floor;
1852 *  
1853 *   quad_viscosities[cell_no][q].reinit(dim+1);
1854 *   for (unsigned int ii=0; ii<dim; ++ii)
1855 *   quad_viscosities[cell_no][q](ii) = fe_values.quadrature_point(q)[ii];
1856 *   quad_viscosities[cell_no][q](dim) = local_quadrature_points_history[q].new_eta;
1857 *   }
1858 *   cell_no++;
1859 *   }
1860 *   smooth_eta_field(failing_cells);
1861 *  
1862 * @endcode
1863 *
1864 * Writes to file the smoothed eta field (which is defined at each quadrature point) for each cell
1865 *
1866 * @code
1867 *   cell_no = 0;
1868 * @endcode
1869 *
1870 * cell_viscosities.resize(triangulation.n_active_cells(), 0);
1871 *
1872 * @code
1873 *   for (typename DoFHandler<dim>::active_cell_iterator cell =
1874 *   dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
1875 *   {
1876 *   if (failing_cells[cell_no])
1877 *   {
1878 *   fe_values.reinit(cell);
1879 * @endcode
1880 *
1881 * Averages across each cell to find mean eta
1882 *
1883 * @code
1884 *   double cell_area = 0;
1885 *   for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
1886 *   {
1887 *   cell_area += fe_values.JxW(q);
1888 *   cell_viscosities[cell_no] += quad_viscosities[cell_no][q][dim] * fe_values.JxW(q);
1889 *   }
1890 *   cell_viscosities[cell_no] /= cell_area;
1891 *  
1892 * @endcode
1893 *
1894 * Writes to file
1895 *
1896 * @code
1897 *   std::ofstream fout_vrnew(plastic_eta_output.str().c_str(), std::ios::app);
1898 *   fout_vrnew << " " << cell_viscosities[cell_no] << "\n";
1899 *   fout_vrnew.close();
1900 *   }
1901 *   else
1902 *   {
1903 *   std::ofstream fout_vrnew(plastic_eta_output.str().c_str(), std::ios::app);
1904 *   fout_vrnew << " " << cell_viscosities[cell_no] << "\n";
1905 *   fout_vrnew.close();
1906 *   }
1907 *   cell_no++;
1908 *   }
1909 *   }
1910 *   }
1911 *  
1912 * @endcode
1913 *
1914 * ====================== SMOOTHES THE VISCOSITY FIELD AT ALL QUADRATURE POINTS ======================
1915 *
1916
1917 *
1918 *
1919 * @code
1920 *   template<int dim>
1921 *   void StokesProblem<dim>::smooth_eta_field(std::vector<bool> failing_cells)
1922 *   {
1923 *   std::cout << " Smoothing viscosity field...\n";
1924 *   unsigned int cell_no = 0;
1925 *   for (typename DoFHandler<dim>::active_cell_iterator cell =
1926 *   dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
1927 *   {
1928 *   if (failing_cells[cell_no])
1929 *   {
1930 *   FEValues<dim> fe_values(fe, quadrature_formula, update_quadrature_points);
1931 *   PointHistory<dim> *local_quadrature_points_history =
1932 *   reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
1933 *  
1934 * @endcode
1935 *
1936 * Currently this algorithm does not permit refinement. To permit refinement, daughter cells of neighbors must be identified
1937 * Find pointers and indices of all cells within certain radius
1938 *
1939 * @code
1940 *   bool find_more_cells = true;
1941 *   std::vector<bool> cell_touched(triangulation.n_active_cells(), false);
1942 *   std::vector< TriaIterator< CellAccessor<dim> > > neighbor_cells;
1943 *   std::vector<int> neighbor_indices;
1944 *   unsigned int start_cell = 0; // Which cell in the neighbor_cells vector to start from
1945 *   int new_cells_found = 0;
1946 *   neighbor_cells.push_back(cell);
1947 *   neighbor_indices.push_back(cell_no);
1948 *   cell_touched[cell_no] = true;
1949 *   while (find_more_cells)
1950 *   {
1951 *   new_cells_found = 0;
1952 *   for (unsigned int i = start_cell; i<neighbor_cells.size(); ++i)
1953 *   {
1954 *   for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
1955 *   {
1956 *   if (!neighbor_cells[i]->face(f)->at_boundary())
1957 *   {
1958 *   int test_cell_no = neighbor_cells[i]->neighbor_index(f);
1959 *   if (!cell_touched[test_cell_no])
1960 *   if (cell->center().distance(neighbor_cells[i]->neighbor(f)->center()) < 2 * system_parameters::smoothing_radius)
1961 *   {
1962 * @endcode
1963 *
1964 * What to do if another nearby cell is found that hasn't been found before
1965 *
1966 * @code
1967 *   neighbor_cells.push_back(neighbor_cells[i]->neighbor(f));
1968 *   neighbor_indices.push_back(test_cell_no);
1969 *   cell_touched[test_cell_no] = true;
1970 *   start_cell++;
1971 *   new_cells_found++;
1972 *   }
1973 *   }
1974 *   }
1975 *   }
1976 *   if (new_cells_found == 0)
1977 *   {
1978 *   find_more_cells = false;
1979 *   }
1980 *   else
1981 *   start_cell = neighbor_cells.size() - new_cells_found;
1982 *   }
1983 *  
1984 *   fe_values.reinit(cell);
1985 * @endcode
1986 *
1987 * Collect the viscosities at nearby quadrature points
1988 *
1989 * @code
1990 *   for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
1991 *   {
1992 *   std::vector<double> nearby_etas_q;
1993 *   for (unsigned int i = 0; i<neighbor_indices.size(); ++i)
1994 *   for (unsigned int j=0; j<quadrature_formula.size(); ++j)
1995 *   {
1996 *   Point<dim> test_q;
1997 *   for (unsigned int d=0; d<dim; ++d)
1998 *   test_q(d) = quad_viscosities[neighbor_indices[i]][j][d];
1999 *   double qq_distance = fe_values.quadrature_point(q).distance(test_q);
2000 *   if (qq_distance < system_parameters::smoothing_radius)
2001 *   nearby_etas_q.push_back(quad_viscosities[neighbor_indices[i]][j][dim]);
2002 *   }
2003 * @endcode
2004 *
2005 * Write smoothed viscosities to quadrature_points_history; simple boxcar function is the smoothing kernel
2006 *
2007 * @code
2008 *   double mean_eta = 0;
2009 *   for (unsigned int l = 0; l<nearby_etas_q.size(); ++l)
2010 *   {
2011 *   mean_eta += nearby_etas_q[l];
2012 *   }
2013 *   mean_eta /= nearby_etas_q.size();
2014 *   local_quadrature_points_history[q].new_eta = mean_eta;
2015 * @endcode
2016 *
2017 * std::cout << local_quadrature_points_history[q].new_eta << " ";
2018 *
2019 * @code
2020 *   }
2021 *   }
2022 *   cell_no++;
2023 *   }
2024 *   }
2025 *  
2026 * @endcode
2027 *
2028 * ====================== SAVE STRESS TENSOR AT QUADRATURE POINTS ======================
2029 *
2030
2031 *
2032 *
2033 * @code
2034 *   template<int dim>
2035 *   void StokesProblem<dim>::update_quadrature_point_history()
2036 *   {
2037 *   std::cout << " Updating stress field...";
2038 *  
2039 *   FEValues<dim> fe_values(fe, quadrature_formula,
2040 *   update_values | update_gradients | update_quadrature_points);
2041 *  
2042 * @endcode
2043 *
2044 * Make the object that will hold the velocity gradients
2045 *
2046 * @code
2047 *   std::vector < std::vector<Tensor<1, dim> >> velocity_grads(quadrature_formula.size(),
2048 *   std::vector < Tensor<1, dim> > (dim + 1));
2049 *   std::vector<Vector<double> > velocities(quadrature_formula.size(),
2050 *   Vector<double>(dim + 1));
2051 *  
2052 *   for (typename DoFHandler<dim>::active_cell_iterator cell =
2053 *   dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
2054 *   {
2055 *   PointHistory<dim> *local_quadrature_points_history =
2056 *   reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
2057 *   Assert(
2058 *   local_quadrature_points_history >= &quadrature_point_history.front(),
2059 *   ExcInternalError());
2060 *   Assert(
2061 *   local_quadrature_points_history < &quadrature_point_history.back(),
2062 *   ExcInternalError());
2063 *  
2064 *   fe_values.reinit(cell);
2065 *   fe_values.get_function_gradients(solution, velocity_grads);
2066 *   fe_values.get_function_values(solution, velocities);
2067 *  
2068 *   for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
2069 *   {
2070 * @endcode
2071 *
2072 * Define the local viscoelastic constants
2073 *
2074 * @code
2075 *   double local_eta_ve = 2
2076 *   / ((1 / local_quadrature_points_history[q].new_eta)
2077 *   + (1 / local_quadrature_points_history[q].G
2078 *   / system_parameters::current_time_interval));
2079 *   double local_chi_ve =
2080 *   1
2081 *   / (1
2082 *   + (local_quadrature_points_history[q].G
2083 *   * system_parameters::current_time_interval
2084 *   / local_quadrature_points_history[q].new_eta));
2085 *  
2086 * @endcode
2087 *
2088 * Compute new stress at each quadrature point
2089 *
2090 * @code
2091 *   SymmetricTensor<2, dim> new_stress;
2092 *   for (unsigned int i = 0; i < dim; ++i)
2093 *   new_stress[i][i] =
2094 *   local_eta_ve * velocity_grads[q][i][i]
2095 *   + local_chi_ve
2096 *   * local_quadrature_points_history[q].old_stress[i][i];
2097 *  
2098 *   for (unsigned int i = 0; i < dim; ++i)
2099 *   for (unsigned int j = i + 1; j < dim; ++j)
2100 *   new_stress[i][j] =
2101 *   local_eta_ve
2102 *   * (velocity_grads[q][i][j]
2103 *   + velocity_grads[q][j][i]) / 2
2104 *   + local_chi_ve
2105 *   * local_quadrature_points_history[q].old_stress[i][j];
2106 *  
2107 * @endcode
2108 *
2109 * Rotate new stress
2110 *
2111 * @code
2112 *   AuxFunctions<dim> rotation_object;
2113 *   const Tensor<2, dim> rotation = rotation_object.get_rotation_matrix(
2114 *   velocity_grads[q]);
2115 *   const SymmetricTensor<2, dim> rotated_new_stress = symmetrize(
2116 *   transpose(rotation)
2117 *   * static_cast<Tensor<2, dim> >(new_stress)
2118 *   * rotation);
2119 *   local_quadrature_points_history[q].old_stress = rotated_new_stress;
2120 *  
2121 * @endcode
2122 *
2123 * For axisymmetric case, make the phi-phi element of stress tensor
2124 *
2125 * @code
2126 *   local_quadrature_points_history[q].old_phiphi_stress =
2127 *   (2 * local_eta_ve * velocities[q](0)
2128 *   / fe_values.quadrature_point(q)[0]
2129 *   + local_chi_ve
2130 *   * local_quadrature_points_history[q].old_phiphi_stress);
2131 *   }
2132 *   }
2133 *   }
2134 *  
2135 * @endcode
2136 *
2137 * ====================== REDEFINE THE TIME INTERVAL FOR THE VISCOUS STEPS ======================
2138 *
2139 * @code
2140 *   template<int dim>
2141 *   void StokesProblem<dim>::update_time_interval()
2142 *   {
2143 *   double move_goal_per_step = system_parameters::initial_disp_target;
2144 *   if (system_parameters::present_timestep > system_parameters::initial_elastic_iterations)
2145 *   {
2146 *   move_goal_per_step = system_parameters::initial_disp_target -
2147 *   ((system_parameters::initial_disp_target - system_parameters::final_disp_target) /
2148 *   system_parameters::total_viscous_steps *
2149 *   (system_parameters::present_timestep - system_parameters::initial_elastic_iterations));
2150 *   }
2151 *  
2152 *   double zero_tolerance = 1e-3;
2153 *   double max_velocity = 0;
2154 *   for (typename DoFHandler<dim>::active_cell_iterator cell =
2155 *   dof_handler.begin_active(); cell != dof_handler.end(); ++cell)// loop over all cells
2156 *   {
2157 *   if (cell->at_boundary())
2158 *   {
2159 *   int zero_faces = 0;
2160 *   for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
2161 *   for (unsigned int i=0; i<dim; ++i)
2162 *   if (fabs(cell->face(f)->center()[i]) < zero_tolerance)
2163 *   zero_faces++;
2164 *   if (zero_faces==0)
2165 *   {
2166 *   for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
2167 *   {
2168 *   Point<dim> vertex_velocity;
2169 *   Point<dim> vertex_position;
2170 *   for (unsigned int d = 0; d < dim; ++d)
2171 *   {
2172 *   vertex_velocity[d] = solution(cell->vertex_dof_index(v, d));
2173 *   vertex_position[d] = cell->vertex(v)[d];
2174 *   }
2175 * @endcode
2176 *
2177 * velocity to be evaluated is the radial component of a surface vertex
2178 *
2179 * @code
2180 *   double local_velocity = 0;
2181 *   for (unsigned int d = 0; d < dim; ++d)
2182 *   {
2183 *   local_velocity += vertex_velocity[d] * vertex_position [d];
2184 *   }
2185 *   local_velocity /= std::sqrt( vertex_position.square() );
2186 *   if (local_velocity < 0)
2187 *   local_velocity *= -1;
2188 *   if (local_velocity > max_velocity)
2189 *   {
2190 *   max_velocity = local_velocity;
2191 *   }
2192 *   }
2193 *   }
2194 *   }
2195 *   }
2196 * @endcode
2197 *
2198 * NOTE: It is possible for this time interval to be very different from that used in the viscoelasticity calculation.
2199 *
2200 * @code
2201 *   system_parameters::current_time_interval = move_goal_per_step / max_velocity;
2202 *   double step_time_yr = system_parameters::current_time_interval / SECSINYEAR;
2203 *   std::cout << "Timestep interval changed to: "
2204 *   << step_time_yr
2205 *   << " years\n";
2206 *   }
2207 *  
2208 * @endcode
2209 *
2210 * ====================== MOVE MESH ======================
2211 *
2212
2213 *
2214 *
2215 * @code
2216 *   template<int dim>
2217 *   void StokesProblem<dim>::move_mesh()
2218 *   {
2219 *  
2220 *   std::cout << "\n" << " Moving mesh...\n";
2221 *   std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
2222 *   for (typename DoFHandler<dim>::active_cell_iterator cell =
2223 *   dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
2224 *   for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
2225 *   if (vertex_touched[cell->vertex_index(v)] == false)
2226 *   {
2227 *   vertex_touched[cell->vertex_index(v)] = true;
2228 *  
2229 *   Point<dim> vertex_displacement;
2230 *   for (unsigned int d = 0; d < dim; ++d)
2231 *   vertex_displacement[d] = solution(
2232 *   cell->vertex_dof_index(v, d));
2233 *   cell->vertex(v) += vertex_displacement
2234 *   * system_parameters::current_time_interval;
2235 *   }
2236 *   }
2237 *  
2238 * @endcode
2239 *
2240 * ====================== WRITE MESH TO FILE ======================
2241 *
2242
2243 *
2244 *
2245 * @code
2246 *   template<int dim>
2247 *   void StokesProblem<dim>::write_mesh()
2248 *   {
2249 * @endcode
2250 *
2251 * output mesh in ucd
2252 *
2253 * @code
2254 *   std::ostringstream initial_mesh_file;
2255 *   initial_mesh_file << system_parameters::output_folder << "/time" <<
2256 *   Utilities::int_to_string(system_parameters::present_timestep, 2) <<
2257 *   "_mesh.inp";
2258 *   std::ofstream out_ucd (initial_mesh_file.str().c_str());
2259 *   GridOut grid_out;
2260 *   grid_out.write_ucd (triangulation, out_ucd);
2261 *   }
2262 *  
2263 * @endcode
2264 *
2265 * ====================== FIT ELLIPSE TO SURFACE AND WRITE RADII TO FILE ======================
2266 *
2267
2268 *
2269 *
2270 * @code
2271 *   template<int dim>
2272 *   void StokesProblem<dim>::do_ellipse_fits()
2273 *   {
2274 *   std::ostringstream ellipses_filename;
2275 *   ellipses_filename << system_parameters::output_folder << "/ellipse_fits.txt";
2276 * @endcode
2277 *
2278 * Find ellipsoidal axes for all layers
2279 *
2280 * @code
2281 *   std::vector<double> ellipse_axes(0);
2282 * @endcode
2283 *
2284 * compute fit to boundary 0, 1, 2 ...
2285 *
2286 * @code
2287 *   std::cout << endl;
2288 *   for (unsigned int i = 0; i<system_parameters::sizeof_material_id; ++i)
2289 *   {
2290 *   ellipsoid.compute_fit(ellipse_axes, system_parameters::material_id[i]);
2291 *   system_parameters::q_axes.push_back(ellipse_axes[0]);
2292 *   system_parameters::p_axes.push_back(ellipse_axes[1]);
2293 *  
2294 *   std::cout << "a_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[0]
2295 *   << " " << " c_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[1] << std::endl;
2296 *   ellipse_axes.clear();
2297 *  
2298 *   std::ofstream fout_ellipses(ellipses_filename.str().c_str(), std::ios::app);
2299 *   fout_ellipses << system_parameters::present_timestep << " a_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[0]
2300 *   << " " << " c_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[1] << endl;
2301 *   fout_ellipses.close();
2302 *   }
2303 *   }
2304 *  
2305 * @endcode
2306 *
2307 * ====================== APPEND LINE TO PHYSICAL_TIMES.TXT FILE WITH STEP NUMBER, PHYSICAL TIME, AND # PLASTIC ITERATIONS ======================
2308 *
2309
2310 *
2311 *
2312 * @code
2313 *   template<int dim>
2314 *   void StokesProblem<dim>::append_physical_times(int max_plastic)
2315 *   {
2316 *   std::ostringstream times_filename;
2317 *   times_filename << system_parameters::output_folder << "/physical_times.txt";
2318 *   std::ofstream fout_times(times_filename.str().c_str(), std::ios::app);
2319 *   fout_times << system_parameters::present_timestep << " "
2320 *   << system_parameters::present_time/SECSINYEAR << " "
2321 *   << max_plastic << "\n";
2322 * @endcode
2323 *
2324 * << system_parameters::q_axes[0] << " " << system_parameters::p_axes[0] << " "
2325 * << system_parameters::q_axes[1] << " " << system_parameters::p_axes[1] << "\n";
2326 *
2327 * @code
2328 *   fout_times.close();
2329 *   }
2330 *  
2331 * @endcode
2332 *
2333 * ====================== WRITE VERTICES TO FILE ======================
2334 *
2335
2336 *
2337 *
2338 * @code
2339 *   template<int dim>
2340 *   void StokesProblem<dim>::write_vertices(unsigned char boundary_that_we_need)
2341 *   {
2342 *   std::ostringstream vertices_output;
2343 *   vertices_output << system_parameters::output_folder << "/time" <<
2344 *   Utilities::int_to_string(system_parameters::present_timestep, 2) << "_" <<
2345 *   Utilities::int_to_string(boundary_that_we_need, 2) <<
2346 *   "_surface.txt";
2347 *   std::ofstream fout_final_vertices(vertices_output.str().c_str());
2348 *   fout_final_vertices.close();
2349 *  
2350 *   std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
2351 *  
2352 *   if (boundary_that_we_need == 0)
2353 *   {
2354 * @endcode
2355 *
2356 * Figure out if the vertex is on the boundary of the domain
2357 *
2358 * @code
2359 *   for (typename Triangulation<dim>::active_cell_iterator cell =
2360 *   triangulation.begin_active(); cell != triangulation.end(); ++cell)
2361 *   for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
2362 *   {
2363 *   unsigned char boundary_ids = cell->face(f)->boundary_id();
2364 *   if (boundary_ids == boundary_that_we_need)
2365 *   {
2366 *   for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
2367 *   if (vertex_touched[cell->face(f)->vertex_index(v)] == false)
2368 *   {
2369 *   vertex_touched[cell->face(f)->vertex_index(v)] = true;
2370 *   std::ofstream fout_final_vertices(vertices_output.str().c_str(), std::ios::app);
2371 *   fout_final_vertices << cell->face(f)->vertex(v) << "\n";
2372 *   fout_final_vertices.close();
2373 *   }
2374 *   }
2375 *   }
2376 *   }
2377 *   else
2378 *   {
2379 * @endcode
2380 *
2381 * Figure out if the vertex is on an internal boundary
2382 *
2383 * @code
2384 *   for (typename Triangulation<dim>::active_cell_iterator cell =
2385 *   triangulation.begin_active(); cell != triangulation.end(); ++cell)
2386 *   for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
2387 *   {
2388 *   if (cell->neighbor(f) != triangulation.end())
2389 *   {
2390 *   if (cell->material_id() != cell->neighbor(f)->material_id()) //finds face is at internal boundary
2391 *   {
2392 *   int high_mat_id = std::max(cell->material_id(),
2393 *   cell->neighbor(f)->material_id());
2394 *   if (high_mat_id == boundary_that_we_need) //finds faces at the correct internal boundary
2395 *   {
2396 *   for (unsigned int v = 0;
2397 *   v < GeometryInfo<dim>::vertices_per_face;
2398 *   ++v)
2399 *   if (vertex_touched[cell->face(f)->vertex_index(
2400 *   v)] == false)
2401 *   {
2402 *   vertex_touched[cell->face(f)->vertex_index(
2403 *   v)] = true;
2404 *   std::ofstream fout_final_vertices(vertices_output.str().c_str(), std::ios::app);
2405 *   fout_final_vertices << cell->face(f)->vertex(v) << "\n";
2406 *   fout_final_vertices.close();
2407 *   }
2408 *   }
2409 *   }
2410 *   }
2411 *   }
2412 *   }
2413 *   }
2414 *  
2415 * @endcode
2416 *
2417 * ====================== SETUP INITIAL MESH ======================
2418 *
2419
2420 *
2421 *
2422 * @code
2423 *   template<int dim>
2424 *   void StokesProblem<dim>::setup_initial_mesh()
2425 *   {
2426 *   GridIn<dim> grid_in;
2427 *   grid_in.attach_triangulation(triangulation);
2428 *   std::ifstream mesh_stream(system_parameters::mesh_filename,
2429 *   std::ifstream::in);
2430 *   grid_in.read_ucd(mesh_stream);
2431 *  
2432 * @endcode
2433 *
2434 * output initial mesh in eps
2435 *
2436 * @code
2437 *   std::ostringstream initial_mesh_file;
2438 *   initial_mesh_file << system_parameters::output_folder << "/initial_mesh.eps";
2439 *   std::ofstream out_eps (initial_mesh_file.str().c_str());
2440 *   GridOut grid_out;
2441 *   grid_out.write_eps (triangulation, out_eps);
2442 *   out_eps.close();
2443 *  
2444 * @endcode
2445 *
2446 * set boundary ids
2447 * boundary indicator 0 is outer free surface; 1, 2, 3 ... is boundary between layers, 99 is flat boundaries
2448 *
2449 * @code
2450 *   typename Triangulation<dim>::active_cell_iterator
2451 *   cell=triangulation.begin_active(), endc=triangulation.end();
2452 *  
2453 *   unsigned int how_many; // how many components away from cardinal planes
2454 *  
2455 *   std::ostringstream boundaries_file;
2456 *   boundaries_file << system_parameters::output_folder << "/boundaries.txt";
2457 *   std::ofstream fout_boundaries(boundaries_file.str().c_str());
2458 *   fout_boundaries.close();
2459 *  
2460 *   double zero_tolerance = 1e-3;
2461 *   for (; cell != endc; ++cell) // loop over all cells
2462 *   {
2463 *   for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f) // loop over all vertices
2464 *   {
2465 *   if (cell->face(f)->at_boundary())
2466 *   {
2467 * @endcode
2468 *
2469 * print boundary
2470 *
2471 * @code
2472 *   std::ofstream fout_boundaries(boundaries_file.str().c_str(), std::ios::app);
2473 *   fout_boundaries << cell->face(f)->center()[0] << " " << cell->face(f)->center()[1]<< "\n";
2474 *   fout_boundaries.close();
2475 *  
2476 *   how_many = 0;
2477 *   for (unsigned int i=0; i<dim; ++i)
2478 *   if (fabs(cell->face(f)->center()[i]) > zero_tolerance)
2479 *   how_many++;
2480 *   if (how_many==dim)
2481 *   cell->face(f)->set_all_boundary_ids(0); // if face center coordinates > zero_tol, set bnry indicators to 0
2482 *   else
2483 *   cell->face(f)->set_all_boundary_ids(99);
2484 *   }
2485 *   }
2486 *   }
2487 *  
2488 *   std::ostringstream ellipses_filename;
2489 *   ellipses_filename << system_parameters::output_folder << "/ellipse_fits.txt";
2490 *   std::ofstream fout_ellipses(ellipses_filename.str().c_str());
2491 *   fout_ellipses.close();
2492 *  
2493 * @endcode
2494 *
2495 * Find ellipsoidal axes for all layers
2496 *
2497 * @code
2498 *   std::vector<double> ellipse_axes(0);
2499 * @endcode
2500 *
2501 * compute fit to boundary 0, 1, 2 ...
2502 *
2503 * @code
2504 *   std::cout << endl;
2505 *   for (unsigned int i = 0; i<system_parameters::sizeof_material_id; ++i)
2506 *   {
2507 *   ellipsoid.compute_fit(ellipse_axes, system_parameters::material_id[i]);
2508 *   system_parameters::q_axes.push_back(ellipse_axes[0]);
2509 *   system_parameters::p_axes.push_back(ellipse_axes[1]);
2510 *  
2511 *   std::cout << "a_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[0]
2512 *   << " " << " c_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[1] << std::endl;
2513 *   ellipse_axes.clear();
2514 *  
2515 *   std::ofstream fout_ellipses(ellipses_filename.str().c_str(), std::ios::app);
2516 *   fout_ellipses << system_parameters::present_timestep << " a_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[0]
2517 *   << " " << " c_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[1] << endl;
2518 *   fout_ellipses.close();
2519 *   }
2520 *  
2521 *   triangulation.refine_global(system_parameters::global_refinement);
2522 *  
2523 *  
2524 * @endcode
2525 *
2526 * refines crustal region
2527 *
2528 * @code
2529 *   if (system_parameters::crustal_refinement != 0)
2530 *   {
2531 *   double a = system_parameters::q_axes[0] - system_parameters::crust_refine_region;
2532 *   double b = system_parameters::p_axes[0] - system_parameters::crust_refine_region;
2533 *  
2534 *  
2535 *   for (unsigned int step = 0;
2536 *   step < system_parameters::crustal_refinement; ++step)
2537 *   {
2538 *   typename ::Triangulation<dim>::active_cell_iterator cell =
2539 *   triangulation.begin_active(), endc = triangulation.end();
2540 *   for (; cell != endc; ++cell)
2541 *   for (unsigned int v = 0;
2542 *   v < GeometryInfo<dim>::vertices_per_cell; ++v)
2543 *   {
2544 *   Point<dim> current_vertex = cell->vertex(v);
2545 *  
2546 *   const double x_coord = current_vertex.operator()(0);
2547 *   const double y_coord = current_vertex.operator()(1);
2548 *   double expected_z = -1;
2549 *  
2550 *   if ((x_coord - a) < -1e-10)
2551 *   expected_z = b
2552 *   * std::sqrt(1 - (x_coord * x_coord / a / a));
2553 *  
2554 *   if (y_coord >= expected_z)
2555 *   {
2556 *   cell->set_refine_flag();
2557 *   break;
2558 *   }
2559 *   }
2560 *   triangulation.execute_coarsening_and_refinement();
2561 *   }
2562 *   }
2563 *  
2564 *  
2565 * @endcode
2566 *
2567 * output initial mesh in eps
2568 *
2569 * @code
2570 *   std::ostringstream refined_mesh_file;
2571 *   refined_mesh_file << system_parameters::output_folder << "/refined_mesh.eps";
2572 *   std::ofstream out_eps_refined (refined_mesh_file.str().c_str());
2573 *   GridOut grid_out_refined;
2574 *   grid_out_refined.write_eps (triangulation, out_eps_refined);
2575 *   out_eps_refined.close();
2576 *   write_vertices(0);
2577 *   write_vertices(1);
2578 *   write_mesh();
2579 *   }
2580 *  
2581 * @endcode
2582 *
2583 * ====================== REFINE MESH ======================
2584 *
2585
2586 *
2587 *
2588 * @code
2589 *   template<int dim>
2590 *   void StokesProblem<dim>::refine_mesh()
2591 *   {
2592 *   using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
2593 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
2594 *  
2595 *   std::vector<bool> component_mask(dim + 1, false);
2596 *   component_mask[dim] = true;
2597 *   KellyErrorEstimator<dim>::estimate(dof_handler, QGauss<dim - 1>(degree + 1),
2598 *   FunctionMap(), solution,
2599 *   estimated_error_per_cell, component_mask);
2600 *  
2601 *   GridRefinement::refine_and_coarsen_fixed_number(triangulation,
2602 *   estimated_error_per_cell, 0.3, 0.0);
2603 *   triangulation.execute_coarsening_and_refinement();
2604 *   }
2605 *  
2606 * @endcode
2607 *
2608 * ====================== SET UP THE DATA STRUCTURES TO REMEMBER STRESS FIELD ======================
2609 *
2610 * @code
2611 *   template<int dim>
2612 *   void StokesProblem<dim>::setup_quadrature_point_history()
2613 *   {
2614 *   unsigned int our_cells = 0;
2615 *   for (typename Triangulation<dim>::active_cell_iterator cell =
2616 *   triangulation.begin_active(); cell != triangulation.end(); ++cell)
2617 *   ++our_cells;
2618 *  
2619 *   triangulation.clear_user_data();
2620 *  
2621 *   quadrature_point_history.resize(our_cells * quadrature_formula.size());
2622 *  
2623 *   unsigned int history_index = 0;
2624 *   for (typename Triangulation<dim>::active_cell_iterator cell =
2625 *   triangulation.begin_active(); cell != triangulation.end(); ++cell)
2626 *   {
2627 *   cell->set_user_pointer(&quadrature_point_history[history_index]);
2628 *   history_index += quadrature_formula.size();
2629 *   }
2630 *  
2631 *   Assert(history_index == quadrature_point_history.size(), ExcInternalError());
2632 *   }
2633 *  
2634 * @endcode
2635 *
2636 * ====================== DOES ELASTIC STEPS ======================
2637 *
2638 * @code
2639 *   template<int dim>
2640 *   void StokesProblem<dim>::do_elastic_steps()
2641 *   {
2642 *   unsigned int elastic_iteration = 0;
2643 *  
2644 *   while (elastic_iteration < system_parameters::initial_elastic_iterations)
2645 *   {
2646 *  
2647 *   std::cout << "\n\nElastic iteration " << elastic_iteration
2648 *   << "\n";
2649 *   setup_dofs();
2650 *  
2651 *   if (system_parameters::present_timestep == 0)
2652 *   initialize_eta_and_G();
2653 *  
2654 *   if (elastic_iteration == 0)
2655 *   system_parameters::current_time_interval =
2656 *   system_parameters::viscous_time; //This is the time interval needed in assembling the problem
2657 *  
2658 *   std::cout << " Assembling..." << std::endl << std::flush;
2659 *   assemble_system();
2660 *  
2661 *   std::cout << " Solving..." << std::flush;
2662 *   solve();
2663 *  
2664 *   output_results();
2665 *   update_quadrature_point_history();
2666 *  
2667 *   append_physical_times(0);
2668 *   elastic_iteration++;
2669 *   system_parameters::present_timestep++;
2670 *   do_ellipse_fits();
2671 *   write_vertices(0);
2672 *   write_vertices(1);
2673 *   write_mesh();
2674 *   update_time_interval();
2675 *   }
2676 *   }
2677 *  
2678 * @endcode
2679 *
2680 * ====================== DO A SINGLE VISCOELASTOPLASTIC TIMESTEP ======================
2681 *
2682 * @code
2683 *   template<int dim>
2684 *   void StokesProblem<dim>::do_flow_step()
2685 *   {
2686 *   plastic_iteration = 0;
2687 *   while (plastic_iteration < system_parameters::max_plastic_iterations)
2688 *   {
2689 *   if (system_parameters::continue_plastic_iterations == true)
2690 *   {
2691 *   std::cout << "Plasticity iteration " << plastic_iteration << "\n";
2692 *   setup_dofs();
2693 *  
2694 *   std::cout << " Assembling..." << std::endl << std::flush;
2695 *   assemble_system();
2696 *  
2697 *   std::cout << " Solving..." << std::flush;
2698 *   solve();
2699 *  
2700 *   output_results();
2701 *   solution_stesses();
2702 *  
2703 *   if (system_parameters::continue_plastic_iterations == false)
2704 *   break;
2705 *  
2706 *   plastic_iteration++;
2707 *   }
2708 *   }
2709 *   }
2710 *  
2711 * @endcode
2712 *
2713 * ====================== RUN ======================
2714 *
2715
2716 *
2717 *
2718 * @code
2719 *   template<int dim>
2720 *   void StokesProblem<dim>::run()
2721 *   {
2722 * @endcode
2723 *
2724 * Sets up mesh and data structure for viscosity and stress at quadrature points
2725 *
2726 * @code
2727 *   setup_initial_mesh();
2728 *   setup_quadrature_point_history();
2729 *  
2730 * @endcode
2731 *
2732 * Makes the physical_times.txt file
2733 *
2734 * @code
2735 *   std::ostringstream times_filename;
2736 *   times_filename << system_parameters::output_folder << "/physical_times.txt";
2737 *   std::ofstream fout_times(times_filename.str().c_str());
2738 *   fout_times.close();
2739 *  
2740 * @endcode
2741 *
2742 * Computes elastic timesteps
2743 *
2744 * @code
2745 *   do_elastic_steps();
2746 * @endcode
2747 *
2748 * Computes viscous timesteps
2749 *
2750 * @code
2751 *   unsigned int VEPstep = 0;
2752 *   while (system_parameters::present_timestep
2753 *   < (system_parameters::initial_elastic_iterations
2754 *   + system_parameters::total_viscous_steps))
2755 *   {
2756 *  
2757 *   if (system_parameters::continue_plastic_iterations == false)
2758 *   system_parameters::continue_plastic_iterations = true;
2759 *   std::cout << "\n\nViscoelastoplastic iteration " << VEPstep << "\n\n";
2760 * @endcode
2761 *
2762 * Computes plasticity
2763 *
2764 * @code
2765 *   do_flow_step();
2766 *   update_quadrature_point_history();
2767 *   move_mesh();
2768 *   append_physical_times(plastic_iteration);
2769 *   system_parameters::present_timestep++;
2770 *   system_parameters::present_time = system_parameters::present_time + system_parameters::current_time_interval;
2771 *   do_ellipse_fits();
2772 *   write_vertices(0);
2773 *   write_vertices(1);
2774 *   write_mesh();
2775 *   VEPstep++;
2776 *   }
2777 *   append_physical_times(-1);
2778 *   }
2779 *  
2780 *   }
2781 *  
2782 * @endcode
2783 *
2784 * ====================== MAIN ======================
2785 *
2786
2787 *
2788 *
2789 * @code
2790 *   int main(int argc, char *argv[])
2791 *   {
2792 *  
2793 * @endcode
2794 *
2795 * output program name
2796 *
2797 * @code
2798 *   std::cout << "Running: " << argv[0] << std::endl;
2799 *  
2800 *   char *cfg_filename = new char[120];
2801 *  
2802 *   if (argc == 1) // if no input parameters (as if launched from eclipse)
2803 *   {
2804 *   std::strcpy(cfg_filename,"config/sample_CeresFE_config.cfg");
2805 *   }
2806 *   else
2807 *   std::strcpy(cfg_filename,argv[1]);
2808 *  
2809 *   try
2810 *   {
2811 *   using namespace dealii;
2812 *   using namespace Step22;
2813 *   config_in cfg(cfg_filename);
2814 *  
2815 *   std::clock_t t1;
2816 *   std::clock_t t2;
2817 *   t1 = std::clock();
2818 *  
2819 *   deallog.depth_console(0);
2820 *  
2821 *   StokesProblem<2> flow_problem(1);
2822 *   flow_problem.run();
2823 *  
2824 *   std::cout << std::endl << "\a";
2825 *  
2826 *   t2 = std::clock();
2827 *   float diff (((float)t2 - (float)t1) / (float)CLOCKS_PER_SEC);
2828 *   std::cout << "\n Program run in: " << diff << " seconds" << endl;
2829 *   }
2830 *   catch (std::exception &exc)
2831 *   {
2832 *   std::cerr << std::endl << std::endl
2833 *   << "----------------------------------------------------"
2834 *   << std::endl;
2835 *   std::cerr << "Exception on processing: " << std::endl << exc.what()
2836 *   << std::endl << "Aborting!" << std::endl
2837 *   << "----------------------------------------------------"
2838 *   << std::endl;
2839 *  
2840 *   return 1;
2841 *   }
2842 *   catch (...)
2843 *   {
2844 *   std::cerr << std::endl << std::endl
2845 *   << "----------------------------------------------------"
2846 *   << std::endl;
2847 *   std::cerr << "Unknown exception!" << std::endl << "Aborting!"
2848 *   << std::endl
2849 *   << "----------------------------------------------------"
2850 *   << std::endl;
2851 *   return 1;
2852 *   }
2853 *  
2854 *   return 0;
2855 *   }
2856 * @endcode
2857
2858
2859<a name="ann-support_code/config_in.h"></a>
2860<h1>Annotated version of support_code/config_in.h</h1>
2861 *
2862 *
2863 *
2864 *
2865 * @code
2866 *   /*
2867 *   * config_in.h
2868 *   *
2869 *   * Created on: Aug 17, 2015
2870 *   * Author: antonermakov
2871 *   */
2872 *  
2873 *  
2874 *   #include <iostream>
2875 *   #include <iomanip>
2876 *   #include <cstdlib>
2877 *   #include <string>
2878 *   #include <libconfig.h++>
2879 *  
2880 *   #include "local_math.h"
2881 *  
2882 *   using namespace std;
2883 *   using namespace libconfig;
2884 *  
2885 *   namespace Step22
2886 *   {
2887 *   namespace system_parameters
2888 *   {
2889 *  
2890 * @endcode
2891 *
2892 * Mesh file name
2893 *
2894 * @code
2895 *   string mesh_filename;
2896 *   string output_folder;
2897 *  
2898 * @endcode
2899 *
2900 * Body parameters
2901 *
2902 * @code
2903 *   double r_mean;
2904 *   double period;
2905 *   double omegasquared;
2906 *   double beta;
2907 *   double intercept;
2908 *  
2909 * @endcode
2910 *
2911 * Rheology parameters
2912 *
2913 * @code
2914 *   vector<double> depths_eta;
2915 *   vector<double> eta_kinks;
2916 *   vector<double> depths_rho;
2917 *   vector<double> rho;
2918 *   vector<int> material_id;
2919 *   vector<double> G;
2920 *  
2921 *   double eta_ceiling;
2922 *   double eta_floor;
2923 *   double eta_Ea;
2924 *   bool lat_dependence;
2925 *  
2926 *   unsigned int sizeof_depths_eta;
2927 *   unsigned int sizeof_depths_rho;
2928 *   unsigned int sizeof_rho;
2929 *   unsigned int sizeof_eta_kinks;
2930 *   unsigned int sizeof_material_id;
2931 *   unsigned int sizeof_G;
2932 *  
2933 *   double pressure_scale;
2934 *   double q;
2935 *   bool cylindrical;
2936 *   bool continue_plastic_iterations;
2937 *  
2938 * @endcode
2939 *
2940 * plasticity variables
2941 *
2942 * @code
2943 *   bool plasticity_on;
2944 *   unsigned int failure_criterion;
2945 *   unsigned int max_plastic_iterations;
2946 *   double smoothing_radius;
2947 *  
2948 * @endcode
2949 *
2950 * viscoelasticity variables
2951 *
2952 * @code
2953 *   unsigned int initial_elastic_iterations;
2954 *   double elastic_time;
2955 *   double viscous_time;
2956 *   double initial_disp_target;
2957 *   double final_disp_target;
2958 *   double current_time_interval;
2959 *  
2960 * @endcode
2961 *
2962 * mesh refinement variables
2963 *
2964 * @code
2965 *   unsigned int global_refinement;
2966 *   unsigned int small_r_refinement;
2967 *   unsigned int crustal_refinement;
2968 *   double crust_refine_region;
2969 *   unsigned int surface_refinement;
2970 *  
2971 * @endcode
2972 *
2973 * solver variables
2974 *
2975 * @code
2976 *   int iteration_coefficient;
2977 *   double tolerance_coefficient;
2978 *  
2979 * @endcode
2980 *
2981 * time step variables
2982 *
2983 * @code
2984 *   double present_time;
2985 *   unsigned int present_timestep;
2986 *   unsigned int total_viscous_steps;
2987 *  
2988 *  
2989 * @endcode
2990 *
2991 * ellipse axes
2992 *
2993 * @code
2994 *   vector<double> q_axes;
2995 *   vector<double> p_axes;
2996 *  
2997 *   }
2998 *  
2999 *   class config_in
3000 *   {
3001 *   public:
3002 *   config_in(char *);
3003 *  
3004 *   private:
3005 *   void write_config();
3006 *   };
3007 *  
3008 *   void config_in::write_config()
3009 *   {
3010 *   std::ostringstream config_parameters;
3011 *   config_parameters << system_parameters::output_folder << "/run_parameters.txt";
3012 *   std::ofstream fout_config(config_parameters.str().c_str());
3013 *  
3014 * @endcode
3015 *
3016 * mesh filename
3017 *
3018 * @code
3019 *   fout_config << "mesh filename: " << system_parameters::mesh_filename << endl << endl;
3020 *  
3021 * @endcode
3022 *
3023 * body parameters
3024 *
3025 * @code
3026 *   fout_config << "r_mean = " << system_parameters::r_mean << endl;
3027 *   fout_config << "period = " << system_parameters::period << endl;
3028 *   fout_config << "omegasquared = " << system_parameters::omegasquared << endl;
3029 *   fout_config << "beta = " << system_parameters::beta << endl;
3030 *   fout_config << "intercept = " << system_parameters::intercept << endl;
3031 *  
3032 * @endcode
3033 *
3034 * rheology parameters
3035 *
3036
3037 *
3038 *
3039 * @code
3040 *   for (unsigned int i=0; i<system_parameters::sizeof_depths_eta; ++i)
3041 *   fout_config << "depths_eta[" << i << "] = " << system_parameters::depths_eta[i] << endl;
3042 *  
3043 *   for (unsigned int i=0; i<system_parameters::sizeof_eta_kinks; ++i)
3044 *   fout_config << "eta_kinks[" << i << "] = " << system_parameters::eta_kinks[i] << endl;
3045 *  
3046 *   for (unsigned int i=0; i<system_parameters::sizeof_depths_rho; ++i)
3047 *   fout_config << "depths_rho[" << i << "] = " << system_parameters::depths_rho[i] << endl;
3048 *  
3049 *   for (unsigned int i=0; i<system_parameters::sizeof_rho; ++i)
3050 *   fout_config << "rho[" << i << "] = " << system_parameters::rho[i] << endl;
3051 *  
3052 *   for (unsigned int i=0; i<system_parameters::sizeof_material_id; ++i)
3053 *   fout_config << "material_id[" << i << "] = " << system_parameters::material_id[i] << endl;
3054 *  
3055 *   for (unsigned int i=0; i<system_parameters::sizeof_G; ++i)
3056 *   fout_config << "G[" << i << "] = " << system_parameters::G[i] << endl;
3057 *  
3058 *   fout_config << "eta_ceiling = " << system_parameters::eta_ceiling << endl;
3059 *   fout_config << "eta_floor = " << system_parameters::eta_floor << endl;
3060 *   fout_config << "eta_Ea = " << system_parameters::eta_Ea << endl;
3061 *   fout_config << "lat_dependence = " << system_parameters::lat_dependence << endl;
3062 *   fout_config << "pressure_scale = " << system_parameters::pressure_scale << endl;
3063 *   fout_config << "q = " << system_parameters::q << endl;
3064 *   fout_config << "cylindrical = " << system_parameters::cylindrical << endl;
3065 *   fout_config << "continue_plastic_iterations = " << system_parameters::continue_plastic_iterations << endl;
3066 *  
3067 * @endcode
3068 *
3069 * Plasticity parameters
3070 *
3071 * @code
3072 *   fout_config << "plasticity_on = " << system_parameters::plasticity_on << endl;
3073 *   fout_config << "failure_criterion = " << system_parameters::failure_criterion << endl;
3074 *   fout_config << "max_plastic_iterations = " << system_parameters::max_plastic_iterations << endl;
3075 *   fout_config << "smoothing_radius = " << system_parameters::smoothing_radius << endl;
3076 *  
3077 * @endcode
3078 *
3079 * Viscoelasticity parameters
3080 *
3081 * @code
3082 *   fout_config << "initial_elastic_iterations = " << system_parameters::initial_elastic_iterations << endl;
3083 *   fout_config << "elastic_time = " << system_parameters::elastic_time << endl;
3084 *   fout_config << "viscous_time = " << system_parameters::viscous_time << endl;
3085 *   fout_config << "initial_disp_target = " << system_parameters::initial_disp_target << endl;
3086 *   fout_config << "final_disp_target = " << system_parameters::final_disp_target << endl;
3087 *   fout_config << "current_time_interval = " << system_parameters::current_time_interval << endl;
3088 *  
3089 * @endcode
3090 *
3091 * Mesh refinement parameters
3092 *
3093 * @code
3094 *   fout_config << "global_refinement = " << system_parameters::global_refinement << endl;
3095 *   fout_config << "small_r_refinement = " << system_parameters::small_r_refinement << endl;
3096 *   fout_config << "crustal_refinement = " << system_parameters::crustal_refinement << endl;
3097 *   fout_config << "crust_refine_region = " << system_parameters::crust_refine_region << endl;
3098 *   fout_config << "surface_refinement = " << system_parameters::surface_refinement << endl;
3099 *  
3100 * @endcode
3101 *
3102 * Solver parameters
3103 *
3104 * @code
3105 *   fout_config << "iteration_coefficient = " << system_parameters::iteration_coefficient << endl;
3106 *   fout_config << "tolerance_coefficient = " << system_parameters::tolerance_coefficient << endl;
3107 *  
3108 * @endcode
3109 *
3110 * Time step parameters
3111 *
3112 * @code
3113 *   fout_config << "present_time = " << system_parameters::present_time << endl;
3114 *   fout_config << "present_timestep = " << system_parameters::present_timestep << endl;
3115 *   fout_config << "total_viscous_steps = " << system_parameters::total_viscous_steps << endl;
3116 *  
3117 *   fout_config.close();
3118 *   }
3119 *  
3120 *   config_in::config_in(char *filename)
3121 *   {
3122 *  
3123 * @endcode
3124 *
3125 * This example reads the configuration file 'example.cfg' and displays
3126 * some of its contents.
3127 *
3128
3129 *
3130 *
3131 * @code
3132 *   Config cfg;
3133 *  
3134 * @endcode
3135 *
3136 * Read the file. If there is an error, report it and exit.
3137 *
3138 * @code
3139 *   try
3140 *   {
3141 *   cfg.readFile(filename);
3142 *   }
3143 *   catch (const FileIOException &fioex)
3144 *   {
3145 *   std::cerr << "I/O error while reading file:" << filename << std::endl;
3146 *   }
3147 *   catch (const ParseException &pex)
3148 *   {
3149 *   std::cerr << "Parse error at " << pex.getFile() << ":" << pex.getLine()
3150 *   << " - " << pex.getError() << std::endl;
3151 *   }
3152 *  
3153 * @endcode
3154 *
3155 * Get mesh name.
3156 *
3157 * @code
3158 *   try
3159 *   {
3160 *   string msh = cfg.lookup("mesh_filename");
3161 *   system_parameters::mesh_filename = msh;
3162 *   }
3163 *   catch (const SettingNotFoundException &nfex)
3164 *   {
3165 *   cerr << "No 'mesh_filename' setting in configuration file." << endl;
3166 *   }
3167 *  
3168 * @endcode
3169 *
3170 * get output folder
3171 *
3172
3173 *
3174 *
3175 * @code
3176 *   try
3177 *   {
3178 *   string output = cfg.lookup("output_folder");
3179 *   system_parameters::output_folder = output;
3180 *  
3181 *   std::cout << "Writing to folder: " << output << endl;
3182 *   }
3183 *   catch (const SettingNotFoundException &nfex)
3184 *   {
3185 *   cerr << "No 'output_folder' setting in configuration file." << endl;
3186 *   }
3187 *  
3188 * @endcode
3189 *
3190 * get radii
3191 *
3192
3193 *
3194 *
3195 * @code
3196 *   const Setting &root = cfg.getRoot();
3197 *  
3198 *  
3199 * @endcode
3200 *
3201 * get body parameters
3202 *
3203 * @code
3204 *   try
3205 *   {
3206 *   const Setting &body_parameters = root["body_parameters"];
3207 *  
3208 *   body_parameters.lookupValue("period", system_parameters::period);
3209 *   system_parameters::omegasquared = pow(TWOPI / 3600.0 / system_parameters::period, 2.0);
3210 *   body_parameters.lookupValue("r_mean", system_parameters::r_mean);
3211 *   body_parameters.lookupValue("beta", system_parameters::beta);
3212 *   body_parameters.lookupValue("intercept", system_parameters::intercept);
3213 *   }
3214 *   catch (const SettingNotFoundException &nfex)
3215 *   {
3216 *   cerr << "We've got a problem in the body parameters block" << endl;
3217 *  
3218 *   }
3219 *  
3220 * @endcode
3221 *
3222 * Rheology parameters
3223 *
3224
3225 *
3226 *
3227 * @code
3228 *   try
3229 *   {
3230 * @endcode
3231 *
3232 * get depths_eta ---------------------
3233 *
3234 * @code
3235 *   const Setting &set_depths_eta = cfg.lookup("rheology_parameters.depths_eta");
3236 *  
3237 *   unsigned int ndepths_eta = set_depths_eta.getLength();
3238 *   system_parameters::sizeof_depths_eta = ndepths_eta;
3239 *  
3240 *   for (unsigned int i=0; i<ndepths_eta; ++i)
3241 *   {
3242 *   system_parameters::depths_eta.push_back(set_depths_eta[i]);
3243 *   cout << "depth_eta[" << i << "] = " << system_parameters::depths_eta[i] << endl;
3244 *   }
3245 *  
3246 * @endcode
3247 *
3248 * get eta_kinks -------------------------
3249 *
3250 * @code
3251 *   const Setting &set_eta_kinks = cfg.lookup("rheology_parameters.eta_kinks");
3252 *  
3253 *   unsigned int neta_kinks = set_eta_kinks.getLength();
3254 *   system_parameters::sizeof_eta_kinks = neta_kinks;
3255 *  
3256 * @endcode
3257 *
3258 * cout << "Number of depth = " << ndepths << endl;
3259 *
3260
3261 *
3262 *
3263 * @code
3264 *   for (unsigned int i=0; i<neta_kinks; ++i)
3265 *   {
3266 *   system_parameters::eta_kinks.push_back(set_eta_kinks[i]);
3267 *   cout << "eta_kinks[" << i << "] = " << system_parameters::eta_kinks[i] << endl;
3268 *   }
3269 *  
3270 * @endcode
3271 *
3272 * get depths_rho -------------------------
3273 *
3274 * @code
3275 *   const Setting &set_depths_rho = cfg.lookup("rheology_parameters.depths_rho");
3276 *  
3277 *   unsigned int ndepths_rho = set_depths_rho.getLength();
3278 *   system_parameters::sizeof_depths_rho = ndepths_rho;
3279 *  
3280 * @endcode
3281 *
3282 * cout << "Number of depth = " << ndepths << endl;
3283 *
3284
3285 *
3286 *
3287 * @code
3288 *   for (unsigned int i=0; i<ndepths_rho; ++i)
3289 *   {
3290 *   system_parameters::depths_rho.push_back(set_depths_rho[i]);
3291 *   cout << "depths_rho[" << i << "] = " << system_parameters::depths_rho[i] << endl;
3292 *   }
3293 *  
3294 * @endcode
3295 *
3296 * get rho -------------------------
3297 *
3298 * @code
3299 *   const Setting &set_rho = cfg.lookup("rheology_parameters.rho");
3300 *  
3301 *   unsigned int nrho = set_rho.getLength();
3302 *   system_parameters::sizeof_rho = nrho;
3303 *  
3304 * @endcode
3305 *
3306 * cout << "Number of depth = " << ndepths << endl;
3307 *
3308
3309 *
3310 *
3311 * @code
3312 *   for (unsigned int i=0; i<nrho; ++i)
3313 *   {
3314 *   system_parameters::rho.push_back(set_rho[i]);
3315 *   cout << "rho[" << i << "] = " << system_parameters::rho[i] << endl;
3316 *   }
3317 *  
3318 * @endcode
3319 *
3320 * get material_id -------------------------
3321 *
3322 * @code
3323 *   const Setting &set_material_id = cfg.lookup("rheology_parameters.material_id");
3324 *  
3325 *   unsigned int nmaterial_id = set_material_id.getLength();
3326 *   system_parameters::sizeof_material_id = nmaterial_id;
3327 *  
3328 * @endcode
3329 *
3330 * cout << "Number of depth = " << ndepths << endl;
3331 *
3332
3333 *
3334 *
3335 * @code
3336 *   for (unsigned int i=0; i<nmaterial_id; ++i)
3337 *   {
3338 *   system_parameters::material_id.push_back(set_material_id[i]);
3339 *   cout << "material_id[" << i << "] = " << system_parameters::material_id[i] << endl;
3340 *   }
3341 *  
3342 * @endcode
3343 *
3344 * get G -------------------------
3345 *
3346 * @code
3347 *   const Setting &set_G = cfg.lookup("rheology_parameters.G");
3348 *  
3349 *   unsigned int nG = set_G.getLength();
3350 *   system_parameters::sizeof_G = nG;
3351 *  
3352 * @endcode
3353 *
3354 * cout << "Number of depth = " << ndepths << endl;
3355 *
3356
3357 *
3358 *
3359 * @code
3360 *   for (unsigned int i=0; i<nG; ++i)
3361 *   {
3362 *   system_parameters::G.push_back(set_G[i]);
3363 *   cout << "G[" << i << "] = " << system_parameters::G[i] << endl;
3364 *   }
3365 *  
3366 *   const Setting &rheology_parameters = root["rheology_parameters"];
3367 *   rheology_parameters.lookupValue("eta_ceiling", system_parameters::eta_ceiling);
3368 *   rheology_parameters.lookupValue("eta_floor", system_parameters::eta_floor);
3369 *   rheology_parameters.lookupValue("eta_Ea", system_parameters::eta_Ea);
3370 *   rheology_parameters.lookupValue("lat_dependence", system_parameters::lat_dependence);
3371 *   rheology_parameters.lookupValue("pressure_scale", system_parameters::pressure_scale);
3372 *   rheology_parameters.lookupValue("q", system_parameters::q);
3373 *   rheology_parameters.lookupValue("cylindrical", system_parameters::cylindrical);
3374 *   rheology_parameters.lookupValue("continue_plastic_iterations", system_parameters::continue_plastic_iterations);
3375 *   }
3376 *   catch (const SettingNotFoundException &nfex)
3377 *   {
3378 *   cerr << "We've got a problem in the rheology parameters block" << endl;
3379 *   }
3380 *  
3381 * @endcode
3382 *
3383 * Plasticity parameters
3384 *
3385 * @code
3386 *   try
3387 *   {
3388 *  
3389 *   const Setting &plasticity_parameters = root["plasticity_parameters"];
3390 *   plasticity_parameters.lookupValue("plasticity_on", system_parameters::plasticity_on);
3391 *   plasticity_parameters.lookupValue("failure_criterion", system_parameters::failure_criterion);
3392 *   plasticity_parameters.lookupValue("max_plastic_iterations", system_parameters::max_plastic_iterations);
3393 *   plasticity_parameters.lookupValue("smoothing_radius", system_parameters::smoothing_radius);
3394 *   }
3395 *   catch (const SettingNotFoundException &nfex)
3396 *   {
3397 *   cerr << "We've got a problem in the plasticity parameters block" << endl;
3398 *   }
3399 *  
3400 * @endcode
3401 *
3402 * Viscoelasticity parameters
3403 *
3404
3405 *
3406 *
3407 * @code
3408 *   try
3409 *   {
3410 *  
3411 *   const Setting &viscoelasticity_parameters = root["viscoelasticity_parameters"];
3412 *   viscoelasticity_parameters.lookupValue("initial_elastic_iterations", system_parameters::initial_elastic_iterations);
3413 *   viscoelasticity_parameters.lookupValue("elastic_time", system_parameters::elastic_time);
3414 *   viscoelasticity_parameters.lookupValue("viscous_time", system_parameters::viscous_time);
3415 *   viscoelasticity_parameters.lookupValue("initial_disp_target", system_parameters::initial_disp_target);
3416 *   viscoelasticity_parameters.lookupValue("final_disp_target", system_parameters::final_disp_target);
3417 *   viscoelasticity_parameters.lookupValue("current_time_interval", system_parameters::current_time_interval);
3418 *  
3419 *   system_parameters::viscous_time *= SECSINYEAR;
3420 *   }
3421 *   catch (const SettingNotFoundException &nfex)
3422 *   {
3423 *   cerr << "We've got a problem in the viscoelasticity parameters block" << endl;
3424 *   }
3425 *  
3426 * @endcode
3427 *
3428 * Mesh refinement parameters
3429 *
3430 * @code
3431 *   try
3432 *   {
3433 *  
3434 *   const Setting &mesh_refinement_parameters = root["mesh_refinement_parameters"];
3435 *   mesh_refinement_parameters.lookupValue("global_refinement", system_parameters::global_refinement);
3436 *   mesh_refinement_parameters.lookupValue("small_r_refinement", system_parameters::small_r_refinement);
3437 *   mesh_refinement_parameters.lookupValue("crustal_refinement", system_parameters::crustal_refinement);
3438 *   mesh_refinement_parameters.lookupValue("crust_refine_region", system_parameters::crust_refine_region);
3439 *   mesh_refinement_parameters.lookupValue("surface_refinement", system_parameters::surface_refinement);
3440 *   }
3441 *   catch (const SettingNotFoundException &nfex)
3442 *   {
3443 *   cerr << "We've got a problem in the mesh refinement parameters block" << endl;
3444 *   }
3445 *  
3446 * @endcode
3447 *
3448 * Solver parameters
3449 *
3450 * @code
3451 *   try
3452 *   {
3453 *   const Setting &solve_parameters = root["solve_parameters"];
3454 *   solve_parameters.lookupValue("iteration_coefficient", system_parameters::iteration_coefficient);
3455 *   solve_parameters.lookupValue("tolerance_coefficient", system_parameters::tolerance_coefficient);
3456 *  
3457 *  
3458 *   }
3459 *   catch (const SettingNotFoundException &nfex)
3460 *   {
3461 *   cerr << "We've got a problem in the solver parameters block" << endl;
3462 *   }
3463 *  
3464 * @endcode
3465 *
3466 * Time step parameters
3467 *
3468 * @code
3469 *   try
3470 *   {
3471 *   const Setting &time_step_parameters = root["time_step_parameters"];
3472 *   time_step_parameters.lookupValue("present_time", system_parameters::present_time);
3473 *   time_step_parameters.lookupValue("present_timestep", system_parameters::present_timestep);
3474 *   time_step_parameters.lookupValue("total_viscous_steps", system_parameters::total_viscous_steps);
3475 *   }
3476 *   catch (const SettingNotFoundException &nfex)
3477 *   {
3478 *   cerr << "We've got a problem in the time step parameters block" << endl;
3479 *   }
3480 *  
3481 *   write_config();
3482 *   }
3483 *   }
3484 *  
3485 *  
3486 *  
3487 *  
3488 *  
3489 * @endcode
3490
3491
3492<a name="ann-support_code/ellipsoid_fit.h"></a>
3493<h1>Annotated version of support_code/ellipsoid_fit.h</h1>
3494 *
3495 *
3496 *
3497 *
3498 * @code
3499 *   /*
3500 *   * ellipsoid_fit.h
3501 *   *
3502 *   * Created on: Jul 24, 2015
3503 *   * Author: antonermakov
3504 *   */
3505 *  
3506 *   #include <deal.II/base/quadrature_lib.h>
3507 *   #include <deal.II/base/function.h>
3508 *   #include <deal.II/base/logstream.h>
3509 *   #include <deal.II/lac/vector.h>
3510 *   #include <deal.II/lac/full_matrix.h>
3511 *   #include <deal.II/lac/sparse_matrix.h>
3512 *   #include <deal.II/lac/solver_cg.h>
3513 *   #include <deal.II/lac/precondition.h>
3514 *   #include <deal.II/grid/tria.h>
3515 *   #include <deal.II/dofs/dof_handler.h>
3516 *   #include <deal.II/grid/tria_accessor.h>
3517 *   #include <deal.II/grid/tria_iterator.h>
3518 *   #include <deal.II/dofs/dof_accessor.h>
3519 *   #include <deal.II/dofs/dof_tools.h>
3520 *   #include <deal.II/fe/fe_q.h>
3521 *   #include <deal.II/fe/fe_values.h>
3522 *   #include <deal.II/numerics/vector_tools.h>
3523 *   #include <deal.II/numerics/matrix_tools.h>
3524 *   #include <deal.II/numerics/data_out.h>
3525 *   #include <deal.II/grid/grid_in.h>
3526 *   #include <deal.II/grid/grid_out.h>
3527 *   #include <deal.II/base/point.h>
3528 *   #include <deal.II/grid/grid_generator.h>
3529 *  
3530 *   #include <fstream>
3531 *   #include <sstream>
3532 *   #include <iostream>
3533 *   #include <iomanip>
3534 *   #include <cstdlib>
3535 *  
3536 *  
3537 *   #include "local_math.h"
3538 *  
3539 *   using namespace dealii;
3540 *  
3541 *   template <int dim>
3542 *   class ellipsoid_fit
3543 *   {
3544 *   public:
3545 *   inline ellipsoid_fit (Triangulation<dim,dim> *pi)
3546 *   {
3547 *   p_triangulation = pi;
3548 *   };
3549 *   void compute_fit(std::vector<double> &ell, unsigned char bndry);
3550 *  
3551 *  
3552 *   private:
3553 *   Triangulation<dim,dim> *p_triangulation;
3554 *  
3555 *   };
3556 *  
3557 *  
3558 * @endcode
3559 *
3560 * This function computes ellipsoid fit to a set of vertices that lie on the
3561 * boundary_that_we_need
3562 *
3563 * @code
3564 *   template <int dim>
3565 *   void ellipsoid_fit<dim>::compute_fit(std::vector<double> &ell, unsigned char boundary_that_we_need)
3566 *   {
3567 *   typename Triangulation<dim>::active_cell_iterator cell = p_triangulation->begin_active();
3568 *   typename Triangulation<dim>::active_cell_iterator endc = p_triangulation->end();
3569 *  
3570 *   FullMatrix<double> A(p_triangulation->n_vertices(),dim);
3571 *   Vector<double> x(dim);
3572 *   Vector<double> b(p_triangulation->n_vertices());
3573 *  
3574 *   std::vector<bool> vertex_touched (p_triangulation->n_vertices(),
3575 *   false);
3576 *  
3577 *   unsigned int j = 0;
3578 *   unsigned char boundary_ids;
3579 *   std::vector<unsigned int> ind_bnry_row;
3580 *   std::vector<unsigned int> ind_bnry_col;
3581 *  
3582 * @endcode
3583 *
3584 * assemble the sensitivity matrix and r.h.s.
3585 *
3586 * @code
3587 *   for (; cell != endc; ++cell)
3588 *   {
3589 *   if (boundary_that_we_need != 0)
3590 *   cell->set_manifold_id(cell->material_id());
3591 *   for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
3592 *   {
3593 *   if (boundary_that_we_need == 0) //if this is the outer surface, then look for boundary ID 0; otherwise look for material ID change.
3594 *   {
3595 *   boundary_ids = cell->face(f)->boundary_id();
3596 *   if (boundary_ids == boundary_that_we_need)
3597 *   {
3598 *   for (unsigned int v = 0;
3599 *   v < GeometryInfo<dim>::vertices_per_face; ++v)
3600 *   if (vertex_touched[cell->face(f)->vertex_index(v)]
3601 *   == false)
3602 *   {
3603 *   vertex_touched[cell->face(f)->vertex_index(v)] =
3604 *   true;
3605 *   for (unsigned int i = 0; i < dim; ++i)
3606 *   {
3607 * @endcode
3608 *
3609 * stiffness matrix entry
3610 *
3611 * @code
3612 *   A(j, i) = pow(cell->face(f)->vertex(v)[i], 2);
3613 * @endcode
3614 *
3615 * r.h.s. entry
3616 *
3617 * @code
3618 *   b[j] = 1.0;
3619 * @endcode
3620 *
3621 * if mesh if not full: set the indicator
3622 *
3623 * @code
3624 *   }
3625 *   ind_bnry_row.push_back(j);
3626 *   j++;
3627 *   }
3628 *   }
3629 *   }
3630 *   else //find the faces that are at the boundary between materials, get the vertices, and write them into the stiffness matrix
3631 *   {
3632 *   if (cell->neighbor(f) != endc)
3633 *   {
3634 *   if (cell->material_id() != cell->neighbor(f)->material_id()) //finds face is at internal boundary
3635 *   {
3636 *   int high_mat_id = std::max(cell->material_id(),
3637 *   cell->neighbor(f)->material_id());
3638 *   if (high_mat_id == boundary_that_we_need) //finds faces at the correct internal boundary
3639 *   {
3640 *   for (unsigned int v = 0;
3641 *   v < GeometryInfo<dim>::vertices_per_face;
3642 *   ++v)
3643 *   if (vertex_touched[cell->face(f)->vertex_index(
3644 *   v)] == false)
3645 *   {
3646 *   vertex_touched[cell->face(f)->vertex_index(
3647 *   v)] = true;
3648 *   for (unsigned int i = 0; i < dim; ++i)
3649 *   {
3650 * @endcode
3651 *
3652 * stiffness matrix entry
3653 *
3654 * @code
3655 *   A(j, i) = pow(
3656 *   cell->face(f)->vertex(v)[i], 2);
3657 * @endcode
3658 *
3659 * r.h.s. entry
3660 *
3661 * @code
3662 *   b[j] = 1.0;
3663 * @endcode
3664 *
3665 * if mesh if not full: set the indicator
3666 *
3667 * @code
3668 *   }
3669 *   ind_bnry_row.push_back(j);
3670 *   j++;
3671 *   }
3672 *   }
3673 *   }
3674 *   }
3675 *   }
3676 *   }
3677 *   }
3678 *   if (ind_bnry_row.size()>0)
3679 *   {
3680 *  
3681 * @endcode
3682 *
3683 * maxtrix A'*A and vector A'*b; A'*A*x = A'*b -- normal system of equations
3684 *
3685 * @code
3686 *   FullMatrix<double> AtA(dim,dim);
3687 *   Vector<double> Atb(dim);
3688 *  
3689 *   FullMatrix<double> A_out(ind_bnry_row.size(),dim);
3690 *   Vector<double> b_out(ind_bnry_row.size());
3691 *  
3692 *   for (unsigned int i=0; i<dim; ++i)
3693 *   ind_bnry_col.push_back(i);
3694 *  
3695 *   for (unsigned int i=0; i<ind_bnry_row.size(); ++i)
3696 *   b_out(i) = 1;
3697 *  
3698 *   A_out.extract_submatrix_from(A, ind_bnry_row, ind_bnry_col);
3699 *   A_out.Tmmult(AtA,A_out,true);
3700 *   A_out.Tvmult(Atb,b_out,true);
3701 *  
3702 * @endcode
3703 *
3704 * solve normal system of equations
3705 *
3706 * @code
3707 *   SolverControl solver_control (1000, 1e-12);
3708 *   SolverCG<> solver (solver_control);
3709 *   solver.solve (AtA, x, Atb, PreconditionIdentity());
3710 *  
3711 * @endcode
3712 *
3713 * find ellipsoidal axes
3714 *
3715 * @code
3716 *   for (unsigned int i=0; i<dim; ++i)
3717 *   ell.push_back(sqrt(1.0/x[i]));
3718 *   }
3719 *   else
3720 *   std::cerr << "fit_ellipsoid: no points to fit" << std::endl;
3721 *  
3722 *   }
3723 *  
3724 * @endcode
3725
3726
3727<a name="ann-support_code/ellipsoid_grav.h"></a>
3728<h1>Annotated version of support_code/ellipsoid_grav.h</h1>
3729 *
3730 *
3731 *
3732 *
3733 * @code
3734 *   /* -----------------------------------------------------------------------------
3735 *   *
3736 *   * SPDX-License-Identifier: LGPL-2.1-or-later
3737 *   * Copyright (C) 2012 by Roger R. Fu
3738 *   *
3739 *   * This file is part of the deal.II code gallery.
3740 *   *
3741 *   * -----------------------------------------------------------------------------
3742 *   *
3743 *   * Author: Roger R. Fu
3744 *   * Reference Pohanka 2011, Contrib. Geophys. Geodes.
3745 *   */
3746 *  
3747 *   #include <math.h>
3748 *   #include <deal.II/base/point.h>
3749 *   #include <fstream>
3750 *   #include <iostream>
3751 *  
3752 *   namespace A_Grav_namespace
3753 *   {
3754 *   namespace system_parameters
3755 *   {
3756 *   double mantle_rho;
3757 *   double core_rho;
3758 *   double excess_rho;
3759 *   double r_eq;
3760 *   double r_polar;
3761 *  
3762 *   double r_core_eq;
3763 *   double r_core_polar;
3764 *   }
3765 *  
3766 *   template <int dim>
3767 *   class AnalyticGravity
3768 *   {
3769 *   public:
3770 *   void setup_vars (std::vector<double> v);
3771 *   void get_gravity (const ::Point<dim> &p, std::vector<double> &g);
3772 *  
3773 *   private:
3774 *   double ecc;
3775 *   double eV;
3776 *   double ke;
3777 *   double r00;
3778 *   double r01;
3779 *   double r11;
3780 *   double ecc_c;
3781 *   double eV_c;
3782 *   double ke_c;
3783 *   double r00_c;
3784 *   double r01_c;
3785 *   double r11_c;
3786 *   double g_coeff;
3787 *   double g_coeff_c;
3788 *   };
3789 *  
3790 *   template <int dim>
3791 *   void AnalyticGravity<dim>::get_gravity (const ::Point<dim> &p, std::vector<double> &g)
3792 *   {
3793 *   double rsph = std::sqrt(p[0] * p[0] + p[1] * p[1]);
3794 *   double thetasph = std::atan2(p[0], p[1]);
3795 *   double costhetasph = std::cos(thetasph);
3796 *  
3797 * @endcode
3798 *
3799 * convert to elliptical coordinates for silicates
3800 *
3801 * @code
3802 *   double stemp = std::sqrt((rsph * rsph - eV * eV + std::sqrt((rsph * rsph - eV * eV) * (rsph * rsph - eV * eV)
3803 *   + 4 * eV * eV * rsph * rsph * costhetasph *costhetasph)) / 2);
3804 *   double vout = stemp / system_parameters::r_eq / std::sqrt(1 - ecc * ecc);
3805 *   double eout = std::acos(rsph * costhetasph / stemp);
3806 *  
3807 * @endcode
3808 *
3809 * convert to elliptical coordinates for core correction
3810 *
3811 * @code
3812 *   double stemp_c = std::sqrt((rsph * rsph - eV_c * eV_c + std::sqrt((rsph * rsph - eV_c * eV_c) * (rsph * rsph - eV_c * eV_c)
3813 *   + 4 * eV_c * eV_c * rsph * rsph * costhetasph *costhetasph)) / 2);
3814 *   double vout_c = stemp_c / system_parameters::r_core_eq / std::sqrt(1 - ecc_c * ecc_c);
3815 *   double eout_c = std::acos(rsph * costhetasph / stemp_c);
3816 *  
3817 * @endcode
3818 *
3819 * shell contribution
3820 *
3821 * @code
3822 *   g[0] = g_coeff * r11 * std::sqrt((1 - ecc * ecc) * vout * vout + ecc * ecc) * std::sin(eout);
3823 *   g[1] = g_coeff * r01 * vout * std::cos(eout) / std::sqrt(1 - ecc * ecc);
3824 *  
3825 * @endcode
3826 *
3827 * core contribution
3828 *
3829 * @code
3830 *   double expected_y = system_parameters::r_core_polar * std::sqrt(1 -
3831 *   (p[0] * p[0] / system_parameters::r_core_eq / system_parameters::r_core_eq));
3832 *  
3833 *  
3834 *   if (p[1] <= expected_y)
3835 *   {
3836 *   g[0] += g_coeff_c * r11_c * std::sqrt((1 - ecc_c * ecc_c) * vout_c * vout_c + ecc_c * ecc_c) * std::sin(eout_c);
3837 *   g[1] += g_coeff_c * r01_c * vout_c * std::cos(eout_c) / std::sqrt(1 - ecc_c * ecc_c);
3838 *   }
3839 *   else
3840 *   {
3841 *   double g_coeff_co = - 2.795007963255562e-10 * system_parameters::excess_rho * system_parameters::r_core_eq
3842 *   / vout_c / vout_c;
3843 *   double r00_co = 0;
3844 *   double r01_co = 0;
3845 *   double r11_co = 0;
3846 *  
3847 *   if (system_parameters::r_core_polar == system_parameters::r_core_eq)
3848 *   {
3849 *   r00_co = 1;
3850 *   r01_co = 1;
3851 *   r11_co = 1;
3852 *   }
3853 *   else
3854 *   {
3855 *   r00_co = ke_c * vout_c * std::atan2(1, ke_c * vout_c);
3856 *   double ke_co2 = ke_c * ke_c * vout_c * vout_c;
3857 *   r01_co = 3 * ke_co2 * (1 - r00_co);
3858 *   r11_co = 3 * ((ke_co2 + 1) * r00_co - ke_co2) / 2;
3859 *   }
3860 *   g[0] += g_coeff_co * vout_c * r11_co / std::sqrt((1 - ecc_c* ecc_c) * vout_c * vout_c + ecc_c * ecc_c) * std::sin(eout_c);
3861 *   g[1] += g_coeff_co * r01_co * std::cos(eout_c) / std::sqrt(1 - ecc_c * ecc_c);
3862 *   }
3863 *   }
3864 *  
3865 *   template <int dim>
3866 *   void AnalyticGravity<dim>::setup_vars (std::vector<double> v)
3867 *   {
3868 *   system_parameters::r_eq = v[0];
3869 *   system_parameters::r_polar = v[1];
3870 *   system_parameters::r_core_eq = v[2];
3871 *   system_parameters::r_core_polar = v[3];
3872 *   system_parameters::mantle_rho = v[4];
3873 *   system_parameters::core_rho = v[5];
3874 *   system_parameters::excess_rho = system_parameters::core_rho - system_parameters::mantle_rho;
3875 *  
3876 *  
3877 * @endcode
3878 *
3879 * Shell
3880 *
3881 * @code
3882 *   if (system_parameters::r_polar > system_parameters::r_eq)
3883 *   {
3884 * @endcode
3885 *
3886 * This makes the gravity field nearly that of a sphere in case the body becomes prolate
3887 *
3888 * @code
3889 *   std::cout << "\nWarning: The model body has become prolate. \n";
3890 *   ecc = 0.001;
3891 *   }
3892 *   else
3893 *   {
3894 *   ecc = std::sqrt(1 - (system_parameters::r_polar * system_parameters::r_polar / system_parameters::r_eq / system_parameters::r_eq));
3895 *   }
3896 *  
3897 *   eV = ecc * system_parameters::r_eq;
3898 *   ke = std::sqrt(1 - (ecc * ecc)) / ecc;
3899 *   r00 = ke * std::atan2(1, ke);
3900 *   double ke2 = ke * ke;
3901 *   r01 = 3 * ke2 * (1 - r00);
3902 *   r11 = 3 * ((ke2 + 1) * r00 - ke2) / 2;
3903 *   g_coeff = - 2.795007963255562e-10 * system_parameters::mantle_rho * system_parameters::r_eq;
3904 *  
3905 * @endcode
3906 *
3907 * Core
3908 *
3909 * @code
3910 *   if (system_parameters::r_core_polar > system_parameters::r_core_eq)
3911 *   {
3912 *   std::cout << "\nWarning: The model core has become prolate. \n";
3913 *   ecc_c = 0.001;
3914 *   }
3915 *   else
3916 *   {
3917 *   ecc_c = std::sqrt(1 - (system_parameters::r_core_polar * system_parameters::r_core_polar / system_parameters::r_core_eq / system_parameters::r_core_eq));
3918 *   }
3919 *   eV_c = ecc_c * system_parameters::r_core_eq;
3920 *   if (system_parameters::r_core_polar == system_parameters::r_core_eq)
3921 *   {
3922 *   ke_c = 1;
3923 *   r00_c = 1;
3924 *   r01_c = 1;
3925 *   r11_c = 1;
3926 *   g_coeff_c = - 2.795007963255562e-10 * system_parameters::excess_rho * system_parameters::r_core_eq;
3927 *   }
3928 *   else
3929 *   {
3930 *   ke_c = std::sqrt(1 - (ecc_c * ecc_c)) / ecc_c;
3931 *   r00_c = ke_c * std::atan2(1, ke_c);
3932 *   double ke2_c = ke_c * ke_c;
3933 *   r01_c = 3 * ke2_c * (1 - r00_c);
3934 *   r11_c = 3 * ((ke2_c + 1) * r00_c - ke2_c) / 2;
3935 *   g_coeff_c = - 2.795007963255562e-10 * system_parameters::excess_rho * system_parameters::r_core_eq;
3936 *   }
3937 * @endcode
3938 *
3939 * std::cout << "Loaded variables: ecc = " << ecc_c << " ke = " << ke_c << " r00 = " << r00_c << " r01 = " << r01_c << " r11 = " << r11_c << "\n";
3940 *
3941 * @code
3942 *   }
3943 *   }
3944 * @endcode
3945
3946
3947<a name="ann-support_code/local_math.h"></a>
3948<h1>Annotated version of support_code/local_math.h</h1>
3949 *
3950 *
3951 *
3952 *
3953 * @code
3954 *   /*
3955 *   * File: localmath.h
3956 *   * Author: antonermakov
3957 *   *
3958 *   * Created on September 21, 2013, 7:14 PM
3959 *   */
3960 *  
3961 *   #ifndef LOCAL_MATH_
3962 *   #define LOCAL_MATH_
3963 *  
3964 *   #define PI 3.14159265358979323846
3965 *   #define TWOPI 6.283185307179586476925287
3966 *   #define SECSINYEAR 3.155692608e+07
3967 * @endcode
3968 *
3969 * #define ABS(a) ((a) < 0 ? -(a) : (a))
3970 *
3971
3972 *
3973 * double factorial(int n)
3974 * {
3975 * if(n == 0) {
3976 * return(1.);
3977 * } else if(n == 1) {
3978 * return(1.);
3979 * } else if(n == 2) {
3980 * return(2.);
3981 * } else if(n == 3) {
3982 * return(6.);
3983 * } else if(n == 4) {
3984 * return(24.);
3985 * } else {
3986 * exit(-1);
3987 * }
3988 * }
3989 *
3990
3991 *
3992 * double fudge(int m)
3993 * {
3994 * if(m == 0) {
3995 * return(1.0);
3996 * } else {
3997 * return(2.0);
3998 * }
3999 * }
4000 *
4001
4002 *
4003 *
4004
4005 *
4006 * double sign(double x)
4007 * {
4008 * if(x > 0) {
4009 * return(1.0);
4010 * } else if(x < 0.0) {
4011 * return(-1.0);
4012 * } else {
4013 * return(0.0);
4014 * }
4015 * }
4016 *
4017
4018 *
4019 * double pv0(double x)
4020 * {
4021 * double ans;
4022 *
4023
4024 *
4025 * ans = x - TWOPI*floor(x/TWOPI);
4026 * if(ans > TWOPI/2.0) {
4027 * ans = ans - TWOPI;
4028 * }
4029 *
4030
4031 *
4032 * return(ans);
4033 * }
4034 *
4035
4036 *
4037 *
4038 * @code
4039 *   /* assumes l=2 */
4040 * @endcode
4041 *
4042 * double System::Plm(int m, double x)
4043 * {
4044 * if(m == 0) {
4045 * return(1.5*x*x - 0.5);
4046 * } else if(m == 1) {
4047 * return(3.0*x*sqrt(1.0 - x*x));
4048 * } else if(m == 2) {
4049 * return(3.0 - 3.0*x*x);
4050 * } else {
4051 * exit(-1);
4052 * }
4053 * }
4054 *
4055
4056 *
4057 * double System::DP(int m, double x)
4058 * {
4059 * if(m == 0) {
4060 * return(3.0*x);
4061 * } else if(m == 1) {
4062 * return((3.0 - 6.0*x*x)/sqrt(1.0 - x*x));
4063 * } else if(m == 2) {
4064 * return(- 6.0*x);
4065 * } else {
4066 * exit(-1);
4067 * }
4068 * }
4069 *
4070
4071 *
4072 *
4073
4074 *
4075 *
4076 * @code
4077 *   #endif /* LOCALMATH_H */
4078 *  
4079 * @endcode
4080
4081
4082*/
DataOut_DoFData::attach_dof_handler
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
FE_Q
Definition fe_q.h:554
Function::vector_value_list
virtual void vector_value_list(const std::vector< Point< dim > > &points, std::vector< Vector< RangeNumberType > > &values) const
Function::value
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
Function::vector_value
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &values) const
Point
Definition point.h:111
SparseILU::AdditionalData
typename SparseLUDecomposition< number >::AdditionalData AdditionalData
Definition sparse_ilu.h:79
SparseILU::initialize
void initialize(const SparseMatrix< somenumber > &matrix, const AdditionalData &parameters=AdditionalData())
DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition config.h:512
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition config.h:555
Assert
#define Assert(cond, exc)
Definition exceptions.h:1659
AssertThrow
#define AssertThrow(cond, exc)
Definition exceptions.h:1780
DoFHandler::active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition dof_handler.h:440
schur_complement
LinearOperator< Range_2, Domain_2, Payload > schur_complement(const LinearOperator< Domain_1, Range_1, Payload > &A_inv, const LinearOperator< Range_1, Domain_2, Payload > &B, const LinearOperator< Range_2, Domain_1, Payload > &C, const LinearOperator< Range_2, Domain_2, Payload > &D)
Definition schur_complement.h:247
MeshWorker::loop
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:564
DoFTools::make_hanging_node_constraints
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
Definition dof_tools_constraints.cc:3073
DoFTools::make_sparsity_pattern
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)
Definition dof_tools_sparsity.cc:56
VectorTools::compute_no_normal_flux_constraints
void compute_no_normal_flux_constraints(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int first_vector_component, const std::set< types::boundary_id > &boundary_ids, AffineConstraints< number > &constraints, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const bool use_manifold_for_normal=true)
update_values
@ update_values
Shape function values.
Definition fe_update_flags.h:75
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition fe_update_flags.h:134
update_gradients
@ update_gradients
Shape function gradients.
Definition fe_update_flags.h:81
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition fe_update_flags.h:127
data
std::vector< index_type > data
Definition mpi.cc:735
size
std::size_t size
Definition mpi.cc:734
AdaptationStrategies::Refinement::l2_norm
std::vector< value_type > l2_norm(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
Differentiation::SD::floor
Expression floor(const Expression &x)
Definition symengine_math.cc:290
DoFRenumbering::component_wise
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
Definition dof_renumbering.cc:666
DoFRenumbering::Cuthill_McKee
void Cuthill_McKee(DoFHandler< dim, spacedim > &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false, const std::vector< types::global_dof_index > &starting_indices=std::vector< types::global_dof_index >())
Definition dof_renumbering.cc:366
DoFTools::count_dofs_per_fe_block
std::vector< types::global_dof_index > count_dofs_per_fe_block(const DoFHandler< dim, spacedim > &dof, const std::vector< unsigned int > &target_block=std::vector< unsigned int >())
Definition dof_tools.cc:2284
IteratorState::valid
@ valid
Iterator points to a valid object.
Definition tria_iterator_base.h:37
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition lapack_support.h:57
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::b
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Physics::VectorRelations::angle
Number angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b)
Utilities::int_to_string
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470
VectorTools::interpolate_boundary_values
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={})
WorkStream::internal::tbb_no_coloring::run
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)
Definition work_stream.h:511
internal::MatrixFreeFunctions::check
bool check(const ConstraintKinds kind_in, const unsigned int dim)
Definition hanging_nodes_internal.h:92
internal::p4est::functions
int(&) functions(const void *v1, const void *v2)
Definition p4est_wrappers.cc:338
std::cos
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6609
std::sin
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6589
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6869
std::pow
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
Definition vectorization.h:6885
std::atan
inline ::VectorizedArray< Number, width > atan(const ::VectorizedArray< Number, width > &x)
Definition vectorization.h:6689
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6927
types::material_id
unsigned int material_id
Definition types.h:167
triangulation
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition p4est_wrappers.cc:68
Generated by doxygen 1.9.8