147 * #include <deal.II/base/quadrature_lib.h>
148 * #include <deal.II/base/logstream.h>
149 * #include <deal.II/base/
function.h>
150 * #include <deal.II/base/utilities.h>
152 * #include <deal.II/lac/block_vector.h>
153 * #include <deal.II/lac/full_matrix.h>
154 * #include <deal.II/lac/block_sparse_matrix.h>
155 * #include <deal.II/lac/solver_cg.h>
156 * #include <deal.II/lac/precondition.h>
157 * #include <deal.II/lac/affine_constraints.h>
159 * #include <deal.II/grid/tria.h>
160 * #include <deal.II/grid/grid_generator.h>
161 * #include <deal.II/grid/tria_accessor.h>
162 * #include <deal.II/grid/tria_iterator.h>
163 * #include <deal.II/grid/grid_tools.h>
164 * #include <deal.II/grid/manifold_lib.h>
165 * #include <deal.II/grid/grid_refinement.h>
166 * #include <deal.II/grid/grid_in.h>
168 * #include <deal.II/dofs/dof_handler.h>
169 * #include <deal.II/dofs/dof_renumbering.h>
170 * #include <deal.II/dofs/dof_accessor.h>
171 * #include <deal.II/dofs/dof_tools.h>
173 * #include <deal.II/fe/fe_q.h>
174 * #include <deal.II/fe/fe_system.h>
175 * #include <deal.II/fe/fe_values.h>
177 * #include <deal.II/numerics/vector_tools.h>
178 * #include <deal.II/numerics/matrix_tools.h>
179 * #include <deal.II/numerics/data_out.h>
180 * #include <deal.II/numerics/error_estimator.h>
181 * #include <deal.II/numerics/derivative_approximation.h>
182 * #include <deal.II/numerics/fe_field_function.h>
184 * #include <deal.II/lac/sparse_direct.h>
185 * #include <deal.II/lac/sparse_ilu.h>
187 * #include <iostream>
193 * #include <armadillo>
195 * #include
"../support_code/ellipsoid_grav.h"
196 * #include
"../support_code/ellipsoid_fit.h"
197 * #include
"../support_code/config_in.h"
201 * As in all programs, the
namespace dealii
209 *
using namespace arma;
212 *
struct InnerPreconditioner;
215 *
struct InnerPreconditioner<2>
221 *
struct InnerPreconditioner<3>
245 *
const double curl = (grad_u[1][0] - grad_u[0][1]);
259 * Class
for remembering material state/properties at each quadrature
point
266 *
struct PointHistory
269 *
double old_phiphi_stress;
277 * Primary
class of this problem
284 *
class StokesProblem
287 * StokesProblem(
const unsigned int degree);
292 *
void assemble_system();
294 *
void output_results()
const;
295 *
void refine_mesh();
296 *
void solution_stesses();
297 *
void smooth_eta_field(std::vector<bool> failing_cells);
299 *
void setup_initial_mesh();
300 *
void do_elastic_steps();
301 *
void do_flow_step();
302 *
void update_time_interval();
303 *
void initialize_eta_and_G();
305 *
void do_ellipse_fits();
306 *
void append_physical_times(
int max_plastic);
307 *
void write_vertices(
unsigned char);
309 *
void setup_quadrature_point_history();
310 *
void update_quadrature_point_history();
312 *
const unsigned int degree;
318 *
unsigned int n_u = 0, n_p = 0;
319 *
unsigned int plastic_iteration = 0;
320 *
unsigned int last_max_plasticity = 0;
323 * std::vector< std::vector <Vector<double> > > quad_viscosities;
324 * std::vector<double> cell_viscosities;
325 * std::vector<PointHistory<dim> > quadrature_point_history;
335 * std::shared_ptr<typename InnerPreconditioner<dim>::type> A_preconditioner;
337 * ellipsoid_fit<dim> ellipsoid;
342 * Class
for boundary conditions and rhs
349 *
class BoundaryValuesP:
public Function<dim>
352 * BoundaryValuesP() :
358 *
const unsigned int component = 0)
const;
365 *
const unsigned int component)
const
376 *
void BoundaryValuesP<dim>::vector_value(
const Point<dim> &p,
384 *
class RightHandSide:
public Function<dim>
387 * RightHandSide () :
Function<dim>(dim+1) {}
393 *
virtual double value(
const Point<dim> &p,
const unsigned int component,
394 * A_Grav_namespace::AnalyticGravity<dim> *aGrav)
const;
397 * A_Grav_namespace::AnalyticGravity<dim> *aGrav)
const;
399 *
virtual void vector_value_list(
const std::vector<
Point<dim> > &points,
401 * A_Grav_namespace::AnalyticGravity<dim> *aGrav)
const;
407 *
const unsigned int component,
408 * A_Grav_namespace::AnalyticGravity<dim> *aGrav)
const
411 * std::vector<double> temp_vector(2);
412 * aGrav->get_gravity(p, temp_vector);
414 *
if (component == 0)
416 *
return temp_vector[0] + system_parameters::omegasquared * p[0];
420 *
if (component == 1)
421 *
return temp_vector[1];
428 *
void RightHandSide<dim>::vector_value(
const Point<dim> &p,
430 * A_Grav_namespace::AnalyticGravity<dim> *aGrav)
const
437 *
void RightHandSide<dim>::vector_value_list(
440 * A_Grav_namespace::AnalyticGravity<dim> *aGrav)
const
444 * check whether component is in
445 * the
valid range is up to the
449 *
Assert(values.size() == points.size(),
452 *
for (
unsigned int i = 0; i < points.size(); ++i)
453 * this->vector_value(points[i], values[i], aGrav);
458 * Class
for linear solvers and preconditioners
464 *
template<
class Matrix,
class Preconditioner>
468 * InverseMatrix(
const Matrix &m,
const Preconditioner &preconditioner);
477 *
template<
class Matrix,
class Preconditioner>
478 * InverseMatrix<Matrix, Preconditioner>::InverseMatrix(
const Matrix &m,
479 *
const Preconditioner &preconditioner) :
480 *
matrix(&m), preconditioner(&preconditioner)
484 *
template<
class Matrix,
class Preconditioner>
485 *
void InverseMatrix<Matrix, Preconditioner>::vmult(
Vector<double> &dst,
488 *
SolverControl solver_control(1000 * src.size(), 1
e-9 * src.l2_norm());
494 * cg.solve(*
matrix, dst, src, *preconditioner);
499 * Class
for the SchurComplement
505 *
template<
class Preconditioner>
521 *
template<
class Preconditioner>
522 * SchurComplement<Preconditioner>::SchurComplement(
525 * system_matrix(&system_matrix), A_inverse(&A_inverse), tmp1(
526 * system_matrix.block(0, 0).m()), tmp2(
527 * system_matrix.block(0, 0).m())
531 *
template<
class Preconditioner>
532 *
void SchurComplement<Preconditioner>::vmult(
Vector<double> &dst,
535 * system_matrix->
block(0, 1).
vmult(tmp1, src);
536 * A_inverse->vmult(tmp2, tmp1);
537 * system_matrix->block(1, 0).vmult(dst, tmp2);
542 * StokesProblem::StokesProblem
549 * StokesProblem<dim>::StokesProblem(
const unsigned int degree) :
553 * fe(
FE_Q<dim>(degree + 1), dim,
FE_Q<dim>(degree), 1),
555 * quadrature_formula(degree + 2),
568 *
void StokesProblem<dim>::setup_dofs()
570 * A_preconditioner.reset();
571 * system_matrix.clear();
573 * dof_handler.distribute_dofs(fe);
576 * std::vector<unsigned int> block_component(dim + 1, 0);
577 * block_component[dim] = 1;
582 * ========================================Apply Boundary Conditions=====================================
586 * constraints.clear();
587 * std::vector<bool> component_maskP(dim + 1,
false);
588 * component_maskP[dim] =
true;
591 * BoundaryValuesP<dim>(), constraints, component_maskP);
594 * std::set<types::boundary_id> no_normal_flux_boundaries;
595 * no_normal_flux_boundaries.insert(99);
597 * no_normal_flux_boundaries, constraints);
600 * constraints.close();
602 * std::vector<types::global_dof_index> dofs_per_block(2);
605 * n_u = dofs_per_block[0];
606 * n_p = dofs_per_block[1];
608 * std::cout <<
" Number of active cells: " <<
triangulation.n_active_cells()
609 * << std::endl <<
" Number of degrees of freedom: "
610 * << dof_handler.n_dofs() <<
" (" << n_u <<
'+' << n_p <<
')'
616 * csp.block(0, 0).reinit(n_u, n_u);
617 * csp.block(1, 0).reinit(n_p, n_u);
618 * csp.block(0, 1).reinit(n_u, n_p);
619 * csp.block(1, 1).reinit(n_p, n_p);
621 * csp.collect_sizes();
624 * sparsity_pattern.copy_from(csp);
627 * system_matrix.reinit(sparsity_pattern);
629 * solution.reinit(2);
630 * solution.block(0).reinit(n_u);
631 * solution.block(1).reinit(n_p);
632 * solution.collect_sizes();
634 * system_rhs.reinit(2);
635 * system_rhs.block(0).reinit(n_u);
636 * system_rhs.block(1).reinit(n_p);
637 * system_rhs.collect_sizes();
653 *
double get_eta(
double &r,
double &z);
654 *
double get_G(
unsigned int mat_id);
657 * std::vector<double> get_manual_eta_profile();
662 * std::vector<double> Rheology<dim>::get_manual_eta_profile()
664 * vector<double> etas;
666 *
for (
unsigned int i=0; i < system_parameters::sizeof_depths_eta; i++)
668 * etas.push_back(system_parameters::depths_eta[i]);
669 * etas.push_back(system_parameters::eta_kinks[i]);
675 *
double Rheology<dim>::get_eta(
double &r,
double &z)
679 * compute local
depth
682 *
double ecc = system_parameters::q_axes[0] / system_parameters::p_axes[0];
683 *
double Rminusr = system_parameters::q_axes[0] - system_parameters::p_axes[0];
684 *
double approx_a =
std::sqrt(r * r + z * z * ecc * ecc);
685 *
double group1 = r * r + z * z - Rminusr * Rminusr;
687 *
double a0 = approx_a;
688 *
double error = 10000;
691 * While
loop finds the a axis of the
"isodepth" ellipse
for which the input
point is on the surface.
692 * An
"isodepth" ellipse is defined as
one whose axes a,
b are related to the global axes
A, B by:
A-h = B-h
698 *
if ((r > system_parameters::q_axes[0] - system_parameters::depths_eta.back()) ||
699 * (z > system_parameters::p_axes[0] - system_parameters::depths_eta.back()))
703 *
while (error >=
eps)
705 *
double a02 = a0 * a0;
706 *
double a03 = a0 * a02;
707 *
double a04 = a0 * a03;
708 *
double fofa = a04 - (2 * Rminusr * a03) - (group1 * a02)
709 * + (2 * r * r * Rminusr * a0) - (r * r * Rminusr * Rminusr);
710 *
double fprimeofa = 4 * a03 - (6 * Rminusr * a02) - (2 * group1 * a0)
711 * + (2 * r * r * Rminusr);
712 *
double deltaa = -fofa / fprimeofa;
717 * cout <<
"error = " << error << endl;
727 *
double local_depth = system_parameters::q_axes[0] - a0;
728 *
if (local_depth < 0)
731 *
if (local_depth > system_parameters::depths_eta.back())
733 *
if (system_parameters::eta_kinks.back() < system_parameters::eta_floor)
734 *
return system_parameters::eta_floor;
735 *
else if (system_parameters::eta_kinks.back() > system_parameters::eta_ceiling)
736 *
return system_parameters::eta_ceiling;
738 *
return system_parameters::eta_kinks.back();
741 * std::vector<double> viscosity_function = get_manual_eta_profile();
743 *
unsigned int n_visc_kinks = viscosity_function.size() / 2;
747 * find the correct interval to
do the interpolation in
750 *
int n_minus_one = -1;
751 *
for (
unsigned int n = 1; n <= n_visc_kinks; n++)
753 *
unsigned int ndeep = 2 * n - 2;
754 *
unsigned int nshallow = 2 * n;
755 *
if (local_depth >= viscosity_function[ndeep] && local_depth <= viscosity_function[nshallow])
756 * n_minus_one = ndeep;
761 * find the viscosity interpolation
764 *
if (n_minus_one == -1)
765 *
return system_parameters::eta_ceiling;
768 *
double visc_exponent =
769 * (viscosity_function[n_minus_one]
771 * / (viscosity_function[n_minus_one]
772 * - viscosity_function[n_minus_one + 2]);
773 *
double visc_base = viscosity_function[n_minus_one + 3]
774 * / viscosity_function[n_minus_one + 1];
777 * This is the
true viscosity given the thermal profile
780 *
double true_eta = viscosity_function[n_minus_one + 1] *
std::pow(visc_base, visc_exponent);
784 * Implement latitude-dependence viscosity
787 *
if (system_parameters::lat_dependence)
793 *
double T_surf = T_eq *
std::sqrt( std::sqrt( std::cos(
PI / 180 * lat ) ) );
794 *
double taper_depth = 40000;
795 *
double surface_taper = (taper_depth - local_depth) / taper_depth;
796 *
if (surface_taper < 0)
798 *
double log_eta_contrast = surface_taper * system_parameters::eta_Ea * 52.5365 * (T_eq - T_surf) / T_eq / T_surf;
799 * true_eta *=
std::pow(10, log_eta_contrast);
802 *
if (true_eta > system_parameters::eta_ceiling)
803 *
return system_parameters::eta_ceiling;
804 *
else if (true_eta < system_parameters::eta_floor)
805 *
return system_parameters::eta_floor;
813 *
double Rheology<dim>::get_G(
unsigned int mat_id)
815 *
return system_parameters::G[mat_id];
821 * Initialize the eta and G parts of the quadrature_point_history
object
828 *
void StokesProblem<dim>::initialize_eta_and_G()
832 *
const unsigned int n_q_points = quadrature_formula.size();
833 * Rheology<dim> rheology;
836 * dof_handler.
begin_active(); cell != dof_handler.end(); ++cell)
838 * PointHistory<dim> *local_quadrature_points_history =
839 *
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
841 * local_quadrature_points_history >= &quadrature_point_history.front(),
844 * local_quadrature_points_history < &quadrature_point_history.back(),
846 * fe_values.reinit(cell);
848 *
for (
unsigned int q = 0; q < n_q_points; ++q)
851 *
double r_value = fe_values.quadrature_point(q)[0];
852 *
double z_value = fe_values.quadrature_point(q)[1];
856 * defines local viscosity
859 *
double local_viscosity = 0;
861 * local_viscosity = rheology.get_eta(r_value, z_value);
863 * local_quadrature_points_history[q].first_eta = local_viscosity;
864 * local_quadrature_points_history[q].new_eta = local_viscosity;
868 * defines local shear modulus
871 *
double local_G = 0;
873 *
unsigned int mat_id = cell->material_id();
875 * local_G = rheology.get_G(mat_id);
876 * local_quadrature_points_history[q].G = local_G;
880 * initializes the phi-phi stress
883 * local_quadrature_points_history[q].old_phiphi_stress = 0;
890 * ====================== ASSEMBLE THE SYSTEM ======================
897 *
void StokesProblem<dim>::assemble_system()
906 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
908 *
const unsigned int n_q_points = quadrature_formula.size();
913 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
917 * runs the gravity script
function
920 *
const RightHandSide<dim> right_hand_side;
922 * A_Grav_namespace::AnalyticGravity<dim> *aGrav =
923 *
new A_Grav_namespace::AnalyticGravity<dim>;
924 * std::vector<double> grav_parameters;
925 * grav_parameters.push_back(system_parameters::q_axes[system_parameters::present_timestep * 2 + 0]);
926 * grav_parameters.push_back(system_parameters::p_axes[system_parameters::present_timestep * 2 + 0]);
927 * grav_parameters.push_back(system_parameters::q_axes[system_parameters::present_timestep * 2 + 1]);
928 * grav_parameters.push_back(system_parameters::p_axes[system_parameters::present_timestep * 2 + 1]);
929 * grav_parameters.push_back(system_parameters::rho[0]);
930 * grav_parameters.push_back(system_parameters::rho[1]);
932 * std::cout <<
"Body parameters are: " ;
933 *
for (
int i=0; i<6; i++)
934 * std::cout << grav_parameters[i] <<
" ";
937 * aGrav->setup_vars(grav_parameters);
939 * std::vector<Vector<double> > rhs_values(n_q_points,
945 * std::vector<SymmetricTensor<2, dim> > phi_grads_u(dofs_per_cell);
946 * std::vector<double> div_phi_u(dofs_per_cell);
947 * std::vector<Tensor<1, dim> > phi_u(dofs_per_cell);
948 * std::vector<double> phi_p(dofs_per_cell);
951 * dof_handler.
begin_active(), first_cell = dof_handler.begin_active(),
952 * endc = dof_handler.end();
954 *
for (; cell != endc; ++cell)
956 * PointHistory<dim> *local_quadrature_points_history =
957 *
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
959 * local_quadrature_points_history >= &quadrature_point_history.front(),
962 * local_quadrature_points_history < &quadrature_point_history.back(),
965 *
double cell_area = cell->measure();
968 * append_physical_times(-1);
974 *
unsigned int m_id = cell->material_id();
978 * initializes the rhs vector to the correct g values
981 * fe_values.reinit(cell);
982 * right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
983 * rhs_values, aGrav);
985 * std::vector<Vector<double> > new_viscosities(quadrature_formula.size(),
Vector<double>(dim + 1));
992 *
bool is_singular =
false;
993 *
for (
unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
994 *
if (cell->face(f)->center()[0] == 0)
995 * is_singular =
true;
997 *
if (is_singular ==
false || system_parameters::cylindrical ==
false)
1004 * ===== outputs the local gravity
1007 * std::vector<Point<dim> > quad_points_list(n_q_points);
1008 * quad_points_list = fe_values.get_quadrature_points();
1010 *
if (plastic_iteration
1011 * == (system_parameters::max_plastic_iterations - 1))
1013 *
if (cell != first_cell)
1015 * std::ofstream fout(
"gravity_field.txt", std::ios::app);
1016 * fout << quad_points_list[0] <<
" " << rhs_values[0];
1021 * std::ofstream fout(
"gravity_field.txt");
1022 * fout << quad_points_list[0] <<
" " << rhs_values[0];
1027 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1030 * local_quadrature_points_history[q].old_stress;
1031 *
double &local_old_phiphi_stress =
1032 * local_quadrature_points_history[q].old_phiphi_stress;
1033 *
double r_value = fe_values.quadrature_point(q)[0];
1037 * get local density based on mat
id
1040 *
double local_density = system_parameters::rho[m_id];
1044 * defines local viscosities
1047 *
double local_viscosity = 0;
1048 *
if (plastic_iteration == 0)
1049 * local_viscosity = local_quadrature_points_history[q].first_eta;
1051 * local_viscosity = local_quadrature_points_history[q].new_eta;
1055 * Define the local viscoelastic constants
1058 *
double local_eta_ve = 2
1059 * / ((1 / local_viscosity)
1060 * + (1 / local_quadrature_points_history[q].G
1061 * / system_parameters::current_time_interval));
1062 *
double local_chi_ve = 1
1064 * + (local_quadrature_points_history[q].G
1065 * * system_parameters::current_time_interval
1066 * / local_viscosity));
1068 *
for (
unsigned int k = 0; k < dofs_per_cell; ++k)
1070 * phi_grads_u[k] = fe_values[velocities].symmetric_gradient(k,
1072 * div_phi_u[k] = (fe_values[velocities].divergence(k, q));
1073 * phi_u[k] = (fe_values[velocities].value(k, q));
1074 *
if (system_parameters::cylindrical ==
true)
1076 * div_phi_u[k] *= (r_value);
1077 * div_phi_u[k] += (phi_u[k][0]);
1079 * phi_p[k] = fe_values[pressure].value(k, q);
1082 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1084 *
for (
unsigned int j = 0; j <= i; ++j)
1086 *
if (system_parameters::cylindrical ==
true)
1088 * local_matrix(i, j) += (phi_grads_u[i]
1089 * * phi_grads_u[j] * 2 * local_eta_ve
1091 * + 2 * phi_u[i][0] * phi_u[j][0]
1092 * * local_eta_ve / r_value
1093 * - div_phi_u[i] * phi_p[j]
1094 * * system_parameters::pressure_scale
1095 * - phi_p[i] * div_phi_u[j]
1096 * * system_parameters::pressure_scale
1097 * + phi_p[i] * phi_p[j] * r_value
1098 * * system_parameters::pressure_scale)
1099 * * fe_values.JxW(q);
1103 * local_matrix(i, j) += (phi_grads_u[i]
1104 * * phi_grads_u[j] * 2 * local_eta_ve
1105 * - div_phi_u[i] * phi_p[j]
1106 * * system_parameters::pressure_scale
1107 * - phi_p[i] * div_phi_u[j]
1108 * * system_parameters::pressure_scale
1109 * + phi_p[i] * phi_p[j]) * fe_values.JxW(q);
1112 *
if (system_parameters::cylindrical ==
true)
1114 *
const unsigned int component_i =
1115 * fe.system_to_component_index(i).first;
1116 * local_rhs(i) += (fe_values.shape_value(i, q)
1117 * * rhs_values[q](component_i) * r_value
1119 * - local_chi_ve * phi_grads_u[i] * old_stress
1121 * - local_chi_ve * phi_u[i][0]
1122 * * local_old_phiphi_stress)
1123 * * fe_values.JxW(q);
1127 *
const unsigned int component_i =
1128 * fe.system_to_component_index(i).first;
1129 * local_rhs(i) += fe_values.shape_value(i, q)
1130 * * rhs_values[q](component_i) * fe_values.JxW(q)
1143 * ===== outputs the local gravity
1146 * std::vector<Point<dim> > quad_points_list(n_q_points);
1147 * quad_points_list = fe_values.get_quadrature_points();
1149 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1152 * local_quadrature_points_history[q].old_stress;
1153 *
double &local_old_phiphi_stress =
1154 * local_quadrature_points_history[q].old_phiphi_stress;
1155 *
double r_value = fe_values.quadrature_point(q)[0];
1157 *
double local_density = system_parameters::rho[m_id];
1161 * defines local viscosities
1164 *
double local_viscosity = 0;
1165 *
if (plastic_iteration == 0)
1167 * local_viscosity = local_quadrature_points_history[q].first_eta;
1170 * local_viscosity = local_quadrature_points_history[q].new_eta;
1174 * Define the local viscoelastic constants
1177 *
double local_eta_ve = 2
1178 * / ((1 / local_viscosity)
1179 * + (1 / local_quadrature_points_history[q].G
1180 * / system_parameters::current_time_interval));
1181 *
double local_chi_ve = 1
1183 * + (local_quadrature_points_history[q].G
1184 * * system_parameters::current_time_interval
1185 * / local_viscosity));
1187 *
for (
unsigned int k = 0; k < dofs_per_cell; ++k)
1189 * phi_grads_u[k] = fe_values[velocities].symmetric_gradient(k,
1191 * div_phi_u[k] = (fe_values[velocities].divergence(k, q));
1192 * phi_u[k] = (fe_values[velocities].value(k, q));
1193 *
if (system_parameters::cylindrical ==
true)
1195 * div_phi_u[k] *= (r_value);
1196 * div_phi_u[k] += (phi_u[k][0]);
1198 * phi_p[k] = fe_values[pressure].value(k, q);
1201 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1203 *
for (
unsigned int j = 0; j <= i; ++j)
1205 *
if (system_parameters::cylindrical ==
true)
1207 * local_matrix(i, j) += (phi_grads_u[i]
1208 * * phi_grads_u[j] * 2 * local_eta_ve
1210 * + 2 * phi_u[i][0] * phi_u[j][0]
1211 * * local_eta_ve / r_value
1212 * - div_phi_u[i] * phi_p[j]
1213 * * system_parameters::pressure_scale
1214 * - phi_p[i] * div_phi_u[j]
1215 * * system_parameters::pressure_scale
1216 * + phi_p[i] * phi_p[j] * r_value
1217 * * system_parameters::pressure_scale)
1218 * * fe_values.JxW(q);
1222 * local_matrix(i, j) += (phi_grads_u[i]
1223 * * phi_grads_u[j] * 2 * local_eta_ve
1224 * - div_phi_u[i] * phi_p[j]
1225 * * system_parameters::pressure_scale
1226 * - phi_p[i] * div_phi_u[j]
1227 * * system_parameters::pressure_scale
1228 * + phi_p[i] * phi_p[j]) * fe_values.JxW(q);
1231 *
if (system_parameters::cylindrical ==
true)
1233 *
const unsigned int component_i =
1234 * fe.system_to_component_index(i).first;
1235 * local_rhs(i) += (fe_values.shape_value(i, q)
1236 * * rhs_values[q](component_i) * r_value
1238 * - local_chi_ve * phi_grads_u[i] * old_stress
1240 * - local_chi_ve * phi_u[i][0]
1241 * * local_old_phiphi_stress)
1242 * * fe_values.JxW(q);
1246 *
const unsigned int component_i =
1247 * fe.system_to_component_index(i).first;
1248 * local_rhs(i) += fe_values.shape_value(i, q)
1249 * * rhs_values[q](component_i) * fe_values.JxW(q)
1256 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1257 *
for (
unsigned int j = i + 1; j < dofs_per_cell; ++j)
1258 * local_matrix(i, j) = local_matrix(j, i);
1260 * cell->get_dof_indices(local_dof_indices);
1261 * constraints.distribute_local_to_global(local_matrix, local_rhs,
1262 * local_dof_indices, system_matrix, system_rhs);
1265 * std::cout <<
" Computing preconditioner..." << std::endl << std::flush;
1267 * A_preconditioner = std::shared_ptr<
1268 *
typename InnerPreconditioner<dim>::type>(
1269 *
new typename InnerPreconditioner<dim>::type());
1270 * A_preconditioner->
initialize(system_matrix.block(0, 0),
1271 *
typename InnerPreconditioner<dim>::type::AdditionalData());
1278 * ====================== SOLVER ======================
1285 *
void StokesProblem<dim>::solve()
1287 *
const InverseMatrix<SparseMatrix<double>,
1288 *
typename InnerPreconditioner<dim>::type> A_inverse(
1289 * system_matrix.block(0, 0), *A_preconditioner);
1294 * A_inverse.vmult(tmp, system_rhs.block(0));
1295 * system_matrix.block(1, 0).vmult(schur_rhs, tmp);
1296 * schur_rhs -= system_rhs.block(1);
1298 * SchurComplement<typename InnerPreconditioner<dim>::type>
schur_complement(
1299 * system_matrix, A_inverse);
1301 *
int n_iterations = system_parameters::iteration_coefficient
1302 * * solution.block(1).size();
1303 *
double tolerance_goal = system_parameters::tolerance_coefficient
1304 * * schur_rhs.l2_norm();
1306 *
SolverControl solver_control(n_iterations, tolerance_goal);
1309 * std::cout <<
"\nMax iterations and tolerance are: " << n_iterations
1310 * <<
" and " << tolerance_goal << std::endl;
1313 * preconditioner.
initialize(system_matrix.block(1, 1),
1317 * system_matrix.block(1, 1), preconditioner);
1321 * constraints.distribute(solution);
1324 * std::cout <<
" " << solver_control.last_step()
1325 * <<
" outer CG Schur complement iterations for pressure"
1330 * system_matrix.block(0, 1).vmult(tmp, solution.block(1));
1332 * tmp += system_rhs.block(0);
1334 * A_inverse.vmult(solution.block(0), tmp);
1335 * constraints.distribute(solution);
1336 * solution.block(1) *= (system_parameters::pressure_scale);
1342 * ====================== OUTPUT RESULTS ======================
1346 *
void StokesProblem<dim>::output_results() const
1348 * std::vector < std::string > solution_names(dim,
"velocity");
1349 * solution_names.push_back(
"pressure");
1351 * std::vector<DataComponentInterpretation::DataComponentInterpretation> data_component_interpretation(
1353 * data_component_interpretation.push_back(
1362 * std::ostringstream filename;
1363 *
if (system_parameters::present_timestep < system_parameters::initial_elastic_iterations)
1365 * filename << system_parameters::output_folder <<
"/time"
1367 * <<
"_elastic_displacements" <<
".txt";
1371 * filename << system_parameters::output_folder <<
"/time"
1376 * std::ofstream output(filename.str().c_str());
1382 * ====================== FIND AND WRITE TO FILE THE STRESS TENSOR; IMPLEMENT PLASTICITY ======================
1389 *
void StokesProblem<dim>::solution_stesses()
1393 * note most of
this section only works with dim=2
1397 * name the output text files
1400 * std::ostringstream stress_output;
1401 * stress_output << system_parameters::output_folder <<
"/time"
1405 * std::ofstream fout_snew(stress_output.str().c_str());
1406 * fout_snew.close();
1408 * std::ostringstream stresstensor_output;
1409 * stresstensor_output << system_parameters::output_folder <<
"/time"
1413 * std::ofstream fout_sfull(stresstensor_output.str().c_str());
1414 * fout_sfull.close();
1416 * std::ostringstream failed_cells_output;
1417 * failed_cells_output << system_parameters::output_folder <<
"/time"
1421 * std::ofstream fout_failed_cells(failed_cells_output.str().c_str());
1422 * fout_failed_cells.close();
1424 * std::ostringstream plastic_eta_output;
1425 * plastic_eta_output << system_parameters::output_folder <<
"/time"
1429 * std::ofstream fout_vrnew(plastic_eta_output.str().c_str());
1430 * fout_vrnew.close();
1432 * std::ostringstream initial_eta_output;
1433 *
if (plastic_iteration == 0)
1435 * initial_eta_output << system_parameters::output_folder <<
"/time"
1437 * <<
"_baseviscosities.txt";
1438 * std::ofstream fout_baseeta(initial_eta_output.str().c_str());
1439 * fout_baseeta.close();
1442 * std::cout <<
"Running stress calculations for plasticity iteration "
1443 * << plastic_iteration <<
"...\n";
1447 * This makes the
set of points at which the stress tensor is calculated
1450 * std::vector<Point<dim> > points_list(0);
1451 * std::vector<unsigned int> material_list(0);
1453 * dof_handler.
begin_active(), endc = dof_handler.end();
1456 * This
loop gets the gradients of the velocity field and saves it in the tensor_gradient_? objects DIM
1459 *
for (; cell != endc; ++cell)
1461 * points_list.push_back(cell->center());
1462 * material_list.push_back(cell->material_id());
1466 * Make the
FEValues object to evaluate values and derivatives at quadrature points
1474 * Make the
object that will hold the velocities and velocity gradients at the quadrature points
1477 * std::vector < std::vector<Tensor<1, dim> >> velocity_grads(quadrature_formula.size(),
1479 * std::vector<Vector<double> > velocities(quadrature_formula.size(),
1483 * Make the
object to find rheology
1486 * Rheology<dim> rheology;
1490 * Write the solution flow velocity and derivative
for each cell
1493 * std::vector<Vector<double> > vector_values(0);
1494 * std::vector < std::vector<Tensor<1, dim> > > gradient_values(0);
1495 * std::vector<bool> failing_cells;
1498 * Write the stresses from the previous step into vectors
1501 * std::vector<SymmetricTensor<2, dim>> old_stress;
1502 * std::vector<double> old_phiphi_stress;
1503 * std::vector<double> cell_Gs;
1505 * dof_handler.
begin_active(); cell != dof_handler.end(); ++cell)
1509 * Makes pointer to data in quadrature_point_history
1512 * PointHistory<dim> *local_quadrature_points_history =
1513 *
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
1515 * fe_values.reinit(cell);
1516 * fe_values.get_function_gradients(solution, velocity_grads);
1517 * fe_values.get_function_values(solution, velocities);
1519 * std::vector<Tensor<1, dim>> current_cell_grads(dim+1);
1521 * current_cell_old_stress = 0;
1522 *
double current_cell_old_phiphi_stress = 0;
1523 *
double cell_area = 0;
1527 * Averages across each cell to find
mean velocities, gradients, and old stresses
1530 *
for (
unsigned int q = 0; q < quadrature_formula.size(); ++q)
1532 * cell_area += fe_values.JxW(q);
1533 * velocities[q] *= fe_values.JxW(q);
1534 * current_cell_velocity += velocities[q];
1535 *
for (
unsigned int i = 0; i < (dim+1); i++)
1537 * velocity_grads[q][i] *= fe_values.JxW(q);
1538 * current_cell_grads[i] += velocity_grads[q][i];
1540 * current_cell_old_stress += local_quadrature_points_history[q].old_stress * fe_values.JxW(q);
1541 * current_cell_old_phiphi_stress += local_quadrature_points_history[q].old_phiphi_stress * fe_values.JxW(q);
1543 * current_cell_velocity /= cell_area;
1544 *
for (
unsigned int i = 0; i < (dim+1); i++)
1545 * current_cell_grads[i] /= cell_area;
1546 * current_cell_old_stress /= cell_area;
1547 * current_cell_old_phiphi_stress /= cell_area;
1549 * vector_values.push_back(current_cell_velocity);
1550 * gradient_values.push_back(current_cell_grads);
1551 * old_stress.push_back(current_cell_old_stress);
1552 * old_phiphi_stress.push_back(current_cell_old_phiphi_stress);
1556 * Get cell shear modulus: assumes it
's constant for the cell
1559 * unsigned int mat_id = cell->material_id();
1560 * double local_G = rheology.get_G(mat_id);
1561 * cell_Gs.push_back(local_G);
1566 * tracks where failure occurred
1569 * std::vector<double> reduction_factor;
1570 * unsigned int total_fails = 0;
1571 * if (plastic_iteration == 0)
1572 * cell_viscosities.resize(0);
1575 * loop across all the cells to find and adjust eta of failing cells
1578 * for (unsigned int i = 0; i < triangulation.n_active_cells(); i++)
1580 * double current_cell_viscosity = 0;
1584 * Fill viscosities vector, analytically if plastic_iteration == 0 and from previous viscosities for later iteration
1587 * if (plastic_iteration == 0)
1589 * double local_viscosity;
1590 * local_viscosity = rheology.get_eta(points_list[i][0], points_list[i][1]);
1591 * current_cell_viscosity = local_viscosity;
1592 * cell_viscosities.push_back(current_cell_viscosity);
1596 * current_cell_viscosity = cell_viscosities[i];
1600 * double cell_eta_ve = 2
1601 * / ((1 / current_cell_viscosity)
1603 * / system_parameters::current_time_interval));
1604 * double cell_chi_ve = 1
1607 * * system_parameters::current_time_interval
1608 * / current_cell_viscosity));
1612 * find local pressure
1615 * double cell_p = vector_values[i].operator()(2);
1618 * find stresses tensor
1619 * makes non-diagonalized local matrix A
1622 * double sigma13 = 0.5
1623 * * (gradient_values[i][0][1] + gradient_values[i][1][0]);
1625 * A << gradient_values[i][0][0] << 0 << sigma13 << endr
1626 * << 0 << vector_values[i].operator()(0) / points_list[i].operator()(0)<< 0 << endr
1627 * << sigma13 << 0 << gradient_values[i][1][1] << endr;
1629 * olddevstress << old_stress[i][0][0] << 0 << old_stress[i][0][1] << endr
1630 * << 0 << old_phiphi_stress[i] << 0 << endr
1631 * << old_stress[i][0][1] << 0 << old_stress[i][1][1] << endr;
1633 * P << cell_p << cell_p << cell_p;
1634 * mat Pmat = diagmat(P);
1636 * B = (cell_eta_ve * A + cell_chi_ve * olddevstress) - Pmat;
1640 * finds principal stresses
1645 * eig_sym(eigval, eigvec, B);
1646 * double sigma1 = -min(eigval);
1647 * double sigma3 = -max(eigval);
1651 * Writes text files for principal stresses, full stress tensor, base viscosities
1654 * std::ofstream fout_snew(stress_output.str().c_str(), std::ios::app);
1655 * fout_snew << " " << sigma1 << " " << sigma3 << "\n";
1656 * fout_snew.close();
1658 * std::ofstream fout_sfull(stresstensor_output.str().c_str(), std::ios::app);
1659 * fout_sfull << A(0,0) << " " << A(1,1) << " " << A(2,2) << " " << A(0,2) << "\n";
1660 * fout_sfull.close();
1662 * if (plastic_iteration == 0)
1664 * std::ofstream fout_baseeta(initial_eta_output.str().c_str(), std::ios::app);
1665 * fout_baseeta << points_list[i]<< " " << current_cell_viscosity << "\n";
1666 * fout_baseeta.close();
1671 * Finds adjusted effective viscosity
1674 * if (system_parameters::plasticity_on)
1676 * if (system_parameters::failure_criterion == 0) //Apply Byerlee's rule
1678 *
if (sigma1 >= 5 * sigma3)
1680 * failing_cells.push_back(
true);
1681 *
double temp_reductionfactor = 1;
1683 * temp_reductionfactor = 100;
1685 * temp_reductionfactor = 1.9 * sigma1 / 5 / sigma3;
1687 * reduction_factor.push_back(temp_reductionfactor);
1692 * Text file of all failure locations
1695 * std::ofstream fout_failed_cells(failed_cells_output.str().c_str(), std::ios::app);
1696 * fout_failed_cells << points_list[i] <<
"\n";
1697 * fout_failed_cells.close();
1701 * reduction_factor.push_back(1);
1702 * failing_cells.push_back(
false);
1707 *
if (system_parameters::failure_criterion == 1)
1709 *
double temp_reductionfactor = 1;
1710 *
if (sigma3 < -114037)
1714 * std::cout <<
" ext ";
1717 * failing_cells.push_back(
true);
1718 * temp_reductionfactor = 10;
1719 * reduction_factor.push_back(temp_reductionfactor);
1724 * Text file of all failure locations
1727 * std::ofstream fout_failed_cells(failed_cells_output.str().c_str(), std::ios::app);
1728 * fout_failed_cells << points_list[i] <<
"\n";
1729 * fout_failed_cells.close();
1733 *
double sigma_c = 160e6;
1734 *
double yield_sigma1 = sigma3 +
std::sqrt( (3.086 * sigma_c * sigma3) + (0.002 * sigma3 * sigma3) );
1735 *
if (sigma1 >= yield_sigma1)
1739 * std::cout <<
" comp ";
1742 * failing_cells.push_back(
true);
1743 * temp_reductionfactor = 1.0 * sigma1 / 5 / sigma3;
1745 * reduction_factor.push_back(temp_reductionfactor);
1750 * Text file of all failure locations
1753 * std::ofstream fout_failed_cells(failed_cells_output.str().c_str(), std::ios::app);
1754 * fout_failed_cells << points_list[i] <<
"\n";
1755 * fout_failed_cells.close();
1759 * reduction_factor.push_back(1);
1760 * failing_cells.push_back(
false);
1766 * std::cout <<
"Specified failure criterion not found\n";
1771 * reduction_factor.push_back(1);
1776 * If there are enough failed cells, update eta at all quadrature points and perform smoothing
1779 * std::cout <<
" Number of failing cells: " << total_fails <<
"\n";
1780 *
double last_max_plasticity_double = last_max_plasticity;
1781 *
double total_fails_double = total_fails;
1782 *
double decrease_in_plasticity = ((last_max_plasticity_double - total_fails_double) / last_max_plasticity_double);
1783 *
if (plastic_iteration == 0)
1784 * decrease_in_plasticity = 1;
1785 * last_max_plasticity = total_fails;
1786 *
if (total_fails <= 100 || decrease_in_plasticity <= 0.2)
1788 * system_parameters::continue_plastic_iterations =
false;
1789 *
for (
unsigned int j=0; j <
triangulation.n_active_cells(); j++)
1793 * Writes to file the undisturbed cell viscosities
1796 * std::ofstream fout_vrnew(plastic_eta_output.str().c_str(), std::ios::app);
1797 * fout_vrnew <<
" " << cell_viscosities[j] <<
"\n";
1798 * fout_vrnew.close();
1806 * Decrease the eta at quadrature points in failing cells
1809 *
unsigned int cell_no = 0;
1811 * dof_handler.
begin_active(); cell != dof_handler.end(); ++cell)
1813 * fe_values.reinit(cell);
1816 * Make local_quadrature_points_history pointer to the cell data
1819 * PointHistory<dim> *local_quadrature_points_history =
1820 *
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
1822 * local_quadrature_points_history >= &quadrature_point_history.front(),
1825 * local_quadrature_points_history < &quadrature_point_history.back(),
1828 * quad_viscosities[cell_no].resize(quadrature_formula.size());
1830 *
for (
unsigned int q = 0; q < quadrature_formula.size(); ++q)
1832 *
if (plastic_iteration == 0)
1833 * local_quadrature_points_history[q].new_eta = local_quadrature_points_history[q].first_eta;
1834 * local_quadrature_points_history[q].new_eta /= reduction_factor[cell_no];
1837 * Prevents viscosities from dropping below the
floor necessary
for numerical stability
1840 *
if (local_quadrature_points_history[q].new_eta < system_parameters::eta_floor)
1841 * local_quadrature_points_history[q].new_eta = system_parameters::eta_floor;
1843 * quad_viscosities[cell_no][q].reinit(dim+1);
1844 *
for (
unsigned int ii=0; ii<dim; ii++)
1845 * quad_viscosities[cell_no][q](ii) = fe_values.quadrature_point(q)[ii];
1846 * quad_viscosities[cell_no][q](dim) = local_quadrature_points_history[q].new_eta;
1850 * smooth_eta_field(failing_cells);
1854 * Writes to file the smoothed eta field (which is defined at each quadrature
point)
for each cell
1860 * cell_viscosities.resize(
triangulation.n_active_cells(), 0);
1864 * dof_handler.
begin_active(); cell != dof_handler.end(); ++cell)
1866 *
if (failing_cells[cell_no])
1868 * fe_values.reinit(cell);
1871 * Averages across each cell to find
mean eta
1874 *
double cell_area = 0;
1875 *
for (
unsigned int q = 0; q < quadrature_formula.size(); ++q)
1877 * cell_area += fe_values.JxW(q);
1878 * cell_viscosities[cell_no] += quad_viscosities[cell_no][q][dim] * fe_values.JxW(q);
1880 * cell_viscosities[cell_no] /= cell_area;
1887 * std::ofstream fout_vrnew(plastic_eta_output.str().c_str(), std::ios::app);
1888 * fout_vrnew <<
" " << cell_viscosities[cell_no] <<
"\n";
1889 * fout_vrnew.close();
1893 * std::ofstream fout_vrnew(plastic_eta_output.str().c_str(), std::ios::app);
1894 * fout_vrnew <<
" " << cell_viscosities[cell_no] <<
"\n";
1895 * fout_vrnew.close();
1904 * ====================== SMOOTHES THE VISCOSITY FIELD AT ALL QUADRATURE POINTS ======================
1911 *
void StokesProblem<dim>::smooth_eta_field(std::vector<bool> failing_cells)
1913 * std::cout <<
" Smoothing viscosity field...\n";
1914 *
unsigned int cell_no = 0;
1916 * dof_handler.
begin_active(); cell != dof_handler.end(); ++cell)
1918 *
if (failing_cells[cell_no])
1921 * PointHistory<dim> *local_quadrature_points_history =
1922 *
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
1926 * Currently
this algorithm does not permit refinement. To permit refinement, daughter cells of neighbors must be identified
1927 * Find pointers and indices of all cells within certain radius
1930 *
bool find_more_cells =
true;
1931 * std::vector<bool> cell_touched(
triangulation.n_active_cells(),
false);
1932 * std::vector< TriaIterator< CellAccessor<dim> > > neighbor_cells;
1933 * std::vector<int> neighbor_indices;
1934 *
unsigned int start_cell = 0;
1935 *
int new_cells_found = 0;
1936 * neighbor_cells.push_back(cell);
1937 * neighbor_indices.push_back(cell_no);
1938 * cell_touched[cell_no] =
true;
1939 *
while (find_more_cells)
1941 * new_cells_found = 0;
1942 *
for (
unsigned int i = start_cell; i<neighbor_cells.size(); i++)
1944 *
for (
unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
1946 *
if (!neighbor_cells[i]->face(f)->at_boundary())
1948 *
int test_cell_no = neighbor_cells[i]->neighbor_index(f);
1949 *
if (!cell_touched[test_cell_no])
1950 *
if (cell->center().distance(neighbor_cells[i]->neighbor(f)->center()) < 2 * system_parameters::smoothing_radius)
1954 * What to
do if another nearby cell is found that hasn
't been found before
1957 * neighbor_cells.push_back(neighbor_cells[i]->neighbor(f));
1958 * neighbor_indices.push_back(test_cell_no);
1959 * cell_touched[test_cell_no] = true;
1961 * new_cells_found++;
1966 * if (new_cells_found == 0)
1968 * find_more_cells = false;
1971 * start_cell = neighbor_cells.size() - new_cells_found;
1974 * fe_values.reinit(cell);
1977 * Collect the viscosities at nearby quadrature points
1980 * for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
1982 * std::vector<double> nearby_etas_q;
1983 * for (unsigned int i = 0; i<neighbor_indices.size(); i++)
1984 * for (unsigned int j=0; j<quadrature_formula.size(); j++)
1986 * Point<dim> test_q;
1987 * for (unsigned int d=0; d<dim; d++)
1988 * test_q(d) = quad_viscosities[neighbor_indices[i]][j][d];
1989 * double qq_distance = fe_values.quadrature_point(q).distance(test_q);
1990 * if (qq_distance < system_parameters::smoothing_radius)
1991 * nearby_etas_q.push_back(quad_viscosities[neighbor_indices[i]][j][dim]);
1995 * Write smoothed viscosities to quadrature_points_history; simple boxcar function is the smoothing kernel
1998 * double mean_eta = 0;
1999 * for (unsigned int l = 0; l<nearby_etas_q.size(); l++)
2001 * mean_eta += nearby_etas_q[l];
2003 * mean_eta /= nearby_etas_q.size();
2004 * local_quadrature_points_history[q].new_eta = mean_eta;
2007 * std::cout << local_quadrature_points_history[q].new_eta << " ";
2018 * ====================== SAVE STRESS TENSOR AT QUADRATURE POINTS ======================
2025 * void StokesProblem<dim>::update_quadrature_point_history()
2027 * std::cout << " Updating stress field...";
2029 * FEValues<dim> fe_values(fe, quadrature_formula,
2030 * update_values | update_gradients | update_quadrature_points);
2034 * Make the object that will hold the velocity gradients
2037 * std::vector < std::vector<Tensor<1, dim> >> velocity_grads(quadrature_formula.size(),
2038 * std::vector < Tensor<1, dim> > (dim + 1));
2039 * std::vector<Vector<double> > velocities(quadrature_formula.size(),
2040 * Vector<double>(dim + 1));
2042 * for (typename DoFHandler<dim>::active_cell_iterator cell =
2043 * dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
2045 * PointHistory<dim> *local_quadrature_points_history =
2046 * reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
2048 * local_quadrature_points_history >= &quadrature_point_history.front(),
2049 * ExcInternalError());
2051 * local_quadrature_points_history < &quadrature_point_history.back(),
2052 * ExcInternalError());
2054 * fe_values.reinit(cell);
2055 * fe_values.get_function_gradients(solution, velocity_grads);
2056 * fe_values.get_function_values(solution, velocities);
2058 * for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
2062 * Define the local viscoelastic constants
2065 * double local_eta_ve = 2
2066 * / ((1 / local_quadrature_points_history[q].new_eta)
2067 * + (1 / local_quadrature_points_history[q].G
2068 * / system_parameters::current_time_interval));
2069 * double local_chi_ve =
2072 * + (local_quadrature_points_history[q].G
2073 * * system_parameters::current_time_interval
2074 * / local_quadrature_points_history[q].new_eta));
2078 * Compute new stress at each quadrature point
2081 * SymmetricTensor<2, dim> new_stress;
2082 * for (unsigned int i = 0; i < dim; ++i)
2083 * new_stress[i][i] =
2084 * local_eta_ve * velocity_grads[q][i][i]
2086 * * local_quadrature_points_history[q].old_stress[i][i];
2088 * for (unsigned int i = 0; i < dim; ++i)
2089 * for (unsigned int j = i + 1; j < dim; ++j)
2090 * new_stress[i][j] =
2092 * * (velocity_grads[q][i][j]
2093 * + velocity_grads[q][j][i]) / 2
2095 * * local_quadrature_points_history[q].old_stress[i][j];
2102 * AuxFunctions<dim> rotation_object;
2103 * const Tensor<2, dim> rotation = rotation_object.get_rotation_matrix(
2104 * velocity_grads[q]);
2105 * const SymmetricTensor<2, dim> rotated_new_stress = symmetrize(
2106 * transpose(rotation)
2107 * * static_cast<Tensor<2, dim> >(new_stress)
2109 * local_quadrature_points_history[q].old_stress = rotated_new_stress;
2113 * For axisymmetric case, make the phi-phi element of stress tensor
2116 * local_quadrature_points_history[q].old_phiphi_stress =
2117 * (2 * local_eta_ve * velocities[q](0)
2118 * / fe_values.quadrature_point(q)[0]
2120 * * local_quadrature_points_history[q].old_phiphi_stress);
2127 * ====================== REDEFINE THE TIME INTERVAL FOR THE VISCOUS STEPS ======================
2131 * void StokesProblem<dim>::update_time_interval()
2133 * double move_goal_per_step = system_parameters::initial_disp_target;
2134 * if (system_parameters::present_timestep > system_parameters::initial_elastic_iterations)
2136 * move_goal_per_step = system_parameters::initial_disp_target -
2137 * ((system_parameters::initial_disp_target - system_parameters::final_disp_target) /
2138 * system_parameters::total_viscous_steps *
2139 * (system_parameters::present_timestep - system_parameters::initial_elastic_iterations));
2142 * double zero_tolerance = 1e-3;
2143 * double max_velocity = 0;
2144 * for (typename DoFHandler<dim>::active_cell_iterator cell =
2145 * dof_handler.begin_active(); cell != dof_handler.end(); ++cell)// loop over all cells
2147 * if (cell->at_boundary())
2149 * int zero_faces = 0;
2150 * for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; f++)
2151 * for (unsigned int i=0; i<dim; i++)
2152 * if (fabs(cell->face(f)->center()[i]) < zero_tolerance)
2154 * if (zero_faces==0)
2156 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
2158 * Point<dim> vertex_velocity;
2159 * Point<dim> vertex_position;
2160 * for (unsigned int d = 0; d < dim; ++d)
2162 * vertex_velocity[d] = solution(cell->vertex_dof_index(v, d));
2163 * vertex_position[d] = cell->vertex(v)[d];
2167 * velocity to be evaluated is the radial component of a surface vertex
2170 * double local_velocity = 0;
2171 * for (unsigned int d = 0; d < dim; ++d)
2173 * local_velocity += vertex_velocity[d] * vertex_position [d];
2175 * local_velocity /= std::sqrt( vertex_position.square() );
2176 * if (local_velocity < 0)
2177 * local_velocity *= -1;
2178 * if (local_velocity > max_velocity)
2180 * max_velocity = local_velocity;
2188 * NOTE: It is possible for this time interval to be very different from that used in the viscoelasticity calculation.
2191 * system_parameters::current_time_interval = move_goal_per_step / max_velocity;
2192 * double step_time_yr = system_parameters::current_time_interval / SECSINYEAR;
2193 * std::cout << "Timestep interval changed to: "
2200 * ====================== MOVE MESH ======================
2207 * void StokesProblem<dim>::move_mesh()
2210 * std::cout << "\n" << " Moving mesh...\n";
2211 * std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
2212 * for (typename DoFHandler<dim>::active_cell_iterator cell =
2213 * dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
2214 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
2215 * if (vertex_touched[cell->vertex_index(v)] == false)
2217 * vertex_touched[cell->vertex_index(v)] = true;
2219 * Point<dim> vertex_displacement;
2220 * for (unsigned int d = 0; d < dim; ++d)
2221 * vertex_displacement[d] = solution(
2222 * cell->vertex_dof_index(v, d));
2223 * cell->vertex(v) += vertex_displacement
2224 * * system_parameters::current_time_interval;
2230 * ====================== WRITE MESH TO FILE ======================
2237 * void StokesProblem<dim>::write_mesh()
2241 * output mesh in ucd
2244 * std::ostringstream initial_mesh_file;
2245 * initial_mesh_file << system_parameters::output_folder << "/time" <<
2246 * Utilities::int_to_string(system_parameters::present_timestep, 2) <<
2248 * std::ofstream out_ucd (initial_mesh_file.str().c_str());
2250 * grid_out.write_ucd (triangulation, out_ucd);
2255 * ====================== FIT ELLIPSE TO SURFACE AND WRITE RADII TO FILE ======================
2262 * void StokesProblem<dim>::do_ellipse_fits()
2264 * std::ostringstream ellipses_filename;
2265 * ellipses_filename << system_parameters::output_folder << "/ellipse_fits.txt";
2268 * Find ellipsoidal axes for all layers
2271 * std::vector<double> ellipse_axes(0);
2274 * compute fit to boundary 0, 1, 2 ...
2277 * std::cout << endl;
2278 * for (unsigned int i = 0; i<system_parameters::sizeof_material_id; i++)
2280 * ellipsoid.compute_fit(ellipse_axes, system_parameters::material_id[i]);
2281 * system_parameters::q_axes.push_back(ellipse_axes[0]);
2282 * system_parameters::p_axes.push_back(ellipse_axes[1]);
2284 * std::cout << "a_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[0]
2285 * << " " << " c_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[1] << std::endl;
2286 * ellipse_axes.clear();
2288 * std::ofstream fout_ellipses(ellipses_filename.str().c_str(), std::ios::app);
2289 * fout_ellipses << system_parameters::present_timestep << " a_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[0]
2290 * << " " << " c_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[1] << endl;
2291 * fout_ellipses.close();
2297 * ====================== APPEND LINE TO PHYSICAL_TIMES.TXT FILE WITH STEP NUMBER, PHYSICAL TIME, AND # PLASTIC ITERATIONS ======================
2304 * void StokesProblem<dim>::append_physical_times(int max_plastic)
2306 * std::ostringstream times_filename;
2307 * times_filename << system_parameters::output_folder << "/physical_times.txt";
2308 * std::ofstream fout_times(times_filename.str().c_str(), std::ios::app);
2309 * fout_times << system_parameters::present_timestep << " "
2310 * << system_parameters::present_time/SECSINYEAR << " "
2311 * << max_plastic << "\n";
2314 * << system_parameters::q_axes[0] << " " << system_parameters::p_axes[0] << " "
2315 * << system_parameters::q_axes[1] << " " << system_parameters::p_axes[1] << "\n";
2318 * fout_times.close();
2323 * ====================== WRITE VERTICES TO FILE ======================
2330 * void StokesProblem<dim>::write_vertices(unsigned char boundary_that_we_need)
2332 * std::ostringstream vertices_output;
2333 * vertices_output << system_parameters::output_folder << "/time" <<
2334 * Utilities::int_to_string(system_parameters::present_timestep, 2) << "_" <<
2335 * Utilities::int_to_string(boundary_that_we_need, 2) <<
2337 * std::ofstream fout_final_vertices(vertices_output.str().c_str());
2338 * fout_final_vertices.close();
2340 * std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
2342 * if (boundary_that_we_need == 0)
2346 * Figure out if the vertex is on the boundary of the domain
2349 * for (typename Triangulation<dim>::active_cell_iterator cell =
2350 * triangulation.begin_active(); cell != triangulation.end(); ++cell)
2351 * for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
2353 * unsigned char boundary_ids = cell->face(f)->boundary_id();
2354 * if (boundary_ids == boundary_that_we_need)
2356 * for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
2357 * if (vertex_touched[cell->face(f)->vertex_index(v)] == false)
2359 * vertex_touched[cell->face(f)->vertex_index(v)] = true;
2360 * std::ofstream fout_final_vertices(vertices_output.str().c_str(), std::ios::app);
2361 * fout_final_vertices << cell->face(f)->vertex(v) << "\n";
2362 * fout_final_vertices.close();
2371 * Figure out if the vertex is on an internal boundary
2374 * for (typename Triangulation<dim>::active_cell_iterator cell =
2375 * triangulation.begin_active(); cell != triangulation.end(); ++cell)
2376 * for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
2378 * if (cell->neighbor(f) != triangulation.end())
2380 * if (cell->material_id() != cell->neighbor(f)->material_id()) //finds face is at internal boundary
2382 * int high_mat_id = std::max(cell->material_id(),
2383 * cell->neighbor(f)->material_id());
2384 * if (high_mat_id == boundary_that_we_need) //finds faces at the correct internal boundary
2386 * for (unsigned int v = 0;
2387 * v < GeometryInfo<dim>::vertices_per_face;
2389 * if (vertex_touched[cell->face(f)->vertex_index(
2392 * vertex_touched[cell->face(f)->vertex_index(
2394 * std::ofstream fout_final_vertices(vertices_output.str().c_str(), std::ios::app);
2395 * fout_final_vertices << cell->face(f)->vertex(v) << "\n";
2396 * fout_final_vertices.close();
2407 * ====================== SETUP INITIAL MESH ======================
2414 * void StokesProblem<dim>::setup_initial_mesh()
2416 * GridIn<dim> grid_in;
2417 * grid_in.attach_triangulation(triangulation);
2418 * std::ifstream mesh_stream(system_parameters::mesh_filename,
2419 * std::ifstream::in);
2420 * grid_in.read_ucd(mesh_stream);
2424 * output initial mesh in eps
2427 * std::ostringstream initial_mesh_file;
2428 * initial_mesh_file << system_parameters::output_folder << "/initial_mesh.eps";
2429 * std::ofstream out_eps (initial_mesh_file.str().c_str());
2431 * grid_out.write_eps (triangulation, out_eps);
2437 * boundary indicator 0 is outer free surface; 1, 2, 3 ... is boundary between layers, 99 is flat boundaries
2440 * typename Triangulation<dim>::active_cell_iterator
2441 * cell=triangulation.begin_active(), endc=triangulation.end();
2443 * unsigned int how_many; // how many components away from cardinal planes
2445 * std::ostringstream boundaries_file;
2446 * boundaries_file << system_parameters::output_folder << "/boundaries.txt";
2447 * std::ofstream fout_boundaries(boundaries_file.str().c_str());
2448 * fout_boundaries.close();
2450 * double zero_tolerance = 1e-3;
2451 * for (; cell != endc; ++cell) // loop over all cells
2453 * for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f) // loop over all vertices
2455 * if (cell->face(f)->at_boundary())
2462 * std::ofstream fout_boundaries(boundaries_file.str().c_str(), std::ios::app);
2463 * fout_boundaries << cell->face(f)->center()[0] << " " << cell->face(f)->center()[1]<< "\n";
2464 * fout_boundaries.close();
2467 * for (unsigned int i=0; i<dim; i++)
2468 * if (fabs(cell->face(f)->center()[i]) > zero_tolerance)
2470 * if (how_many==dim)
2471 * cell->face(f)->set_all_boundary_ids(0); // if face center coordinates > zero_tol, set bnry indicators to 0
2473 * cell->face(f)->set_all_boundary_ids(99);
2478 * std::ostringstream ellipses_filename;
2479 * ellipses_filename << system_parameters::output_folder << "/ellipse_fits.txt";
2480 * std::ofstream fout_ellipses(ellipses_filename.str().c_str());
2481 * fout_ellipses.close();
2485 * Find ellipsoidal axes for all layers
2488 * std::vector<double> ellipse_axes(0);
2491 * compute fit to boundary 0, 1, 2 ...
2494 * std::cout << endl;
2495 * for (unsigned int i = 0; i<system_parameters::sizeof_material_id; i++)
2497 * ellipsoid.compute_fit(ellipse_axes, system_parameters::material_id[i]);
2498 * system_parameters::q_axes.push_back(ellipse_axes[0]);
2499 * system_parameters::p_axes.push_back(ellipse_axes[1]);
2501 * std::cout << "a_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[0]
2502 * << " " << " c_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[1] << std::endl;
2503 * ellipse_axes.clear();
2505 * std::ofstream fout_ellipses(ellipses_filename.str().c_str(), std::ios::app);
2506 * fout_ellipses << system_parameters::present_timestep << " a_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[0]
2507 * << " " << " c_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[1] << endl;
2508 * fout_ellipses.close();
2511 * triangulation.refine_global(system_parameters::global_refinement);
2516 * refines crustal region
2519 * if (system_parameters::crustal_refinement != 0)
2521 * double a = system_parameters::q_axes[0] - system_parameters::crust_refine_region;
2522 * double b = system_parameters::p_axes[0] - system_parameters::crust_refine_region;
2525 * for (unsigned int step = 0;
2526 * step < system_parameters::crustal_refinement; ++step)
2528 * typename ::Triangulation<dim>::active_cell_iterator cell =
2529 * triangulation.begin_active(), endc = triangulation.end();
2530 * for (; cell != endc; ++cell)
2531 * for (unsigned int v = 0;
2532 * v < GeometryInfo<dim>::vertices_per_cell; ++v)
2534 * Point<dim> current_vertex = cell->vertex(v);
2536 * const double x_coord = current_vertex.operator()(0);
2537 * const double y_coord = current_vertex.operator()(1);
2538 * double expected_z = -1;
2540 * if ((x_coord - a) < -1e-10)
2542 * * std::sqrt(1 - (x_coord * x_coord / a / a));
2544 * if (y_coord >= expected_z)
2546 * cell->set_refine_flag();
2550 * triangulation.execute_coarsening_and_refinement();
2557 * output initial mesh in eps
2560 * std::ostringstream refined_mesh_file;
2561 * refined_mesh_file << system_parameters::output_folder << "/refined_mesh.eps";
2562 * std::ofstream out_eps_refined (refined_mesh_file.str().c_str());
2563 * GridOut grid_out_refined;
2564 * grid_out_refined.write_eps (triangulation, out_eps_refined);
2565 * out_eps_refined.close();
2566 * write_vertices(0);
2567 * write_vertices(1);
2573 * ====================== REFINE MESH ======================
2580 * void StokesProblem<dim>::refine_mesh()
2582 * using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
2583 * Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
2585 * std::vector<bool> component_mask(dim + 1, false);
2586 * component_mask[dim] = true;
2587 * KellyErrorEstimator<dim>::estimate(dof_handler, QGauss<dim - 1>(degree + 1),
2588 * FunctionMap(), solution,
2589 * estimated_error_per_cell, component_mask);
2591 * GridRefinement::refine_and_coarsen_fixed_number(triangulation,
2592 * estimated_error_per_cell, 0.3, 0.0);
2593 * triangulation.execute_coarsening_and_refinement();
2598 * ====================== SET UP THE DATA STRUCTURES TO REMEMBER STRESS FIELD ======================
2602 * void StokesProblem<dim>::setup_quadrature_point_history()
2604 * unsigned int our_cells = 0;
2605 * for (typename Triangulation<dim>::active_cell_iterator cell =
2606 * triangulation.begin_active(); cell != triangulation.end(); ++cell)
2609 * triangulation.clear_user_data();
2611 * quadrature_point_history.resize(our_cells * quadrature_formula.size());
2613 * unsigned int history_index = 0;
2614 * for (typename Triangulation<dim>::active_cell_iterator cell =
2615 * triangulation.begin_active(); cell != triangulation.end(); ++cell)
2617 * cell->set_user_pointer(&quadrature_point_history[history_index]);
2618 * history_index += quadrature_formula.size();
2621 * Assert(history_index == quadrature_point_history.size(), ExcInternalError());
2626 * ====================== DOES ELASTIC STEPS ======================
2630 * void StokesProblem<dim>::do_elastic_steps()
2632 * unsigned int elastic_iteration = 0;
2634 * while (elastic_iteration < system_parameters::initial_elastic_iterations)
2637 * std::cout << "\n\nElastic iteration " << elastic_iteration
2641 * if (system_parameters::present_timestep == 0)
2642 * initialize_eta_and_G();
2644 * if (elastic_iteration == 0)
2645 * system_parameters::current_time_interval =
2646 * system_parameters::viscous_time; //This is the time interval needed in assembling the problem
2648 * std::cout << " Assembling..." << std::endl << std::flush;
2649 * assemble_system();
2651 * std::cout << " Solving..." << std::flush;
2655 * update_quadrature_point_history();
2657 * append_physical_times(0);
2658 * elastic_iteration++;
2659 * system_parameters::present_timestep++;
2660 * do_ellipse_fits();
2661 * write_vertices(0);
2662 * write_vertices(1);
2664 * update_time_interval();
2670 * ====================== DO A SINGLE VISCOELASTOPLASTIC TIMESTEP ======================
2674 * void StokesProblem<dim>::do_flow_step()
2676 * plastic_iteration = 0;
2677 * while (plastic_iteration < system_parameters::max_plastic_iterations)
2679 * if (system_parameters::continue_plastic_iterations == true)
2681 * std::cout << "Plasticity iteration " << plastic_iteration << "\n";
2684 * std::cout << " Assembling..." << std::endl << std::flush;
2685 * assemble_system();
2687 * std::cout << " Solving..." << std::flush;
2691 * solution_stesses();
2693 * if (system_parameters::continue_plastic_iterations == false)
2696 * plastic_iteration++;
2703 * ====================== RUN ======================
2710 * void StokesProblem<dim>::run()
2714 * Sets up mesh and data structure for viscosity and stress at quadrature points
2717 * setup_initial_mesh();
2718 * setup_quadrature_point_history();
2722 * Makes the physical_times.txt file
2725 * std::ostringstream times_filename;
2726 * times_filename << system_parameters::output_folder << "/physical_times.txt";
2727 * std::ofstream fout_times(times_filename.str().c_str());
2728 * fout_times.close();
2732 * Computes elastic timesteps
2735 * do_elastic_steps();
2738 * Computes viscous timesteps
2741 * unsigned int VEPstep = 0;
2742 * while (system_parameters::present_timestep
2743 * < (system_parameters::initial_elastic_iterations
2744 * + system_parameters::total_viscous_steps))
2747 * if (system_parameters::continue_plastic_iterations == false)
2748 * system_parameters::continue_plastic_iterations = true;
2749 * std::cout << "\n\nViscoelastoplastic iteration " << VEPstep << "\n\n";
2752 * Computes plasticity
2756 * update_quadrature_point_history();
2758 * append_physical_times(plastic_iteration);
2759 * system_parameters::present_timestep++;
2760 * system_parameters::present_time = system_parameters::present_time + system_parameters::current_time_interval;
2761 * do_ellipse_fits();
2762 * write_vertices(0);
2763 * write_vertices(1);
2767 * append_physical_times(-1);
2774 * ====================== MAIN ======================
2780 * int main(int argc, char *argv[])
2785 * output program name
2788 * std::cout << "Running: " << argv[0] << std::endl;
2790 * char *cfg_filename = new char[120];
2792 * if (argc == 1) // if no input parameters (as if launched from eclipse)
2794 * std::strcpy(cfg_filename,"config/ConfigurationV2.cfg");
2797 * std::strcpy(cfg_filename,argv[1]);
2801 * using namespace dealii;
2802 * using namespace Step22;
2803 * config_in cfg(cfg_filename);
2807 * t1 = std::clock();
2809 * deallog.depth_console(0);
2811 * StokesProblem<2> flow_problem(1);
2812 * flow_problem.run();
2814 * std::cout << std::endl << "\a";
2816 * t2 = std::clock();
2817 * float diff (((float)t2 - (float)t1) / (float)CLOCKS_PER_SEC);
2818 * std::cout << "\n Program run in: " << diff << " seconds" << endl;
2820 * catch (std::exception &exc)
2822 * std::cerr << std::endl << std::endl
2823 * << "----------------------------------------------------"
2825 * std::cerr << "Exception on processing: " << std::endl << exc.what()
2826 * << std::endl << "Aborting!" << std::endl
2827 * << "----------------------------------------------------"
2834 * std::cerr << std::endl << std::endl
2835 * << "----------------------------------------------------"
2837 * std::cerr << "Unknown exception!" << std::endl << "Aborting!"
2839 * << "----------------------------------------------------"
2849 <a name="ann-support_code/config_in.h"></a>
2850 <h1>Annotated version of support_code/config_in.h</h1>
2859 * * Created on: Aug 17, 2015
2860 * * Author: antonermakov
2864 * #include <iostream>
2865 * #include <iomanip>
2866 * #include <cstdlib>
2868 * #include <libconfig.h++>
2870 * #include "local_math.h"
2872 * using namespace std;
2873 * using namespace libconfig;
2877 * namespace system_parameters
2885 * string mesh_filename;
2886 * string output_folder;
2895 * double omegasquared;
2901 * Rheology parameters
2904 * vector<double> depths_eta;
2905 * vector<double> eta_kinks;
2906 * vector<double> depths_rho;
2907 * vector<double> rho;
2908 * vector<int> material_id;
2911 * double eta_ceiling;
2914 * bool lat_dependence;
2916 * unsigned int sizeof_depths_eta;
2917 * unsigned int sizeof_depths_rho;
2918 * unsigned int sizeof_rho;
2919 * unsigned int sizeof_eta_kinks;
2920 * unsigned int sizeof_material_id;
2921 * unsigned int sizeof_G;
2923 * double pressure_scale;
2926 * bool continue_plastic_iterations;
2930 * plasticity variables
2933 * bool plasticity_on;
2934 * unsigned int failure_criterion;
2935 * unsigned int max_plastic_iterations;
2936 * double smoothing_radius;
2940 * viscoelasticity variables
2943 * unsigned int initial_elastic_iterations;
2944 * double elastic_time;
2945 * double viscous_time;
2946 * double initial_disp_target;
2947 * double final_disp_target;
2948 * double current_time_interval;
2952 * mesh refinement variables
2955 * unsigned int global_refinement;
2956 * unsigned int small_r_refinement;
2957 * unsigned int crustal_refinement;
2958 * double crust_refine_region;
2959 * unsigned int surface_refinement;
2966 * int iteration_coefficient;
2967 * double tolerance_coefficient;
2971 * time step variables
2974 * double present_time;
2975 * unsigned int present_timestep;
2976 * unsigned int total_viscous_steps;
2984 * vector<double> q_axes;
2985 * vector<double> p_axes;
2992 * config_in(char *);
2995 * void write_config();
2998 * void config_in::write_config()
3000 * std::ostringstream config_parameters;
3001 * config_parameters << system_parameters::output_folder << "/run_parameters.txt";
3002 * std::ofstream fout_config(config_parameters.str().c_str());
3009 * fout_config << "mesh filename: " << system_parameters::mesh_filename << endl << endl;
3016 * fout_config << "r_mean = " << system_parameters::r_mean << endl;
3017 * fout_config << "period = " << system_parameters::period << endl;
3018 * fout_config << "omegasquared = " << system_parameters::omegasquared << endl;
3019 * fout_config << "beta = " << system_parameters::beta << endl;
3020 * fout_config << "intercept = " << system_parameters::intercept << endl;
3024 * rheology parameters
3030 * for (unsigned int i=0; i<system_parameters::sizeof_depths_eta; i++)
3031 * fout_config << "depths_eta[" << i << "] = " << system_parameters::depths_eta[i] << endl;
3033 * for (unsigned int i=0; i<system_parameters::sizeof_eta_kinks; i++)
3034 * fout_config << "eta_kinks[" << i << "] = " << system_parameters::eta_kinks[i] << endl;
3036 * for (unsigned int i=0; i<system_parameters::sizeof_depths_rho; i++)
3037 * fout_config << "depths_rho[" << i << "] = " << system_parameters::depths_rho[i] << endl;
3039 * for (unsigned int i=0; i<system_parameters::sizeof_rho; i++)
3040 * fout_config << "rho[" << i << "] = " << system_parameters::rho[i] << endl;
3042 * for (unsigned int i=0; i<system_parameters::sizeof_material_id; i++)
3043 * fout_config << "material_id[" << i << "] = " << system_parameters::material_id[i] << endl;
3045 * for (unsigned int i=0; i<system_parameters::sizeof_G; i++)
3046 * fout_config << "G[" << i << "] = " << system_parameters::G[i] << endl;
3048 * fout_config << "eta_ceiling = " << system_parameters::eta_ceiling << endl;
3049 * fout_config << "eta_floor = " << system_parameters::eta_floor << endl;
3050 * fout_config << "eta_Ea = " << system_parameters::eta_Ea << endl;
3051 * fout_config << "lat_dependence = " << system_parameters::lat_dependence << endl;
3052 * fout_config << "pressure_scale = " << system_parameters::pressure_scale << endl;
3053 * fout_config << "q = " << system_parameters::q << endl;
3054 * fout_config << "cylindrical = " << system_parameters::cylindrical << endl;
3055 * fout_config << "continue_plastic_iterations = " << system_parameters::continue_plastic_iterations << endl;
3059 * Plasticity parameters
3062 * fout_config << "plasticity_on = " << system_parameters::plasticity_on << endl;
3063 * fout_config << "failure_criterion = " << system_parameters::failure_criterion << endl;
3064 * fout_config << "max_plastic_iterations = " << system_parameters::max_plastic_iterations << endl;
3065 * fout_config << "smoothing_radius = " << system_parameters::smoothing_radius << endl;
3069 * Viscoelasticity parameters
3072 * fout_config << "initial_elastic_iterations = " << system_parameters::initial_elastic_iterations << endl;
3073 * fout_config << "elastic_time = " << system_parameters::elastic_time << endl;
3074 * fout_config << "viscous_time = " << system_parameters::viscous_time << endl;
3075 * fout_config << "initial_disp_target = " << system_parameters::initial_disp_target << endl;
3076 * fout_config << "final_disp_target = " << system_parameters::final_disp_target << endl;
3077 * fout_config << "current_time_interval = " << system_parameters::current_time_interval << endl;
3081 * Mesh refinement parameters
3084 * fout_config << "global_refinement = " << system_parameters::global_refinement << endl;
3085 * fout_config << "small_r_refinement = " << system_parameters::small_r_refinement << endl;
3086 * fout_config << "crustal_refinement = " << system_parameters::crustal_refinement << endl;
3087 * fout_config << "crust_refine_region = " << system_parameters::crust_refine_region << endl;
3088 * fout_config << "surface_refinement = " << system_parameters::surface_refinement << endl;
3095 * fout_config << "iteration_coefficient = " << system_parameters::iteration_coefficient << endl;
3096 * fout_config << "tolerance_coefficient = " << system_parameters::tolerance_coefficient << endl;
3100 * Time step parameters
3103 * fout_config << "present_time = " << system_parameters::present_time << endl;
3104 * fout_config << "present_timestep = " << system_parameters::present_timestep << endl;
3105 * fout_config << "total_viscous_steps = " << system_parameters::total_viscous_steps << endl;
3107 * fout_config.close();
3110 * config_in::config_in(char *filename)
3115 * This example reads the configuration file 'example.cfg
' and displays
3116 * some of its contents.
3126 * Read the file. If there is an error, report it and exit.
3131 * cfg.readFile(filename);
3133 * catch (const FileIOException &fioex)
3135 * std::cerr << "I/O error while reading file:" << filename << std::endl;
3137 * catch (const ParseException &pex)
3139 * std::cerr << "Parse error at " << pex.getFile() << ":" << pex.getLine()
3140 * << " - " << pex.getError() << std::endl;
3150 * string msh = cfg.lookup("mesh_filename");
3151 * system_parameters::mesh_filename = msh;
3153 * catch (const SettingNotFoundException &nfex)
3155 * cerr << "No 'mesh_filename
' setting in configuration file." << endl;
3168 * string output = cfg.lookup("output_folder");
3169 * system_parameters::output_folder = output;
3171 * std::cout << "Writing to folder: " << output << endl;
3173 * catch (const SettingNotFoundException &nfex)
3175 * cerr << "No 'output_folder
' setting in configuration file." << endl;
3186 * const Setting &root = cfg.getRoot();
3191 * get body parameters
3196 * const Setting &body_parameters = root["body_parameters"];
3198 * body_parameters.lookupValue("period", system_parameters::period);
3199 * system_parameters::omegasquared = pow(TWOPI / 3600.0 / system_parameters::period, 2.0);
3200 * body_parameters.lookupValue("r_mean", system_parameters::r_mean);
3201 * body_parameters.lookupValue("beta", system_parameters::beta);
3202 * body_parameters.lookupValue("intercept", system_parameters::intercept);
3204 * catch (const SettingNotFoundException &nfex)
3206 * cerr << "We've got a problem in the body parameters block
" << endl;
3212 * Rheology parameters
3222 * get depths_eta ---------------------
3225 * const Setting &set_depths_eta = cfg.lookup("rheology_parameters.depths_eta
");
3227 * unsigned int ndepths_eta = set_depths_eta.getLength();
3228 * system_parameters::sizeof_depths_eta = ndepths_eta;
3230 * for (unsigned int i=0; i<ndepths_eta; i++)
3232 * system_parameters::depths_eta.push_back(set_depths_eta[i]);
3233 * cout << "depth_eta[
" << i << "] =
" << system_parameters::depths_eta[i] << endl;
3238 * get eta_kinks -------------------------
3241 * const Setting &set_eta_kinks = cfg.lookup("rheology_parameters.eta_kinks
");
3243 * unsigned int neta_kinks = set_eta_kinks.getLength();
3244 * system_parameters::sizeof_eta_kinks = neta_kinks;
3248 * cout << "Number of
depth =
" << ndepths << endl;
3254 * for (unsigned int i=0; i<neta_kinks; i++)
3256 * system_parameters::eta_kinks.push_back(set_eta_kinks[i]);
3257 * cout << "eta_kinks[
" << i << "] =
" << system_parameters::eta_kinks[i] << endl;
3262 * get depths_rho -------------------------
3265 * const Setting &set_depths_rho = cfg.lookup("rheology_parameters.depths_rho
");
3267 * unsigned int ndepths_rho = set_depths_rho.getLength();
3268 * system_parameters::sizeof_depths_rho = ndepths_rho;
3272 * cout << "Number of
depth =
" << ndepths << endl;
3278 * for (unsigned int i=0; i<ndepths_rho; i++)
3280 * system_parameters::depths_rho.push_back(set_depths_rho[i]);
3281 * cout << "depths_rho[
" << i << "] =
" << system_parameters::depths_rho[i] << endl;
3286 * get rho -------------------------
3289 * const Setting &set_rho = cfg.lookup("rheology_parameters.rho
");
3291 * unsigned int nrho = set_rho.getLength();
3292 * system_parameters::sizeof_rho = nrho;
3296 * cout << "Number of
depth =
" << ndepths << endl;
3302 * for (unsigned int i=0; i<nrho; i++)
3304 * system_parameters::rho.push_back(set_rho[i]);
3305 * cout << "rho[
" << i << "] =
" << system_parameters::rho[i] << endl;
3310 * get material_id -------------------------
3313 * const Setting &set_material_id = cfg.lookup("rheology_parameters.material_id
");
3315 * unsigned int nmaterial_id = set_material_id.getLength();
3316 * system_parameters::sizeof_material_id = nmaterial_id;
3320 * cout << "Number of
depth =
" << ndepths << endl;
3326 * for (unsigned int i=0; i<nmaterial_id; i++)
3328 * system_parameters::material_id.push_back(set_material_id[i]);
3329 * cout << "material_id[
" << i << "] =
" << system_parameters::material_id[i] << endl;
3334 * get G -------------------------
3337 * const Setting &set_G = cfg.lookup("rheology_parameters.G
");
3339 * unsigned int nG = set_G.getLength();
3340 * system_parameters::sizeof_G = nG;
3344 * cout << "Number of
depth =
" << ndepths << endl;
3350 * for (unsigned int i=0; i<nG; i++)
3352 * system_parameters::G.push_back(set_G[i]);
3353 * cout << "G[
" << i << "] =
" << system_parameters::G[i] << endl;
3356 * const Setting &rheology_parameters = root["rheology_parameters
"];
3357 * rheology_parameters.lookupValue("eta_ceiling
", system_parameters::eta_ceiling);
3358 * rheology_parameters.lookupValue("eta_floor
", system_parameters::eta_floor);
3359 * rheology_parameters.lookupValue("eta_Ea
", system_parameters::eta_Ea);
3360 * rheology_parameters.lookupValue("lat_dependence
", system_parameters::lat_dependence);
3361 * rheology_parameters.lookupValue("pressure_scale
", system_parameters::pressure_scale);
3362 * rheology_parameters.lookupValue("q
", system_parameters::q);
3363 * rheology_parameters.lookupValue("cylindrical
", system_parameters::cylindrical);
3364 * rheology_parameters.lookupValue("continue_plastic_iterations
", system_parameters::continue_plastic_iterations);
3366 * catch (const SettingNotFoundException &nfex)
3368 * cerr << "We
've got a problem in the rheology parameters block" << endl;
3373 * Plasticity parameters
3379 * const Setting &plasticity_parameters = root["plasticity_parameters"];
3380 * plasticity_parameters.lookupValue("plasticity_on", system_parameters::plasticity_on);
3381 * plasticity_parameters.lookupValue("failure_criterion", system_parameters::failure_criterion);
3382 * plasticity_parameters.lookupValue("max_plastic_iterations", system_parameters::max_plastic_iterations);
3383 * plasticity_parameters.lookupValue("smoothing_radius", system_parameters::smoothing_radius);
3385 * catch (const SettingNotFoundException &nfex)
3387 * cerr << "We've got a problem in the plasticity parameters block
" << endl;
3392 * Viscoelasticity parameters
3401 * const Setting &viscoelasticity_parameters = root["viscoelasticity_parameters
"];
3402 * viscoelasticity_parameters.lookupValue("initial_elastic_iterations
", system_parameters::initial_elastic_iterations);
3403 * viscoelasticity_parameters.lookupValue("elastic_time
", system_parameters::elastic_time);
3404 * viscoelasticity_parameters.lookupValue("viscous_time
", system_parameters::viscous_time);
3405 * viscoelasticity_parameters.lookupValue("initial_disp_target
", system_parameters::initial_disp_target);
3406 * viscoelasticity_parameters.lookupValue("final_disp_target
", system_parameters::final_disp_target);
3407 * viscoelasticity_parameters.lookupValue("current_time_interval
", system_parameters::current_time_interval);
3409 * system_parameters::viscous_time *= SECSINYEAR;
3411 * catch (const SettingNotFoundException &nfex)
3413 * cerr << "We
've got a problem in the viscoelasticity parameters block" << endl;
3418 * Mesh refinement parameters
3424 * const Setting &mesh_refinement_parameters = root["mesh_refinement_parameters"];
3425 * mesh_refinement_parameters.lookupValue("global_refinement", system_parameters::global_refinement);
3426 * mesh_refinement_parameters.lookupValue("small_r_refinement", system_parameters::small_r_refinement);
3427 * mesh_refinement_parameters.lookupValue("crustal_refinement", system_parameters::crustal_refinement);
3428 * mesh_refinement_parameters.lookupValue("crust_refine_region", system_parameters::crust_refine_region);
3429 * mesh_refinement_parameters.lookupValue("surface_refinement", system_parameters::surface_refinement);
3431 * catch (const SettingNotFoundException &nfex)
3433 * cerr << "We've got a problem in the mesh refinement parameters block
" << endl;
3443 * const Setting &solve_parameters = root["solve_parameters
"];
3444 * solve_parameters.lookupValue("iteration_coefficient
", system_parameters::iteration_coefficient);
3445 * solve_parameters.lookupValue("tolerance_coefficient
", system_parameters::tolerance_coefficient);
3449 * catch (const SettingNotFoundException &nfex)
3451 * cerr << "We
've got a problem in the solver parameters block" << endl;
3456 * Time step parameters
3461 * const Setting &time_step_parameters = root["time_step_parameters"];
3462 * time_step_parameters.lookupValue("present_time", system_parameters::present_time);
3463 * time_step_parameters.lookupValue("present_timestep", system_parameters::present_timestep);
3464 * time_step_parameters.lookupValue("total_viscous_steps", system_parameters::total_viscous_steps);
3466 * catch (const SettingNotFoundException &nfex)
3468 * cerr << "We've got a problem in the time step parameters block
" << endl;
3482 <a name="ann-support_code/ellipsoid_fit.h
"></a>
3483 <h1>Annotated version of support_code/ellipsoid_fit.h</h1>
3492 * * Created on: Jul 24, 2015
3493 * * Author: antonermakov
3496 * #include <deal.II/base/quadrature_lib.h>
3497 * #include <deal.II/base/function.h>
3498 * #include <deal.II/base/logstream.h>
3499 * #include <deal.II/lac/vector.h>
3500 * #include <deal.II/lac/full_matrix.h>
3501 * #include <deal.II/lac/sparse_matrix.h>
3502 * #include <deal.II/lac/solver_cg.h>
3503 * #include <deal.II/lac/precondition.h>
3504 * #include <deal.II/grid/tria.h>
3505 * #include <deal.II/dofs/dof_handler.h>
3506 * #include <deal.II/grid/tria_accessor.h>
3507 * #include <deal.II/grid/tria_iterator.h>
3508 * #include <deal.II/dofs/dof_accessor.h>
3509 * #include <deal.II/dofs/dof_tools.h>
3510 * #include <deal.II/fe/fe_q.h>
3511 * #include <deal.II/fe/fe_values.h>
3512 * #include <deal.II/numerics/vector_tools.h>
3513 * #include <deal.II/numerics/matrix_tools.h>
3514 * #include <deal.II/numerics/data_out.h>
3515 * #include <deal.II/grid/grid_in.h>
3516 * #include <deal.II/grid/grid_out.h>
3517 * #include <deal.II/grid/tria_boundary_lib.h>
3518 * #include <deal.II/base/point.h>
3519 * #include <deal.II/grid/grid_generator.h>
3521 * #include <fstream>
3522 * #include <sstream>
3523 * #include <iostream>
3524 * #include <iomanip>
3525 * #include <cstdlib>
3528 * #include "local_math.h
"
3530 * using namespace dealii;
3532 * template <int dim>
3533 * class ellipsoid_fit
3536 * inline ellipsoid_fit (Triangulation<dim,dim> *pi)
3538 * p_triangulation = pi;
3540 * void compute_fit(std::vector<double> &ell, unsigned char bndry);
3544 * Triangulation<dim,dim> *p_triangulation;
3551 * This function computes ellipsoid fit to a set of vertices that lie on the
3552 * boundary_that_we_need
3555 * template <int dim>
3556 * void ellipsoid_fit<dim>::compute_fit(std::vector<double> &ell, unsigned char boundary_that_we_need)
3558 * typename Triangulation<dim>::active_cell_iterator cell = p_triangulation->begin_active();
3559 * typename Triangulation<dim>::active_cell_iterator endc = p_triangulation->end();
3561 * FullMatrix<double> A(p_triangulation->n_vertices(),dim);
3562 * Vector<double> x(dim);
3563 * Vector<double> b(p_triangulation->n_vertices());
3565 * std::vector<bool> vertex_touched (p_triangulation->n_vertices(),
3568 * unsigned int j = 0;
3569 * unsigned char boundary_ids;
3570 * std::vector<unsigned int> ind_bnry_row;
3571 * std::vector<unsigned int> ind_bnry_col;
3575 * assemble the sensitivity matrix and r.h.s.
3578 * for (; cell != endc; ++cell)
3580 * if (boundary_that_we_need != 0)
3581 * cell->set_manifold_id(cell->material_id());
3582 * for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
3584 * if (boundary_that_we_need == 0) //if this is the outer surface, then look for boundary ID 0; otherwise look for material ID change.
3586 * boundary_ids = cell->face(f)->boundary_id();
3587 * if (boundary_ids == boundary_that_we_need)
3589 * for (unsigned int v = 0;
3590 * v < GeometryInfo<dim>::vertices_per_face; ++v)
3591 * if (vertex_touched[cell->face(f)->vertex_index(v)]
3594 * vertex_touched[cell->face(f)->vertex_index(v)] =
3596 * for (unsigned int i = 0; i < dim; ++i)
3600 * stiffness matrix entry
3603 * A(j, i) = pow(cell->face(f)->vertex(v)[i], 2);
3612 * if mesh if not full: set the indicator
3616 * ind_bnry_row.push_back(j);
3621 * else //find the faces that are at the boundary between materials, get the vertices, and write them into the stiffness matrix
3623 * if (cell->neighbor(f) != endc)
3625 * if (cell->material_id() != cell->neighbor(f)->material_id()) //finds face is at internal boundary
3627 * int high_mat_id = std::max(cell->material_id(),
3628 * cell->neighbor(f)->material_id());
3629 * if (high_mat_id == boundary_that_we_need) //finds faces at the correct internal boundary
3631 * for (unsigned int v = 0;
3632 * v < GeometryInfo<dim>::vertices_per_face;
3634 * if (vertex_touched[cell->face(f)->vertex_index(
3637 * vertex_touched[cell->face(f)->vertex_index(
3639 * for (unsigned int i = 0; i < dim; ++i)
3643 * stiffness matrix entry
3647 * cell->face(f)->vertex(v)[i], 2);
3656 * if mesh if not full: set the indicator
3660 * ind_bnry_row.push_back(j);
3669 * if (ind_bnry_row.size()>0)
3674 * maxtrix A'*A and vector A'*b; A'*A*x = A'*b -- normal system of equations
3677 * FullMatrix<double> AtA(dim,dim);
3678 * Vector<double> Atb(dim);
3680 * FullMatrix<double> A_out(ind_bnry_row.size(),dim);
3681 * Vector<double> b_out(ind_bnry_row.size());
3683 * for (unsigned int i=0; i<dim; i++)
3684 * ind_bnry_col.push_back(i);
3686 * for (unsigned int i=0; i<ind_bnry_row.size(); i++)
3689 * A_out.extract_submatrix_from(A, ind_bnry_row, ind_bnry_col);
3690 * A_out.Tmmult(AtA,A_out,true);
3691 * A_out.Tvmult(Atb,b_out,true);
3695 * solve normal system of equations
3698 * SolverControl solver_control (1000, 1e-12);
3699 * SolverCG<> solver (solver_control);
3700 * solver.solve (AtA, x, Atb, PreconditionIdentity());
3704 * find ellipsoidal axes
3707 * for (unsigned int i=0; i<dim; i++)
3708 * ell.push_back(sqrt(1.0/x[i]));
3711 * std::cerr << "fit_ellipsoid: no points to fit
" << std::endl;
3718 <a name="ann-support_code/ellipsoid_grav.h
"></a>
3719 <h1>Annotated version of support_code/ellipsoid_grav.h</h1>
3729 * Started 10/8/2012, R. R. Fu
3731 * Reference Pohanka 2011, Contrib. Geophys. Geodes.
3735 * #include <deal.II/base/point.h>
3736 * #include <fstream>
3737 * #include <iostream>
3739 * namespace A_Grav_namespace
3741 * namespace system_parameters
3743 * double mantle_rho;
3745 * double excess_rho;
3750 * double r_core_polar;
3753 * template <int dim>
3754 * class AnalyticGravity
3757 * void setup_vars (std::vector<double> v);
3758 * void get_gravity (const ::Point<dim> &p, std::vector<double> &g);
3777 * template <int dim>
3778 * void AnalyticGravity<dim>::get_gravity (const ::Point<dim> &p, std::vector<double> &g)
3780 * double rsph = std::sqrt(p[0] * p[0] + p[1] * p[1]);
3781 * double thetasph = std::atan2(p[0], p[1]);
3782 * double costhetasph = std::cos(thetasph);
3786 * convert to elliptical coordinates for silicates
3789 * double stemp = std::sqrt((rsph * rsph - eV * eV + std::sqrt((rsph * rsph - eV * eV) * (rsph * rsph - eV * eV)
3790 * + 4 * eV * eV * rsph * rsph * costhetasph *costhetasph)) / 2);
3791 * double vout = stemp / system_parameters::r_eq / std::sqrt(1 - ecc * ecc);
3792 * double eout = std::acos(rsph * costhetasph / stemp);
3796 * convert to elliptical coordinates for core correction
3799 * double stemp_c = std::sqrt((rsph * rsph - eV_c * eV_c + std::sqrt((rsph * rsph - eV_c * eV_c) * (rsph * rsph - eV_c * eV_c)
3800 * + 4 * eV_c * eV_c * rsph * rsph * costhetasph *costhetasph)) / 2);
3801 * double vout_c = stemp_c / system_parameters::r_core_eq / std::sqrt(1 - ecc_c * ecc_c);
3802 * double eout_c = std::acos(rsph * costhetasph / stemp_c);
3806 * shell contribution
3809 * g[0] = g_coeff * r11 * std::sqrt((1 - ecc * ecc) * vout * vout + ecc * ecc) * std::sin(eout);
3810 * g[1] = g_coeff * r01 * vout * std::cos(eout) / std::sqrt(1 - ecc * ecc);
3817 * double expected_y = system_parameters::r_core_polar * std::sqrt(1 -
3818 * (p[0] * p[0] / system_parameters::r_core_eq / system_parameters::r_core_eq));
3821 * if (p[1] <= expected_y)
3823 * g[0] += g_coeff_c * r11_c * std::sqrt((1 - ecc_c * ecc_c) * vout_c * vout_c + ecc_c * ecc_c) * std::sin(eout_c);
3824 * g[1] += g_coeff_c * r01_c * vout_c * std::cos(eout_c) / std::sqrt(1 - ecc_c * ecc_c);
3828 * double g_coeff_co = - 2.795007963255562e-10 * system_parameters::excess_rho * system_parameters::r_core_eq
3829 * / vout_c / vout_c;
3830 * double r00_co = 0;
3831 * double r01_co = 0;
3832 * double r11_co = 0;
3834 * if (system_parameters::r_core_polar == system_parameters::r_core_eq)
3842 * r00_co = ke_c * vout_c * std::atan2(1, ke_c * vout_c);
3843 * double ke_co2 = ke_c * ke_c * vout_c * vout_c;
3844 * r01_co = 3 * ke_co2 * (1 - r00_co);
3845 * r11_co = 3 * ((ke_co2 + 1) * r00_co - ke_co2) / 2;
3847 * g[0] += g_coeff_co * vout_c * r11_co / std::sqrt((1 - ecc_c* ecc_c) * vout_c * vout_c + ecc_c * ecc_c) * std::sin(eout_c);
3848 * g[1] += g_coeff_co * r01_co * std::cos(eout_c) / std::sqrt(1 - ecc_c * ecc_c);
3852 * template <int dim>
3853 * void AnalyticGravity<dim>::setup_vars (std::vector<double> v)
3855 * system_parameters::r_eq = v[0];
3856 * system_parameters::r_polar = v[1];
3857 * system_parameters::r_core_eq = v[2];
3858 * system_parameters::r_core_polar = v[3];
3859 * system_parameters::mantle_rho = v[4];
3860 * system_parameters::core_rho = v[5];
3861 * system_parameters::excess_rho = system_parameters::core_rho - system_parameters::mantle_rho;
3869 * if (system_parameters::r_polar > system_parameters::r_eq)
3873 * This makes the gravity field nearly that of a sphere in case the body becomes prolate
3876 * std::cout << "\nWarning: The model body has become prolate. \n
";
3881 * ecc = std::sqrt(1 - (system_parameters::r_polar * system_parameters::r_polar / system_parameters::r_eq / system_parameters::r_eq));
3884 * eV = ecc * system_parameters::r_eq;
3885 * ke = std::sqrt(1 - (ecc * ecc)) / ecc;
3886 * r00 = ke * std::atan2(1, ke);
3887 * double ke2 = ke * ke;
3888 * r01 = 3 * ke2 * (1 - r00);
3889 * r11 = 3 * ((ke2 + 1) * r00 - ke2) / 2;
3890 * g_coeff = - 2.795007963255562e-10 * system_parameters::mantle_rho * system_parameters::r_eq;
3897 * if (system_parameters::r_core_polar > system_parameters::r_core_eq)
3899 * std::cout << "\nWarning: The model core has become prolate. \n
";
3904 * ecc_c = std::sqrt(1 - (system_parameters::r_core_polar * system_parameters::r_core_polar / system_parameters::r_core_eq / system_parameters::r_core_eq));
3906 * eV_c = ecc_c * system_parameters::r_core_eq;
3907 * if (system_parameters::r_core_polar == system_parameters::r_core_eq)
3913 * g_coeff_c = - 2.795007963255562e-10 * system_parameters::excess_rho * system_parameters::r_core_eq;
3917 * ke_c = std::sqrt(1 - (ecc_c * ecc_c)) / ecc_c;
3918 * r00_c = ke_c * std::atan2(1, ke_c);
3919 * double ke2_c = ke_c * ke_c;
3920 * r01_c = 3 * ke2_c * (1 - r00_c);
3921 * r11_c = 3 * ((ke2_c + 1) * r00_c - ke2_c) / 2;
3922 * g_coeff_c = - 2.795007963255562e-10 * system_parameters::excess_rho * system_parameters::r_core_eq;
3926 * std::cout << "Loaded variables: ecc =
" << ecc_c << " ke =
" << ke_c << " r00 =
" << r00_c << " r01 =
" << r01_c << " r11 =
" << r11_c << "\n
";
3934 <a name="ann-support_code/local_math.h
"></a>
3935 <h1>Annotated version of support_code/local_math.h</h1>
3942 * * File: localmath.h
3943 * * Author: antonermakov
3945 * * Created on September 21, 2013, 7:14 PM
3948 * #ifndef LOCAL_MATH_
3949 * #define LOCAL_MATH_
3951 * #define PI 3.14159265358979323846
3952 * #define TWOPI 6.283185307179586476925287
3953 * #define SECSINYEAR 3.155692608e+07
3956 * #define ABS(a) ((a) < 0 ? -(a) : (a))
3960 * double factorial(int n)
3964 * } else if(n == 1) {
3966 * } else if(n == 2) {
3968 * } else if(n == 3) {
3970 * } else if(n == 4) {
3979 * double fudge(int m)
3993 * double sign(double x)
3997 * } else if(x < 0.0) {
4006 * double pv0(double x)
4012 * ans = x - TWOPI*floor(x/TWOPI);
4013 * if(ans > TWOPI/2.0) {
4014 * ans = ans - TWOPI;
4029 * double System::Plm(int m, double x)
4032 * return(1.5*x*x - 0.5);
4033 * } else if(m == 1) {
4034 * return(3.0*x*sqrt(1.0 - x*x));
4035 * } else if(m == 2) {
4036 * return(3.0 - 3.0*x*x);
4044 * double System::DP(int m, double x)
4048 * } else if(m == 1) {
4049 * return((3.0 - 6.0*x*x)/sqrt(1.0 - x*x));
4050 * } else if(m == 2) {
4064 * #endif /* LOCALMATH_H */