Reference documentation for deal.II version GIT relicensing-245-g36f19064f7 2024-03-29 07:20:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
The 'Viscoelastoplastic topography evolution' code gallery program

This program was contributed by Roger Fu <rogerfu@fas.harvard.edu>.
It comes without any warranty or support by its authors or the authors of deal.II.

This program is part of the deal.II code gallery and consists of the following files (click to inspect):

Annotated version of Readme.md

Readme file for CeresFE

Motivation for project

This code was made to simulate the evolution of global-scale topography on planetary bodies. Specifically, it is designed to compute the rates of topography relaxation on the dwarf planet Ceres. The NASA Dawn mission, in orbit around Ceres since March, 2015, has produced a high resolution shape model of its surface. As on other planets including the Earth, topography on Ceres is subject to decay over time due to processes such as viscous flow and brittle failure. Because the efficiency of these processes is dependent on the material properties of the body at depth, simulating the decay of topography and comparing it to the observed shape model permits insights into Ceres' internal stucture.

Some previous applications of this basic idea- using topography to constrain internal structure- may be found in the following references:

  1. Takeuchi, H. and Hasegawa, Y. (1965) Viscosity distribution within the Earth. Geophys. J. R. astr. Soc. 9, 503-508.
  2. Anderson, D. L. and O'Connell, R. (1967) Viscosity of the Earth. Geophys. J. R. astr. Soc. 14, 287-295.
  3. Solomon, S. C., Comer, R. P., Head, J. W. (1982) The Evolution of impact basins: Viscous relaxation of topographic relief.
  4. Zuber, M. T. et al. (2000) Internal structure and early thermal evolution of Mars from Mars Global Surveyor topography and gravity. Science 287, 1788-1793.
  5. Fu, R. R. et al. (2014) Efficient early global relaxation of asteroid Vesta. Icarus 240, 133-145.

The code included here is a development of a simpler code for the asteroid Vesta, published as reference 5 above. Because both versions of the code were written specifically to model long wavelength topography on these small bodies, the code is rather specific. We hope certain components of it may be useful to the reader even if the problem of topographic relaxation on asteroid belt bodies is not on everyone's radar.

Quick facts about the code

Viscoelastoplastic Asymmetric Lagrangian Uses analytical self-gravity One sentence purpose: Simulates evolution of topography due to self-gravity on axisymmetry planetary body.

More detailed properties of the code in CeresFE

Viscoelastoplasticity

The code is viscoelastoplastic: it solves the Stokes equations modified to include elasticity and iteratively uses the stress solution to account for displacement due to brittle failure The implementation of viscoelasticity follows mainly section 2.2.3 of Keller, T. et al. (2013) Numerical modelling of magma dynamics coupled to tectonic deformation of lithosphere and crust. Geophys. J. Int. 195, 1406-1442. At the end of each FE calculation, the principle stresses (in 3D) are computed in all cells. Each cell is evaluated according to either Byerlee's Rule or a damaged rock brittle failure criterion to determine if failure occurred. See Byerlee, J. (1978) Friction of rocks, Pageoph. 116, 615-626. and Schultz, R. A. (1993) Brittle strength of basalitc rock masses with applications to Venus. J. Geophys. Res. 98, 10,883-10,895. If a cell failed, its viscosity is lowered by a computed amount to simulate plastic yielding. The viscosity fields is smoothed and the FE model run again. This is repeated until the number of failed cells falls below a prescribed number. The final viscosity field (i.e., the effective viscosity) is then used to compute velocities and advance the mesh.

Domain and boundary conditions

The domain of the model is 2D, but the Stokes equations are cast in axisymmetric form. The domain consists of approximately a quarter ellipse, with two straight edges corresponding to the rotation axis and equator of the body. No normal flux boundary conditions are applied to these edges. The remaining curved edge that corresponds to the surface of the body is assigned a zero pressure boundary condition With respect to self-gravity, an ellipse is fitted to the outer surface and any internal density surfaces at each time step and a gravity field is computed analytically following Pohanka, V. (2011) Gravitational field of the homogeneous rotational ellipsoidal body: a simple derivation and applications. Contrib. Geophys. Geodesy 41, 117-157.

Description of files in repo

src/ceres.cc Main code support_code/config_in.h Reads config file and intializes system parameters support_code/ellipsoid_fit.h Finds best-fit ellipse for surface and internal density boundaries. Also uses deal.II support_code/ellipsoid_grav.h Analytically computes self gravity of layered ellipsoids structure support_code/local_math.h Defines some constants for convenience meshes/sample_CeresFE_mesh.inp Sample input mesh config/sample_CeresFE_config.cfg ample configurations file with simulation parameters

Other dependencies

Two more code packages are necessary to run CeresFE:

  1. config++: https://sourceforge.net/projects/config/
  2. Armadillo: http://arma.sourceforge.net

To run the code

After running cmake and compiling, run the executable with one argument, which is the config file :

\(ceres \){SOURCE_DIR}/config/sample_CeresFE_config.cfg

Annotated version of src/ceres.cc

  /* -----------------------------------------------------------------------------
  *
  * SPDX-License-Identifier: LGPL-2.1-or-later
  * Copyright (C) 2012 by Roger R. Fu
  *
  * This file is part of the deal.II code gallery.
  *
  * -----------------------------------------------------------------------------
  *
  * Viscoelastoplastic relaxation of Ceres
  * Author: Roger R. Fu
  * Adapted from Fu et al. 2014 Icarus 240, 133-145 starting Oct. 19, 2014
  */
 
  /*
  Summary of output files:
 
  One per run:
 
  - initial_mesh.eps : Visualization of initially imported mesh
  - 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
 
  One per time step:
 
  - 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().
  - timeXX_baseviscosities.txt : Columns (1) cell x, (2) cell y, (3) base viscosity in Pa s. Written in solution_stresses().
  - 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()
 
  One per plasticity step:
 
  - timeXX_flowYY.txt : Same as timeXX_elastic_displacements.txt above
  - timeXX_principalstressesYY.txt : Columns with sigma1 and sigma3 at each cell. Same order as timeXX_baseviscosities.txt. Written in solution_stresses().
  - timeXX_stresstensorYY.txt : Columns with components 11, 22, 33, and 13 of stress tensor at each cell. Written in solution_stresses().
  - timeXX_failurelocationsYY.txt : Gives x,y coordinates of all cells where failure occurred. Written in solution_stresses().
  - timeXX_viscositiesregYY.txt : Gives smoothed and regularized (i.e., floor and ceiling-filtered) effective viscosities. Written at end of solution_stresses().
 
  */
 
  #include <deal.II/base/quadrature_lib.h>
  #include <deal.II/base/logstream.h>
  #include <deal.II/base/function.h>
  #include <deal.II/base/utilities.h>
 
  #include <deal.II/lac/block_vector.h>
  #include <deal.II/lac/full_matrix.h>
  #include <deal.II/lac/block_sparse_matrix.h>
  #include <deal.II/lac/solver_cg.h>
  #include <deal.II/lac/precondition.h>
  #include <deal.II/lac/affine_constraints.h>
 
  #include <deal.II/grid/tria.h>
  #include <deal.II/grid/grid_generator.h>
  #include <deal.II/grid/tria_accessor.h>
  #include <deal.II/grid/tria_iterator.h>
  #include <deal.II/grid/grid_tools.h>
  #include <deal.II/grid/manifold_lib.h>
  #include <deal.II/grid/grid_refinement.h>
  #include <deal.II/grid/grid_in.h>
 
  #include <deal.II/dofs/dof_handler.h>
  #include <deal.II/dofs/dof_renumbering.h>
  #include <deal.II/dofs/dof_accessor.h>
  #include <deal.II/dofs/dof_tools.h>
 
  #include <deal.II/fe/fe_q.h>
  #include <deal.II/fe/fe_system.h>
  #include <deal.II/fe/fe_values.h>
 
  #include <deal.II/numerics/vector_tools.h>
  #include <deal.II/numerics/matrix_tools.h>
  #include <deal.II/numerics/data_out.h>
  #include <deal.II/numerics/error_estimator.h>
  #include <deal.II/numerics/derivative_approximation.h>
  #include <deal.II/numerics/fe_field_function.h>
 
  #include <deal.II/lac/sparse_direct.h>
  #include <deal.II/lac/sparse_ilu.h>
 
  #include <iostream>
  #include <fstream>
  #include <sstream>
  #include <string>
  #include <time.h>
  #include <math.h>
 
  #include <armadillo>
 
  #include "../support_code/ellipsoid_grav.h"
  #include "../support_code/ellipsoid_fit.h"
  #include "../support_code/config_in.h"
 
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition config.h:516
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition config.h:559

As in all programs, the namespace dealii is included:

  namespace Step22
  {
 
  using namespace dealii;
  using namespace arma;
 
  template<int dim>
  struct InnerPreconditioner;
 
  template<>
  struct InnerPreconditioner<2>
  {
  typedef SparseDirectUMFPACK type;
  };
 
  template<>
  struct InnerPreconditioner<3>
  {
  typedef SparseILU<double> type;
  };
 

Auxiliary functions

  template<int dim>
  class AuxFunctions
  {
  public:
  Tensor<2, 2> get_rotation_matrix(const std::vector<Tensor<1, 2> > &grad_u);
  };
 
  template<int dim>
  Tensor<2, 2> AuxFunctions<dim>::get_rotation_matrix(
  const std::vector<Tensor<1, 2> > &grad_u)
  {
  const double curl = (grad_u[1][0] - grad_u[0][1]);
 
  const double angle = std::atan(curl);
 
  const double t[2][2] = { { cos(angle), sin(angle) }, {
  -sin(angle), cos(
  angle)
  }
  };
  return Tensor<2, 2>(t);
  }
 
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)

Class for remembering material state/properties at each quadrature point

  template<int dim>
  struct PointHistory
  {
  double old_phiphi_stress;
  double first_eta;
  double new_eta;
  double G;
  };
 

Primary class of this problem

  template<int dim>
  class StokesProblem
  {
  public:
  StokesProblem(const unsigned int degree);
  void run();
 
  private:
  void setup_dofs();
  void assemble_system();
  void solve();
  void output_results() const;
  void refine_mesh();
  void solution_stesses();
  void smooth_eta_field(std::vector<bool> failing_cells);
 
  void setup_initial_mesh();
  void do_elastic_steps();
  void do_flow_step();
  void update_time_interval();
  void initialize_eta_and_G();
  void move_mesh();
  void do_ellipse_fits();
  void append_physical_times(int max_plastic);
  void write_vertices(unsigned char);
  void write_mesh();
  void setup_quadrature_point_history();
  void update_quadrature_point_history();
 
  const unsigned int degree;
 
  const MappingQ1<dim> mapping;
  DoFHandler<dim> dof_handler;
  unsigned int n_u = 0, n_p = 0;
  unsigned int plastic_iteration = 0;
  unsigned int last_max_plasticity = 0;
 
  QGauss<dim> quadrature_formula;
  std::vector< std::vector <Vector<double> > > quad_viscosities; // Indices for this object are [cell][q][q coords, eta]
  std::vector<double> cell_viscosities; // This vector is only used for output, not FE computations
  std::vector<PointHistory<dim> > quadrature_point_history;
 
 
  BlockSparsityPattern sparsity_pattern;
  BlockSparseMatrix<double> system_matrix;
 
  BlockVector<double> solution;
  BlockVector<double> system_rhs;
 
  std::shared_ptr<typename InnerPreconditioner<dim>::type> A_preconditioner;
 
  ellipsoid_fit<dim> ellipsoid;
  };
 
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation

Class for boundary conditions and rhs

  template<int dim>
  class BoundaryValuesP: public Function<dim>
  {
  public:
  BoundaryValuesP() :
  Function<dim>(dim + 1)
  {
  }
 
  virtual double value(const Point<dim> &p,
  const unsigned int component = 0) const override;
 
  virtual void vector_value(const Point<dim> &p, Vector<double> &value) const override;
  };
 
  template<int dim>
  double BoundaryValuesP<dim>::value(const Point<dim> &p,
  const unsigned int component) const
  {
  Assert(component < this->n_components,
  ExcIndexRange (component, 0, this->n_components));
  (void)p;
  (void)component;
 
  Assert(p[0] >= -10.0, ExcLowerRangeType<double>(p[0], -10.0)); //value of -10 is to permit some small numerical error moving nodes left of x=0; a << value is in fact sufficient
 
  return 0;
  }
 
  template<int dim>
  void BoundaryValuesP<dim>::vector_value(const Point<dim> &p,
  Vector<double> &values) const
  {
  for (unsigned int c = 0; c < this->n_components; ++c)
  values(c) = BoundaryValuesP<dim>::value(p, c);
  }
 
  template<int dim>
  class RightHandSide: public Function<dim>
  {
  public:
  RightHandSide () : Function<dim>(dim+1) {}
 
  using Function<dim>::value;
  using Function<dim>::vector_value;
 
  virtual double value(const Point<dim> &p, const unsigned int component,
  A_Grav_namespace::AnalyticGravity<dim> *aGrav) const;
 
  virtual void vector_value(const Point<dim> &p, Vector<double> &value,
  A_Grav_namespace::AnalyticGravity<dim> *aGrav) const;
 
  virtual void vector_value_list(const std::vector<Point<dim> > &points,
  std::vector<Vector<double> > &values,
  A_Grav_namespace::AnalyticGravity<dim> *aGrav) const;
 
  };
 
  template<int dim>
  double RightHandSide<dim>::value(const Point<dim> &p,
  const unsigned int component,
  A_Grav_namespace::AnalyticGravity<dim> *aGrav) const
  {
 
  std::vector<double> temp_vector(2);
  aGrav->get_gravity(p, temp_vector);
 
  if (component == 0)
  {
  return temp_vector[0] + system_parameters::omegasquared * p[0]; // * 1.2805;
  }
  else
  {
  if (component == 1)
  return temp_vector[1];
  else
  return 0;
  }
  }
 
  template<int dim>
  void RightHandSide<dim>::vector_value(const Point<dim> &p,
  Vector<double> &values,
  A_Grav_namespace::AnalyticGravity<dim> *aGrav) const
  {
  for (unsigned int c = 0; c < this->n_components; ++c)
  values(c) = RightHandSide<dim>::value(p, c, aGrav);
  }
 
  template<int dim>
  void RightHandSide<dim>::vector_value_list(
  const std::vector<Point<dim> > &points,
  std::vector<Vector<double> > &values,
  A_Grav_namespace::AnalyticGravity<dim> *aGrav) const
  {
virtual void vector_value_list(const std::vector< Point< dim > > &points, std::vector< Vector< RangeNumberType > > &values) const
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &values) const
Definition point.h:111
#define Assert(cond, exc)
spacedim const Point< spacedim > & p
Definition grid_tools.h:980

check whether component is in the valid range is up to the derived class

  Assert(values.size() == points.size(),
  ExcDimensionMismatch(values.size(), points.size()));
 
  for (unsigned int i = 0; i < points.size(); ++i)
  this->vector_value(points[i], values[i], aGrav);
  }
 

Class for linear solvers and preconditioners

  template<class Matrix, class Preconditioner>
  class InverseMatrix: public Subscriptor
  {
  public:
  InverseMatrix(const Matrix &m, const Preconditioner &preconditioner);
 
  void vmult(Vector<double> &dst, const Vector<double> &src) const;
 
  private:
  const SmartPointer<const Preconditioner> preconditioner;
  };
 
  template<class Matrix, class Preconditioner>
  InverseMatrix<Matrix, Preconditioner>::InverseMatrix(const Matrix &m,
  const Preconditioner &preconditioner) :
  matrix(&m), preconditioner(&preconditioner)
  {
  }
 
  template<class Matrix, class Preconditioner>
  void InverseMatrix<Matrix, Preconditioner>::vmult(Vector<double> &dst,
  const Vector<double> &src) const
  {
  SolverControl solver_control(1000 * src.size(), 1e-9 * src.l2_norm());
 
  SolverCG<> cg(solver_control);
 
  dst = 0;
 
  cg.solve(*matrix, dst, src, *preconditioner);
  }
 

Class for the SchurComplement

  template<class Preconditioner>
  class SchurComplement: public Subscriptor
  {
  public:
  SchurComplement(const BlockSparseMatrix<double> &system_matrix,
  const InverseMatrix<SparseMatrix<double>, Preconditioner> &A_inverse);
 
  void vmult(Vector<double> &dst, const Vector<double> &src) const;
 
  private:
  const SmartPointer<const InverseMatrix<SparseMatrix<double>, Preconditioner> > A_inverse;
 
  mutable Vector<double> tmp1, tmp2;
  };
 
  template<class Preconditioner>
  SchurComplement<Preconditioner>::SchurComplement(
  const BlockSparseMatrix<double> &system_matrix,
  const InverseMatrix<SparseMatrix<double>, Preconditioner> &A_inverse) :
  system_matrix(&system_matrix), A_inverse(&A_inverse), tmp1(
  system_matrix.block(0, 0).m()), tmp2(
  system_matrix.block(0, 0).m())
  {
  }
 
  template<class Preconditioner>
  void SchurComplement<Preconditioner>::vmult(Vector<double> &dst,
  const Vector<double> &src) const
  {
  system_matrix->block(0, 1).vmult(tmp1, src);
  A_inverse->vmult(tmp2, tmp1);
  system_matrix->block(1, 0).vmult(dst, tmp2);
  }
 

StokesProblem::StokesProblem

  template<int dim>
  StokesProblem<dim>::StokesProblem(const unsigned int degree) :
  degree(degree),
  triangulation(Triangulation<dim>::maximum_smoothing),
  mapping(),
  fe(FE_Q<dim>(degree + 1), dim, FE_Q<dim>(degree), 1),
  dof_handler(triangulation),
  quadrature_formula(degree + 2),
  ellipsoid(&triangulation)
  {}
 
Definition fe_q.h:550

Set up dofs

  template<int dim>
  void StokesProblem<dim>::setup_dofs()
  {
  A_preconditioner.reset();
  system_matrix.clear();
 
  dof_handler.distribute_dofs(fe);
 
  std::vector<unsigned int> block_component(dim + 1, 0);
  block_component[dim] = 1;
  DoFRenumbering::component_wise(dof_handler, block_component);
 
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
void Cuthill_McKee(DoFHandler< dim, spacedim > &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false, const std::vector< types::global_dof_index > &starting_indices=std::vector< types::global_dof_index >())

========================================Apply Boundary Conditions=====================================

  {
  constraints.clear();
  std::vector<bool> component_maskP(dim + 1, false);
  component_maskP[dim] = true;
  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
  BoundaryValuesP<dim>(), constraints, component_maskP);
  }
  {
  std::set<types::boundary_id> no_normal_flux_boundaries;
  no_normal_flux_boundaries.insert(99);
  no_normal_flux_boundaries, constraints);
  }
 
  constraints.close();
 
  std::vector<types::global_dof_index> dofs_per_block = DoFTools::count_dofs_per_fe_block(dof_handler, block_component);
  n_u = dofs_per_block[0];
  n_p = dofs_per_block[1];
 
  std::cout << " Number of active cells: " << triangulation.n_active_cells()
  << std::endl << " Number of degrees of freedom: "
  << dof_handler.n_dofs() << " (" << n_u << '+' << n_p << ')'
  << std::endl;
 
  {
 
  csp.block(0, 0).reinit(n_u, n_u);
  csp.block(1, 0).reinit(n_p, n_u);
  csp.block(0, 1).reinit(n_u, n_p);
  csp.block(1, 1).reinit(n_p, n_p);
 
  csp.collect_sizes();
 
  DoFTools::make_sparsity_pattern(dof_handler, csp, constraints, false);
  sparsity_pattern.copy_from(csp);
  }
 
  system_matrix.reinit(sparsity_pattern);
 
  solution.reinit(2);
  solution.block(0).reinit(n_u);
  solution.block(1).reinit(n_p);
  solution.collect_sizes();
 
  system_rhs.reinit(2);
  system_rhs.block(0).reinit(n_u);
  system_rhs.block(1).reinit(n_p);
  system_rhs.collect_sizes();
 
  }
 
void compute_no_normal_flux_constraints(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int first_vector_component, const std::set< types::boundary_id > &boundary_ids, AffineConstraints< double > &constraints, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
std::vector< types::global_dof_index > count_dofs_per_fe_block(const DoFHandler< dim, spacedim > &dof, const std::vector< unsigned int > &target_block=std::vector< unsigned int >())
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask={})

Viscosity and Shear modulus functions

  template<int dim>
  class Rheology
  {
  public:
  double get_eta(double &r, double &z);
  double get_G(unsigned int mat_id);
 
  private:
  std::vector<double> get_manual_eta_profile();
 
  };
 
  template<int dim>
  std::vector<double> Rheology<dim>::get_manual_eta_profile()
  {
  vector<double> etas;
 
  for (unsigned int i=0; i < system_parameters::sizeof_depths_eta; ++i)
  {
  etas.push_back(system_parameters::depths_eta[i]);
  etas.push_back(system_parameters::eta_kinks[i]);
  }
  return etas;
  }
 
  template<int dim>
  double Rheology<dim>::get_eta(double &r, double &z)
  {

compute local depth

  double ecc = system_parameters::q_axes[0] / system_parameters::p_axes[0];
  double Rminusr = system_parameters::q_axes[0] - system_parameters::p_axes[0];
  double approx_a = std::sqrt(r * r + z * z * ecc * ecc);
  double group1 = r * r + z * z - Rminusr * Rminusr;
 
  double a0 = approx_a;
  double error = 10000;
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)

While loop finds the a axis of the "isodepth" ellipse for which the input point is on the surface. An "isodepth" ellipse is defined as one whose axes a,b are related to the global axes A, B by: A-h = B-h

  if ((r > system_parameters::q_axes[0] - system_parameters::depths_eta.back()) ||
  (z > system_parameters::p_axes[0] - system_parameters::depths_eta.back()))
  {
 
  double eps = 10.0;
  while (error >= eps)
  {
  double a02 = a0 * a0;
  double a03 = a0 * a02;
  double a04 = a0 * a03;
  double fofa = a04 - (2 * Rminusr * a03) - (group1 * a02)
  + (2 * r * r * Rminusr * a0) - (r * r * Rminusr * Rminusr);
  double fprimeofa = 4 * a03 - (6 * Rminusr * a02) - (2 * group1 * a0)
  + (2 * r * r * Rminusr);
  double deltaa = -fofa / fprimeofa;
  a0 += deltaa;
  error = std::abs(deltaa);
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)

cout << "error = " << error << endl;

  }
  }
  else
  {
  a0 = 0.0;
  }
 
  double local_depth = system_parameters::q_axes[0] - a0;
  if (local_depth < 0)
  local_depth = 0;
 
  if (local_depth > system_parameters::depths_eta.back())
  {
  if (system_parameters::eta_kinks.back() < system_parameters::eta_floor)
  return system_parameters::eta_floor;
  else if (system_parameters::eta_kinks.back() > system_parameters::eta_ceiling)
  return system_parameters::eta_ceiling;
  else
  return system_parameters::eta_kinks.back();
  }
 
  std::vector<double> viscosity_function = get_manual_eta_profile();
 
  unsigned int n_visc_kinks = viscosity_function.size() / 2;
 

find the correct interval to do the interpolation in

  int n_minus_one = -1;
  for (unsigned int n = 1; n <= n_visc_kinks; ++n)
  {
  unsigned int ndeep = 2 * n - 2;
  unsigned int nshallow = 2 * n;
  if (local_depth >= viscosity_function[ndeep] && local_depth <= viscosity_function[nshallow])
  n_minus_one = ndeep;
  }
 

find the viscosity interpolation

  if (n_minus_one == -1)
  return system_parameters::eta_ceiling;
  else
  {
  double visc_exponent =
  (viscosity_function[n_minus_one]
  - local_depth)
  / (viscosity_function[n_minus_one]
  - viscosity_function[n_minus_one + 2]);
  double visc_base = viscosity_function[n_minus_one + 3]
  / viscosity_function[n_minus_one + 1];

This is the true viscosity given the thermal profile

  double true_eta = viscosity_function[n_minus_one + 1] * std::pow(visc_base, visc_exponent);
 
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)

Implement latitude-dependence viscosity

  if (system_parameters::lat_dependence)
  {
  double lat = 180 / PI * std::atan(z / r);
  if (lat > 80)
  lat = 80;
  double T_eq = 155;
  double T_surf = T_eq * std::sqrt( std::sqrt( std::cos( PI / 180 * lat ) ) );
  double taper_depth = 40000;
  double surface_taper = (taper_depth - local_depth) / taper_depth;
  if (surface_taper < 0)
  surface_taper = 0;
  double log_eta_contrast = surface_taper * system_parameters::eta_Ea * 52.5365 * (T_eq - T_surf) / T_eq / T_surf;
  true_eta *= std::pow(10, log_eta_contrast);
  }
 
  if (true_eta > system_parameters::eta_ceiling)
  return system_parameters::eta_ceiling;
  else if (true_eta < system_parameters::eta_floor)
  return system_parameters::eta_floor;
  else
  return true_eta;
  }
  }
 
 
  template<int dim>
  double Rheology<dim>::get_G(unsigned int mat_id)
  {
  return system_parameters::G[mat_id];
  }
 
 

Initialize the eta and G parts of the quadrature_point_history object

  template<int dim>
  void StokesProblem<dim>::initialize_eta_and_G()
  {
  FEValues<dim> fe_values(fe, quadrature_formula, update_quadrature_points);
 
  const unsigned int n_q_points = quadrature_formula.size();
  Rheology<dim> rheology;
 
  dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
  {
  PointHistory<dim> *local_quadrature_points_history =
  reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
  Assert(
  local_quadrature_points_history >= &quadrature_point_history.front(),
  ExcInternalError());
  Assert(
  local_quadrature_points_history < &quadrature_point_history.back(),
  ExcInternalError());
  fe_values.reinit(cell);
 
  for (unsigned int q = 0; q < n_q_points; ++q)
  {
 
  double r_value = fe_values.quadrature_point(q)[0];
  double z_value = fe_values.quadrature_point(q)[1];
 
typename ActiveSelector::active_cell_iterator active_cell_iterator
@ update_quadrature_points
Transformed quadrature points.

defines local viscosity

  double local_viscosity = 0;
 
  local_viscosity = rheology.get_eta(r_value, z_value);
 
  local_quadrature_points_history[q].first_eta = local_viscosity;
  local_quadrature_points_history[q].new_eta = local_viscosity;
 

defines local shear modulus

  double local_G = 0;
 
  unsigned int mat_id = cell->material_id();
 
  local_G = rheology.get_G(mat_id);
  local_quadrature_points_history[q].G = local_G;
 

initializes the phi-phi stress

  local_quadrature_points_history[q].old_phiphi_stress = 0;
  }
  }
  }
 

====================== ASSEMBLE THE SYSTEM ======================

  template<int dim>
  void StokesProblem<dim>::assemble_system()
  {
  system_matrix = 0;
  system_rhs = 0;
 
  FEValues<dim> fe_values(fe, quadrature_formula,
 
  const unsigned int dofs_per_cell = fe.dofs_per_cell;
 
  const unsigned int n_q_points = quadrature_formula.size();
 
  FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
  Vector<double> local_rhs(dofs_per_cell);
 
  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.

runs the gravity script function

  const RightHandSide<dim> right_hand_side;
 
  A_Grav_namespace::AnalyticGravity<dim> *aGrav =
  new A_Grav_namespace::AnalyticGravity<dim>;
  std::vector<double> grav_parameters;
  grav_parameters.push_back(system_parameters::q_axes[system_parameters::present_timestep * 2 + 0]);
  grav_parameters.push_back(system_parameters::p_axes[system_parameters::present_timestep * 2 + 0]);
  grav_parameters.push_back(system_parameters::q_axes[system_parameters::present_timestep * 2 + 1]);
  grav_parameters.push_back(system_parameters::p_axes[system_parameters::present_timestep * 2 + 1]);
  grav_parameters.push_back(system_parameters::rho[0]);
  grav_parameters.push_back(system_parameters::rho[1]);
 
  std::cout << "Body parameters are: " ;
  for (int i=0; i<6; ++i)
  std::cout << grav_parameters[i] << " ";
  std::cout << endl;
 
  aGrav->setup_vars(grav_parameters);
 
  std::vector<Vector<double> > rhs_values(n_q_points,
  Vector<double>(dim + 1));
 
  const FEValuesExtractors::Vector velocities(0);
  const FEValuesExtractors::Scalar pressure(dim);
 
  std::vector<SymmetricTensor<2, dim> > phi_grads_u(dofs_per_cell);
  std::vector<double> div_phi_u(dofs_per_cell);
  std::vector<Tensor<1, dim> > phi_u(dofs_per_cell);
  std::vector<double> phi_p(dofs_per_cell);
 
  dof_handler.begin_active(), first_cell = dof_handler.begin_active(),
  endc = dof_handler.end();
 
  for (; cell != endc; ++cell)
  {
  PointHistory<dim> *local_quadrature_points_history =
  reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
  Assert(
  local_quadrature_points_history >= &quadrature_point_history.front(),
  ExcInternalError());
  Assert(
  local_quadrature_points_history < &quadrature_point_history.back(),
  ExcInternalError());
 
  double cell_area = cell->measure();
 
  if (cell_area<0)
  append_physical_times(-1); // This writes final line to physical_times.txt if step terminates prematurely
  AssertThrow(cell_area > 0
  ,
  ExcInternalError());
 
 
  unsigned int m_id = cell->material_id();
 
#define AssertThrow(cond, exc)

initializes the rhs vector to the correct g values

  fe_values.reinit(cell);
  right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
  rhs_values, aGrav);
 
  std::vector<Vector<double> > new_viscosities(quadrature_formula.size(), Vector<double>(dim + 1));
 

Finds vertices where the radius is zero DIM

  bool is_singular = false;
  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
  if (cell->face(f)->center()[0] == 0)
  is_singular = true;
 
  if (is_singular == false || system_parameters::cylindrical == false)
  {
  local_matrix = 0;
  local_rhs = 0;
 

===== outputs the local gravity

  std::vector<Point<dim> > quad_points_list(n_q_points);
  quad_points_list = fe_values.get_quadrature_points();
 
  if (plastic_iteration
  == (system_parameters::max_plastic_iterations - 1))
  {
  if (cell != first_cell)
  {
  std::ofstream fout("gravity_field.txt", std::ios::app);
  fout << quad_points_list[0] << " " << rhs_values[0];
  fout.close();
  }
  else
  {
  std::ofstream fout("gravity_field.txt");
  fout << quad_points_list[0] << " " << rhs_values[0];
  fout.close();
  }
  }
 
  for (unsigned int q = 0; q < n_q_points; ++q)
  {
  SymmetricTensor<2, dim> &old_stress =
  local_quadrature_points_history[q].old_stress;
  double &local_old_phiphi_stress =
  local_quadrature_points_history[q].old_phiphi_stress;
  double r_value = fe_values.quadrature_point(q)[0];
 

get local density based on mat id

  double local_density = system_parameters::rho[m_id];
 

defines local viscosities

  double local_viscosity = 0;
  if (plastic_iteration == 0)
  local_viscosity = local_quadrature_points_history[q].first_eta;
  else
  local_viscosity = local_quadrature_points_history[q].new_eta;
 

Define the local viscoelastic constants

  double local_eta_ve = 2
  / ((1 / local_viscosity)
  + (1 / local_quadrature_points_history[q].G
  / system_parameters::current_time_interval));
  double local_chi_ve = 1
  / (1
  + (local_quadrature_points_history[q].G
  * system_parameters::current_time_interval
  / local_viscosity));
 
  for (unsigned int k = 0; k < dofs_per_cell; ++k)
  {
  phi_grads_u[k] = fe_values[velocities].symmetric_gradient(k,
  q);
  div_phi_u[k] = (fe_values[velocities].divergence(k, q));
  phi_u[k] = (fe_values[velocities].value(k, q));
  if (system_parameters::cylindrical == true)
  {
  div_phi_u[k] *= (r_value);
  div_phi_u[k] += (phi_u[k][0]);
  }
  phi_p[k] = fe_values[pressure].value(k, q);
  }
 
  for (unsigned int i = 0; i < dofs_per_cell; ++i)
  {
  for (unsigned int j = 0; j <= i; ++j)
  {
  if (system_parameters::cylindrical == true)
  {
  local_matrix(i, j) += (phi_grads_u[i]
  * phi_grads_u[j] * 2 * local_eta_ve
  * r_value
  + 2 * phi_u[i][0] * phi_u[j][0]
  * local_eta_ve / r_value
  - div_phi_u[i] * phi_p[j]
  * system_parameters::pressure_scale
  - phi_p[i] * div_phi_u[j]
  * system_parameters::pressure_scale
  + phi_p[i] * phi_p[j] * r_value
  * system_parameters::pressure_scale)
  * fe_values.JxW(q);
  }
  else
  {
  local_matrix(i, j) += (phi_grads_u[i]
  * phi_grads_u[j] * 2 * local_eta_ve
  - div_phi_u[i] * phi_p[j]
  * system_parameters::pressure_scale
  - phi_p[i] * div_phi_u[j]
  * system_parameters::pressure_scale
  + phi_p[i] * phi_p[j]) * fe_values.JxW(q);
  }
  }
  if (system_parameters::cylindrical == true)
  {
  const unsigned int component_i =
  fe.system_to_component_index(i).first;
  local_rhs(i) += (fe_values.shape_value(i, q)
  * rhs_values[q](component_i) * r_value
  * local_density
  - local_chi_ve * phi_grads_u[i] * old_stress
  * r_value
  - local_chi_ve * phi_u[i][0]
  * local_old_phiphi_stress)
  * fe_values.JxW(q);
  }
  else
  {
  const unsigned int component_i =
  fe.system_to_component_index(i).first;
  local_rhs(i) += fe_values.shape_value(i, q)
  * rhs_values[q](component_i) * fe_values.JxW(q)
  * local_density;
  }
  }
  }
  } // end of non-singular
  else
  {
  local_matrix = 0;
  local_rhs = 0;
 

===== outputs the local gravity

  std::vector<Point<dim> > quad_points_list(n_q_points);
  quad_points_list = fe_values.get_quadrature_points();
 
  for (unsigned int q = 0; q < n_q_points; ++q)
  {
  const SymmetricTensor<2, dim> &old_stress =
  local_quadrature_points_history[q].old_stress;
  double &local_old_phiphi_stress =
  local_quadrature_points_history[q].old_phiphi_stress;
  double r_value = fe_values.quadrature_point(q)[0];
 
  double local_density = system_parameters::rho[m_id];
 

defines local viscosities

  double local_viscosity = 0;
  if (plastic_iteration == 0)
  {
  local_viscosity = local_quadrature_points_history[q].first_eta;
  }
  else
  local_viscosity = local_quadrature_points_history[q].new_eta;
 

Define the local viscoelastic constants

  double local_eta_ve = 2
  / ((1 / local_viscosity)
  + (1 / local_quadrature_points_history[q].G
  / system_parameters::current_time_interval));
  double local_chi_ve = 1
  / (1
  + (local_quadrature_points_history[q].G
  * system_parameters::current_time_interval
  / local_viscosity));
 
  for (unsigned int k = 0; k < dofs_per_cell; ++k)
  {
  phi_grads_u[k] = fe_values[velocities].symmetric_gradient(k,
  q);
  div_phi_u[k] = (fe_values[velocities].divergence(k, q));
  phi_u[k] = (fe_values[velocities].value(k, q));
  if (system_parameters::cylindrical == true)
  {
  div_phi_u[k] *= (r_value);
  div_phi_u[k] += (phi_u[k][0]);
  }
  phi_p[k] = fe_values[pressure].value(k, q);
  }
 
  for (unsigned int i = 0; i < dofs_per_cell; ++i)
  {
  for (unsigned int j = 0; j <= i; ++j)
  {
  if (system_parameters::cylindrical == true)
  {
  local_matrix(i, j) += (phi_grads_u[i]
  * phi_grads_u[j] * 2 * local_eta_ve
  * r_value
  + 2 * phi_u[i][0] * phi_u[j][0]
  * local_eta_ve / r_value
  - div_phi_u[i] * phi_p[j]
  * system_parameters::pressure_scale
  - phi_p[i] * div_phi_u[j]
  * system_parameters::pressure_scale
  + phi_p[i] * phi_p[j] * r_value
  * system_parameters::pressure_scale)
  * fe_values.JxW(q);
  }
  else
  {
  local_matrix(i, j) += (phi_grads_u[i]
  * phi_grads_u[j] * 2 * local_eta_ve
  - div_phi_u[i] * phi_p[j]
  * system_parameters::pressure_scale
  - phi_p[i] * div_phi_u[j]
  * system_parameters::pressure_scale
  + phi_p[i] * phi_p[j]) * fe_values.JxW(q);
  }
  }
  if (system_parameters::cylindrical == true)
  {
  const unsigned int component_i =
  fe.system_to_component_index(i).first;
  local_rhs(i) += (fe_values.shape_value(i, q)
  * rhs_values[q](component_i) * r_value
  * local_density
  - local_chi_ve * phi_grads_u[i] * old_stress
  * r_value
  - local_chi_ve * phi_u[i][0]
  * local_old_phiphi_stress)
  * fe_values.JxW(q);
  }
  else
  {
  const unsigned int component_i =
  fe.system_to_component_index(i).first;
  local_rhs(i) += fe_values.shape_value(i, q)
  * rhs_values[q](component_i) * fe_values.JxW(q)
  * local_density;
  }
  }
  }
  } // end of singular
 
  for (unsigned int i = 0; i < dofs_per_cell; ++i)
  for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
  local_matrix(i, j) = local_matrix(j, i);
 
  cell->get_dof_indices(local_dof_indices);
  constraints.distribute_local_to_global(local_matrix, local_rhs,
  local_dof_indices, system_matrix, system_rhs);
  }
 
  std::cout << " Computing preconditioner..." << std::endl << std::flush;
 
  A_preconditioner = std::shared_ptr<
  typename InnerPreconditioner<dim>::type>(
  new typename InnerPreconditioner<dim>::type());
  A_preconditioner->initialize(system_matrix.block(0, 0),
  typename InnerPreconditioner<dim>::type::AdditionalData());
 
  delete aGrav;
  }
 

====================== SOLVER ======================

  template<int dim>
  void StokesProblem<dim>::solve()
  {
  const InverseMatrix<SparseMatrix<double>,
  typename InnerPreconditioner<dim>::type> A_inverse(
  system_matrix.block(0, 0), *A_preconditioner);
  Vector<double> tmp(solution.block(0).size());
 
  {
  Vector<double> schur_rhs(solution.block(1).size());
  A_inverse.vmult(tmp, system_rhs.block(0));
  system_matrix.block(1, 0).vmult(schur_rhs, tmp);
  schur_rhs -= system_rhs.block(1);
 
  SchurComplement<typename InnerPreconditioner<dim>::type> schur_complement(
  system_matrix, A_inverse);
 
  int n_iterations = system_parameters::iteration_coefficient
  * solution.block(1).size();
  double tolerance_goal = system_parameters::tolerance_coefficient
  * schur_rhs.l2_norm();
 
  SolverControl solver_control(n_iterations, tolerance_goal);
  SolverCG<> cg(solver_control);
 
  std::cout << "\nMax iterations and tolerance are: " << n_iterations
  << " and " << tolerance_goal << std::endl;
 
  SparseILU<double> preconditioner;
  preconditioner.initialize(system_matrix.block(1, 1),
 
  InverseMatrix<SparseMatrix<double>, SparseILU<double> > m_inverse(
  system_matrix.block(1, 1), preconditioner);
 
  cg.solve(schur_complement, solution.block(1), schur_rhs, m_inverse);
 
  constraints.distribute(solution);
 
 
  std::cout << " " << solver_control.last_step()
  << " outer CG Schur complement iterations for pressure"
  << std::endl;
  }
 
  {
  system_matrix.block(0, 1).vmult(tmp, solution.block(1));
  tmp *= -1;
  tmp += system_rhs.block(0);
 
  A_inverse.vmult(solution.block(0), tmp);
  constraints.distribute(solution);
  solution.block(1) *= (system_parameters::pressure_scale);
  }
  }
 
typename SparseLUDecomposition< number >::AdditionalData AdditionalData
Definition sparse_ilu.h:79
void initialize(const SparseMatrix< somenumber > &matrix, const AdditionalData &parameters=AdditionalData())
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)

====================== OUTPUT RESULTS ======================

  template<int dim>
  void StokesProblem<dim>::output_results() const
  {
  std::vector < std::string > solution_names(dim, "velocity");
  solution_names.push_back("pressure");
 
  std::vector<DataComponentInterpretation::DataComponentInterpretation> data_component_interpretation(
  data_component_interpretation.push_back(
 
  DataOut<dim> data_out;
  data_out.attach_dof_handler(dof_handler);
  data_out.add_data_vector(solution, solution_names,
  DataOut<dim>::type_dof_data, data_component_interpretation);
  data_out.build_patches();
 
  std::ostringstream filename;
  if (system_parameters::present_timestep < system_parameters::initial_elastic_iterations)
  {
  filename << system_parameters::output_folder << "/time"
  << Utilities::int_to_string(system_parameters::present_timestep, 2)
  << "_elastic_displacements" << ".txt";
  }
  else
  {
  filename << system_parameters::output_folder << "/time"
  << Utilities::int_to_string(system_parameters::present_timestep, 2)
  << "_flow" << Utilities::int_to_string(plastic_iteration, 2) << ".txt";
  }
 
  std::ofstream output(filename.str().c_str());
  data_out.write_gnuplot(output);
  }
 
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470
STL namespace.

====================== FIND AND WRITE TO FILE THE STRESS TENSOR; IMPLEMENT PLASTICITY ======================

  template<int dim>
  void StokesProblem<dim>::solution_stesses()
  {

note most of this section only works with dim=2

name the output text files

  std::ostringstream stress_output;
  stress_output << system_parameters::output_folder << "/time"
  << Utilities::int_to_string(system_parameters::present_timestep, 2)
  << "_principalstresses" << Utilities::int_to_string(plastic_iteration, 2)
  << ".txt";
  std::ofstream fout_snew(stress_output.str().c_str());
  fout_snew.close();
 
  std::ostringstream stresstensor_output;
  stresstensor_output << system_parameters::output_folder << "/time"
  << Utilities::int_to_string(system_parameters::present_timestep, 2)
  << "_stresstensor" << Utilities::int_to_string(plastic_iteration, 2)
  << ".txt";
  std::ofstream fout_sfull(stresstensor_output.str().c_str());
  fout_sfull.close();
 
  std::ostringstream failed_cells_output;
  failed_cells_output << system_parameters::output_folder << "/time"
  << Utilities::int_to_string(system_parameters::present_timestep, 2)
  << "_failurelocations" << Utilities::int_to_string(plastic_iteration, 2)
  << ".txt";
  std::ofstream fout_failed_cells(failed_cells_output.str().c_str());
  fout_failed_cells.close();
 
  std::ostringstream plastic_eta_output;
  plastic_eta_output << system_parameters::output_folder << "/time"
  << Utilities::int_to_string(system_parameters::present_timestep, 2)
  << "_viscositiesreg" << Utilities::int_to_string(plastic_iteration, 2)
  << ".txt";
  std::ofstream fout_vrnew(plastic_eta_output.str().c_str());
  fout_vrnew.close();
 
  std::ostringstream initial_eta_output;
  if (plastic_iteration == 0)
  {
  initial_eta_output << system_parameters::output_folder << "/time"
  << Utilities::int_to_string(system_parameters::present_timestep, 2)
  << "_baseviscosities.txt";
  std::ofstream fout_baseeta(initial_eta_output.str().c_str());
  fout_baseeta.close();
  }
 
  std::cout << "Running stress calculations for plasticity iteration "
  << plastic_iteration << "...\n";
 

This makes the set of points at which the stress tensor is calculated

  std::vector<Point<dim> > points_list(0);
  std::vector<unsigned int> material_list(0);
  dof_handler.begin_active(), endc = dof_handler.end();

This loop gets the gradients of the velocity field and saves it in the tensor_gradient_? objects DIM

  for (; cell != endc; ++cell)
  {
  points_list.push_back(cell->center());
  material_list.push_back(cell->material_id());
  }

Make the FEValues object to evaluate values and derivatives at quadrature points

  FEValues<dim> fe_values(fe, quadrature_formula,
 

Make the object that will hold the velocities and velocity gradients at the quadrature points

  std::vector < std::vector<Tensor<1, dim> >> velocity_grads(quadrature_formula.size(),
  std::vector < Tensor<1, dim> > (dim + 1));
  std::vector<Vector<double> > velocities(quadrature_formula.size(),
  Vector<double>(dim + 1));

Make the object to find rheology

  Rheology<dim> rheology;
 

Write the solution flow velocity and derivative for each cell

  std::vector<Vector<double> > vector_values(0);
  std::vector < std::vector<Tensor<1, dim> > > gradient_values(0);
  std::vector<bool> failing_cells;

Write the stresses from the previous step into vectors

  std::vector<SymmetricTensor<2, dim>> old_stress;
  std::vector<double> old_phiphi_stress;
  std::vector<double> cell_Gs;
  dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
  {

Makes pointer to data in quadrature_point_history

  PointHistory<dim> *local_quadrature_points_history =
  reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
 
  fe_values.reinit(cell);
  fe_values.get_function_gradients(solution, velocity_grads);
  fe_values.get_function_values(solution, velocities);
  Vector<double> current_cell_velocity(dim+1);
  std::vector<Tensor<1, dim>> current_cell_grads(dim+1);
  SymmetricTensor<2, dim> current_cell_old_stress;
  current_cell_old_stress = 0;
  double current_cell_old_phiphi_stress = 0;
  double cell_area = 0;
 

Averages across each cell to find mean velocities, gradients, and old stresses

  for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
  {
  cell_area += fe_values.JxW(q);
  velocities[q] *= fe_values.JxW(q);
  current_cell_velocity += velocities[q];
  for (unsigned int i = 0; i < (dim+1); ++i)
  {
  velocity_grads[q][i] *= fe_values.JxW(q);
  current_cell_grads[i] += velocity_grads[q][i];
  }
  current_cell_old_stress += local_quadrature_points_history[q].old_stress * fe_values.JxW(q);
  current_cell_old_phiphi_stress += local_quadrature_points_history[q].old_phiphi_stress * fe_values.JxW(q);
  }
  current_cell_velocity /= cell_area;
  for (unsigned int i = 0; i < (dim+1); ++i)
  current_cell_grads[i] /= cell_area;
  current_cell_old_stress /= cell_area;
  current_cell_old_phiphi_stress /= cell_area;
 
  vector_values.push_back(current_cell_velocity);
  gradient_values.push_back(current_cell_grads);
  old_stress.push_back(current_cell_old_stress);
  old_phiphi_stress.push_back(current_cell_old_phiphi_stress);
 

Get cell shear modulus: assumes it's constant for the cell

  unsigned int mat_id = cell->material_id();
  double local_G = rheology.get_G(mat_id);
  cell_Gs.push_back(local_G);
  }
 

tracks where failure occurred

  std::vector<double> reduction_factor;
  unsigned int total_fails = 0;
  if (plastic_iteration == 0)
  cell_viscosities.resize(0);

loop across all the cells to find and adjust eta of failing cells

  for (unsigned int i = 0; i < triangulation.n_active_cells(); ++i)
  {
  double current_cell_viscosity = 0;
 

Fill viscosities vector, analytically if plastic_iteration == 0 and from previous viscosities for later iteration

  if (plastic_iteration == 0)
  {
  double local_viscosity;
  local_viscosity = rheology.get_eta(points_list[i][0], points_list[i][1]);
  current_cell_viscosity = local_viscosity;
  cell_viscosities.push_back(current_cell_viscosity);
  }
  else
  {
  current_cell_viscosity = cell_viscosities[i];
  }
 
 
  double cell_eta_ve = 2
  / ((1 / current_cell_viscosity)
  + (1 / cell_Gs[i]
  / system_parameters::current_time_interval));
  double cell_chi_ve = 1
  / (1
  + (cell_Gs[i]
  * system_parameters::current_time_interval
  / current_cell_viscosity));
 

find local pressure

  double cell_p = vector_values[i].operator()(2);

find stresses tensor makes non-diagonalized local matrix A

  double sigma13 = 0.5
  * (gradient_values[i][0][1] + gradient_values[i][1][0]);
  mat A;
  A << gradient_values[i][0][0] << 0 << sigma13 << endr
  << 0 << vector_values[i].operator()(0) / points_list[i].operator()(0)<< 0 << endr
  << sigma13 << 0 << gradient_values[i][1][1] << endr;
  mat olddevstress;
  olddevstress << old_stress[i][0][0] << 0 << old_stress[i][0][1] << endr
  << 0 << old_phiphi_stress[i] << 0 << endr
  << old_stress[i][0][1] << 0 << old_stress[i][1][1] << endr;
  vec P;
  P << cell_p << cell_p << cell_p;
  mat Pmat = diagmat(P);
  mat B;
  B = (cell_eta_ve * A + cell_chi_ve * olddevstress) - Pmat;
 

finds principal stresses

  vec eigval;
  mat eigvec;
  eig_sym(eigval, eigvec, B);
  double sigma1 = -min(eigval);
  double sigma3 = -max(eigval);
 

Writes text files for principal stresses, full stress tensor, base viscosities

  std::ofstream fout_snew(stress_output.str().c_str(), std::ios::app);
  fout_snew << " " << sigma1 << " " << sigma3 << "\n";
  fout_snew.close();
 
  std::ofstream fout_sfull(stresstensor_output.str().c_str(), std::ios::app);
  fout_sfull << A(0,0) << " " << A(1,1) << " " << A(2,2) << " " << A(0,2) << "\n";
  fout_sfull.close();
 
  if (plastic_iteration == 0)
  {
  std::ofstream fout_baseeta(initial_eta_output.str().c_str(), std::ios::app);
  fout_baseeta << points_list[i]<< " " << current_cell_viscosity << "\n";
  fout_baseeta.close();
  }
 

Finds adjusted effective viscosity

  if (system_parameters::plasticity_on)
  {
  if (system_parameters::failure_criterion == 0) //Apply Byerlee's rule
  {
  if (sigma1 >= 5 * sigma3) // this guarantees that viscosities only go down, never up
  {
  failing_cells.push_back(true);
  double temp_reductionfactor = 1;
  if (sigma3 < 0)
  temp_reductionfactor = 100;
  else
  temp_reductionfactor = 1.9 * sigma1 / 5 / sigma3;
 
  reduction_factor.push_back(temp_reductionfactor);
  total_fails++;
 

Text file of all failure locations

  std::ofstream fout_failed_cells(failed_cells_output.str().c_str(), std::ios::app);
  fout_failed_cells << points_list[i] << "\n";
  fout_failed_cells.close();
  }
  else
  {
  reduction_factor.push_back(1);
  failing_cells.push_back(false);
  }
  }
  else
  {
  if (system_parameters::failure_criterion == 1) //Apply Schultz criterion for frozen sand, RMR=45
  {
  double temp_reductionfactor = 1;
  if (sigma3 < -114037)
  {

std::cout << " ext ";

  failing_cells.push_back(true);
  temp_reductionfactor = 10;
  reduction_factor.push_back(temp_reductionfactor);
  total_fails++;
 

Text file of all failure locations

  std::ofstream fout_failed_cells(failed_cells_output.str().c_str(), std::ios::app);
  fout_failed_cells << points_list[i] << "\n";
  fout_failed_cells.close();
  }
  else
  {
  double sigma_c = 160e6; //Unconfined compressive strength
  double yield_sigma1 = sigma3 + std::sqrt( (3.086 * sigma_c * sigma3) + (0.002 * sigma3 * sigma3) );
  if (sigma1 >= yield_sigma1)
  {

std::cout << " comp ";

  failing_cells.push_back(true);
  temp_reductionfactor = 1.0 * sigma1 / 5 / sigma3;
 
  reduction_factor.push_back(temp_reductionfactor);
  total_fails++;
 

Text file of all failure locations

  std::ofstream fout_failed_cells(failed_cells_output.str().c_str(), std::ios::app);
  fout_failed_cells << points_list[i] << "\n";
  fout_failed_cells.close();
  }
  else
  {
  reduction_factor.push_back(1);
  failing_cells.push_back(false);
  }
  }
  }
  else
  {
  std::cout << "Specified failure criterion not found\n";
  }
  }
  }
  else
  reduction_factor.push_back(1);
  }
 

If there are enough failed cells, update eta at all quadrature points and perform smoothing

  std::cout << " Number of failing cells: " << total_fails << "\n";
  double last_max_plasticity_double = last_max_plasticity;
  double total_fails_double = total_fails;
  double decrease_in_plasticity = ((last_max_plasticity_double - total_fails_double) / last_max_plasticity_double);
  if (plastic_iteration == 0)
  decrease_in_plasticity = 1;
  last_max_plasticity = total_fails;
  if (total_fails <= 100 || decrease_in_plasticity <= 0.2)
  {
  system_parameters::continue_plastic_iterations = false;
  for (unsigned int j=0; j < triangulation.n_active_cells(); ++j)
  {

Writes to file the undisturbed cell viscosities

  std::ofstream fout_vrnew(plastic_eta_output.str().c_str(), std::ios::app);
  fout_vrnew << " " << cell_viscosities[j] << "\n";
  fout_vrnew.close();
  }
  }
  else
  {
  quad_viscosities.resize(triangulation.n_active_cells());

Decrease the eta at quadrature points in failing cells

  unsigned int cell_no = 0;
  dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
  {
  fe_values.reinit(cell);

Make local_quadrature_points_history pointer to the cell data

  PointHistory<dim> *local_quadrature_points_history =
  reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
  Assert(
  local_quadrature_points_history >= &quadrature_point_history.front(),
  ExcInternalError());
  Assert(
  local_quadrature_points_history < &quadrature_point_history.back(),
  ExcInternalError());
 
  quad_viscosities[cell_no].resize(quadrature_formula.size());
 
  for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
  {
  if (plastic_iteration == 0)
  local_quadrature_points_history[q].new_eta = local_quadrature_points_history[q].first_eta;
  local_quadrature_points_history[q].new_eta /= reduction_factor[cell_no];

Prevents viscosities from dropping below the floor necessary for numerical stability

  if (local_quadrature_points_history[q].new_eta < system_parameters::eta_floor)
  local_quadrature_points_history[q].new_eta = system_parameters::eta_floor;
 
  quad_viscosities[cell_no][q].reinit(dim+1);
  for (unsigned int ii=0; ii<dim; ++ii)
  quad_viscosities[cell_no][q](ii) = fe_values.quadrature_point(q)[ii];
  quad_viscosities[cell_no][q](dim) = local_quadrature_points_history[q].new_eta;
  }
  cell_no++;
  }
  smooth_eta_field(failing_cells);
 

Writes to file the smoothed eta field (which is defined at each quadrature point) for each cell

  cell_no = 0;

cell_viscosities.resize(triangulation.n_active_cells(), 0);

  dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
  {
  if (failing_cells[cell_no])
  {
  fe_values.reinit(cell);

Averages across each cell to find mean eta

  double cell_area = 0;
  for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
  {
  cell_area += fe_values.JxW(q);
  cell_viscosities[cell_no] += quad_viscosities[cell_no][q][dim] * fe_values.JxW(q);
  }
  cell_viscosities[cell_no] /= cell_area;
 

Writes to file

  std::ofstream fout_vrnew(plastic_eta_output.str().c_str(), std::ios::app);
  fout_vrnew << " " << cell_viscosities[cell_no] << "\n";
  fout_vrnew.close();
  }
  else
  {
  std::ofstream fout_vrnew(plastic_eta_output.str().c_str(), std::ios::app);
  fout_vrnew << " " << cell_viscosities[cell_no] << "\n";
  fout_vrnew.close();
  }
  cell_no++;
  }
  }
  }
 

====================== SMOOTHES THE VISCOSITY FIELD AT ALL QUADRATURE POINTS ======================

  template<int dim>
  void StokesProblem<dim>::smooth_eta_field(std::vector<bool> failing_cells)
  {
  std::cout << " Smoothing viscosity field...\n";
  unsigned int cell_no = 0;
  dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
  {
  if (failing_cells[cell_no])
  {
  FEValues<dim> fe_values(fe, quadrature_formula, update_quadrature_points);
  PointHistory<dim> *local_quadrature_points_history =
  reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
 

Currently this algorithm does not permit refinement. To permit refinement, daughter cells of neighbors must be identified Find pointers and indices of all cells within certain radius

  bool find_more_cells = true;
  std::vector<bool> cell_touched(triangulation.n_active_cells(), false);
  std::vector< TriaIterator< CellAccessor<dim> > > neighbor_cells;
  std::vector<int> neighbor_indices;
  unsigned int start_cell = 0; // Which cell in the neighbor_cells vector to start from
  int new_cells_found = 0;
  neighbor_cells.push_back(cell);
  neighbor_indices.push_back(cell_no);
  cell_touched[cell_no] = true;
  while (find_more_cells)
  {
  new_cells_found = 0;
  for (unsigned int i = start_cell; i<neighbor_cells.size(); ++i)
  {
  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
  {
  if (!neighbor_cells[i]->face(f)->at_boundary())
  {
  int test_cell_no = neighbor_cells[i]->neighbor_index(f);
  if (!cell_touched[test_cell_no])
  if (cell->center().distance(neighbor_cells[i]->neighbor(f)->center()) < 2 * system_parameters::smoothing_radius)
  {

What to do if another nearby cell is found that hasn't been found before

  neighbor_cells.push_back(neighbor_cells[i]->neighbor(f));
  neighbor_indices.push_back(test_cell_no);
  cell_touched[test_cell_no] = true;
  start_cell++;
  new_cells_found++;
  }
  }
  }
  }
  if (new_cells_found == 0)
  {
  find_more_cells = false;
  }
  else
  start_cell = neighbor_cells.size() - new_cells_found;
  }
 
  fe_values.reinit(cell);

Collect the viscosities at nearby quadrature points

  for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
  {
  std::vector<double> nearby_etas_q;
  for (unsigned int i = 0; i<neighbor_indices.size(); ++i)
  for (unsigned int j=0; j<quadrature_formula.size(); ++j)
  {
  Point<dim> test_q;
  for (unsigned int d=0; d<dim; ++d)
  test_q(d) = quad_viscosities[neighbor_indices[i]][j][d];
  double qq_distance = fe_values.quadrature_point(q).distance(test_q);
  if (qq_distance < system_parameters::smoothing_radius)
  nearby_etas_q.push_back(quad_viscosities[neighbor_indices[i]][j][dim]);
  }
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const

Write smoothed viscosities to quadrature_points_history; simple boxcar function is the smoothing kernel

  double mean_eta = 0;
  for (unsigned int l = 0; l<nearby_etas_q.size(); ++l)
  {
  mean_eta += nearby_etas_q[l];
  }
  mean_eta /= nearby_etas_q.size();
  local_quadrature_points_history[q].new_eta = mean_eta;

std::cout << local_quadrature_points_history[q].new_eta << " ";

  }
  }
  cell_no++;
  }
  }
 

====================== SAVE STRESS TENSOR AT QUADRATURE POINTS ======================

  template<int dim>
  void StokesProblem<dim>::update_quadrature_point_history()
  {
  std::cout << " Updating stress field...";
 
  FEValues<dim> fe_values(fe, quadrature_formula,
 

Make the object that will hold the velocity gradients

  std::vector < std::vector<Tensor<1, dim> >> velocity_grads(quadrature_formula.size(),
  std::vector < Tensor<1, dim> > (dim + 1));
  std::vector<Vector<double> > velocities(quadrature_formula.size(),
  Vector<double>(dim + 1));
 
  dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
  {
  PointHistory<dim> *local_quadrature_points_history =
  reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
  Assert(
  local_quadrature_points_history >= &quadrature_point_history.front(),
  ExcInternalError());
  Assert(
  local_quadrature_points_history < &quadrature_point_history.back(),
  ExcInternalError());
 
  fe_values.reinit(cell);
  fe_values.get_function_gradients(solution, velocity_grads);
  fe_values.get_function_values(solution, velocities);
 
  for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
  {

Define the local viscoelastic constants

  double local_eta_ve = 2
  / ((1 / local_quadrature_points_history[q].new_eta)
  + (1 / local_quadrature_points_history[q].G
  / system_parameters::current_time_interval));
  double local_chi_ve =
  1
  / (1
  + (local_quadrature_points_history[q].G
  * system_parameters::current_time_interval
  / local_quadrature_points_history[q].new_eta));
 

Compute new stress at each quadrature point

  for (unsigned int i = 0; i < dim; ++i)
  new_stress[i][i] =
  local_eta_ve * velocity_grads[q][i][i]
  + local_chi_ve
  * local_quadrature_points_history[q].old_stress[i][i];
 
  for (unsigned int i = 0; i < dim; ++i)
  for (unsigned int j = i + 1; j < dim; ++j)
  new_stress[i][j] =
  local_eta_ve
  * (velocity_grads[q][i][j]
  + velocity_grads[q][j][i]) / 2
  + local_chi_ve
  * local_quadrature_points_history[q].old_stress[i][j];
 

Rotate new stress

  AuxFunctions<dim> rotation_object;
  const Tensor<2, dim> rotation = rotation_object.get_rotation_matrix(
  velocity_grads[q]);
  const SymmetricTensor<2, dim> rotated_new_stress = symmetrize(
  transpose(rotation)
  * static_cast<Tensor<2, dim> >(new_stress)
  * rotation);
  local_quadrature_points_history[q].old_stress = rotated_new_stress;
 
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)

For axisymmetric case, make the phi-phi element of stress tensor

  local_quadrature_points_history[q].old_phiphi_stress =
  (2 * local_eta_ve * velocities[q](0)
  / fe_values.quadrature_point(q)[0]
  + local_chi_ve
  * local_quadrature_points_history[q].old_phiphi_stress);
  }
  }
  }
 

====================== REDEFINE THE TIME INTERVAL FOR THE VISCOUS STEPS ======================

  template<int dim>
  void StokesProblem<dim>::update_time_interval()
  {
  double move_goal_per_step = system_parameters::initial_disp_target;
  if (system_parameters::present_timestep > system_parameters::initial_elastic_iterations)
  {
  move_goal_per_step = system_parameters::initial_disp_target -
  ((system_parameters::initial_disp_target - system_parameters::final_disp_target) /
  system_parameters::total_viscous_steps *
  (system_parameters::present_timestep - system_parameters::initial_elastic_iterations));
  }
 
  double zero_tolerance = 1e-3;
  double max_velocity = 0;
  dof_handler.begin_active(); cell != dof_handler.end(); ++cell)// loop over all cells
  {
  if (cell->at_boundary())
  {
  int zero_faces = 0;
  for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
  for (unsigned int i=0; i<dim; ++i)
  if (fabs(cell->face(f)->center()[i]) < zero_tolerance)
  zero_faces++;
  if (zero_faces==0)
  {
  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
  {
  Point<dim> vertex_velocity;
  Point<dim> vertex_position;
  for (unsigned int d = 0; d < dim; ++d)
  {
  vertex_velocity[d] = solution(cell->vertex_dof_index(v, d));
  vertex_position[d] = cell->vertex(v)[d];
  }

velocity to be evaluated is the radial component of a surface vertex

  double local_velocity = 0;
  for (unsigned int d = 0; d < dim; ++d)
  {
  local_velocity += vertex_velocity[d] * vertex_position [d];
  }
  local_velocity /= std::sqrt( vertex_position.square() );
  if (local_velocity < 0)
  local_velocity *= -1;
  if (local_velocity > max_velocity)
  {
  max_velocity = local_velocity;
  }
  }
  }
  }
  }

NOTE: It is possible for this time interval to be very different from that used in the viscoelasticity calculation.

  system_parameters::current_time_interval = move_goal_per_step / max_velocity;
  double step_time_yr = system_parameters::current_time_interval / SECSINYEAR;
  std::cout << "Timestep interval changed to: "
  << step_time_yr
  << " years\n";
  }
 

====================== MOVE MESH ======================

  template<int dim>
  void StokesProblem<dim>::move_mesh()
  {
 
  std::cout << "\n" << " Moving mesh...\n";
  std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
  dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
  if (vertex_touched[cell->vertex_index(v)] == false)
  {
  vertex_touched[cell->vertex_index(v)] = true;
 
  Point<dim> vertex_displacement;
  for (unsigned int d = 0; d < dim; ++d)
  vertex_displacement[d] = solution(
  cell->vertex_dof_index(v, d));
  cell->vertex(v) += vertex_displacement
  * system_parameters::current_time_interval;
  }
  }
 

====================== WRITE MESH TO FILE ======================

  template<int dim>
  void StokesProblem<dim>::write_mesh()
  {

output mesh in ucd

  std::ostringstream initial_mesh_file;
  initial_mesh_file << system_parameters::output_folder << "/time" <<
  Utilities::int_to_string(system_parameters::present_timestep, 2) <<
  "_mesh.inp";
  std::ofstream out_ucd (initial_mesh_file.str().c_str());
  GridOut grid_out;
  grid_out.write_ucd (triangulation, out_ucd);
  }
 
void write_ucd(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition grid_out.cc:1126

====================== FIT ELLIPSE TO SURFACE AND WRITE RADII TO FILE ======================

  template<int dim>
  void StokesProblem<dim>::do_ellipse_fits()
  {
  std::ostringstream ellipses_filename;
  ellipses_filename << system_parameters::output_folder << "/ellipse_fits.txt";

Find ellipsoidal axes for all layers

  std::vector<double> ellipse_axes(0);

compute fit to boundary 0, 1, 2 ...

  std::cout << endl;
  for (unsigned int i = 0; i<system_parameters::sizeof_material_id; ++i)
  {
  ellipsoid.compute_fit(ellipse_axes, system_parameters::material_id[i]);
  system_parameters::q_axes.push_back(ellipse_axes[0]);
  system_parameters::p_axes.push_back(ellipse_axes[1]);
 
  std::cout << "a_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[0]
  << " " << " c_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[1] << std::endl;
  ellipse_axes.clear();
 
  std::ofstream fout_ellipses(ellipses_filename.str().c_str(), std::ios::app);
  fout_ellipses << system_parameters::present_timestep << " a_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[0]
  << " " << " c_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[1] << endl;
  fout_ellipses.close();
  }
  }
 

====================== APPEND LINE TO PHYSICAL_TIMES.TXT FILE WITH STEP NUMBER, PHYSICAL TIME, AND # PLASTIC ITERATIONS ======================

  template<int dim>
  void StokesProblem<dim>::append_physical_times(int max_plastic)
  {
  std::ostringstream times_filename;
  times_filename << system_parameters::output_folder << "/physical_times.txt";
  std::ofstream fout_times(times_filename.str().c_str(), std::ios::app);
  fout_times << system_parameters::present_timestep << " "
  << system_parameters::present_time/SECSINYEAR << " "
  << max_plastic << "\n";

<< system_parameters::q_axes[0] << " " << system_parameters::p_axes[0] << " " << system_parameters::q_axes[1] << " " << system_parameters::p_axes[1] << "\n";

  fout_times.close();
  }
 

====================== WRITE VERTICES TO FILE ======================

  template<int dim>
  void StokesProblem<dim>::write_vertices(unsigned char boundary_that_we_need)
  {
  std::ostringstream vertices_output;
  vertices_output << system_parameters::output_folder << "/time" <<
  Utilities::int_to_string(system_parameters::present_timestep, 2) << "_" <<
  Utilities::int_to_string(boundary_that_we_need, 2) <<
  "_surface.txt";
  std::ofstream fout_final_vertices(vertices_output.str().c_str());
  fout_final_vertices.close();
 
  std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
 
  if (boundary_that_we_need == 0)
  {

Figure out if the vertex is on the boundary of the domain

  triangulation.begin_active(); cell != triangulation.end(); ++cell)
  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
  {
  unsigned char boundary_ids = cell->face(f)->boundary_id();
  if (boundary_ids == boundary_that_we_need)
  {
  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
  if (vertex_touched[cell->face(f)->vertex_index(v)] == false)
  {
  vertex_touched[cell->face(f)->vertex_index(v)] = true;
  std::ofstream fout_final_vertices(vertices_output.str().c_str(), std::ios::app);
  fout_final_vertices << cell->face(f)->vertex(v) << "\n";
  fout_final_vertices.close();
  }
  }
  }
  }
  else
  {

Figure out if the vertex is on an internal boundary

  triangulation.begin_active(); cell != triangulation.end(); ++cell)
  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
  {
  if (cell->neighbor(f) != triangulation.end())
  {
  if (cell->material_id() != cell->neighbor(f)->material_id()) //finds face is at internal boundary
  {
  int high_mat_id = std::max(cell->material_id(),
  cell->neighbor(f)->material_id());
  if (high_mat_id == boundary_that_we_need) //finds faces at the correct internal boundary
  {
  for (unsigned int v = 0;
  v < GeometryInfo<dim>::vertices_per_face;
  ++v)
  if (vertex_touched[cell->face(f)->vertex_index(
  v)] == false)
  {
  vertex_touched[cell->face(f)->vertex_index(
  v)] = true;
  std::ofstream fout_final_vertices(vertices_output.str().c_str(), std::ios::app);
  fout_final_vertices << cell->face(f)->vertex(v) << "\n";
  fout_final_vertices.close();
  }
  }
  }
  }
  }
  }
  }
 
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)

====================== SETUP INITIAL MESH ======================

  template<int dim>
  void StokesProblem<dim>::setup_initial_mesh()
  {
  GridIn<dim> grid_in;
  std::ifstream mesh_stream(system_parameters::mesh_filename,
  std::ifstream::in);
  grid_in.read_ucd(mesh_stream);
 
void attach_triangulation(Triangulation< dim, spacedim > &tria)
Definition grid_in.cc:154

output initial mesh in eps

  std::ostringstream initial_mesh_file;
  initial_mesh_file << system_parameters::output_folder << "/initial_mesh.eps";
  std::ofstream out_eps (initial_mesh_file.str().c_str());
  GridOut grid_out;
  grid_out.write_eps (triangulation, out_eps);
  out_eps.close();
 
void write_eps(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
Definition grid_out.cc:5109

set boundary ids boundary indicator 0 is outer free surface; 1, 2, 3 ... is boundary between layers, 99 is flat boundaries

  cell=triangulation.begin_active(), endc=triangulation.end();
 
  unsigned int how_many; // how many components away from cardinal planes
 
  std::ostringstream boundaries_file;
  boundaries_file << system_parameters::output_folder << "/boundaries.txt";
  std::ofstream fout_boundaries(boundaries_file.str().c_str());
  fout_boundaries.close();
 
  double zero_tolerance = 1e-3;
  for (; cell != endc; ++cell) // loop over all cells
  {
  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f) // loop over all vertices
  {
  if (cell->face(f)->at_boundary())
  {

print boundary

  std::ofstream fout_boundaries(boundaries_file.str().c_str(), std::ios::app);
  fout_boundaries << cell->face(f)->center()[0] << " " << cell->face(f)->center()[1]<< "\n";
  fout_boundaries.close();
 
  how_many = 0;
  for (unsigned int i=0; i<dim; ++i)
  if (fabs(cell->face(f)->center()[i]) > zero_tolerance)
  how_many++;
  if (how_many==dim)
  cell->face(f)->set_all_boundary_ids(0); // if face center coordinates > zero_tol, set bnry indicators to 0
  else
  cell->face(f)->set_all_boundary_ids(99);
  }
  }
  }
 
  std::ostringstream ellipses_filename;
  ellipses_filename << system_parameters::output_folder << "/ellipse_fits.txt";
  std::ofstream fout_ellipses(ellipses_filename.str().c_str());
  fout_ellipses.close();
 

Find ellipsoidal axes for all layers

  std::vector<double> ellipse_axes(0);

compute fit to boundary 0, 1, 2 ...

  std::cout << endl;
  for (unsigned int i = 0; i<system_parameters::sizeof_material_id; ++i)
  {
  ellipsoid.compute_fit(ellipse_axes, system_parameters::material_id[i]);
  system_parameters::q_axes.push_back(ellipse_axes[0]);
  system_parameters::p_axes.push_back(ellipse_axes[1]);
 
  std::cout << "a_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[0]
  << " " << " c_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[1] << std::endl;
  ellipse_axes.clear();
 
  std::ofstream fout_ellipses(ellipses_filename.str().c_str(), std::ios::app);
  fout_ellipses << system_parameters::present_timestep << " a_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[0]
  << " " << " c_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[1] << endl;
  fout_ellipses.close();
  }
 
  triangulation.refine_global(system_parameters::global_refinement);
 
 

refines crustal region

  if (system_parameters::crustal_refinement != 0)
  {
  double a = system_parameters::q_axes[0] - system_parameters::crust_refine_region;
  double b = system_parameters::p_axes[0] - system_parameters::crust_refine_region;
 
 
  for (unsigned int step = 0;
  step < system_parameters::crustal_refinement; ++step)
  {
  typename ::Triangulation<dim>::active_cell_iterator cell =
  triangulation.begin_active(), endc = triangulation.end();
  for (; cell != endc; ++cell)
  for (unsigned int v = 0;
  v < GeometryInfo<dim>::vertices_per_cell; ++v)
  {
  Point<dim> current_vertex = cell->vertex(v);
 
  const double x_coord = current_vertex.operator()(0);
  const double y_coord = current_vertex.operator()(1);
  double expected_z = -1;
 
  if ((x_coord - a) < -1e-10)
  expected_z = b
  * std::sqrt(1 - (x_coord * x_coord / a / a));
 
  if (y_coord >= expected_z)
  {
  cell->set_refine_flag();
  break;
  }
  }
  triangulation.execute_coarsening_and_refinement();
  }
  }
 
 

output initial mesh in eps

  std::ostringstream refined_mesh_file;
  refined_mesh_file << system_parameters::output_folder << "/refined_mesh.eps";
  std::ofstream out_eps_refined (refined_mesh_file.str().c_str());
  GridOut grid_out_refined;
  grid_out_refined.write_eps (triangulation, out_eps_refined);
  out_eps_refined.close();
  write_vertices(0);
  write_vertices(1);
  write_mesh();
  }
 

====================== REFINE MESH ======================

  template<int dim>
  void StokesProblem<dim>::refine_mesh()
  {
  using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
  Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
 
  std::vector<bool> component_mask(dim + 1, false);
  component_mask[dim] = true;
  FunctionMap(), solution,
  estimated_error_per_cell, component_mask);
 
  estimated_error_per_cell, 0.3, 0.0);
  triangulation.execute_coarsening_and_refinement();
  }
 
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())

====================== SET UP THE DATA STRUCTURES TO REMEMBER STRESS FIELD ======================

  template<int dim>
  void StokesProblem<dim>::setup_quadrature_point_history()
  {
  unsigned int our_cells = 0;
  triangulation.begin_active(); cell != triangulation.end(); ++cell)
  ++our_cells;
 
  triangulation.clear_user_data();
 
  quadrature_point_history.resize(our_cells * quadrature_formula.size());
 
  unsigned int history_index = 0;
  triangulation.begin_active(); cell != triangulation.end(); ++cell)
  {
  cell->set_user_pointer(&quadrature_point_history[history_index]);
  history_index += quadrature_formula.size();
  }
 
  Assert(history_index == quadrature_point_history.size(), ExcInternalError());
  }
 

====================== DOES ELASTIC STEPS ======================

  template<int dim>
  void StokesProblem<dim>::do_elastic_steps()
  {
  unsigned int elastic_iteration = 0;
 
  while (elastic_iteration < system_parameters::initial_elastic_iterations)
  {
 
  std::cout << "\n\nElastic iteration " << elastic_iteration
  << "\n";
  setup_dofs();
 
  if (system_parameters::present_timestep == 0)
  initialize_eta_and_G();
 
  if (elastic_iteration == 0)
  system_parameters::current_time_interval =
  system_parameters::viscous_time; //This is the time interval needed in assembling the problem
 
  std::cout << " Assembling..." << std::endl << std::flush;
  assemble_system();
 
  std::cout << " Solving..." << std::flush;
  solve();
 
  output_results();
  update_quadrature_point_history();
 
  append_physical_times(0);
  elastic_iteration++;
  system_parameters::present_timestep++;
  do_ellipse_fits();
  write_vertices(0);
  write_vertices(1);
  write_mesh();
  update_time_interval();
  }
  }
 

====================== DO A SINGLE VISCOELASTOPLASTIC TIMESTEP ======================

  template<int dim>
  void StokesProblem<dim>::do_flow_step()
  {
  plastic_iteration = 0;
  while (plastic_iteration < system_parameters::max_plastic_iterations)
  {
  if (system_parameters::continue_plastic_iterations == true)
  {
  std::cout << "Plasticity iteration " << plastic_iteration << "\n";
  setup_dofs();
 
  std::cout << " Assembling..." << std::endl << std::flush;
  assemble_system();
 
  std::cout << " Solving..." << std::flush;
  solve();
 
  output_results();
  solution_stesses();
 
  if (system_parameters::continue_plastic_iterations == false)
  break;
 
  plastic_iteration++;
  }
  }
  }
 

====================== RUN ======================

  template<int dim>
  void StokesProblem<dim>::run()
  {

Sets up mesh and data structure for viscosity and stress at quadrature points

  setup_initial_mesh();
  setup_quadrature_point_history();
 

Makes the physical_times.txt file

  std::ostringstream times_filename;
  times_filename << system_parameters::output_folder << "/physical_times.txt";
  std::ofstream fout_times(times_filename.str().c_str());
  fout_times.close();
 

Computes elastic timesteps

  do_elastic_steps();

Computes viscous timesteps

  unsigned int VEPstep = 0;
  while (system_parameters::present_timestep
  < (system_parameters::initial_elastic_iterations
  + system_parameters::total_viscous_steps))
  {
 
  if (system_parameters::continue_plastic_iterations == false)
  system_parameters::continue_plastic_iterations = true;
  std::cout << "\n\nViscoelastoplastic iteration " << VEPstep << "\n\n";

Computes plasticity

  do_flow_step();
  update_quadrature_point_history();
  move_mesh();
  append_physical_times(plastic_iteration);
  system_parameters::present_timestep++;
  system_parameters::present_time = system_parameters::present_time + system_parameters::current_time_interval;
  do_ellipse_fits();
  write_vertices(0);
  write_vertices(1);
  write_mesh();
  VEPstep++;
  }
  append_physical_times(-1);
  }
 
  }
 

====================== MAIN ======================

  int main(int argc, char *argv[])
  {
 

output program name

  std::cout << "Running: " << argv[0] << std::endl;
 
  char *cfg_filename = new char[120];
 
  if (argc == 1) // if no input parameters (as if launched from eclipse)
  {
  std::strcpy(cfg_filename,"config/sample_CeresFE_config.cfg");
  }
  else
  std::strcpy(cfg_filename,argv[1]);
 
  try
  {
  using namespace dealii;
  using namespace Step22;
  config_in cfg(cfg_filename);
 
  std::clock_t t1;
  std::clock_t t2;
  t1 = std::clock();
 
 
  StokesProblem<2> flow_problem(1);
  flow_problem.run();
 
  std::cout << std::endl << "\a";
 
  t2 = std::clock();
  float diff (((float)t2 - (float)t1) / (float)CLOCKS_PER_SEC);
  std::cout << "\n Program run in: " << diff << " seconds" << endl;
  }
  catch (std::exception &exc)
  {
  std::cerr << std::endl << std::endl
  << "----------------------------------------------------"
  << std::endl;
  std::cerr << "Exception on processing: " << std::endl << exc.what()
  << std::endl << "Aborting!" << std::endl
  << "----------------------------------------------------"
  << std::endl;
 
  return 1;
  }
  catch (...)
  {
  std::cerr << std::endl << std::endl
  << "----------------------------------------------------"
  << std::endl;
  std::cerr << "Unknown exception!" << std::endl << "Aborting!"
  << std::endl
  << "----------------------------------------------------"
  << std::endl;
  return 1;
  }
 
  return 0;
  }
unsigned int depth_console(const unsigned int n)
Definition logstream.cc:349
LogStream deallog
Definition logstream.cc:36

Annotated version of support_code/config_in.h

  /*
  * config_in.h
  *
  * Created on: Aug 17, 2015
  * Author: antonermakov
  */
 
 
  #include <iostream>
  #include <iomanip>
  #include <cstdlib>
  #include <string>
  #include <libconfig.h++>
 
  #include "local_math.h"
 
  using namespace std;
  using namespace libconfig;
 
  namespace Step22
  {
  namespace system_parameters
  {
 

Mesh file name

  string mesh_filename;
  string output_folder;
 

Body parameters

  double r_mean;
  double period;
  double omegasquared;
  double beta;
  double intercept;
 

Rheology parameters

  vector<double> depths_eta;
  vector<double> eta_kinks;
  vector<double> depths_rho;
  vector<double> rho;
  vector<int> material_id;
  vector<double> G;
 
  double eta_ceiling;
  double eta_floor;
  double eta_Ea;
  bool lat_dependence;
 
  unsigned int sizeof_depths_eta;
  unsigned int sizeof_depths_rho;
  unsigned int sizeof_rho;
  unsigned int sizeof_eta_kinks;
  unsigned int sizeof_material_id;
  unsigned int sizeof_G;
 
  double pressure_scale;
  double q;
  bool cylindrical;
  bool continue_plastic_iterations;
 

plasticity variables

  bool plasticity_on;
  unsigned int failure_criterion;
  unsigned int max_plastic_iterations;
  double smoothing_radius;
 

viscoelasticity variables

  unsigned int initial_elastic_iterations;
  double elastic_time;
  double viscous_time;
  double initial_disp_target;
  double final_disp_target;
  double current_time_interval;
 

mesh refinement variables

  unsigned int global_refinement;
  unsigned int small_r_refinement;
  unsigned int crustal_refinement;
  double crust_refine_region;
  unsigned int surface_refinement;
 

solver variables

  int iteration_coefficient;
  double tolerance_coefficient;
 

time step variables

  double present_time;
  unsigned int present_timestep;
  unsigned int total_viscous_steps;
 
 

ellipse axes

  vector<double> q_axes;
  vector<double> p_axes;
 
  }
 
  class config_in
  {
  public:
  config_in(char *);
 
  private:
  void write_config();
  };
 
  void config_in::write_config()
  {
  std::ostringstream config_parameters;
  config_parameters << system_parameters::output_folder << "/run_parameters.txt";
  std::ofstream fout_config(config_parameters.str().c_str());
 

mesh filename

  fout_config << "mesh filename: " << system_parameters::mesh_filename << endl << endl;
 

body parameters

  fout_config << "r_mean = " << system_parameters::r_mean << endl;
  fout_config << "period = " << system_parameters::period << endl;
  fout_config << "omegasquared = " << system_parameters::omegasquared << endl;
  fout_config << "beta = " << system_parameters::beta << endl;
  fout_config << "intercept = " << system_parameters::intercept << endl;
 

rheology parameters

  for (unsigned int i=0; i<system_parameters::sizeof_depths_eta; ++i)
  fout_config << "depths_eta[" << i << "] = " << system_parameters::depths_eta[i] << endl;
 
  for (unsigned int i=0; i<system_parameters::sizeof_eta_kinks; ++i)
  fout_config << "eta_kinks[" << i << "] = " << system_parameters::eta_kinks[i] << endl;
 
  for (unsigned int i=0; i<system_parameters::sizeof_depths_rho; ++i)
  fout_config << "depths_rho[" << i << "] = " << system_parameters::depths_rho[i] << endl;
 
  for (unsigned int i=0; i<system_parameters::sizeof_rho; ++i)
  fout_config << "rho[" << i << "] = " << system_parameters::rho[i] << endl;
 
  for (unsigned int i=0; i<system_parameters::sizeof_material_id; ++i)
  fout_config << "material_id[" << i << "] = " << system_parameters::material_id[i] << endl;
 
  for (unsigned int i=0; i<system_parameters::sizeof_G; ++i)
  fout_config << "G[" << i << "] = " << system_parameters::G[i] << endl;
 
  fout_config << "eta_ceiling = " << system_parameters::eta_ceiling << endl;
  fout_config << "eta_floor = " << system_parameters::eta_floor << endl;
  fout_config << "eta_Ea = " << system_parameters::eta_Ea << endl;
  fout_config << "lat_dependence = " << system_parameters::lat_dependence << endl;
  fout_config << "pressure_scale = " << system_parameters::pressure_scale << endl;
  fout_config << "q = " << system_parameters::q << endl;
  fout_config << "cylindrical = " << system_parameters::cylindrical << endl;
  fout_config << "continue_plastic_iterations = " << system_parameters::continue_plastic_iterations << endl;
 

Plasticity parameters

  fout_config << "plasticity_on = " << system_parameters::plasticity_on << endl;
  fout_config << "failure_criterion = " << system_parameters::failure_criterion << endl;
  fout_config << "max_plastic_iterations = " << system_parameters::max_plastic_iterations << endl;
  fout_config << "smoothing_radius = " << system_parameters::smoothing_radius << endl;
 

Viscoelasticity parameters

  fout_config << "initial_elastic_iterations = " << system_parameters::initial_elastic_iterations << endl;
  fout_config << "elastic_time = " << system_parameters::elastic_time << endl;
  fout_config << "viscous_time = " << system_parameters::viscous_time << endl;
  fout_config << "initial_disp_target = " << system_parameters::initial_disp_target << endl;
  fout_config << "final_disp_target = " << system_parameters::final_disp_target << endl;
  fout_config << "current_time_interval = " << system_parameters::current_time_interval << endl;
 

Mesh refinement parameters

  fout_config << "global_refinement = " << system_parameters::global_refinement << endl;
  fout_config << "small_r_refinement = " << system_parameters::small_r_refinement << endl;
  fout_config << "crustal_refinement = " << system_parameters::crustal_refinement << endl;
  fout_config << "crust_refine_region = " << system_parameters::crust_refine_region << endl;
  fout_config << "surface_refinement = " << system_parameters::surface_refinement << endl;
 

Solver parameters

  fout_config << "iteration_coefficient = " << system_parameters::iteration_coefficient << endl;
  fout_config << "tolerance_coefficient = " << system_parameters::tolerance_coefficient << endl;
 

Time step parameters

  fout_config << "present_time = " << system_parameters::present_time << endl;
  fout_config << "present_timestep = " << system_parameters::present_timestep << endl;
  fout_config << "total_viscous_steps = " << system_parameters::total_viscous_steps << endl;
 
  fout_config.close();
  }
 
  config_in::config_in(char *filename)
  {
 

This example reads the configuration file 'example.cfg' and displays some of its contents.

  Config cfg;
 

Read the file. If there is an error, report it and exit.

  try
  {
  cfg.readFile(filename);
  }
  catch (const FileIOException &fioex)
  {
  std::cerr << "I/O error while reading file:" << filename << std::endl;
  }
  catch (const ParseException &pex)
  {
  std::cerr << "Parse error at " << pex.getFile() << ":" << pex.getLine()
  << " - " << pex.getError() << std::endl;
  }
 
std::vector< unsigned int > pex(const std::vector< unsigned int > &targets, const std::function< RequestType(const unsigned int)> &create_request, const std::function< AnswerType(const unsigned int, const RequestType &)> &answer_request, const std::function< void(const unsigned int, const AnswerType &)> &process_answer, const MPI_Comm comm)

Get mesh name.

  try
  {
  string msh = cfg.lookup("mesh_filename");
  system_parameters::mesh_filename = msh;
  }
  catch (const SettingNotFoundException &nfex)
  {
  cerr << "No 'mesh_filename' setting in configuration file." << endl;
  }
 

get output folder

  try
  {
  string output = cfg.lookup("output_folder");
  system_parameters::output_folder = output;
 
  std::cout << "Writing to folder: " << output << endl;
  }
  catch (const SettingNotFoundException &nfex)
  {
  cerr << "No 'output_folder' setting in configuration file." << endl;
  }
 

get radii

  const Setting &root = cfg.getRoot();
 
 

get body parameters

  try
  {
  const Setting &body_parameters = root["body_parameters"];
 
  body_parameters.lookupValue("period", system_parameters::period);
  system_parameters::omegasquared = pow(TWOPI / 3600.0 / system_parameters::period, 2.0);
  body_parameters.lookupValue("r_mean", system_parameters::r_mean);
  body_parameters.lookupValue("beta", system_parameters::beta);
  body_parameters.lookupValue("intercept", system_parameters::intercept);
  }
  catch (const SettingNotFoundException &nfex)
  {
  cerr << "We've got a problem in the body parameters block" << endl;
 
  }
 

Rheology parameters

  try
  {

get depths_eta ------------------—

  const Setting &set_depths_eta = cfg.lookup("rheology_parameters.depths_eta");
 
  unsigned int ndepths_eta = set_depths_eta.getLength();
  system_parameters::sizeof_depths_eta = ndepths_eta;
 
  for (unsigned int i=0; i<ndepths_eta; ++i)
  {
  system_parameters::depths_eta.push_back(set_depths_eta[i]);
  cout << "depth_eta[" << i << "] = " << system_parameters::depths_eta[i] << endl;
  }
 

get eta_kinks ----------------------—

  const Setting &set_eta_kinks = cfg.lookup("rheology_parameters.eta_kinks");
 
  unsigned int neta_kinks = set_eta_kinks.getLength();
  system_parameters::sizeof_eta_kinks = neta_kinks;
 

cout << "Number of depth = " << ndepths << endl;

  for (unsigned int i=0; i<neta_kinks; ++i)
  {
  system_parameters::eta_kinks.push_back(set_eta_kinks[i]);
  cout << "eta_kinks[" << i << "] = " << system_parameters::eta_kinks[i] << endl;
  }
 

get depths_rho ----------------------—

  const Setting &set_depths_rho = cfg.lookup("rheology_parameters.depths_rho");
 
  unsigned int ndepths_rho = set_depths_rho.getLength();
  system_parameters::sizeof_depths_rho = ndepths_rho;
 

cout << "Number of depth = " << ndepths << endl;

  for (unsigned int i=0; i<ndepths_rho; ++i)
  {
  system_parameters::depths_rho.push_back(set_depths_rho[i]);
  cout << "depths_rho[" << i << "] = " << system_parameters::depths_rho[i] << endl;
  }
 

get rho ----------------------—

  const Setting &set_rho = cfg.lookup("rheology_parameters.rho");
 
  unsigned int nrho = set_rho.getLength();
  system_parameters::sizeof_rho = nrho;
 

cout << "Number of depth = " << ndepths << endl;

  for (unsigned int i=0; i<nrho; ++i)
  {
  system_parameters::rho.push_back(set_rho[i]);
  cout << "rho[" << i << "] = " << system_parameters::rho[i] << endl;
  }
 

get material_id ----------------------—

  const Setting &set_material_id = cfg.lookup("rheology_parameters.material_id");
 
  unsigned int nmaterial_id = set_material_id.getLength();
  system_parameters::sizeof_material_id = nmaterial_id;
 

cout << "Number of depth = " << ndepths << endl;

  for (unsigned int i=0; i<nmaterial_id; ++i)
  {
  system_parameters::material_id.push_back(set_material_id[i]);
  cout << "material_id[" << i << "] = " << system_parameters::material_id[i] << endl;
  }
 

get G ----------------------—

  const Setting &set_G = cfg.lookup("rheology_parameters.G");
 
  unsigned int nG = set_G.getLength();
  system_parameters::sizeof_G = nG;
 

cout << "Number of depth = " << ndepths << endl;

  for (unsigned int i=0; i<nG; ++i)
  {
  system_parameters::G.push_back(set_G[i]);
  cout << "G[" << i << "] = " << system_parameters::G[i] << endl;
  }
 
  const Setting &rheology_parameters = root["rheology_parameters"];
  rheology_parameters.lookupValue("eta_ceiling", system_parameters::eta_ceiling);
  rheology_parameters.lookupValue("eta_floor", system_parameters::eta_floor);
  rheology_parameters.lookupValue("eta_Ea", system_parameters::eta_Ea);
  rheology_parameters.lookupValue("lat_dependence", system_parameters::lat_dependence);
  rheology_parameters.lookupValue("pressure_scale", system_parameters::pressure_scale);
  rheology_parameters.lookupValue("q", system_parameters::q);
  rheology_parameters.lookupValue("cylindrical", system_parameters::cylindrical);
  rheology_parameters.lookupValue("continue_plastic_iterations", system_parameters::continue_plastic_iterations);
  }
  catch (const SettingNotFoundException &nfex)
  {
  cerr << "We've got a problem in the rheology parameters block" << endl;
  }
 

Plasticity parameters

  try
  {
 
  const Setting &plasticity_parameters = root["plasticity_parameters"];
  plasticity_parameters.lookupValue("plasticity_on", system_parameters::plasticity_on);
  plasticity_parameters.lookupValue("failure_criterion", system_parameters::failure_criterion);
  plasticity_parameters.lookupValue("max_plastic_iterations", system_parameters::max_plastic_iterations);
  plasticity_parameters.lookupValue("smoothing_radius", system_parameters::smoothing_radius);
  }
  catch (const SettingNotFoundException &nfex)
  {
  cerr << "We've got a problem in the plasticity parameters block" << endl;
  }
 

Viscoelasticity parameters

  try
  {
 
  const Setting &viscoelasticity_parameters = root["viscoelasticity_parameters"];
  viscoelasticity_parameters.lookupValue("initial_elastic_iterations", system_parameters::initial_elastic_iterations);
  viscoelasticity_parameters.lookupValue("elastic_time", system_parameters::elastic_time);
  viscoelasticity_parameters.lookupValue("viscous_time", system_parameters::viscous_time);
  viscoelasticity_parameters.lookupValue("initial_disp_target", system_parameters::initial_disp_target);
  viscoelasticity_parameters.lookupValue("final_disp_target", system_parameters::final_disp_target);
  viscoelasticity_parameters.lookupValue("current_time_interval", system_parameters::current_time_interval);
 
  system_parameters::viscous_time *= SECSINYEAR;
  }
  catch (const SettingNotFoundException &nfex)
  {
  cerr << "We've got a problem in the viscoelasticity parameters block" << endl;
  }
 

Mesh refinement parameters

  try
  {
 
  const Setting &mesh_refinement_parameters = root["mesh_refinement_parameters"];
  mesh_refinement_parameters.lookupValue("global_refinement", system_parameters::global_refinement);
  mesh_refinement_parameters.lookupValue("small_r_refinement", system_parameters::small_r_refinement);
  mesh_refinement_parameters.lookupValue("crustal_refinement", system_parameters::crustal_refinement);
  mesh_refinement_parameters.lookupValue("crust_refine_region", system_parameters::crust_refine_region);
  mesh_refinement_parameters.lookupValue("surface_refinement", system_parameters::surface_refinement);
  }
  catch (const SettingNotFoundException &nfex)
  {
  cerr << "We've got a problem in the mesh refinement parameters block" << endl;
  }
 

Solver parameters

  try
  {
  const Setting &solve_parameters = root["solve_parameters"];
  solve_parameters.lookupValue("iteration_coefficient", system_parameters::iteration_coefficient);
  solve_parameters.lookupValue("tolerance_coefficient", system_parameters::tolerance_coefficient);
 
 
  }
  catch (const SettingNotFoundException &nfex)
  {
  cerr << "We've got a problem in the solver parameters block" << endl;
  }
 

Time step parameters

  try
  {
  const Setting &time_step_parameters = root["time_step_parameters"];
  time_step_parameters.lookupValue("present_time", system_parameters::present_time);
  time_step_parameters.lookupValue("present_timestep", system_parameters::present_timestep);
  time_step_parameters.lookupValue("total_viscous_steps", system_parameters::total_viscous_steps);
  }
  catch (const SettingNotFoundException &nfex)
  {
  cerr << "We've got a problem in the time step parameters block" << endl;
  }
 
  write_config();
  }
  }
 
 
 
 
 

Annotated version of support_code/ellipsoid_fit.h

  /*
  * ellipsoid_fit.h
  *
  * Created on: Jul 24, 2015
  * Author: antonermakov
  */
 
  #include <deal.II/base/quadrature_lib.h>
  #include <deal.II/base/function.h>
  #include <deal.II/base/logstream.h>
  #include <deal.II/lac/vector.h>
  #include <deal.II/lac/full_matrix.h>
  #include <deal.II/lac/sparse_matrix.h>
  #include <deal.II/lac/solver_cg.h>
  #include <deal.II/lac/precondition.h>
  #include <deal.II/grid/tria.h>
  #include <deal.II/dofs/dof_handler.h>
  #include <deal.II/grid/tria_accessor.h>
  #include <deal.II/grid/tria_iterator.h>
  #include <deal.II/dofs/dof_accessor.h>
  #include <deal.II/dofs/dof_tools.h>
  #include <deal.II/fe/fe_q.h>
  #include <deal.II/fe/fe_values.h>
  #include <deal.II/numerics/vector_tools.h>
  #include <deal.II/numerics/matrix_tools.h>
  #include <deal.II/numerics/data_out.h>
  #include <deal.II/grid/grid_in.h>
  #include <deal.II/grid/grid_out.h>
  #include <deal.II/base/point.h>
  #include <deal.II/grid/grid_generator.h>
 
  #include <fstream>
  #include <sstream>
  #include <iostream>
  #include <iomanip>
  #include <cstdlib>
 
 
  #include "local_math.h"
 
  using namespace dealii;
 
  template <int dim>
  class ellipsoid_fit
  {
  public:
  inline ellipsoid_fit (Triangulation<dim,dim> *pi)
  {
  p_triangulation = pi;
  };
  void compute_fit(std::vector<double> &ell, unsigned char bndry);
 
 
  private:
  Triangulation<dim,dim> *p_triangulation;
 
  };
 
 

This function computes ellipsoid fit to a set of vertices that lie on the boundary_that_we_need

  template <int dim>
  void ellipsoid_fit<dim>::compute_fit(std::vector<double> &ell, unsigned char boundary_that_we_need)
  {
  typename Triangulation<dim>::active_cell_iterator cell = p_triangulation->begin_active();
  typename Triangulation<dim>::active_cell_iterator endc = p_triangulation->end();
 
  FullMatrix<double> A(p_triangulation->n_vertices(),dim);
  Vector<double> x(dim);
  Vector<double> b(p_triangulation->n_vertices());
 
  std::vector<bool> vertex_touched (p_triangulation->n_vertices(),
  false);
 
  unsigned int j = 0;
  unsigned char boundary_ids;
  std::vector<unsigned int> ind_bnry_row;
  std::vector<unsigned int> ind_bnry_col;
 

assemble the sensitivity matrix and r.h.s.

  for (; cell != endc; ++cell)
  {
  if (boundary_that_we_need != 0)
  cell->set_manifold_id(cell->material_id());
  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
  {
  if (boundary_that_we_need == 0) //if this is the outer surface, then look for boundary ID 0; otherwise look for material ID change.
  {
  boundary_ids = cell->face(f)->boundary_id();
  if (boundary_ids == boundary_that_we_need)
  {
  for (unsigned int v = 0;
  v < GeometryInfo<dim>::vertices_per_face; ++v)
  if (vertex_touched[cell->face(f)->vertex_index(v)]
  == false)
  {
  vertex_touched[cell->face(f)->vertex_index(v)] =
  true;
  for (unsigned int i = 0; i < dim; ++i)
  {

stiffness matrix entry

  A(j, i) = pow(cell->face(f)->vertex(v)[i], 2);

r.h.s. entry

  b[j] = 1.0;

if mesh if not full: set the indicator

  }
  ind_bnry_row.push_back(j);
  j++;
  }
  }
  }
  else //find the faces that are at the boundary between materials, get the vertices, and write them into the stiffness matrix
  {
  if (cell->neighbor(f) != endc)
  {
  if (cell->material_id() != cell->neighbor(f)->material_id()) //finds face is at internal boundary
  {
  int high_mat_id = std::max(cell->material_id(),
  cell->neighbor(f)->material_id());
  if (high_mat_id == boundary_that_we_need) //finds faces at the correct internal boundary
  {
  for (unsigned int v = 0;
  v < GeometryInfo<dim>::vertices_per_face;
  ++v)
  if (vertex_touched[cell->face(f)->vertex_index(
  v)] == false)
  {
  vertex_touched[cell->face(f)->vertex_index(
  v)] = true;
  for (unsigned int i = 0; i < dim; ++i)
  {

stiffness matrix entry

  A(j, i) = pow(
  cell->face(f)->vertex(v)[i], 2);

r.h.s. entry

  b[j] = 1.0;

if mesh if not full: set the indicator

  }
  ind_bnry_row.push_back(j);
  j++;
  }
  }
  }
  }
  }
  }
  }
  if (ind_bnry_row.size()>0)
  {
 

maxtrix A'*A and vector A'*b; A'*A*x = A'*b – normal system of equations

  FullMatrix<double> AtA(dim,dim);
  Vector<double> Atb(dim);
 
  FullMatrix<double> A_out(ind_bnry_row.size(),dim);
  Vector<double> b_out(ind_bnry_row.size());
 
  for (unsigned int i=0; i<dim; ++i)
  ind_bnry_col.push_back(i);
 
  for (unsigned int i=0; i<ind_bnry_row.size(); ++i)
  b_out(i) = 1;
 
  A_out.extract_submatrix_from(A, ind_bnry_row, ind_bnry_col);
  A_out.Tmmult(AtA,A_out,true);
  A_out.Tvmult(Atb,b_out,true);
 

solve normal system of equations

  SolverControl solver_control (1000, 1e-12);
  SolverCG<> solver (solver_control);
  solver.solve (AtA, x, Atb, PreconditionIdentity());
 

find ellipsoidal axes

  for (unsigned int i=0; i<dim; ++i)
  ell.push_back(sqrt(1.0/x[i]));
  }
  else
  std::cerr << "fit_ellipsoid: no points to fit" << std::endl;
 
  }
 

Annotated version of support_code/ellipsoid_grav.h

  /* -----------------------------------------------------------------------------
  *
  * SPDX-License-Identifier: LGPL-2.1-or-later
  * Copyright (C) 2012 by Roger R. Fu
  *
  * This file is part of the deal.II code gallery.
  *
  * -----------------------------------------------------------------------------
  *
  * Author: Roger R. Fu
  * Reference Pohanka 2011, Contrib. Geophys. Geodes.
  */
 
  #include <math.h>
  #include <deal.II/base/point.h>
  #include <fstream>
  #include <iostream>
 
  namespace A_Grav_namespace
  {
  namespace system_parameters
  {
  double mantle_rho;
  double core_rho;
  double excess_rho;
  double r_eq;
  double r_polar;
 
  double r_core_eq;
  double r_core_polar;
  }
 
  template <int dim>
  class AnalyticGravity
  {
  public:
  void setup_vars (std::vector<double> v);
  void get_gravity (const ::Point<dim> &p, std::vector<double> &g);
 
  private:
  double ecc;
  double eV;
  double ke;
  double r00;
  double r01;
  double r11;
  double ecc_c;
  double eV_c;
  double ke_c;
  double r00_c;
  double r01_c;
  double r11_c;
  double g_coeff;
  double g_coeff_c;
  };
 
  template <int dim>
  void AnalyticGravity<dim>::get_gravity (const ::Point<dim> &p, std::vector<double> &g)
  {
  double rsph = std::sqrt(p[0] * p[0] + p[1] * p[1]);
  double thetasph = std::atan2(p[0], p[1]);
  double costhetasph = std::cos(thetasph);
 

convert to elliptical coordinates for silicates

  double stemp = std::sqrt((rsph * rsph - eV * eV + std::sqrt((rsph * rsph - eV * eV) * (rsph * rsph - eV * eV)
  + 4 * eV * eV * rsph * rsph * costhetasph *costhetasph)) / 2);
  double vout = stemp / system_parameters::r_eq / std::sqrt(1 - ecc * ecc);
  double eout = std::acos(rsph * costhetasph / stemp);
 

convert to elliptical coordinates for core correction

  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)
  + 4 * eV_c * eV_c * rsph * rsph * costhetasph *costhetasph)) / 2);
  double vout_c = stemp_c / system_parameters::r_core_eq / std::sqrt(1 - ecc_c * ecc_c);
  double eout_c = std::acos(rsph * costhetasph / stemp_c);
 

shell contribution

  g[0] = g_coeff * r11 * std::sqrt((1 - ecc * ecc) * vout * vout + ecc * ecc) * std::sin(eout);
  g[1] = g_coeff * r01 * vout * std::cos(eout) / std::sqrt(1 - ecc * ecc);
 

core contribution

  double expected_y = system_parameters::r_core_polar * std::sqrt(1 -
  (p[0] * p[0] / system_parameters::r_core_eq / system_parameters::r_core_eq));
 
 
  if (p[1] <= expected_y)
  {
  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);
  g[1] += g_coeff_c * r01_c * vout_c * std::cos(eout_c) / std::sqrt(1 - ecc_c * ecc_c);
  }
  else
  {
  double g_coeff_co = - 2.795007963255562e-10 * system_parameters::excess_rho * system_parameters::r_core_eq
  / vout_c / vout_c;
  double r00_co = 0;
  double r01_co = 0;
  double r11_co = 0;
 
  if (system_parameters::r_core_polar == system_parameters::r_core_eq)
  {
  r00_co = 1;
  r01_co = 1;
  r11_co = 1;
  }
  else
  {
  r00_co = ke_c * vout_c * std::atan2(1, ke_c * vout_c);
  double ke_co2 = ke_c * ke_c * vout_c * vout_c;
  r01_co = 3 * ke_co2 * (1 - r00_co);
  r11_co = 3 * ((ke_co2 + 1) * r00_co - ke_co2) / 2;
  }
  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);
  g[1] += g_coeff_co * r01_co * std::cos(eout_c) / std::sqrt(1 - ecc_c * ecc_c);
  }
  }
 
  template <int dim>
  void AnalyticGravity<dim>::setup_vars (std::vector<double> v)
  {
  system_parameters::r_eq = v[0];
  system_parameters::r_polar = v[1];
  system_parameters::r_core_eq = v[2];
  system_parameters::r_core_polar = v[3];
  system_parameters::mantle_rho = v[4];
  system_parameters::core_rho = v[5];
  system_parameters::excess_rho = system_parameters::core_rho - system_parameters::mantle_rho;
 
 

Shell

  if (system_parameters::r_polar > system_parameters::r_eq)
  {

This makes the gravity field nearly that of a sphere in case the body becomes prolate

  std::cout << "\nWarning: The model body has become prolate. \n";
  ecc = 0.001;
  }
  else
  {
  ecc = std::sqrt(1 - (system_parameters::r_polar * system_parameters::r_polar / system_parameters::r_eq / system_parameters::r_eq));
  }
 
  eV = ecc * system_parameters::r_eq;
  ke = std::sqrt(1 - (ecc * ecc)) / ecc;
  r00 = ke * std::atan2(1, ke);
  double ke2 = ke * ke;
  r01 = 3 * ke2 * (1 - r00);
  r11 = 3 * ((ke2 + 1) * r00 - ke2) / 2;
  g_coeff = - 2.795007963255562e-10 * system_parameters::mantle_rho * system_parameters::r_eq;
 

Core

  if (system_parameters::r_core_polar > system_parameters::r_core_eq)
  {
  std::cout << "\nWarning: The model core has become prolate. \n";
  ecc_c = 0.001;
  }
  else
  {
  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));
  }
  eV_c = ecc_c * system_parameters::r_core_eq;
  if (system_parameters::r_core_polar == system_parameters::r_core_eq)
  {
  ke_c = 1;
  r00_c = 1;
  r01_c = 1;
  r11_c = 1;
  g_coeff_c = - 2.795007963255562e-10 * system_parameters::excess_rho * system_parameters::r_core_eq;
  }
  else
  {
  ke_c = std::sqrt(1 - (ecc_c * ecc_c)) / ecc_c;
  r00_c = ke_c * std::atan2(1, ke_c);
  double ke2_c = ke_c * ke_c;
  r01_c = 3 * ke2_c * (1 - r00_c);
  r11_c = 3 * ((ke2_c + 1) * r00_c - ke2_c) / 2;
  g_coeff_c = - 2.795007963255562e-10 * system_parameters::excess_rho * system_parameters::r_core_eq;
  }

std::cout << "Loaded variables: ecc = " << ecc_c << " ke = " << ke_c << " r00 = " << r00_c << " r01 = " << r01_c << " r11 = " << r11_c << "\n";

  }
  }

Annotated version of support_code/local_math.h

  /*
  * File: localmath.h
  * Author: antonermakov
  *
  * Created on September 21, 2013, 7:14 PM
  */
 
  #ifndef LOCAL_MATH_
  #define LOCAL_MATH_
 
  #define PI 3.14159265358979323846
  #define TWOPI 6.283185307179586476925287
  #define SECSINYEAR 3.155692608e+07

#define ABS(a) ((a) < 0 ? -(a) : (a))

double factorial(int n) { if(n == 0) { return(1.); } else if(n == 1) { return(1.); } else if(n == 2) { return(2.); } else if(n == 3) { return(6.); } else if(n == 4) { return(24.); } else { exit(-1); } }

double fudge(int m) { if(m == 0) { return(1.0); } else { return(2.0); } }

double sign(double x) { if(x > 0) { return(1.0); } else if(x < 0.0) { return(-1.0); } else { return(0.0); } }

double pv0(double x) { double ans;

ans = x - TWOPI*floor(x/TWOPI); if(ans > TWOPI/2.0) { ans = ans - TWOPI; }

return(ans); }

  /* assumes l=2 */

double System::Plm(int m, double x) { if(m == 0) { return(1.5*x*x - 0.5); } else if(m == 1) { return(3.0*x*sqrt(1.0 - x*x)); } else if(m == 2) { return(3.0 - 3.0*x*x); } else { exit(-1); } }

double System::DP(int m, double x) { if(m == 0) { return(3.0*x); } else if(m == 1) { return((3.0 - 6.0*x*x)/sqrt(1.0 - x*x)); } else if(m == 2) { return(- 6.0*x); } else { exit(-1); } }

  #endif /* LOCALMATH_H */