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