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