DoFHandler<dim> dof_handler;
unsigned int n_u = 0, n_p = 0;
unsigned int plastic_iteration = 0;
+ unsigned int last_max_plasticity = 0;
QGauss<dim> quadrature_formula;
std::vector< std::vector <Vector<double> > > quad_viscosities; // Indices for this object are [cell][q][q coords, eta]
Assert(
local_quadrature_points_history < &quadrature_point_history.back(),
ExcInternalError());
+
+ double cell_area = cell->measure();
+ if(cell_area<0)
+ append_physical_times(-1);
+ AssertThrow(cell_area > 0
+ ,
+ ExcInternalError());
+
unsigned int m_id = cell->material_id();
}
for (unsigned int q = 0; q < n_q_points; ++q) {
- const SymmetricTensor<2, dim> &old_stress =
+ SymmetricTensor<2, dim> &old_stress =
local_quadrature_points_history[q].old_stress;
double &local_old_phiphi_stress =
local_quadrature_points_history[q].old_phiphi_stress;
double r_value = fe_values.quadrature_point(q)[0];
double z_value = fe_values.quadrature_point(q)[1];
+
+ // if(system_parameters::present_timestep == system_parameters::initial_elastic_iterations)
+ // {
+ // old_stress *= 0;
+ // local_old_phiphi_stress = 0;
+ // }
// get local density based on mat id
double local_density = system_parameters::rho[m_id];
std::vector<Vector<double> > vector_values(0);
std::vector < std::vector<Tensor<1, dim> > > gradient_values(0);
std::vector<bool> failing_cells;
+ // Write the stresses from the previous step into vectors
+ std::vector<SymmetricTensor<2, dim>> old_stress;
+ std::vector<double> old_phiphi_stress;
+ std::vector<double> cell_Gs;
for (typename DoFHandler<dim>::active_cell_iterator cell =
dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
{
+ // Makes pointer to data in quadrature_point_history
+ PointHistory<dim> *local_quadrature_points_history =
+ reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
+
fe_values.reinit(cell);
fe_values.get_function_gradients(solution, velocity_grads);
fe_values.get_function_values(solution, velocities);
Vector<double> current_cell_velocity(dim+1);
std::vector<Tensor<1, dim>> current_cell_grads(dim+1);
+ SymmetricTensor<2, dim> current_cell_old_stress;
+ current_cell_old_stress = 0;
+ double current_cell_old_phiphi_stress = 0;
double cell_area = 0;
- // Averages across each cell to find mean velocities and gradients
+ // Averages across each cell to find mean velocities, gradients, and old stresses
for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
{
cell_area += fe_values.JxW(q);
velocity_grads[q][i] *= fe_values.JxW(q);
current_cell_grads[i] += velocity_grads[q][i];
}
+ current_cell_old_stress += local_quadrature_points_history[q].old_stress * fe_values.JxW(q);
+ current_cell_old_phiphi_stress += local_quadrature_points_history[q].old_phiphi_stress * fe_values.JxW(q);
}
current_cell_velocity /= cell_area;
for (unsigned int i = 0; i < (dim+1); i++)
current_cell_grads[i] /= cell_area;
+ current_cell_old_stress /= cell_area;
+ current_cell_old_phiphi_stress /= cell_area;
vector_values.push_back(current_cell_velocity);
gradient_values.push_back(current_cell_grads);
+ old_stress.push_back(current_cell_old_stress);
+ old_phiphi_stress.push_back(current_cell_old_phiphi_stress);
+
+ // Get cell shear modulus: assumes it's constant for the cell
+ unsigned int mat_id = cell->material_id();
+ double local_G = rheology.get_G(mat_id);
+ cell_Gs.push_back(local_G);
}
//tracks where failure occurred
current_cell_viscosity = cell_viscosities[i];
}
+
+ double cell_eta_ve = 2
+ / ((1 / current_cell_viscosity)
+ + (1 / cell_Gs[i]
+ / system_parameters::current_time_interval));
+ double cell_chi_ve = 1
+ / (1
+ + (cell_Gs[i]
+ * system_parameters::current_time_interval
+ / current_cell_viscosity));
+
//find local pressure
double cell_p = vector_values[i].operator()(2);
//find stresses tensor
A << gradient_values[i][0][0] << 0 << sigma13 << endr
<< 0 << vector_values[i].operator()(0) / points_list[i].operator()(0)<< 0 << endr
<< sigma13 << 0 << gradient_values[i][1][1] << endr;
+ mat olddevstress;
+ olddevstress << old_stress[i][0][0] << 0 << old_stress[i][0][1] << endr
+ << 0 << old_phiphi_stress[i] << 0 << endr
+ << old_stress[i][0][1] << 0 << old_stress[i][1][1] << endr;
vec P;
P << cell_p << cell_p << cell_p;
mat Pmat = diagmat(P);
mat B;
- B = (2 * current_cell_viscosity * A) - Pmat;
+ B = (cell_eta_ve * A + cell_chi_ve * olddevstress) - Pmat;
//finds principal stresses
vec eigval;
if (sigma3 < 0)
temp_reductionfactor = 100;
else
- temp_reductionfactor = 1.5 * sigma1 / 5 / sigma3;
+ temp_reductionfactor = 1.9 * sigma1 / 5 / sigma3;
reduction_factor.push_back(temp_reductionfactor);
total_fails++;
// If there are enough failed cells, update eta at all quadrature points and perform smoothing
std::cout << " Number of failing cells: " << total_fails << "\n";
- if (total_fails <= 120)
+ double last_max_plasticity_double = last_max_plasticity;
+ double total_fails_double = total_fails;
+ double decrease_in_plasticity = ((last_max_plasticity_double - total_fails_double) / last_max_plasticity_double);
+ if(plastic_iteration == 0)
+ decrease_in_plasticity = 1;
+ last_max_plasticity = total_fails;
+ if (total_fails <= 100 || decrease_in_plasticity <= 0.2)
{
system_parameters::continue_plastic_iterations = false;
for(unsigned int j=0; j < triangulation.n_active_cells(); j++)
template<int dim>
void StokesProblem<dim>::update_time_interval()
{
- double move_goal_per_step = system_parameters::initial_disp_target -
+ double move_goal_per_step = system_parameters::initial_disp_target;
+ if(system_parameters::present_timestep > system_parameters::initial_elastic_iterations)
+ {
+ move_goal_per_step = system_parameters::initial_disp_target -
((system_parameters::initial_disp_target - system_parameters::final_disp_target) /
- system_parameters::total_viscous_steps *
- (system_parameters::present_timestep - system_parameters::initial_elastic_iterations));
+ system_parameters::total_viscous_steps *
+ (system_parameters::present_timestep - system_parameters::initial_elastic_iterations));
+ }
double zero_tolerance = 1e-3;
double max_velocity = 0;
}
}
}
- std:: cout << "Breaker " << max_velocity;
// NOTE: It is possible for this time interval to be very different from that used in the viscoelasticity calculation.
system_parameters::current_time_interval = move_goal_per_step / max_velocity;
double step_time_yr = system_parameters::current_time_interval / SECSINYEAR;
- std::cout << "\n Viscous time for moving mesh: " << step_time_yr << " yr";
+ std::cout << "Timestep interval changed to: "
+ << step_time_yr
+ << " years\n";
}
//====================== MOVE MESH ======================
for (unsigned int d = 0; d < dim; ++d)
vertex_displacement[d] = solution(
cell->vertex_dof_index(v, d));
-
cell->vertex(v) += vertex_displacement
* system_parameters::current_time_interval;
}
std::ofstream fout_final_vertices(vertices_output.str().c_str());
fout_final_vertices.close();
- // Figure out if the vertex is on the boundary of the domain
std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
+
+ if (boundary_that_we_need == 0)
+ {
+ // Figure out if the vertex is on the boundary of the domain
for (typename Triangulation<dim>::active_cell_iterator cell =
triangulation.begin_active(); cell != triangulation.end(); ++cell)
for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
}
}
}
+ }
+ else
+ {
+ // Figure out if the vertex is on an internal boundary
+ for (typename Triangulation<dim>::active_cell_iterator cell =
+ triangulation.begin_active(); cell != triangulation.end(); ++cell)
+ for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+ {
+ if (cell->neighbor(f) != triangulation.end()) {
+ if (cell->material_id() != cell->neighbor(f)->material_id()) //finds face is at internal boundary
+ {
+ int high_mat_id = std::max(cell->material_id(),
+ cell->neighbor(f)->material_id());
+ if (high_mat_id == boundary_that_we_need) //finds faces at the correct internal boundary
+ {
+ for (unsigned int v = 0;
+ v < GeometryInfo<dim>::vertices_per_face;
+ ++v)
+ if (vertex_touched[cell->face(f)->vertex_index(
+ v)] == false) {
+ vertex_touched[cell->face(f)->vertex_index(
+ v)] = true;
+ std::ofstream fout_final_vertices(vertices_output.str().c_str(), std::ios::app);
+ fout_final_vertices << cell->face(f)->vertex(v) << "\n";
+ fout_final_vertices.close();
+ }
+ }
+ }
+ }
+ }
+ }
}
//====================== SETUP INITIAL MESH ======================
triangulation.refine_global(system_parameters::global_refinement);
-//
-////refines region near r=0
-// if (system_parameters::small_r_refinement != 0) {
-// for (unsigned int step = 0;
-// step < system_parameters::small_r_refinement; ++step) {
-// // std::cout << "iteration " << step + 1 << "\n";
-// typename dealii::Triangulation<dim>::active_cell_iterator cell =
-// triangulation.begin_active(), endc = triangulation.end();
-//
-// for (; cell != endc; ++cell)
-// for (unsigned int v = 0;
-// v < GeometryInfo<dim>::vertices_per_cell; ++v) {
-// Point<dim> current_vertex = cell->vertex(v);
-//
-// const double x_coord = current_vertex.operator()(0);
-//
-// if (std::fabs(x_coord) < 1e-10) {
-// cell->set_refine_flag();
-// break;
-// }
-//
-// }
-// triangulation.execute_coarsening_and_refinement();
-// }
-// }
-//
//refines crustal region
if (system_parameters::crustal_refinement != 0) {
}
}
-//
-// //refines surface region
-// if (system_parameters::surface_refinement != 0) {
-// for (unsigned int step = 0;
-// step < system_parameters::surface_refinement; ++step) {
-// // std::cout << "iteration " << step + 1 << "\n";
-// typename dealii::Triangulation<dim>::active_cell_iterator cell =
-// triangulation.begin_active(), endc = triangulation.end();
-// for (; cell != endc; ++cell)
-// for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell;
-// ++f) {
-// if (cell->face(f)->at_boundary()) {
-// if (cell->face(f)->center()[0] != 0) {
-// if (cell->face(f)->center()[1] != 0)
-// if (cell->face(f)->center()[0] > 1) {
-// cell->set_refine_flag();
-// break;
-// }
-// }
-// }
-// }
-// triangulation.execute_coarsening_and_refinement();
-// }
-// }
// output initial mesh in eps
std::ostringstream refined_mesh_file;
if (system_parameters::present_timestep == 0)
initialize_eta_and_G();
-
- system_parameters::current_time_interval =
- system_parameters::elastic_time; //This is the time interval needed in assembling the problem
+
+ if(elastic_iteration == 0)
+ system_parameters::current_time_interval =
+ system_parameters::viscous_time; //This is the time interval needed in assembling the problem
std::cout << " Assembling..." << std::endl << std::flush;
assemble_system();
output_results();
update_quadrature_point_history();
- // std::cout << std::endl << "\a";
+ append_physical_times(0);
elastic_iteration++;
system_parameters::present_timestep++;
- system_parameters::present_time = system_parameters::present_time + system_parameters::current_time_interval;
- move_mesh();
do_ellipse_fits();
write_vertices(0);
- write_vertices(1);
+ write_vertices(1);
write_mesh();
+ update_time_interval();
}
}
times_filename << system_parameters::output_folder << "/physical_times.txt";
std::ofstream fout_times(times_filename.str().c_str());
fout_times.close();
- append_physical_times(0);
// Computes elastic timesteps
do_elastic_steps();
-
// Computes viscous timesteps
- system_parameters::current_time_interval = system_parameters::viscous_time;
unsigned int VEPstep = 0;
while (system_parameters::present_timestep
< (system_parameters::initial_elastic_iterations
// Computes plasticity
do_flow_step();
update_quadrature_point_history();
- update_time_interval();
+ move_mesh();
append_physical_times(plastic_iteration);
system_parameters::present_timestep++;
system_parameters::present_time = system_parameters::present_time + system_parameters::current_time_interval;
- move_mesh();
do_ellipse_fits();
write_vertices(0);
write_vertices(1);
write_mesh();
VEPstep++;
}
- append_physical_times(0);
+ append_physical_times(-1);
}
}
}
return 0;
-}
\ No newline at end of file
+}
projects / code-gallery.git / commitdiff