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