154 * #include <deal.II/base/quadrature_lib.h>
155 * #include <deal.II/base/logstream.h>
156 * #include <deal.II/base/function.h>
157 * #include <deal.II/base/utilities.h>
159 * #include <deal.II/lac/block_vector.h>
160 * #include <deal.II/lac/full_matrix.h>
161 * #include <deal.II/lac/block_sparse_matrix.h>
162 * #include <deal.II/lac/solver_cg.h>
163 * #include <deal.II/lac/precondition.h>
164 * #include <deal.II/lac/affine_constraints.h>
166 * #include <deal.II/grid/tria.h>
167 * #include <deal.II/grid/grid_generator.h>
168 * #include <deal.II/grid/tria_accessor.h>
169 * #include <deal.II/grid/tria_iterator.h>
170 * #include <deal.II/grid/grid_tools.h>
171 * #include <deal.II/grid/manifold_lib.h>
172 * #include <deal.II/grid/grid_refinement.h>
173 * #include <deal.II/grid/grid_in.h>
175 * #include <deal.II/dofs/dof_handler.h>
176 * #include <deal.II/dofs/dof_renumbering.h>
177 * #include <deal.II/dofs/dof_accessor.h>
178 * #include <deal.II/dofs/dof_tools.h>
180 * #include <deal.II/fe/fe_q.h>
181 * #include <deal.II/fe/fe_system.h>
182 * #include <deal.II/fe/fe_values.h>
184 * #include <deal.II/numerics/vector_tools.h>
185 * #include <deal.II/numerics/matrix_tools.h>
186 * #include <deal.II/numerics/data_out.h>
187 * #include <deal.II/numerics/error_estimator.h>
188 * #include <deal.II/numerics/derivative_approximation.h>
189 * #include <deal.II/numerics/fe_field_function.h>
191 * #include <deal.II/lac/sparse_direct.h>
192 * #include <deal.II/lac/sparse_ilu.h>
194 * #include <iostream>
202 * #include <armadillo>
205 * #include
"../support_code/ellipsoid_grav.h"
206 * #include
"../support_code/ellipsoid_fit.h"
207 * #include
"../support_code/config_in.h"
211 * As in all programs, the
namespace dealii
219 *
using namespace arma;
222 *
struct InnerPreconditioner;
225 *
struct InnerPreconditioner<2>
231 *
struct InnerPreconditioner<3>
255 *
const double curl = (grad_u[1][0] - grad_u[0][1]);
259 *
const double t[2][2] = { {
cos(angle),
sin(angle) }, {
322 *
const unsigned int degree;
328 *
unsigned int n_u = 0,
n_p = 0;
345 *
std::shared_ptr<typename InnerPreconditioner<dim>::type>
A_preconditioner;
375 *
const unsigned int component)
const
378 *
ExcIndexRange (component, 0, this->n_components));
391 *
for (
unsigned int c = 0; c < this->n_components; ++c)
402 *
using Function<dim>::vector_value;
403 *
using Function<dim>::vector_value_list;
405 *
virtual double value(
const Point<dim> &p,
const unsigned int component,
411 *
virtual void vector_value_list(
const std::vector<
Point<dim> > &points,
419 *
const unsigned int component,
426 *
if (component == 0)
432 *
if (component == 1)
444 *
for (
unsigned int c = 0; c < this->n_components; ++c)
462 *
ExcDimensionMismatch(
values.size(), points.size()));
464 *
for (
unsigned int i = 0; i < points.size(); ++i)
465 *
this->vector_value(points[i], values[i],
aGrav);
476 *
template<
class Matrix,
class Preconditioner>
489 *
template<
class Matrix,
class Preconditioner>
492 *
matrix(&m), preconditioner(&preconditioner)
506 *
cg.solve(*matrix, dst, src, *preconditioner);
517 *
template<
class Preconditioner>
533 *
template<
class Preconditioner>
538 *
system_matrix.block(0, 0).m()),
tmp2(
539 *
system_matrix.block(0, 0).m())
547 *
system_matrix->block(0, 1).vmult(
tmp1, src);
549 *
system_matrix->block(1, 0).vmult(dst,
tmp2);
565 *
fe(
FE_Q<dim>(degree + 1), dim,
FE_Q<dim>(degree), 1),
583 *
system_matrix.
clear();
585 *
dof_handler.distribute_dofs(fe);
594 * ========================================Apply
Boundary Conditions=====================================
598 *
constraints.
clear();
612 *
constraints.close();
618 *
std::cout <<
" Number of active cells: " <<
triangulation.n_active_cells()
619 *
<< std::endl <<
" Number of degrees of freedom: "
620 *
<< dof_handler.n_dofs() <<
" (" <<
n_u <<
'+' <<
n_p <<
')'
631 *
csp.collect_sizes();
634 *
sparsity_pattern.copy_from(
csp);
637 *
system_matrix.reinit(sparsity_pattern);
639 *
solution.reinit(2);
640 *
solution.block(0).reinit(
n_u);
641 *
solution.block(1).reinit(
n_p);
642 *
solution.collect_sizes();
663 *
double get_eta(
double &
r,
double &z);
689 * compute local
depth
698 *
double error = 10000;
713 *
while (
error >= eps)
763 *
unsigned int ndeep = 2 * n - 2;
825 *
return system_parameters::G[
mat_id];
846 *
dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
856 *
fe_values.reinit(cell);
858 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
861 *
double r_value = fe_values.quadrature_point(
q)[0];
862 *
double z_value = fe_values.quadrature_point(
q)[1];
883 *
unsigned int mat_id = cell->material_id();
916 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
923 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
942 *
std::cout <<
"Body parameters are: " ;
943 *
for (
int i=0; i<6; ++i)
949 *
std::vector<Vector<double> >
rhs_values(n_q_points,
955 *
std::vector<SymmetricTensor<2, dim> >
phi_grads_u(dofs_per_cell);
956 *
std::vector<double>
div_phi_u(dofs_per_cell);
957 *
std::vector<Tensor<1, dim> >
phi_u(dofs_per_cell);
958 *
std::vector<double>
phi_p(dofs_per_cell);
961 *
dof_handler.begin_active(), first_cell = dof_handler.begin_active(),
962 *
endc = dof_handler.end();
964 *
for (; cell !=
endc; ++cell)
984 *
unsigned int m_id = cell->material_id();
991 *
fe_values.reinit(cell);
1004 *
if (cell->face(f)->center()[0] == 0)
1023 *
if (cell != first_cell)
1025 *
std::ofstream
fout(
"gravity_field.txt", std::ios::app);
1031 *
std::ofstream
fout(
"gravity_field.txt");
1037 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
1043 *
double r_value = fe_values.quadrature_point(
q)[0];
1078 *
for (
unsigned int k = 0;
k < dofs_per_cell; ++
k)
1092 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1094 *
for (
unsigned int j = 0;
j <= i; ++
j)
1109 *
* fe_values.JxW(
q);
1125 *
fe.system_to_component_index(i).first;
1126 *
local_rhs(i) += (fe_values.shape_value(i,
q)
1133 *
* fe_values.JxW(
q);
1138 *
fe.system_to_component_index(i).first;
1139 *
local_rhs(i) += fe_values.shape_value(i,
q)
1159 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
1165 *
double r_value = fe_values.quadrature_point(
q)[0];
1197 *
for (
unsigned int k = 0;
k < dofs_per_cell; ++
k)
1211 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1213 *
for (
unsigned int j = 0;
j <= i; ++
j)
1228 *
* fe_values.JxW(
q);
1244 *
fe.system_to_component_index(i).first;
1245 *
local_rhs(i) += (fe_values.shape_value(i,
q)
1252 *
* fe_values.JxW(
q);
1257 *
fe.system_to_component_index(i).first;
1258 *
local_rhs(i) += fe_values.shape_value(i,
q)
1266 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1267 *
for (
unsigned int j = i + 1;
j < dofs_per_cell; ++
j)
1270 *
cell->get_dof_indices(local_dof_indices);
1272 *
local_dof_indices, system_matrix,
system_rhs);
1275 *
std::cout <<
" Computing preconditioner..." << std::endl << std::flush;
1288 * ======================
SOLVER ======================
1305 *
system_matrix.block(1, 0).vmult(
schur_rhs, tmp);
1312 *
* solution.block(1).size();
1319 *
std::cout <<
"\nMax iterations and tolerance are: " << n_iterations
1323 *
preconditioner.initialize(system_matrix.block(1, 1),
1327 *
system_matrix.block(1, 1), preconditioner);
1331 *
constraints.distribute(solution);
1334 *
std::cout <<
" " << solver_control.last_step()
1335 *
<<
" outer CG Schur complement iterations for pressure"
1340 *
system_matrix.block(0, 1).vmult(tmp, solution.block(1));
1344 *
A_inverse.vmult(solution.block(0), tmp);
1345 *
constraints.distribute(solution);
1352 * ======================
OUTPUT RESULTS ======================
1361 *
std::vector<DataComponentInterpretation::DataComponentInterpretation> data_component_interpretation(
1363 *
data_component_interpretation.push_back(
1367 *
data_out.attach_dof_handler(dof_handler);
1370 *
data_out.build_patches();
1377 *
<<
"_elastic_displacements" <<
".txt";
1386 *
std::ofstream output(
filename.str().c_str());
1387 *
data_out.write_gnuplot(output);
1447 *
<<
"_baseviscosities.txt";
1452 *
std::cout <<
"Running stress calculations for plasticity iteration "
1463 *
dof_handler.begin_active(),
endc = dof_handler.end();
1469 *
for (; cell !=
endc; ++cell)
1503 *
std::vector<Vector<double> > vector_values(0);
1511 *
std::vector<SymmetricTensor<2, dim>>
old_stress;
1513 *
std::vector<double>
cell_Gs;
1515 *
dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
1525 *
fe_values.reinit(cell);
1527 *
fe_values.get_function_values(solution,
velocities);
1545 *
for (
unsigned int i = 0; i < (dim+1); ++i)
1554 *
for (
unsigned int i = 0; i < (dim+1); ++i)
1569 * unsigned int mat_id = cell->material_id();
1570 * double local_G = rheology.get_G(mat_id);
1571 * cell_Gs.push_back(local_G);
1576 * tracks where failure occurred
1579 * std::vector<double> reduction_factor;
1580 * unsigned int total_fails = 0;
1581 * if (plastic_iteration == 0)
1582 * cell_viscosities.resize(0);
1585 * loop across all the cells to find and adjust eta of failing cells
1588 * for (unsigned int i = 0; i < triangulation.n_active_cells(); ++i)
1590 * double current_cell_viscosity = 0;
1594 * Fill viscosities vector, analytically if plastic_iteration == 0 and from previous viscosities for later iteration
1597 * if (plastic_iteration == 0)
1599 * double local_viscosity;
1600 * local_viscosity = rheology.get_eta(points_list[i][0], points_list[i][1]);
1601 * current_cell_viscosity = local_viscosity;
1602 * cell_viscosities.push_back(current_cell_viscosity);
1606 * current_cell_viscosity = cell_viscosities[i];
1610 * double cell_eta_ve = 2
1611 * / ((1 / current_cell_viscosity)
1613 * / system_parameters::current_time_interval));
1614 * double cell_chi_ve = 1
1617 * * system_parameters::current_time_interval
1618 * / current_cell_viscosity));
1622 * find local pressure
1625 * double cell_p = vector_values[i].operator()(2);
1628 * find stresses tensor
1629 * makes non-diagonalized local matrix A
1632 * double sigma13 = 0.5
1633 * * (gradient_values[i][0][1] + gradient_values[i][1][0]);
1635 * A << gradient_values[i][0][0] << 0 << sigma13 << endr
1636 * << 0 << vector_values[i].operator()(0) / points_list[i].operator()(0)<< 0 << endr
1637 * << sigma13 << 0 << gradient_values[i][1][1] << endr;
1639 * olddevstress << old_stress[i][0][0] << 0 << old_stress[i][0][1] << endr
1640 * << 0 << old_phiphi_stress[i] << 0 << endr
1641 * << old_stress[i][0][1] << 0 << old_stress[i][1][1] << endr;
1643 * P << cell_p << cell_p << cell_p;
1644 * mat Pmat = diagmat(P);
1646 * B = (cell_eta_ve * A + cell_chi_ve * olddevstress) - Pmat;
1650 * finds principal stresses
1655 * eig_sym(eigval, eigvec, B);
1656 * double sigma1 = -min(eigval);
1657 * double sigma3 = -max(eigval);
1661 * Writes text files for principal stresses, full stress tensor, base viscosities
1664 * std::ofstream fout_snew(stress_output.str().c_str(), std::ios::app);
1665 * fout_snew << " " << sigma1 << " " << sigma3 << "\n";
1666 * fout_snew.close();
1668 * std::ofstream fout_sfull(stresstensor_output.str().c_str(), std::ios::app);
1669 * fout_sfull << A(0,0) << " " << A(1,1) << " " << A(2,2) << " " << A(0,2) << "\n";
1670 * fout_sfull.close();
1672 * if (plastic_iteration == 0)
1674 * std::ofstream fout_baseeta(initial_eta_output.str().c_str(), std::ios::app);
1675 * fout_baseeta << points_list[i]<< " " << current_cell_viscosity << "\n";
1676 * fout_baseeta.close();
1681 * Finds adjusted effective viscosity
1684 * if (system_parameters::plasticity_on)
1686 * if (system_parameters::failure_criterion == 0) //Apply Byerlee's
rule
1702 * Text file
of all failure locations
1724 * std::cout <<
" ext ";
1734 * Text file
of all failure locations
1749 * std::cout <<
" comp ";
1760 * Text file
of all failure locations
1776 *
std::cout <<
"Specified failure criterion not found\n";
1789 *
std::cout <<
" Number of failing cells: " <<
total_fails <<
"\n";
1821 *
dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
1823 *
fe_values.reinit(cell);
1854 *
for (
unsigned int ii=0;
ii<dim; ++
ii)
1874 *
dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
1878 *
fe_values.reinit(cell);
1923 *
std::cout <<
" Smoothing viscosity field...\n";
1926 *
dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
1942 *
std::vector< TriaIterator< CellAccessor<dim> > >
neighbor_cells;
1967 * neighbor_cells.push_back(neighbor_cells[i]->neighbor(f));
1968 * neighbor_indices.push_back(test_cell_no);
1969 * cell_touched[test_cell_no] = true;
1971 * new_cells_found++;
1976 * if (new_cells_found == 0)
1978 * find_more_cells = false;
1981 * start_cell = neighbor_cells.size() - new_cells_found;
1984 * fe_values.reinit(cell);
1987 * Collect the viscosities at nearby quadrature points
1990 * for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
1992 * std::vector<double> nearby_etas_q;
1993 * for (unsigned int i = 0; i<neighbor_indices.size(); ++i)
1994 * for (unsigned int j=0; j<quadrature_formula.size(); ++j)
1996 * Point<dim> test_q;
1997 * for (unsigned int d=0; d<dim; ++d)
1998 * test_q(d) = quad_viscosities[neighbor_indices[i]][j][d];
1999 * double qq_distance = fe_values.quadrature_point(q).distance(test_q);
2000 * if (qq_distance < system_parameters::smoothing_radius)
2001 * nearby_etas_q.push_back(quad_viscosities[neighbor_indices[i]][j][dim]);
2005 * Write smoothed viscosities to quadrature_points_history; simple boxcar function is the smoothing kernel
2008 * double mean_eta = 0;
2009 * for (unsigned int l = 0; l<nearby_etas_q.size(); ++l)
2011 * mean_eta += nearby_etas_q[l];
2013 * mean_eta /= nearby_etas_q.size();
2014 * local_quadrature_points_history[q].new_eta = mean_eta;
2017 * std::cout << local_quadrature_points_history[q].new_eta << " ";
2028 * ====================== SAVE STRESS TENSOR AT QUADRATURE POINTS ======================
2035 * void StokesProblem<dim>::update_quadrature_point_history()
2037 * std::cout << " Updating stress field...";
2039 * FEValues<dim> fe_values(fe, quadrature_formula,
2040 * update_values | update_gradients | update_quadrature_points);
2044 * Make the object that will hold the velocity gradients
2047 * std::vector < std::vector<Tensor<1, dim> >> velocity_grads(quadrature_formula.size(),
2048 * std::vector < Tensor<1, dim> > (dim + 1));
2049 * std::vector<Vector<double> > velocities(quadrature_formula.size(),
2050 * Vector<double>(dim + 1));
2052 * for (typename DoFHandler<dim>::active_cell_iterator cell =
2053 * dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
2055 * PointHistory<dim> *local_quadrature_points_history =
2056 * reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
2058 * local_quadrature_points_history >= &quadrature_point_history.front(),
2059 * ExcInternalError());
2061 * local_quadrature_points_history < &quadrature_point_history.back(),
2062 * ExcInternalError());
2064 * fe_values.reinit(cell);
2065 * fe_values.get_function_gradients(solution, velocity_grads);
2066 * fe_values.get_function_values(solution, velocities);
2068 * for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
2072 * Define the local viscoelastic constants
2075 * double local_eta_ve = 2
2076 * / ((1 / local_quadrature_points_history[q].new_eta)
2077 * + (1 / local_quadrature_points_history[q].G
2078 * / system_parameters::current_time_interval));
2079 * double local_chi_ve =
2082 * + (local_quadrature_points_history[q].G
2083 * * system_parameters::current_time_interval
2084 * / local_quadrature_points_history[q].new_eta));
2088 * Compute new stress at each quadrature point
2091 * SymmetricTensor<2, dim> new_stress;
2092 * for (unsigned int i = 0; i < dim; ++i)
2093 * new_stress[i][i] =
2094 * local_eta_ve * velocity_grads[q][i][i]
2096 * * local_quadrature_points_history[q].old_stress[i][i];
2098 * for (unsigned int i = 0; i < dim; ++i)
2099 * for (unsigned int j = i + 1; j < dim; ++j)
2100 * new_stress[i][j] =
2102 * * (velocity_grads[q][i][j]
2103 * + velocity_grads[q][j][i]) / 2
2105 * * local_quadrature_points_history[q].old_stress[i][j];
2112 * AuxFunctions<dim> rotation_object;
2113 * const Tensor<2, dim> rotation = rotation_object.get_rotation_matrix(
2114 * velocity_grads[q]);
2115 * const SymmetricTensor<2, dim> rotated_new_stress = symmetrize(
2116 * transpose(rotation)
2117 * * static_cast<Tensor<2, dim> >(new_stress)
2119 * local_quadrature_points_history[q].old_stress = rotated_new_stress;
2123 * For axisymmetric case, make the phi-phi element of stress tensor
2126 * local_quadrature_points_history[q].old_phiphi_stress =
2127 * (2 * local_eta_ve * velocities[q](0)
2128 * / fe_values.quadrature_point(q)[0]
2130 * * local_quadrature_points_history[q].old_phiphi_stress);
2137 * ====================== REDEFINE THE TIME INTERVAL FOR THE VISCOUS STEPS ======================
2141 * void StokesProblem<dim>::update_time_interval()
2143 * double move_goal_per_step = system_parameters::initial_disp_target;
2144 * if (system_parameters::present_timestep > system_parameters::initial_elastic_iterations)
2146 * move_goal_per_step = system_parameters::initial_disp_target -
2147 * ((system_parameters::initial_disp_target - system_parameters::final_disp_target) /
2148 * system_parameters::total_viscous_steps *
2149 * (system_parameters::present_timestep - system_parameters::initial_elastic_iterations));
2152 * double zero_tolerance = 1e-3;
2153 * double max_velocity = 0;
2154 * for (typename DoFHandler<dim>::active_cell_iterator cell =
2155 * dof_handler.begin_active(); cell != dof_handler.end(); ++cell)// loop over all cells
2157 * if (cell->at_boundary())
2159 * int zero_faces = 0;
2160 * for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
2161 * for (unsigned int i=0; i<dim; ++i)
2162 * if (fabs(cell->face(f)->center()[i]) < zero_tolerance)
2164 * if (zero_faces==0)
2166 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
2168 * Point<dim> vertex_velocity;
2169 * Point<dim> vertex_position;
2170 * for (unsigned int d = 0; d < dim; ++d)
2172 * vertex_velocity[d] = solution(cell->vertex_dof_index(v, d));
2173 * vertex_position[d] = cell->vertex(v)[d];
2177 * velocity to be evaluated is the radial component of a surface vertex
2180 * double local_velocity = 0;
2181 * for (unsigned int d = 0; d < dim; ++d)
2183 * local_velocity += vertex_velocity[d] * vertex_position [d];
2185 * local_velocity /= std::sqrt( vertex_position.square() );
2186 * if (local_velocity < 0)
2187 * local_velocity *= -1;
2188 * if (local_velocity > max_velocity)
2190 * max_velocity = local_velocity;
2198 * NOTE: It is possible for this time interval to be very different from that used in the viscoelasticity calculation.
2201 * system_parameters::current_time_interval = move_goal_per_step / max_velocity;
2202 * double step_time_yr = system_parameters::current_time_interval / SECSINYEAR;
2203 * std::cout << "Timestep interval changed to: "
2210 * ====================== MOVE MESH ======================
2217 * void StokesProblem<dim>::move_mesh()
2220 * std::cout << "\n" << " Moving mesh...\n";
2221 * std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
2222 * for (typename DoFHandler<dim>::active_cell_iterator cell =
2223 * dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
2224 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
2225 * if (vertex_touched[cell->vertex_index(v)] == false)
2227 * vertex_touched[cell->vertex_index(v)] = true;
2229 * Point<dim> vertex_displacement;
2230 * for (unsigned int d = 0; d < dim; ++d)
2231 * vertex_displacement[d] = solution(
2232 * cell->vertex_dof_index(v, d));
2233 * cell->vertex(v) += vertex_displacement
2234 * * system_parameters::current_time_interval;
2240 * ====================== WRITE MESH TO FILE ======================
2247 * void StokesProblem<dim>::write_mesh()
2251 * output mesh in ucd
2254 * std::ostringstream initial_mesh_file;
2255 * initial_mesh_file << system_parameters::output_folder << "/time" <<
2256 * Utilities::int_to_string(system_parameters::present_timestep, 2) <<
2258 * std::ofstream out_ucd (initial_mesh_file.str().c_str());
2260 * grid_out.write_ucd (triangulation, out_ucd);
2265 * ====================== FIT ELLIPSE TO SURFACE AND WRITE RADII TO FILE ======================
2272 * void StokesProblem<dim>::do_ellipse_fits()
2274 * std::ostringstream ellipses_filename;
2275 * ellipses_filename << system_parameters::output_folder << "/ellipse_fits.txt";
2278 * Find ellipsoidal axes for all layers
2281 * std::vector<double> ellipse_axes(0);
2284 * compute fit to boundary 0, 1, 2 ...
2287 * std::cout << endl;
2288 * for (unsigned int i = 0; i<system_parameters::sizeof_material_id; ++i)
2290 * ellipsoid.compute_fit(ellipse_axes, system_parameters::material_id[i]);
2291 * system_parameters::q_axes.push_back(ellipse_axes[0]);
2292 * system_parameters::p_axes.push_back(ellipse_axes[1]);
2294 * std::cout << "a_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[0]
2295 * << " " << " c_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[1] << std::endl;
2296 * ellipse_axes.clear();
2298 * std::ofstream fout_ellipses(ellipses_filename.str().c_str(), std::ios::app);
2299 * fout_ellipses << system_parameters::present_timestep << " a_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[0]
2300 * << " " << " c_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[1] << endl;
2301 * fout_ellipses.close();
2307 * ====================== APPEND LINE TO PHYSICAL_TIMES.TXT FILE WITH STEP NUMBER, PHYSICAL TIME, AND # PLASTIC ITERATIONS ======================
2314 * void StokesProblem<dim>::append_physical_times(int max_plastic)
2316 * std::ostringstream times_filename;
2317 * times_filename << system_parameters::output_folder << "/physical_times.txt";
2318 * std::ofstream fout_times(times_filename.str().c_str(), std::ios::app);
2319 * fout_times << system_parameters::present_timestep << " "
2320 * << system_parameters::present_time/SECSINYEAR << " "
2321 * << max_plastic << "\n";
2324 * << system_parameters::q_axes[0] << " " << system_parameters::p_axes[0] << " "
2325 * << system_parameters::q_axes[1] << " " << system_parameters::p_axes[1] << "\n";
2328 * fout_times.close();
2333 * ====================== WRITE VERTICES TO FILE ======================
2340 * void StokesProblem<dim>::write_vertices(unsigned char boundary_that_we_need)
2342 * std::ostringstream vertices_output;
2343 * vertices_output << system_parameters::output_folder << "/time" <<
2344 * Utilities::int_to_string(system_parameters::present_timestep, 2) << "_" <<
2345 * Utilities::int_to_string(boundary_that_we_need, 2) <<
2347 * std::ofstream fout_final_vertices(vertices_output.str().c_str());
2348 * fout_final_vertices.close();
2350 * std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
2352 * if (boundary_that_we_need == 0)
2356 * Figure out if the vertex is on the boundary of the domain
2359 * for (typename Triangulation<dim>::active_cell_iterator cell =
2360 * triangulation.begin_active(); cell != triangulation.end(); ++cell)
2361 * for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
2363 * unsigned char boundary_ids = cell->face(f)->boundary_id();
2364 * if (boundary_ids == boundary_that_we_need)
2366 * for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
2367 * if (vertex_touched[cell->face(f)->vertex_index(v)] == false)
2369 * vertex_touched[cell->face(f)->vertex_index(v)] = true;
2370 * std::ofstream fout_final_vertices(vertices_output.str().c_str(), std::ios::app);
2371 * fout_final_vertices << cell->face(f)->vertex(v) << "\n";
2372 * fout_final_vertices.close();
2381 * Figure out if the vertex is on an internal boundary
2384 * for (typename Triangulation<dim>::active_cell_iterator cell =
2385 * triangulation.begin_active(); cell != triangulation.end(); ++cell)
2386 * for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
2388 * if (cell->neighbor(f) != triangulation.end())
2390 * if (cell->material_id() != cell->neighbor(f)->material_id()) //finds face is at internal boundary
2392 * int high_mat_id = std::max(cell->material_id(),
2393 * cell->neighbor(f)->material_id());
2394 * if (high_mat_id == boundary_that_we_need) //finds faces at the correct internal boundary
2396 * for (unsigned int v = 0;
2397 * v < GeometryInfo<dim>::vertices_per_face;
2399 * if (vertex_touched[cell->face(f)->vertex_index(
2402 * vertex_touched[cell->face(f)->vertex_index(
2404 * std::ofstream fout_final_vertices(vertices_output.str().c_str(), std::ios::app);
2405 * fout_final_vertices << cell->face(f)->vertex(v) << "\n";
2406 * fout_final_vertices.close();
2417 * ====================== SETUP INITIAL MESH ======================
2424 * void StokesProblem<dim>::setup_initial_mesh()
2426 * GridIn<dim> grid_in;
2427 * grid_in.attach_triangulation(triangulation);
2428 * std::ifstream mesh_stream(system_parameters::mesh_filename,
2429 * std::ifstream::in);
2430 * grid_in.read_ucd(mesh_stream);
2434 * output initial mesh in eps
2437 * std::ostringstream initial_mesh_file;
2438 * initial_mesh_file << system_parameters::output_folder << "/initial_mesh.eps";
2439 * std::ofstream out_eps (initial_mesh_file.str().c_str());
2441 * grid_out.write_eps (triangulation, out_eps);
2447 * boundary indicator 0 is outer free surface; 1, 2, 3 ... is boundary between layers, 99 is flat boundaries
2450 * typename Triangulation<dim>::active_cell_iterator
2451 * cell=triangulation.begin_active(), endc=triangulation.end();
2453 * unsigned int how_many; // how many components away from cardinal planes
2455 * std::ostringstream boundaries_file;
2456 * boundaries_file << system_parameters::output_folder << "/boundaries.txt";
2457 * std::ofstream fout_boundaries(boundaries_file.str().c_str());
2458 * fout_boundaries.close();
2460 * double zero_tolerance = 1e-3;
2461 * for (; cell != endc; ++cell) // loop over all cells
2463 * for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f) // loop over all vertices
2465 * if (cell->face(f)->at_boundary())
2472 * std::ofstream fout_boundaries(boundaries_file.str().c_str(), std::ios::app);
2473 * fout_boundaries << cell->face(f)->center()[0] << " " << cell->face(f)->center()[1]<< "\n";
2474 * fout_boundaries.close();
2477 * for (unsigned int i=0; i<dim; ++i)
2478 * if (fabs(cell->face(f)->center()[i]) > zero_tolerance)
2480 * if (how_many==dim)
2481 * cell->face(f)->set_all_boundary_ids(0); // if face center coordinates > zero_tol, set bnry indicators to 0
2483 * cell->face(f)->set_all_boundary_ids(99);
2488 * std::ostringstream ellipses_filename;
2489 * ellipses_filename << system_parameters::output_folder << "/ellipse_fits.txt";
2490 * std::ofstream fout_ellipses(ellipses_filename.str().c_str());
2491 * fout_ellipses.close();
2495 * Find ellipsoidal axes for all layers
2498 * std::vector<double> ellipse_axes(0);
2501 * compute fit to boundary 0, 1, 2 ...
2504 * std::cout << endl;
2505 * for (unsigned int i = 0; i<system_parameters::sizeof_material_id; ++i)
2507 * ellipsoid.compute_fit(ellipse_axes, system_parameters::material_id[i]);
2508 * system_parameters::q_axes.push_back(ellipse_axes[0]);
2509 * system_parameters::p_axes.push_back(ellipse_axes[1]);
2511 * std::cout << "a_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[0]
2512 * << " " << " c_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[1] << std::endl;
2513 * ellipse_axes.clear();
2515 * std::ofstream fout_ellipses(ellipses_filename.str().c_str(), std::ios::app);
2516 * fout_ellipses << system_parameters::present_timestep << " a_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[0]
2517 * << " " << " c_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[1] << endl;
2518 * fout_ellipses.close();
2521 * triangulation.refine_global(system_parameters::global_refinement);
2526 * refines crustal region
2529 * if (system_parameters::crustal_refinement != 0)
2531 * double a = system_parameters::q_axes[0] - system_parameters::crust_refine_region;
2532 * double b = system_parameters::p_axes[0] - system_parameters::crust_refine_region;
2535 * for (unsigned int step = 0;
2536 * step < system_parameters::crustal_refinement; ++step)
2538 * typename ::Triangulation<dim>::active_cell_iterator cell =
2539 * triangulation.begin_active(), endc = triangulation.end();
2540 * for (; cell != endc; ++cell)
2541 * for (unsigned int v = 0;
2542 * v < GeometryInfo<dim>::vertices_per_cell; ++v)
2544 * Point<dim> current_vertex = cell->vertex(v);
2546 * const double x_coord = current_vertex.operator()(0);
2547 * const double y_coord = current_vertex.operator()(1);
2548 * double expected_z = -1;
2550 * if ((x_coord - a) < -1e-10)
2552 * * std::sqrt(1 - (x_coord * x_coord / a / a));
2554 * if (y_coord >= expected_z)
2556 * cell->set_refine_flag();
2560 * triangulation.execute_coarsening_and_refinement();
2567 * output initial mesh in eps
2570 * std::ostringstream refined_mesh_file;
2571 * refined_mesh_file << system_parameters::output_folder << "/refined_mesh.eps";
2572 * std::ofstream out_eps_refined (refined_mesh_file.str().c_str());
2573 * GridOut grid_out_refined;
2574 * grid_out_refined.write_eps (triangulation, out_eps_refined);
2575 * out_eps_refined.close();
2576 * write_vertices(0);
2577 * write_vertices(1);
2583 * ====================== REFINE MESH ======================
2590 * void StokesProblem<dim>::refine_mesh()
2592 * using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
2593 * Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
2595 * std::vector<bool> component_mask(dim + 1, false);
2596 * component_mask[dim] = true;
2597 * KellyErrorEstimator<dim>::estimate(dof_handler, QGauss<dim - 1>(degree + 1),
2598 * FunctionMap(), solution,
2599 * estimated_error_per_cell, component_mask);
2601 * GridRefinement::refine_and_coarsen_fixed_number(triangulation,
2602 * estimated_error_per_cell, 0.3, 0.0);
2603 * triangulation.execute_coarsening_and_refinement();
2608 * ====================== SET UP THE DATA STRUCTURES TO REMEMBER STRESS FIELD ======================
2612 * void StokesProblem<dim>::setup_quadrature_point_history()
2614 * unsigned int our_cells = 0;
2615 * for (typename Triangulation<dim>::active_cell_iterator cell =
2616 * triangulation.begin_active(); cell != triangulation.end(); ++cell)
2619 * triangulation.clear_user_data();
2621 * quadrature_point_history.resize(our_cells * quadrature_formula.size());
2623 * unsigned int history_index = 0;
2624 * for (typename Triangulation<dim>::active_cell_iterator cell =
2625 * triangulation.begin_active(); cell != triangulation.end(); ++cell)
2627 * cell->set_user_pointer(&quadrature_point_history[history_index]);
2628 * history_index += quadrature_formula.size();
2631 * Assert(history_index == quadrature_point_history.size(), ExcInternalError());
2636 * ====================== DOES ELASTIC STEPS ======================
2640 * void StokesProblem<dim>::do_elastic_steps()
2642 * unsigned int elastic_iteration = 0;
2644 * while (elastic_iteration < system_parameters::initial_elastic_iterations)
2647 * std::cout << "\n\nElastic iteration " << elastic_iteration
2651 * if (system_parameters::present_timestep == 0)
2652 * initialize_eta_and_G();
2654 * if (elastic_iteration == 0)
2655 * system_parameters::current_time_interval =
2656 * system_parameters::viscous_time; //This is the time interval needed in assembling the problem
2658 * std::cout << " Assembling..." << std::endl << std::flush;
2659 * assemble_system();
2661 * std::cout << " Solving..." << std::flush;
2665 * update_quadrature_point_history();
2667 * append_physical_times(0);
2668 * elastic_iteration++;
2669 * system_parameters::present_timestep++;
2670 * do_ellipse_fits();
2671 * write_vertices(0);
2672 * write_vertices(1);
2674 * update_time_interval();
2680 * ====================== DO A SINGLE VISCOELASTOPLASTIC TIMESTEP ======================
2684 * void StokesProblem<dim>::do_flow_step()
2686 * plastic_iteration = 0;
2687 * while (plastic_iteration < system_parameters::max_plastic_iterations)
2689 * if (system_parameters::continue_plastic_iterations == true)
2691 * std::cout << "Plasticity iteration " << plastic_iteration << "\n";
2694 * std::cout << " Assembling..." << std::endl << std::flush;
2695 * assemble_system();
2697 * std::cout << " Solving..." << std::flush;
2701 * solution_stesses();
2703 * if (system_parameters::continue_plastic_iterations == false)
2706 * plastic_iteration++;
2713 * ====================== RUN ======================
2720 * void StokesProblem<dim>::run()
2724 * Sets up mesh and data structure for viscosity and stress at quadrature points
2727 * setup_initial_mesh();
2728 * setup_quadrature_point_history();
2732 * Makes the physical_times.txt file
2735 * std::ostringstream times_filename;
2736 * times_filename << system_parameters::output_folder << "/physical_times.txt";
2737 * std::ofstream fout_times(times_filename.str().c_str());
2738 * fout_times.close();
2742 * Computes elastic timesteps
2745 * do_elastic_steps();
2748 * Computes viscous timesteps
2751 * unsigned int VEPstep = 0;
2752 * while (system_parameters::present_timestep
2753 * < (system_parameters::initial_elastic_iterations
2754 * + system_parameters::total_viscous_steps))
2757 * if (system_parameters::continue_plastic_iterations == false)
2758 * system_parameters::continue_plastic_iterations = true;
2759 * std::cout << "\n\nViscoelastoplastic iteration " << VEPstep << "\n\n";
2762 * Computes plasticity
2766 * update_quadrature_point_history();
2768 * append_physical_times(plastic_iteration);
2769 * system_parameters::present_timestep++;
2770 * system_parameters::present_time = system_parameters::present_time + system_parameters::current_time_interval;
2771 * do_ellipse_fits();
2772 * write_vertices(0);
2773 * write_vertices(1);
2777 * append_physical_times(-1);
2784 * ====================== MAIN ======================
2790 * int main(int argc, char *argv[])
2795 * output program name
2798 * std::cout << "Running: " << argv[0] << std::endl;
2800 * char *cfg_filename = new char[120];
2802 * if (argc == 1) // if no input parameters (as if launched from eclipse)
2804 * std::strcpy(cfg_filename,"config/sample_CeresFE_config.cfg");
2807 * std::strcpy(cfg_filename,argv[1]);
2811 * using namespace dealii;
2812 * using namespace Step22;
2813 * config_in cfg(cfg_filename);
2817 * t1 = std::clock();
2819 * deallog.depth_console(0);
2821 * StokesProblem<2> flow_problem(1);
2822 * flow_problem.run();
2824 * std::cout << std::endl << "\a";
2826 * t2 = std::clock();
2827 * float diff (((float)t2 - (float)t1) / (float)CLOCKS_PER_SEC);
2828 * std::cout << "\n Program run in: " << diff << " seconds" << endl;
2830 * catch (std::exception &exc)
2832 * std::cerr << std::endl << std::endl
2833 * << "----------------------------------------------------"
2835 * std::cerr << "Exception on processing: " << std::endl << exc.what()
2836 * << std::endl << "Aborting!" << std::endl
2837 * << "----------------------------------------------------"
2844 * std::cerr << std::endl << std::endl
2845 * << "----------------------------------------------------"
2847 * std::cerr << "Unknown exception!" << std::endl << "Aborting!"
2849 * << "----------------------------------------------------"
2859<a name="ann-support_code/config_in.h"></a>
2860<h1>Annotated version of support_code/config_in.h</h1>
2866 * /* -----------------------------------------------------------------------------
2868 * * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
2869 * * Copyright (C) 2015 by Anton Ermakov
2871 * * This file is part of the deal.II code gallery.
2873 * * -----------------------------------------------------------------------------
2875 * * Author: antonermakov
2879 * #include <iostream>
2880 * #include <iomanip>
2881 * #include <cstdlib>
2883 * #include <libconfig.h++>
2885 * #include "local_math.h"
2887 * using namespace std;
2888 * using namespace libconfig;
2892 * namespace system_parameters
2900 * string mesh_filename;
2901 * string output_folder;
2910 * double omegasquared;
2916 * Rheology parameters
2919 * vector<double> depths_eta;
2920 * vector<double> eta_kinks;
2921 * vector<double> depths_rho;
2922 * vector<double> rho;
2923 * vector<int> material_id;
2926 * double eta_ceiling;
2929 * bool lat_dependence;
2931 * unsigned int sizeof_depths_eta;
2932 * unsigned int sizeof_depths_rho;
2933 * unsigned int sizeof_rho;
2934 * unsigned int sizeof_eta_kinks;
2935 * unsigned int sizeof_material_id;
2936 * unsigned int sizeof_G;
2938 * double pressure_scale;
2941 * bool continue_plastic_iterations;
2945 * plasticity variables
2948 * bool plasticity_on;
2949 * unsigned int failure_criterion;
2950 * unsigned int max_plastic_iterations;
2951 * double smoothing_radius;
2955 * viscoelasticity variables
2958 * unsigned int initial_elastic_iterations;
2959 * double elastic_time;
2960 * double viscous_time;
2961 * double initial_disp_target;
2962 * double final_disp_target;
2963 * double current_time_interval;
2967 * mesh refinement variables
2970 * unsigned int global_refinement;
2971 * unsigned int small_r_refinement;
2972 * unsigned int crustal_refinement;
2973 * double crust_refine_region;
2974 * unsigned int surface_refinement;
2981 * int iteration_coefficient;
2982 * double tolerance_coefficient;
2986 * time step variables
2989 * double present_time;
2990 * unsigned int present_timestep;
2991 * unsigned int total_viscous_steps;
2999 * vector<double> q_axes;
3000 * vector<double> p_axes;
3007 * config_in(char *);
3010 * void write_config();
3013 * void config_in::write_config()
3015 * std::ostringstream config_parameters;
3016 * config_parameters << system_parameters::output_folder << "/run_parameters.txt";
3017 * std::ofstream fout_config(config_parameters.str().c_str());
3024 * fout_config << "mesh filename: " << system_parameters::mesh_filename << endl << endl;
3031 * fout_config << "r_mean = " << system_parameters::r_mean << endl;
3032 * fout_config << "period = " << system_parameters::period << endl;
3033 * fout_config << "omegasquared = " << system_parameters::omegasquared << endl;
3034 * fout_config << "beta = " << system_parameters::beta << endl;
3035 * fout_config << "intercept = " << system_parameters::intercept << endl;
3039 * rheology parameters
3045 * for (unsigned int i=0; i<system_parameters::sizeof_depths_eta; ++i)
3046 * fout_config << "depths_eta[" << i << "] = " << system_parameters::depths_eta[i] << endl;
3048 * for (unsigned int i=0; i<system_parameters::sizeof_eta_kinks; ++i)
3049 * fout_config << "eta_kinks[" << i << "] = " << system_parameters::eta_kinks[i] << endl;
3051 * for (unsigned int i=0; i<system_parameters::sizeof_depths_rho; ++i)
3052 * fout_config << "depths_rho[" << i << "] = " << system_parameters::depths_rho[i] << endl;
3054 * for (unsigned int i=0; i<system_parameters::sizeof_rho; ++i)
3055 * fout_config << "rho[" << i << "] = " << system_parameters::rho[i] << endl;
3057 * for (unsigned int i=0; i<system_parameters::sizeof_material_id; ++i)
3058 * fout_config << "material_id[" << i << "] = " << system_parameters::material_id[i] << endl;
3060 * for (unsigned int i=0; i<system_parameters::sizeof_G; ++i)
3061 * fout_config << "G[" << i << "] = " << system_parameters::G[i] << endl;
3063 * fout_config << "eta_ceiling = " << system_parameters::eta_ceiling << endl;
3064 * fout_config << "eta_floor = " << system_parameters::eta_floor << endl;
3065 * fout_config << "eta_Ea = " << system_parameters::eta_Ea << endl;
3066 * fout_config << "lat_dependence = " << system_parameters::lat_dependence << endl;
3067 * fout_config << "pressure_scale = " << system_parameters::pressure_scale << endl;
3068 * fout_config << "q = " << system_parameters::q << endl;
3069 * fout_config << "cylindrical = " << system_parameters::cylindrical << endl;
3070 * fout_config << "continue_plastic_iterations = " << system_parameters::continue_plastic_iterations << endl;
3074 * Plasticity parameters
3077 * fout_config << "plasticity_on = " << system_parameters::plasticity_on << endl;
3078 * fout_config << "failure_criterion = " << system_parameters::failure_criterion << endl;
3079 * fout_config << "max_plastic_iterations = " << system_parameters::max_plastic_iterations << endl;
3080 * fout_config << "smoothing_radius = " << system_parameters::smoothing_radius << endl;
3084 * Viscoelasticity parameters
3087 * fout_config << "initial_elastic_iterations = " << system_parameters::initial_elastic_iterations << endl;
3088 * fout_config << "elastic_time = " << system_parameters::elastic_time << endl;
3089 * fout_config << "viscous_time = " << system_parameters::viscous_time << endl;
3090 * fout_config << "initial_disp_target = " << system_parameters::initial_disp_target << endl;
3091 * fout_config << "final_disp_target = " << system_parameters::final_disp_target << endl;
3092 * fout_config << "current_time_interval = " << system_parameters::current_time_interval << endl;
3096 * Mesh refinement parameters
3099 * fout_config << "global_refinement = " << system_parameters::global_refinement << endl;
3100 * fout_config << "small_r_refinement = " << system_parameters::small_r_refinement << endl;
3101 * fout_config << "crustal_refinement = " << system_parameters::crustal_refinement << endl;
3102 * fout_config << "crust_refine_region = " << system_parameters::crust_refine_region << endl;
3103 * fout_config << "surface_refinement = " << system_parameters::surface_refinement << endl;
3110 * fout_config << "iteration_coefficient = " << system_parameters::iteration_coefficient << endl;
3111 * fout_config << "tolerance_coefficient = " << system_parameters::tolerance_coefficient << endl;
3115 * Time step parameters
3118 * fout_config << "present_time = " << system_parameters::present_time << endl;
3119 * fout_config << "present_timestep = " << system_parameters::present_timestep << endl;
3120 * fout_config << "total_viscous_steps = " << system_parameters::total_viscous_steps << endl;
3122 * fout_config.close();
3125 * config_in::config_in(char *filename)
3130 * This example reads the configuration file 'example.cfg
' and displays
3131 * some of its contents.
3141 * Read the file. If there is an error, report it and exit.
3146 * cfg.readFile(filename);
3148 * catch (const FileIOException &fioex)
3150 * std::cerr << "I/O error while reading file:" << filename << std::endl;
3152 * catch (const ParseException &pex)
3154 * std::cerr << "Parse error at " << pex.getFile() << ":" << pex.getLine()
3155 * << " - " << pex.getError() << std::endl;
3165 * string msh = cfg.lookup("mesh_filename");
3166 * system_parameters::mesh_filename = msh;
3168 * catch (const SettingNotFoundException &nfex)
3170 * cerr << "No 'mesh_filename' setting in configuration file." << endl;
3183 * string output = cfg.lookup("output_folder");
3184 * system_parameters::output_folder = output;
3186 * std::cout << "Writing to folder: " << output << endl;
3188 * catch (const SettingNotFoundException &nfex)
3190 * cerr << "No 'output_folder' setting in configuration file." << endl;
3201 * const Setting &root = cfg.getRoot();
3206 * get body parameters
3211 * const Setting &body_parameters = root["body_parameters"];
3213 * body_parameters.lookupValue("period", system_parameters::period);
3214 * system_parameters::omegasquared = pow(TWOPI / 3600.0 / system_parameters::period, 2.0);
3215 * body_parameters.lookupValue("r_mean", system_parameters::r_mean);
3216 * body_parameters.lookupValue("beta", system_parameters::beta);
3217 * body_parameters.lookupValue("intercept", system_parameters::intercept);
3219 * catch (const SettingNotFoundException &nfex)
3227 * Rheology parameters
3237 * get depths_eta ---------------------
3242 * unsigned int ndepths_eta = set_depths_eta.getLength();
3243 * system_parameters::sizeof_depths_eta = ndepths_eta;
3245 * for (unsigned int i=0; i<ndepths_eta; ++i)
3247 * system_parameters::depths_eta.push_back(set_depths_eta[i]);
3248 * cout << "depth_eta[
" << i << "] =
" << system_parameters::depths_eta[i] << endl;
3253 * get eta_kinks -------------------------
3258 * unsigned int neta_kinks = set_eta_kinks.getLength();
3259 * system_parameters::sizeof_eta_kinks = neta_kinks;
3263 * cout << "Number
of depth =
" << ndepths << endl;
3269 * for (unsigned int i=0; i<neta_kinks; ++i)
3271 * system_parameters::eta_kinks.push_back(set_eta_kinks[i]);
3272 * cout << "eta_kinks[
" << i << "] =
" << system_parameters::eta_kinks[i] << endl;
3277 * get depths_rho -------------------------
3282 * unsigned int ndepths_rho = set_depths_rho.getLength();
3283 * system_parameters::sizeof_depths_rho = ndepths_rho;
3287 * cout << "Number
of depth =
" << ndepths << endl;
3293 * for (unsigned int i=0; i<ndepths_rho; ++i)
3295 * system_parameters::depths_rho.push_back(set_depths_rho[i]);
3296 * cout << "depths_rho[
" << i << "] =
" << system_parameters::depths_rho[i] << endl;
3301 * get rho -------------------------
3306 * unsigned int nrho = set_rho.getLength();
3307 * system_parameters::sizeof_rho = nrho;
3311 * cout << "Number
of depth =
" << ndepths << endl;
3317 * for (unsigned int i=0; i<nrho; ++i)
3319 * system_parameters::rho.push_back(set_rho[i]);
3320 * cout << "rho[
" << i << "] =
" << system_parameters::rho[i] << endl;
3325 * get material_id -------------------------
3330 * unsigned int nmaterial_id = set_material_id.getLength();
3331 * system_parameters::sizeof_material_id = nmaterial_id;
3335 * cout << "Number
of depth =
" << ndepths << endl;
3341 * for (unsigned int i=0; i<nmaterial_id; ++i)
3343 * system_parameters::material_id.push_back(set_material_id[i]);
3344 * cout << "material_id[
" << i << "] =
" << system_parameters::material_id[i] << endl;
3349 * get G -------------------------
3354 * unsigned int nG = set_G.getLength();
3355 * system_parameters::sizeof_G = nG;
3359 * cout << "Number
of depth =
" << ndepths << endl;
3365 * for (unsigned int i=0; i<nG; ++i)
3367 * system_parameters::G.push_back(set_G[i]);
3368 * cout << "G[
" << i << "] =
" << system_parameters::G[i] << endl;
3372 * rheology_parameters.lookupValue("eta_ceiling", system_parameters::eta_ceiling);
3373 * rheology_parameters.lookupValue("eta_floor", system_parameters::eta_floor);
3374 * rheology_parameters.lookupValue("eta_Ea", system_parameters::eta_Ea);
3375 * rheology_parameters.lookupValue("lat_dependence", system_parameters::lat_dependence);
3376 * rheology_parameters.lookupValue("pressure_scale", system_parameters::pressure_scale);
3377 * rheology_parameters.lookupValue("q", system_parameters::q);
3378 * rheology_parameters.lookupValue("cylindrical", system_parameters::cylindrical);
3381 * catch (const SettingNotFoundException &nfex)
3383 * cerr << "We've got a problem in the rheology parameters block" << endl;
3388 * Plasticity parameters
3394 * const Setting &plasticity_parameters = root["plasticity_parameters"];
3395 * plasticity_parameters.lookupValue("plasticity_on", system_parameters::plasticity_on);
3396 * plasticity_parameters.lookupValue("failure_criterion", system_parameters::failure_criterion);
3397 * plasticity_parameters.lookupValue("max_plastic_iterations", system_parameters::max_plastic_iterations);
3398 * plasticity_parameters.lookupValue("smoothing_radius", system_parameters::smoothing_radius);
3400 * catch (const SettingNotFoundException &nfex)
3407 * Viscoelasticity parameters
3418 * viscoelasticity_parameters.lookupValue("elastic_time", system_parameters::elastic_time);
3419 * viscoelasticity_parameters.lookupValue("viscous_time", system_parameters::viscous_time);
3420 * viscoelasticity_parameters.lookupValue("initial_disp_target", system_parameters::initial_disp_target);
3421 * viscoelasticity_parameters.lookupValue("final_disp_target", system_parameters::final_disp_target);
3422 * viscoelasticity_parameters.lookupValue("current_time_interval", system_parameters::current_time_interval);
3424 * system_parameters::viscous_time *= SECSINYEAR;
3426 * catch (const SettingNotFoundException &nfex)
3428 * cerr << "We've got a problem in the viscoelasticity parameters block" << endl;
3433 * Mesh refinement parameters
3439 * const Setting &mesh_refinement_parameters = root["mesh_refinement_parameters"];
3440 * mesh_refinement_parameters.lookupValue("global_refinement", system_parameters::global_refinement);
3441 * mesh_refinement_parameters.lookupValue("small_r_refinement", system_parameters::small_r_refinement);
3442 * mesh_refinement_parameters.lookupValue("crustal_refinement", system_parameters::crustal_refinement);
3443 * mesh_refinement_parameters.lookupValue("crust_refine_region", system_parameters::crust_refine_region);
3444 * mesh_refinement_parameters.lookupValue("surface_refinement", system_parameters::surface_refinement);
3446 * catch (const SettingNotFoundException &nfex)
3464 * catch (const SettingNotFoundException &nfex)
3466 * cerr << "We've got a problem in the solver parameters block" << endl;
3471 * Time step parameters
3476 * const Setting &time_step_parameters = root["time_step_parameters"];
3477 * time_step_parameters.lookupValue("present_time", system_parameters::present_time);
3478 * time_step_parameters.lookupValue("present_timestep", system_parameters::present_timestep);
3479 * time_step_parameters.lookupValue("total_viscous_steps", system_parameters::total_viscous_steps);
3481 * catch (const SettingNotFoundException &nfex)
3483 * cerr << "We've got a
problem in
the time step parameters block
" << endl;
3498<h1>Annotated version of support_code/ellipsoid_fit.h</h1>
3504 * /* -----------------------------------------------------------------------------
3506 * * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
3507 * * Copyright (C) 2015 by Anton Ermakov
3509 * * This file is part of the deal.II code gallery.
3511 * * -----------------------------------------------------------------------------
3513 * * Author: antonermakov
3516 * #include <deal.II/base/quadrature_lib.h>
3517 * #include <deal.II/base/function.h>
3518 * #include <deal.II/base/logstream.h>
3519 * #include <deal.II/lac/vector.h>
3520 * #include <deal.II/lac/full_matrix.h>
3521 * #include <deal.II/lac/sparse_matrix.h>
3522 * #include <deal.II/lac/solver_cg.h>
3523 * #include <deal.II/lac/precondition.h>
3524 * #include <deal.II/grid/tria.h>
3525 * #include <deal.II/dofs/dof_handler.h>
3526 * #include <deal.II/grid/tria_accessor.h>
3527 * #include <deal.II/grid/tria_iterator.h>
3528 * #include <deal.II/dofs/dof_accessor.h>
3529 * #include <deal.II/dofs/dof_tools.h>
3530 * #include <deal.II/fe/fe_q.h>
3531 * #include <deal.II/fe/fe_values.h>
3532 * #include <deal.II/numerics/vector_tools.h>
3533 * #include <deal.II/numerics/matrix_tools.h>
3534 * #include <deal.II/numerics/data_out.h>
3535 * #include <deal.II/grid/grid_in.h>
3536 * #include <deal.II/grid/grid_out.h>
3537 * #include <deal.II/base/point.h>
3538 * #include <deal.II/grid/grid_generator.h>
3540 * #include <fstream>
3541 * #include <sstream>
3542 * #include <iostream>
3543 * #include <iomanip>
3544 * #include <cstdlib>
3549 * using namespace dealii;
3551 * template <int dim>
3552 * class ellipsoid_fit
3555 * inline ellipsoid_fit (Triangulation<dim,dim> *pi)
3557 * p_triangulation = pi;
3559 * void compute_fit(std::vector<double> &ell, unsigned char bndry);
3563 * Triangulation<dim,dim> *p_triangulation;
3570 * This function computes ellipsoid fit to a set of vertices that lie on the
3571 * boundary_that_we_need
3574 * template <int dim>
3575 * void ellipsoid_fit<dim>::compute_fit(std::vector<double> &ell, unsigned char boundary_that_we_need)
3577 * typename Triangulation<dim>::active_cell_iterator cell = p_triangulation->begin_active();
3578 * typename Triangulation<dim>::active_cell_iterator endc = p_triangulation->end();
3580 * FullMatrix<double> A(p_triangulation->n_vertices(),dim);
3581 * Vector<double> x(dim);
3582 * Vector<double> b(p_triangulation->n_vertices());
3584 * std::vector<bool> vertex_touched (p_triangulation->n_vertices(),
3587 * unsigned int j = 0;
3588 * unsigned char boundary_ids;
3589 * std::vector<unsigned int> ind_bnry_row;
3590 * std::vector<unsigned int> ind_bnry_col;
3594 * assemble the sensitivity matrix and r.h.s.
3597 * for (; cell != endc; ++cell)
3599 * if (boundary_that_we_need != 0)
3600 * cell->set_manifold_id(cell->material_id());
3601 * for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
3603 * if (boundary_that_we_need == 0) //if this is the outer surface, then look for boundary ID 0; otherwise look for material ID change.
3605 * boundary_ids = cell->face(f)->boundary_id();
3606 * if (boundary_ids == boundary_that_we_need)
3608 * for (unsigned int v = 0;
3609 * v < GeometryInfo<dim>::vertices_per_face; ++v)
3610 * if (vertex_touched[cell->face(f)->vertex_index(v)]
3613 * vertex_touched[cell->face(f)->vertex_index(v)] =
3615 * for (unsigned int i = 0; i < dim; ++i)
3619 * stiffness matrix entry
3622 * A(j, i) = pow(cell->face(f)->vertex(v)[i], 2);
3631 * if mesh if not full: set the indicator
3635 * ind_bnry_row.push_back(j);
3640 * else //find the faces that are at the boundary between materials, get the vertices, and write them into the stiffness matrix
3642 * if (cell->neighbor(f) != endc)
3644 * if (cell->material_id() != cell->neighbor(f)->material_id()) //finds face is at internal boundary
3646 * int high_mat_id = std::max(cell->material_id(),
3647 * cell->neighbor(f)->material_id());
3648 * if (high_mat_id == boundary_that_we_need) //finds faces at the correct internal boundary
3650 * for (unsigned int v = 0;
3651 * v < GeometryInfo<dim>::vertices_per_face;
3653 * if (vertex_touched[cell->face(f)->vertex_index(
3656 * vertex_touched[cell->face(f)->vertex_index(
3658 * for (unsigned int i = 0; i < dim; ++i)
3662 * stiffness matrix entry
3666 * cell->face(f)->vertex(v)[i], 2);
3675 * if mesh if not full: set the indicator
3679 * ind_bnry_row.push_back(j);
3688 * if (ind_bnry_row.size()>0)
3693 * maxtrix A'*A and vector A'*b; A'*A*x = A'*b -- normal system of equations
3696 * FullMatrix<double> AtA(dim,dim);
3697 * Vector<double> Atb(dim);
3699 * FullMatrix<double> A_out(ind_bnry_row.size(),dim);
3700 * Vector<double> b_out(ind_bnry_row.size());
3702 * for (unsigned int i=0; i<dim; ++i)
3703 * ind_bnry_col.push_back(i);
3705 * for (unsigned int i=0; i<ind_bnry_row.size(); ++i)
3708 * A_out.extract_submatrix_from(A, ind_bnry_row, ind_bnry_col);
3709 * A_out.Tmmult(AtA,A_out,true);
3710 * A_out.Tvmult(Atb,b_out,true);
3714 * solve normal system of equations
3717 * SolverControl solver_control (1000, 1e-12);
3718 * SolverCG<> solver (solver_control);
3719 * solver.solve (AtA, x, Atb, PreconditionIdentity());
3723 * find ellipsoidal axes
3726 * for (unsigned int i=0; i<dim; ++i)
3727 * ell.push_back(sqrt(1.0/x[i]));
3738<h1>Annotated version of support_code/ellipsoid_grav.h</h1>
3744 * /* -----------------------------------------------------------------------------
3746 * * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
3747 * * Copyright (C) 2012 by Roger R. Fu
3749 * * This file is part of the deal.II code gallery.
3751 * * -----------------------------------------------------------------------------
3753 * * Author: Roger R. Fu
3754 * * Reference Pohanka 2011, Contrib. Geophys. Geodes.
3758 * #include <deal.II/base/point.h>
3759 * #include <fstream>
3760 * #include <iostream>
3762 * namespace A_Grav_namespace
3764 * namespace system_parameters
3766 * double mantle_rho;
3768 * double excess_rho;
3773 * double r_core_polar;
3776 * template <int dim>
3777 * class AnalyticGravity
3780 * void setup_vars (std::vector<double> v);
3781 * void get_gravity (const ::Point<dim> &p, std::vector<double> &g);
3800 * template <int dim>
3801 * void AnalyticGravity<dim>::get_gravity (const ::Point<dim> &p, std::vector<double> &g)
3803 * double rsph = std::sqrt(p[0] * p[0] + p[1] * p[1]);
3804 * double thetasph = std::atan2(p[0], p[1]);
3805 * double costhetasph = std::cos(thetasph);
3809 * convert to elliptical coordinates for silicates
3812 * double stemp = std::sqrt((rsph * rsph - eV * eV + std::sqrt((rsph * rsph - eV * eV) * (rsph * rsph - eV * eV)
3813 * + 4 * eV * eV * rsph * rsph * costhetasph *costhetasph)) / 2);
3814 * double vout = stemp / system_parameters::r_eq / std::sqrt(1 - ecc * ecc);
3815 * double eout = std::acos(rsph * costhetasph / stemp);
3819 * convert to elliptical coordinates for core correction
3822 * double stemp_c = std::sqrt((rsph * rsph - eV_c * eV_c + std::sqrt((rsph * rsph - eV_c * eV_c) * (rsph * rsph - eV_c * eV_c)
3823 * + 4 * eV_c * eV_c * rsph * rsph * costhetasph *costhetasph)) / 2);
3824 * double vout_c = stemp_c / system_parameters::r_core_eq / std::sqrt(1 - ecc_c * ecc_c);
3825 * double eout_c = std::acos(rsph * costhetasph / stemp_c);
3829 * shell contribution
3832 * g[0] = g_coeff * r11 * std::sqrt((1 - ecc * ecc) * vout * vout + ecc * ecc) * std::sin(eout);
3833 * g[1] = g_coeff * r01 * vout * std::cos(eout) / std::sqrt(1 - ecc * ecc);
3840 * double expected_y = system_parameters::r_core_polar * std::sqrt(1 -
3841 * (p[0] * p[0] / system_parameters::r_core_eq / system_parameters::r_core_eq));
3844 * if (p[1] <= expected_y)
3846 * g[0] += g_coeff_c * r11_c * std::sqrt((1 - ecc_c * ecc_c) * vout_c * vout_c + ecc_c * ecc_c) * std::sin(eout_c);
3847 * g[1] += g_coeff_c * r01_c * vout_c * std::cos(eout_c) / std::sqrt(1 - ecc_c * ecc_c);
3851 * double g_coeff_co = - 2.795007963255562e-10 * system_parameters::excess_rho * system_parameters::r_core_eq
3852 * / vout_c / vout_c;
3853 * double r00_co = 0;
3854 * double r01_co = 0;
3855 * double r11_co = 0;
3857 * if (system_parameters::r_core_polar == system_parameters::r_core_eq)
3865 * r00_co = ke_c * vout_c * std::atan2(1, ke_c * vout_c);
3866 * double ke_co2 = ke_c * ke_c * vout_c * vout_c;
3867 * r01_co = 3 * ke_co2 * (1 - r00_co);
3868 * r11_co = 3 * ((ke_co2 + 1) * r00_co - ke_co2) / 2;
3870 * g[0] += g_coeff_co * vout_c * r11_co / std::sqrt((1 - ecc_c* ecc_c) * vout_c * vout_c + ecc_c * ecc_c) * std::sin(eout_c);
3871 * g[1] += g_coeff_co * r01_co * std::cos(eout_c) / std::sqrt(1 - ecc_c * ecc_c);
3875 * template <int dim>
3876 * void AnalyticGravity<dim>::setup_vars (std::vector<double> v)
3878 * system_parameters::r_eq = v[0];
3879 * system_parameters::r_polar = v[1];
3880 * system_parameters::r_core_eq = v[2];
3881 * system_parameters::r_core_polar = v[3];
3882 * system_parameters::mantle_rho = v[4];
3883 * system_parameters::core_rho = v[5];
3884 * system_parameters::excess_rho = system_parameters::core_rho - system_parameters::mantle_rho;
3892 * if (system_parameters::r_polar > system_parameters::r_eq)
3896 * This makes the gravity field nearly that of a sphere in case the body becomes prolate
3904 * ecc = std::sqrt(1 - (system_parameters::r_polar * system_parameters::r_polar / system_parameters::r_eq / system_parameters::r_eq));
3907 * eV = ecc * system_parameters::r_eq;
3908 * ke = std::sqrt(1 - (ecc * ecc)) / ecc;
3909 * r00 = ke * std::atan2(1, ke);
3910 * double ke2 = ke * ke;
3911 * r01 = 3 * ke2 * (1 - r00);
3912 * r11 = 3 * ((ke2 + 1) * r00 - ke2) / 2;
3913 * g_coeff = - 2.795007963255562e-10 * system_parameters::mantle_rho * system_parameters::r_eq;
3920 * if (system_parameters::r_core_polar > system_parameters::r_core_eq)
3927 * ecc_c = std::sqrt(1 - (system_parameters::r_core_polar * system_parameters::r_core_polar / system_parameters::r_core_eq / system_parameters::r_core_eq));
3929 * eV_c = ecc_c * system_parameters::r_core_eq;
3930 * if (system_parameters::r_core_polar == system_parameters::r_core_eq)
3936 * g_coeff_c = - 2.795007963255562e-10 * system_parameters::excess_rho * system_parameters::r_core_eq;
3940 * ke_c = std::sqrt(1 - (ecc_c * ecc_c)) / ecc_c;
3941 * r00_c = ke_c * std::atan2(1, ke_c);
3942 * double ke2_c = ke_c * ke_c;
3943 * r01_c = 3 * ke2_c * (1 - r00_c);
3944 * r11_c = 3 * ((ke2_c + 1) * r00_c - ke2_c) / 2;
3945 * g_coeff_c = - 2.795007963255562e-10 * system_parameters::excess_rho * system_parameters::r_core_eq;
3949 * std::cout << "Loaded variables:
ecc =
" << ecc_c << " ke =
" << ke_c << " r00 =
" << r00_c << " r01 =
" << r01_c << " r11 =
" << r11_c << "\n
";
3958<h1>Annotated version of support_code/local_math.h</h1>
3964 * /* -----------------------------------------------------------------------------
3966 * * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
3967 * * Copyright (C) 2013 by Anton Ermakov
3969 * * This file is part of the deal.II code gallery.
3971 * * -----------------------------------------------------------------------------
3973 * * Author: antonermakov
3976 * #ifndef LOCAL_MATH_
3977 * #define LOCAL_MATH_
3979 * #define PI 3.14159265358979323846
3980 * #define TWOPI 6.283185307179586476925287
3981 * #define SECSINYEAR 3.155692608e+07
3985 * double factorial(int n)
3989 * } else if(n == 1) {
3991 * } else if(n == 2) {
3993 * } else if(n == 3) {
3995 * } else if(n == 4) {
4004 * double fudge(int m)
4018 * double sign(double x)
4022 * } else if(x < 0.0) {
4031 * double pv0(double x)
4037 * ans = x - TWOPI*floor(x/TWOPI);
4038 * if(ans > TWOPI/2.0) {
4039 * ans = ans - TWOPI;
4054 * double System::Plm(int m, double x)
4057 * return(1.5*x*x - 0.5);
4058 * } else if(m == 1) {
4059 * return(3.0*x*sqrt(1.0 - x*x));
4060 * } else if(m == 2) {
4061 * return(3.0 - 3.0*x*x);
4069 * double System::DP(int m, double x)
4073 * } else if(m == 1) {
4074 * return((3.0 - 6.0*x*x)/sqrt(1.0 - x*x));
4075 * } else if(m == 2) {
4089 * #endif /* LOCALMATH_H */
typename SparseLUDecomposition< number >::AdditionalData AdditionalData
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
#define AssertThrow(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
LinearOperator< Range_2, Domain_2, Payload > schur_complement(const LinearOperator< Domain_1, Range_1, Payload > &A_inv, const LinearOperator< Range_1, Domain_2, Payload > &B, const LinearOperator< Range_2, Domain_1, Payload > &C, const LinearOperator< Range_2, Domain_2, Payload > &D)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void compute_no_normal_flux_constraints(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int first_vector_component, const std::set< types::boundary_id > &boundary_ids, AffineConstraints< number > &constraints, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const bool use_manifold_for_normal=true)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
std::vector< value_type > l2_norm(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
@ component_is_part_of_vector
Expression floor(const Expression &x)
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
void Cuthill_McKee(DoFHandler< dim, spacedim > &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false, const std::vector< types::global_dof_index > &starting_indices=std::vector< types::global_dof_index >())
@ valid
Iterator points to a valid object.
@ matrix
Contents is actually a matrix.
constexpr types::blas_int zero
constexpr types::blas_int one
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Number angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
int(&) functions(const void *v1, const void *v2)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
inline ::VectorizedArray< Number, width > atan(const ::VectorizedArray< Number, width > &x)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation