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