233 * #include <deal.II/base/mpi.h>
234 * #include <deal.II/base/
function.h>
235 * #include <deal.II/base/parameter_handler.h>
236 * #include <deal.II/base/
point.h>
237 * #include <deal.II/base/quadrature_lib.h>
238 * #include <deal.II/base/symmetric_tensor.h>
239 * #include <deal.II/base/tensor.h>
240 * #include <deal.II/base/timer.h>
241 * #include <deal.II/base/work_stream.h>
242 * #include <deal.II/dofs/dof_renumbering.h>
243 * #include <deal.II/dofs/dof_tools.h>
244 * #include <deal.II/base/quadrature_point_data.h>
245 * #include <deal.II/grid/filtered_iterator.h>
246 * #include <deal.II/grid/grid_generator.h>
247 * #include <deal.II/grid/grid_tools.h>
248 * #include <deal.II/grid/grid_in.h>
249 * #include <deal.II/grid/manifold_lib.h>
250 * #include <deal.II/grid/tria.h>
251 * #include <deal.II/grid/tria_boundary_lib.h>
252 * #include <deal.II/fe/fe_dgp_monomial.h>
253 * #include <deal.II/fe/fe_q.h>
254 * #include <deal.II/fe/fe_system.h>
255 * #include <deal.II/fe/fe_tools.h>
256 * #include <deal.II/fe/fe_values.h>
257 * #include <deal.II/fe/mapping_q_eulerian.h>
258 * #include <deal.II/lac/block_sparsity_pattern.h>
259 * #include <deal.II/lac/dynamic_sparsity_pattern.h>
260 * #include <deal.II/lac/full_matrix.h>
261 * #include <deal.II/lac/constraint_matrix.h>
262 * #include <deal.II/lac/solver_selector.h>
263 * #include <deal.II/lac/trilinos_block_sparse_matrix.h>
264 * #include <deal.II/lac/trilinos_precondition.h>
265 * #include <deal.II/lac/trilinos_sparsity_pattern.h>
266 * #include <deal.II/lac/trilinos_sparse_matrix.h>
267 * #include <deal.II/lac/trilinos_vector.h>
269 * #include <deal.II/lac/packaged_operation.h>
270 * #include <deal.II/lac/trilinos_linear_operator.h>
271 * #include <deal.II/numerics/data_out.h>
272 * #include <deal.II/numerics/vector_tools.h>
273 * #include <deal.II/physics/transformations.h>
274 * #include <deal.II/physics/elasticity/kinematics.h>
275 * #include <deal.II/physics/elasticity/standard_tensors.h>
276 * #include <iostream>
280 * #include <deal.II/grid/grid_out.h>
282 *
namespace ViscoElasStripHole
286 *
namespace Parameters
288 *
struct BoundaryConditions
290 * BoundaryConditions();
292 * std::string driver;
311 * BoundaryConditions::BoundaryConditions()
313 * driver (
"Neumann"),
317 * boundary_id_minus_X (1),
318 * boundary_id_plus_X (2),
319 * boundary_id_minus_Y (3),
320 * boundary_id_plus_Y (4),
321 * boundary_id_minus_Z (5),
322 * boundary_id_plus_Z (6),
323 * boundary_id_hole (10),
324 * manifold_id_hole (10)
328 * prm.enter_subsection(
"Boundary conditions");
330 * prm.declare_entry(
"Driver",
"Dirichlet",
332 *
"Driver boundary condition for the problem");
333 * prm.declare_entry(
"Final stretch",
"2.0",
335 *
"Positive stretch applied length-ways to the strip");
336 * prm.declare_entry(
"Applied pressure",
"0.0",
338 *
"Hydrostatic pressure applied (in the referential configuration) to the interior surface of the hole");
339 * prm.declare_entry(
"Load time",
"2.5",
341 *
"Total time over which the stretch/pressure is ramped up");
343 * prm.leave_subsection();
347 * prm.enter_subsection(
"Boundary conditions");
349 * driver = prm.get(
"Driver");
350 * stretch = prm.get_double(
"Final stretch");
351 * pressure = prm.get_double(
"Applied pressure");
352 * load_time = prm.get_double(
"Load time");
354 * prm.leave_subsection();
358 *
unsigned int poly_degree;
359 *
unsigned int quad_order;
367 * prm.enter_subsection(
"Finite element system");
369 * prm.declare_entry(
"Polynomial degree",
"2",
371 *
"Displacement system polynomial order");
372 * prm.declare_entry(
"Quadrature order",
"3",
374 *
"Gauss quadrature order");
376 * prm.leave_subsection();
380 * prm.enter_subsection(
"Finite element system");
382 * poly_degree = prm.get_integer(
"Polynomial degree");
383 * quad_order = prm.get_integer(
"Quadrature order");
385 * prm.leave_subsection();
392 *
double hole_diameter;
393 *
double hole_division_fraction;
394 *
unsigned int n_repetitions_xy;
395 *
unsigned int n_repetitions_z;
396 *
unsigned int global_refinement;
405 * prm.enter_subsection(
"Geometry");
407 * prm.declare_entry(
"Length",
"100.0",
409 *
"Total sample length");
410 * prm.declare_entry(
"Width",
"50.0",
412 *
"Total sample width");
413 * prm.declare_entry(
"Thickness",
"5.0",
415 *
"Total sample thickness");
416 * prm.declare_entry(
"Hole diameter",
"20.0",
419 * prm.declare_entry(
"Hole division fraction",
"0.5",
421 *
"A geometric factor affecting the discretisation near the hole");
422 * prm.declare_entry(
"Number of subdivisions in cross-section",
"2",
424 *
"A factor defining the number of initial grid subdivisions in the cross-section");
425 * prm.declare_entry(
"Number of subdivisions thickness",
"6",
427 *
"A factor defining the number of initial grid subdivisions through the thickness");
428 * prm.declare_entry(
"Global refinement",
"2",
430 *
"Global refinement level");
431 * prm.declare_entry(
"Grid scale",
"1e-3",
433 *
"Global grid scaling factor");
435 * prm.leave_subsection();
439 * prm.enter_subsection(
"Geometry");
441 * length = prm.get_double(
"Length");
442 * width = prm.get_double(
"Width");
443 * thickness = prm.get_double(
"Thickness");
444 * hole_diameter = prm.get_double(
"Hole diameter");
445 * hole_division_fraction = prm.get_double(
"Hole division fraction");
446 * n_repetitions_xy = prm.get_integer(
"Number of subdivisions in cross-section");
447 * n_repetitions_z = prm.get_integer(
"Number of subdivisions thickness");
448 * global_refinement = prm.get_integer(
"Global refinement");
449 *
scale = prm.get_double(
"Grid scale");
451 * prm.leave_subsection();
466 * prm.enter_subsection(
"Material properties");
468 * prm.declare_entry(
"Poisson's ratio",
"0.4999",
470 *
"Poisson's ratio");
471 * prm.declare_entry(
"Elastic shear modulus",
"80.194e6",
473 *
"Elastic shear modulus");
474 * prm.declare_entry(
"Viscous shear modulus",
"80.194e6",
476 *
"Viscous shear modulus");
477 * prm.declare_entry(
"Viscous relaxation time",
"2.0",
479 *
"Viscous relaxation time");
481 * prm.leave_subsection();
485 * prm.enter_subsection(
"Material properties");
487 * nu_e = prm.get_double(
"Poisson's ratio");
488 * mu_e = prm.get_double(
"Elastic shear modulus");
489 * mu_v = prm.get_double(
"Viscous shear modulus");
490 * tau_v = prm.get_double(
"Viscous relaxation time");
492 * prm.leave_subsection();
494 *
struct LinearSolver
496 * std::string type_lin;
498 *
double max_iterations_lin;
506 * prm.enter_subsection(
"Linear solver");
508 * prm.declare_entry(
"Solver type",
"cg",
510 *
"Type of solver used to solve the linear system");
511 * prm.declare_entry(
"Residual",
"1e-6",
513 *
"Linear solver residual (scaled by residual norm)");
514 * prm.declare_entry(
"Max iteration multiplier",
"1",
516 *
"Linear solver iterations (multiples of the system matrix size)");
518 * prm.leave_subsection();
522 * prm.enter_subsection(
"Linear solver");
524 * type_lin = prm.get(
"Solver type");
525 * tol_lin = prm.get_double(
"Residual");
526 * max_iterations_lin = prm.get_double(
"Max iteration multiplier");
528 * prm.leave_subsection();
530 *
struct NonlinearSolver
532 *
unsigned int max_iterations_NR;
542 * prm.enter_subsection(
"Nonlinear solver");
544 * prm.declare_entry(
"Max iterations Newton-Raphson",
"10",
546 *
"Number of Newton-Raphson iterations allowed");
547 * prm.declare_entry(
"Tolerance displacement",
"1.0e-6",
549 *
"Displacement error tolerance");
550 * prm.declare_entry(
"Tolerance force",
"1.0e-9",
552 *
"Force residual tolerance");
554 * prm.leave_subsection();
558 * prm.enter_subsection(
"Nonlinear solver");
560 * max_iterations_NR = prm.get_integer(
"Max iterations Newton-Raphson");
561 * tol_f = prm.get_double(
"Tolerance force");
562 * tol_u = prm.get_double(
"Tolerance displacement");
564 * prm.leave_subsection();
577 * prm.enter_subsection(
"Time");
579 * prm.declare_entry(
"End time",
"1",
582 * prm.declare_entry(
"Time step size",
"0.1",
586 * prm.leave_subsection();
590 * prm.enter_subsection(
"Time");
592 * end_time = prm.get_double(
"End time");
593 * delta_t = prm.get_double(
"Time step size");
595 * prm.leave_subsection();
597 *
struct AllParameters
598 * :
public BoundaryConditions,
602 *
public LinearSolver,
603 *
public NonlinearSolver,
606 * AllParameters(
const std::string &input_file);
612 * AllParameters::AllParameters(
const std::string &input_file)
615 * declare_parameters(prm);
617 * parse_parameters(prm);
621 * BoundaryConditions::declare_parameters(prm);
622 * FESystem::declare_parameters(prm);
623 * Geometry::declare_parameters(prm);
624 * Materials::declare_parameters(prm);
625 * LinearSolver::declare_parameters(prm);
626 * NonlinearSolver::declare_parameters(prm);
627 * Time::declare_parameters(prm);
631 * BoundaryConditions::parse_parameters(prm);
632 * FESystem::parse_parameters(prm);
633 * Geometry::parse_parameters(prm);
634 * Materials::parse_parameters(prm);
635 * LinearSolver::parse_parameters(prm);
636 * NonlinearSolver::parse_parameters(prm);
637 * Time::parse_parameters(prm);
643 * Time (
const double time_end,
644 *
const double delta_t)
648 * time_end(time_end),
653 *
double current() const
655 *
return time_current;
661 *
double get_delta_t() const
665 *
unsigned int get_timestep() const
671 * time_current += delta_t;
675 *
unsigned int timestep;
676 *
double time_current;
677 *
const double time_end;
678 *
const double delta_t;
681 *
class Material_Compressible_Three_Field_Linear_Viscoelastic
684 * Material_Compressible_Three_Field_Linear_Viscoelastic(
const double mu_e,
687 *
const double tau_v,
690 * kappa((2.0 * mu_e * (1.0 + nu_e)) / (3.0 * (1.0 - 2.0 * nu_e))),
700 * ~Material_Compressible_Three_Field_Linear_Viscoelastic()
705 *
const double &p_tilde)
const
707 *
return get_tau_iso(
F) + get_tau_vol(
F,p_tilde);
710 *
const double &p_tilde)
const
712 *
return get_Jc_iso(
F) + get_Jc_vol(
F,p_tilde);
715 * get_dPsi_vol_dJ(
const double &J_tilde)
const
717 *
return (kappa / 2.0) * (J_tilde - 1.0 / J_tilde);
720 * get_d2Psi_vol_dJ2(
const double &J_tilde)
const
722 *
return ( (kappa / 2.0) * (1.0 + 1.0 / (J_tilde * J_tilde)));
734 * Assumes
first-oder backward Euler time discretisation
737 * Q_n_t = (1.0/(1.0 + time.get_delta_t()/tau_v))*(Q_t1 + (time.get_delta_t()/tau_v)*
invert(C_bar));
740 * update_end_timestep()
746 *
const double kappa;
749 *
const double tau_v;
756 *
const double &p_tilde)
const
775 * Elastic Neo-Hookean + Linder2011 eq 47
778 *
return mu_e * b_bar
782 *
const double &p_tilde)
const
785 *
return p_tilde * det_F
800 *
return (2.0 / dim) *
trace(tau_bar)
802 * - (2.0 / dim) * (tau_iso_x_I + I_x_tau_iso)
810 * Elastic Neo-Hookean + Linder2011 eq 56
822 *
virtual ~PointHistory()
825 * setup_lqp (
const Parameters::AllParameters ¶meters,
828 * material.reset(
new Material_Compressible_Three_Field_Linear_Viscoelastic<dim>(
829 * parameters.mu_e, parameters.nu_e,
830 * parameters.mu_v, parameters.tau_v,
836 *
const double &p_tilde)
const
838 *
return material->get_tau(
F, p_tilde);
842 *
const double &p_tilde)
const
844 *
return material->get_Jc(
F, p_tilde);
847 * get_dPsi_vol_dJ(
const double &J_tilde)
const
849 *
return material->get_dPsi_vol_dJ(J_tilde);
852 * get_d2Psi_vol_dJ2(
const double &J_tilde)
const
854 *
return material->get_d2Psi_vol_dJ2(J_tilde);
858 *
const double &p_tilde,
859 *
const double &J_tilde)
861 * material->update_internal_equilibrium(
F,p_tilde,J_tilde);
864 * update_end_timestep()
866 * material->update_end_timestep();
869 * std::shared_ptr< Material_Compressible_Three_Field_Linear_Viscoelastic<dim> > material;
875 * Solid(
const std::string &input_file);
881 *
struct PerTaskData_ASM;
882 *
struct ScratchData_ASM;
887 *
const double half_length,
888 *
const double half_width,
889 *
const double hole_radius,
890 *
const unsigned int n_repetitions_xy = 1,
891 *
const double hole_division_fraction = 0.25);
895 * determine_component_extractors();
900 * ScratchData_ASM &scratch,
901 * PerTaskData_ASM &data)
const;
903 * copy_local_to_global_system(
const PerTaskData_ASM &data);
905 * make_constraints(
const int &it_nr);
910 * std::pair<unsigned int, double>
915 * update_end_timestep();
917 * output_results(
const unsigned int timestep,
918 *
const double current_time)
const;
920 * compute_vertex_positions(std::vector<double> &real_time,
921 * std::vector<std::vector<
Point<dim> > > &tracked_vertices,
926 * Parallel communication
934 * Parameters::AllParameters parameters;
939 * PointHistory<dim> > quadrature_point_history;
940 *
const unsigned int degree;
943 *
const unsigned int dofs_per_cell;
947 *
static const unsigned int n_blocks = 3;
949 *
static const unsigned int first_u_component = 0;
950 *
static const unsigned int p_component = dim;
951 *
static const unsigned int J_component = dim + 1;
963 * std::vector<unsigned int> block_component;
970 * std::vector<IndexSet> all_locally_owned_dofs;
973 * std::vector<IndexSet> locally_owned_partitioning;
974 * std::vector<IndexSet> locally_relevant_partitioning;
975 * std::vector<types::global_dof_index> dofs_per_block;
976 * std::vector<types::global_dof_index> element_indices_u;
977 * std::vector<types::global_dof_index> element_indices_p;
978 * std::vector<types::global_dof_index> element_indices_J;
980 *
const QGauss<dim - 1> qf_face;
981 *
const unsigned int n_q_points;
982 *
const unsigned int n_q_points_f;
991 *
norm(1.0), u(1.0), p(1.0), J(1.0)
1000 *
void normalise(
const Errors &rhs)
1002 *
if (rhs.norm != 0.0)
1011 *
double norm, u, p, J;
1013 * Errors error_residual, error_residual_0, error_residual_norm, error_update,
1014 * error_update_0, error_update_norm;
1016 * get_error_residual(Errors &error_residual);
1019 * Errors &error_update);
1020 * std::pair<double, std::pair<double,double> >
1023 * print_conv_header();
1027 *
template <
int dim>
1028 * Solid<dim>::Solid(
const std::string &input_file)
1030 * mpi_communicator(MPI_COMM_WORLD),
1034 * parameters(input_file),
1036 * time(parameters.end_time, parameters.delta_t),
1037 * timer(mpi_communicator,
1041 * degree(parameters.poly_degree),
1042 * fe(
FE_Q<dim>(parameters.poly_degree), dim,
1046 * dofs_per_cell (fe.dofs_per_cell),
1047 * u_fe(first_u_component),
1048 * p_fe(p_component),
1049 * J_fe(J_component),
1051 * qf_cell(parameters.quad_order),
1052 * qf_face(parameters.quad_order),
1053 * n_q_points (qf_cell.
size()),
1054 * n_q_points_f (qf_face.size())
1056 *
Assert(dim==2 || dim==3,
ExcMessage(
"This problem only works in 2 or 3 space dimensions."));
1057 * determine_component_extractors();
1059 *
template <
int dim>
1060 * Solid<dim>::~Solid()
1062 * dof_handler.
clear();
1064 *
template <
int dim>
1070 * setup_system(solution_delta);
1073 * constraints.
close();
1082 * output_results(time.get_timestep(), time.current());
1087 * Some points
for post-processing
1090 * std::vector<double> real_time;
1091 * real_time.push_back(0);
1092 * std::vector<std::vector<Point<dim> > > tracked_vertices (4);
1095 * p[1] = parameters.length/2.0;
1096 * tracked_vertices[0].push_back(p*parameters.scale);
1100 * p[1] = parameters.hole_diameter/2.0;
1101 * tracked_vertices[1].push_back(p*parameters.scale);
1105 * p[0] = parameters.hole_diameter/2.0;
1106 * tracked_vertices[2].push_back(p*parameters.scale);
1110 * p[0] = parameters.width/2.0;
1111 * tracked_vertices[3].push_back(p*parameters.scale);
1114 *
while (time.current() < time.end()+0.01*time.get_delta_t())
1116 * solve_nonlinear_timestep(solution_delta);
1117 * solution_n += solution_delta;
1118 * solution_delta = 0.0;
1119 * output_results(time.get_timestep(), time.current());
1120 * compute_vertex_positions(real_time,
1122 * get_solution_total(solution_delta));
1123 * update_end_timestep();
1127 * pcout <<
"\n\n*** Spatial position history for tracked vertices ***" << std::endl;
1128 *
for (
unsigned int t=0; t<real_time.size(); ++t)
1133 *
for (
unsigned int p=0; p<tracked_vertices.size(); ++p)
1135 *
for (
unsigned int d=0;
d<dim; ++
d)
1137 * pcout <<
"Point " << p <<
" [" <<
d <<
"]";
1138 *
if (!(p == tracked_vertices.size()-1 &&
d == dim-1))
1142 * pcout << std::endl;
1145 * pcout << std::setprecision(6);
1146 * pcout << real_time[t] <<
",";
1147 *
for (
unsigned int p=0; p<tracked_vertices.size(); ++p)
1149 *
Assert(tracked_vertices[p].size() == real_time.size(),
1150 *
ExcMessage(
"Vertex not tracked at each timestep"));
1151 *
for (
unsigned int d=0;
d<dim; ++
d)
1153 * pcout << tracked_vertices[p][t][
d];
1154 *
if (!(p == tracked_vertices.size()-1 &&
d == dim-1))
1158 * pcout << std::endl;
1161 *
template <
int dim>
1162 *
struct Solid<dim>::PerTaskData_ASM
1166 * std::vector<types::global_dof_index> local_dof_indices;
1167 * PerTaskData_ASM(
const unsigned int dofs_per_cell)
1170 * cell_rhs(dofs_per_cell),
1171 * local_dof_indices(dofs_per_cell)
1179 *
template <
int dim>
1180 *
struct Solid<dim>::ScratchData_ASM
1186 * Integration helper
1197 * std::vector<Tensor<2, dim> > solution_grads_u_total;
1198 * std::vector<double> solution_values_p_total;
1199 * std::vector<double> solution_values_J_total;
1203 * Shape
function values and gradients
1206 * std::vector<std::vector<double> > Nx;
1207 * std::vector<std::vector<Tensor<2, dim> > > grad_Nx;
1208 * std::vector<std::vector<SymmetricTensor<2, dim> > > symm_grad_Nx;
1215 * solution_total (solution_total),
1216 * fe_values_ref(fe_cell, qf_cell, uf_cell),
1217 * fe_face_values_ref(fe_cell, qf_face, uf_face),
1218 * solution_grads_u_total(qf_cell.
size()),
1219 * solution_values_p_total(qf_cell.
size()),
1220 * solution_values_J_total(qf_cell.
size()),
1221 * Nx(qf_cell.
size(),
1222 * std::vector<double>(fe_cell.dofs_per_cell)),
1223 * grad_Nx(qf_cell.
size(),
1225 * symm_grad_Nx(qf_cell.
size(),
1227 * (fe_cell.dofs_per_cell))
1229 * ScratchData_ASM(
const ScratchData_ASM &rhs)
1231 * solution_total (rhs.solution_total),
1232 * fe_values_ref(rhs.fe_values_ref.get_fe(),
1233 * rhs.fe_values_ref.get_quadrature(),
1234 * rhs.fe_values_ref.get_update_flags()),
1235 * fe_face_values_ref(rhs.fe_face_values_ref.get_fe(),
1236 * rhs.fe_face_values_ref.get_quadrature(),
1237 * rhs.fe_face_values_ref.get_update_flags()),
1238 * solution_grads_u_total(rhs.solution_grads_u_total),
1239 * solution_values_p_total(rhs.solution_values_p_total),
1240 * solution_values_J_total(rhs.solution_values_J_total),
1242 * grad_Nx(rhs.grad_Nx),
1243 * symm_grad_Nx(rhs.symm_grad_Nx)
1247 *
const unsigned int n_q_points = solution_grads_u_total.size();
1248 *
const unsigned int n_dofs_per_cell = Nx[0].size();
1250 *
Assert(solution_grads_u_total.size() == n_q_points,
1252 *
Assert(solution_values_p_total.size() == n_q_points,
1254 *
Assert(solution_values_J_total.size() == n_q_points,
1256 *
Assert(Nx.size() == n_q_points,
1258 *
Assert(grad_Nx.size() == n_q_points,
1260 *
Assert(symm_grad_Nx.size() == n_q_points,
1263 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1266 *
Assert( grad_Nx[q_point].size() == n_dofs_per_cell,
1268 *
Assert( symm_grad_Nx[q_point].size() == n_dofs_per_cell,
1271 * solution_grads_u_total[q_point] = 0.0;
1272 * solution_values_p_total[q_point] = 0.0;
1273 * solution_values_J_total[q_point] = 0.0;
1274 *
for (
unsigned int k = 0; k < n_dofs_per_cell; ++k)
1276 * Nx[q_point][k] = 0.0;
1277 * grad_Nx[q_point][k] = 0.0;
1278 * symm_grad_Nx[q_point][k] = 0.0;
1284 *
void Solid<2>::make_grid()
1286 *
const int dim = 2;
1287 *
const double tol = 1
e-12;
1289 * parameters.length/2.0,
1290 * parameters.width/2.0,
1291 * parameters.hole_diameter/2.0,
1292 * parameters.n_repetitions_xy,
1293 * parameters.hole_division_fraction);
1297 * Clear boundary ID
's
1300 * for (typename Triangulation<dim>::active_cell_iterator
1301 * cell = triangulation.begin_active();
1302 * cell != triangulation.end(); ++cell)
1304 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1305 * if (cell->face(face)->at_boundary())
1307 * cell->face(face)->set_all_boundary_ids(0);
1313 * Set boundary IDs and and manifolds
1316 * const Point<dim> centre (0,0);
1317 * for (typename Triangulation<dim>::active_cell_iterator
1318 * cell = triangulation.begin_active();
1319 * cell != triangulation.end(); ++cell)
1321 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1322 * if (cell->face(face)->at_boundary())
1329 * if (std::abs(cell->face(face)->center()[0] - 0.0) < tol)
1331 * cell->face(face)->set_boundary_id(parameters.boundary_id_minus_X);
1333 * else if (std::abs(cell->face(face)->center()[0] - parameters.width/2.0) < tol)
1335 * cell->face(face)->set_boundary_id(parameters.boundary_id_plus_X);
1337 * else if (std::abs(cell->face(face)->center()[1] - 0.0) < tol)
1339 * cell->face(face)->set_boundary_id(parameters.boundary_id_minus_Y);
1341 * else if (std::abs(cell->face(face)->center()[1] - parameters.length/2.0) < tol)
1343 * cell->face(face)->set_boundary_id(parameters.boundary_id_plus_Y);
1347 * for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
1348 * if (std::abs(cell->vertex(vertex).distance(centre) - parameters.hole_diameter/2.0) < tol)
1350 * cell->face(face)->set_boundary_id(parameters.boundary_id_hole);
1360 * for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
1361 * if (std::abs(cell->vertex(vertex).distance(centre) - parameters.hole_diameter/2.0) < tol)
1363 * cell->face(face)->set_manifold_id(parameters.manifold_id_hole);
1368 * static SphericalManifold<dim> spherical_manifold (centre);
1369 * triangulation.set_manifold(parameters.manifold_id_hole,spherical_manifold);
1370 * triangulation.refine_global(parameters.global_refinement);
1371 * GridTools::scale(parameters.scale,triangulation);
1374 * void Solid<3>::make_grid()
1376 * const int dim = 3;
1377 * const double tol = 1e-12;
1378 * Triangulation<2> tria_2d;
1379 * make_2d_quarter_plate_with_hole(tria_2d,
1380 * parameters.length/2.0,
1381 * parameters.width/2.0,
1382 * parameters.hole_diameter/2.0,
1383 * parameters.n_repetitions_xy,
1384 * parameters.hole_division_fraction);
1385 * GridGenerator::extrude_triangulation(tria_2d,
1386 * parameters.n_repetitions_z+1,
1387 * parameters.thickness/2.0,
1392 * Clear boundary ID's
1399 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1400 *
if (cell->face(face)->at_boundary())
1402 * cell->face(face)->set_all_boundary_ids(0);
1408 * Set boundary IDs and and manifolds
1417 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1418 *
if (cell->face(face)->at_boundary())
1425 *
if (std::abs(cell->face(face)->center()[0] - 0.0) < tol)
1427 * cell->face(face)->set_boundary_id(parameters.boundary_id_minus_X);
1429 *
else if (std::abs(cell->face(face)->center()[0] - parameters.width/2.0) < tol)
1431 * cell->face(face)->set_boundary_id(parameters.boundary_id_plus_X);
1433 *
else if (std::abs(cell->face(face)->center()[1] - 0.0) < tol)
1435 * cell->face(face)->set_boundary_id(parameters.boundary_id_minus_Y);
1437 *
else if (std::abs(cell->face(face)->center()[1] - parameters.length/2.0) < tol)
1439 * cell->face(face)->set_boundary_id(parameters.boundary_id_plus_Y);
1441 *
else if (std::abs(cell->face(face)->center()[2] - 0.0) < tol)
1443 * cell->face(face)->set_boundary_id(parameters.boundary_id_minus_Z);
1445 *
else if (std::abs(cell->face(face)->center()[2] - parameters.thickness/2.0) < tol)
1447 * cell->face(face)->set_boundary_id(parameters.boundary_id_plus_Z);
1451 *
for (
unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
1455 * Project the cell vertex to the XY plane and
1456 * test the distance from the
cylinder axis
1459 *
Point<dim> vertex_proj = cell->vertex(vertex);
1460 * vertex_proj[2] = 0.0;
1461 *
if (std::abs(vertex_proj.distance(centre) - parameters.hole_diameter/2.0) < tol)
1463 * cell->face(face)->set_boundary_id(parameters.boundary_id_hole);
1474 *
for (
unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
1478 * Project the cell vertex to the XY plane and
1479 * test the distance from the
cylinder axis
1482 *
Point<dim> vertex_proj = cell->vertex(vertex);
1483 * vertex_proj[2] = 0.0;
1484 *
if (std::abs(vertex_proj.distance(centre) - parameters.hole_diameter/2.0) < 1
e-12)
1488 * Set manifold ID on face and edges
1491 * cell->face(face)->set_all_manifold_ids(parameters.manifold_id_hole);
1498 *
triangulation.set_manifold(parameters.manifold_id_hole,cylindrical_manifold);
1499 *
triangulation.refine_global(parameters.global_refinement);
1502 *
template <
int dim>
1503 *
void Solid<dim>::make_2d_quarter_plate_with_hole(
Triangulation<2> &tria_2d,
1504 *
const double half_length,
1505 *
const double half_width,
1506 *
const double hole_radius,
1507 *
const unsigned int n_repetitions_xy,
1508 *
const double hole_division_fraction)
1510 *
const double length = 2.0*half_length;
1511 *
const double width = 2.0*half_width;
1512 *
const double hole_diameter = 2.0*hole_radius;
1514 *
const double internal_width = hole_diameter + hole_division_fraction*(width - hole_diameter);
1519 * hole_diameter/2.0,
1520 * internal_width/2.0);
1522 * std::set<typename Triangulation<2>::active_cell_iterator > cells_to_remove;
1525 * cell != tria_plate_hole.
end(); ++cell)
1529 * Remove all cells that are not in the
first quadrant
1532 *
if (cell->center()[0] < 0.0 || cell->center()[1] < 0.0)
1533 * cells_to_remove.insert(cell);
1545 * Subdivide the plate so that we
're left one
1546 * cell to remove (we'll replace
this with the
1547 * plate with the hole) and then make the
1548 * rest of the subdivisions so that we
're left
1549 * with cells with a decent aspect ratio
1552 * std::vector<std::vector<double> > step_sizes;
1554 * std::vector<double> subdivision_width;
1555 * subdivision_width.push_back(internal_width/2.0);
1556 * const double width_remaining = (width - internal_width)/2.0;
1557 * const unsigned int n_subs = std::max(1.0,std::ceil(width_remaining/(internal_width/2.0)));
1558 * Assert(n_subs>0, ExcInternalError());
1559 * for (unsigned int s=0; s<n_subs; ++s)
1560 * subdivision_width.push_back(width_remaining/n_subs);
1561 * step_sizes.push_back(subdivision_width);
1563 * const double sum_half_width = std::accumulate(subdivision_width.begin(), subdivision_width.end(), 0.0);
1564 * Assert(std::abs(sum_half_width-width/2.0) < 1e-12, ExcInternalError());
1567 * std::vector<double> subdivision_length;
1568 * subdivision_length.push_back(internal_width/2.0);
1569 * const double length_remaining = (length - internal_width)/2.0;
1570 * const unsigned int n_subs = std::max(1.0,std::ceil(length_remaining/(internal_width/2.0)));
1571 * Assert(n_subs>0, ExcInternalError());
1572 * for (unsigned int s=0; s<n_subs; ++s)
1573 * subdivision_length.push_back(length_remaining/n_subs);
1574 * step_sizes.push_back(subdivision_length);
1576 * const double sum_half_length = std::accumulate(subdivision_length.begin(), subdivision_length.end(), 0.0);
1577 * Assert(std::abs(sum_half_length-length/2.0) < 1e-12, ExcInternalError());
1580 * GridGenerator::subdivided_hyper_rectangle(tria_plate,
1582 * Point<2>(0.0, 0.0),
1583 * Point<2>(width/2.0, length/2.0));
1585 * std::set<typename Triangulation<2>::active_cell_iterator > cells_to_remove;
1586 * for (typename Triangulation<2>::active_cell_iterator
1587 * cell = tria_plate.begin_active();
1588 * cell != tria_plate.end(); ++cell)
1592 * Remove all cells that are in the first quadrant
1595 * if (cell->center()[0] < internal_width/2.0 && cell->center()[1] < internal_width/2.0)
1596 * cells_to_remove.insert(cell);
1598 * Assert(cells_to_remove.size() > 0, ExcInternalError());
1599 * Assert(cells_to_remove.size() != tria_plate.n_active_cells(), ExcInternalError());
1600 * GridGenerator::create_triangulation_with_removed_cells(tria_plate,cells_to_remove,tria_cut_plate);
1603 * Triangulation<2> tria_2d_not_flat;
1604 * GridGenerator::merge_triangulations(tria_quarter_plate_hole,
1606 * tria_2d_not_flat);
1610 * Attach a manifold to the curved boundary and refine
1611 * Note: We can only guarentee that the vertices sit on
1612 * the curve, so we must test with their position instead
1613 * of the cell centre.
1616 * const Point<2> centre_2d (0,0);
1617 * for (typename Triangulation<2>::active_cell_iterator
1618 * cell = tria_2d_not_flat.begin_active();
1619 * cell != tria_2d_not_flat.end(); ++cell)
1621 * for (unsigned int face=0; face<GeometryInfo<2>::faces_per_cell; ++face)
1622 * if (cell->face(face)->at_boundary())
1623 * for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_face; ++vertex)
1624 * if (std::abs(cell->vertex(vertex).distance(centre_2d) - hole_diameter/2.0) < 1e-12)
1626 * cell->face(face)->set_manifold_id(10);
1630 * SphericalManifold<2> spherical_manifold_2d (centre_2d);
1631 * tria_2d_not_flat.set_manifold(10,spherical_manifold_2d);
1632 * tria_2d_not_flat.refine_global(std::max (1U, n_repetitions_xy));
1633 * tria_2d_not_flat.reset_manifold(10); // Clear manifold
1635 * GridGenerator::flatten_triangulation(tria_2d_not_flat,tria_2d);
1637 * template <int dim>
1638 * void Solid<dim>::setup_system(LA::MPI::BlockVector &solution_delta)
1640 * timer.enter_subsection("Setup system");
1641 * pcout << "Setting up linear system..." << std::endl;
1645 * Partition triangulation
1648 * GridTools::partition_triangulation (n_mpi_processes,
1651 * block_component = std::vector<unsigned int> (n_components, u_block); // Displacement
1652 * block_component[p_component] = p_block; // Pressure
1653 * block_component[J_component] = J_block; // Dilatation
1654 * dof_handler.distribute_dofs(fe);
1655 * DoFRenumbering::Cuthill_McKee(dof_handler);
1656 * DoFRenumbering::component_wise(dof_handler, block_component);
1660 * Count DoFs in each block
1663 * dofs_per_block.clear();
1664 * dofs_per_block.resize(n_blocks);
1665 * DoFTools::count_dofs_per_block(dof_handler, dofs_per_block,
1668 * all_locally_owned_dofs = DoFTools::locally_owned_dofs_per_subdomain (dof_handler);
1669 * std::vector<IndexSet> all_locally_relevant_dofs
1670 * = DoFTools::locally_relevant_dofs_per_subdomain (dof_handler);
1672 * locally_owned_dofs.clear();
1673 * locally_owned_partitioning.clear();
1674 * Assert(all_locally_owned_dofs.size() > this_mpi_process, ExcInternalError());
1675 * locally_owned_dofs = all_locally_owned_dofs[this_mpi_process];
1677 * locally_relevant_dofs.clear();
1678 * locally_relevant_partitioning.clear();
1679 * Assert(all_locally_relevant_dofs.size() > this_mpi_process, ExcInternalError());
1680 * locally_relevant_dofs = all_locally_relevant_dofs[this_mpi_process];
1682 * locally_owned_partitioning.reserve(n_blocks);
1683 * locally_relevant_partitioning.reserve(n_blocks);
1684 * for (unsigned int b=0; b<n_blocks; ++b)
1686 * const types::global_dof_index idx_begin
1687 * = std::accumulate(dofs_per_block.begin(),
1688 * std::next(dofs_per_block.begin(),b), 0);
1689 * const types::global_dof_index idx_end
1690 * = std::accumulate(dofs_per_block.begin(),
1691 * std::next(dofs_per_block.begin(),b+1), 0);
1692 * locally_owned_partitioning.push_back(locally_owned_dofs.get_view(idx_begin, idx_end));
1693 * locally_relevant_partitioning.push_back(locally_relevant_dofs.get_view(idx_begin, idx_end));
1697 * << " Number of active cells: " << triangulation.n_active_cells()
1698 * << " (by partition:";
1699 * for (unsigned int p=0; p<n_mpi_processes; ++p)
1701 * << (p==0 ? ' ' : '+
')
1702 * << (GridTools::count_cells_with_subdomain_association (triangulation,p));
1703 * pcout << ")" << std::endl;
1706 * << " Number of degrees of freedom: " << dof_handler.n_dofs()
1707 * << " (by partition:";
1708 * for (unsigned int p=0; p<n_mpi_processes; ++p)
1710 * << (p==0 ? ' ' : '+
')
1711 * << (DoFTools::count_dofs_with_subdomain_association (dof_handler,p));
1712 * pcout << ")" << std::endl;
1714 * << " Number of degrees of freedom per block: "
1715 * << "[n_u, n_p, n_J] = ["
1716 * << dofs_per_block[u_block] << ", "
1717 * << dofs_per_block[p_block] << ", "
1718 * << dofs_per_block[J_block] << "]"
1722 * Table<2, DoFTools::Coupling> coupling(n_components, n_components);
1723 * for (unsigned int ii = 0; ii < n_components; ++ii)
1724 * for (unsigned int jj = 0; jj < n_components; ++jj)
1725 * if (((ii < p_component) && (jj == J_component))
1726 * || ((ii == J_component) && (jj < p_component))
1727 * || ((ii == p_component) && (jj == p_component)))
1728 * coupling[ii][jj] = DoFTools::none;
1730 * coupling[ii][jj] = DoFTools::always;
1732 * TrilinosWrappers::BlockSparsityPattern bsp (locally_owned_partitioning,
1733 * locally_owned_partitioning,
1734 * locally_relevant_partitioning,
1735 * mpi_communicator);
1736 * DoFTools::make_sparsity_pattern (dof_handler, bsp,
1737 * constraints, false,
1738 * this_mpi_process);
1740 * tangent_matrix.reinit (bsp);
1744 * We then set up storage vectors
1747 * system_rhs.reinit(locally_owned_partitioning,
1748 * mpi_communicator);
1749 * solution_n.reinit(locally_owned_partitioning,
1750 * mpi_communicator);
1751 * solution_delta.reinit(locally_owned_partitioning,
1752 * mpi_communicator);
1754 * timer.leave_subsection();
1756 * template <int dim>
1758 * Solid<dim>::determine_component_extractors()
1760 * element_indices_u.clear();
1761 * element_indices_p.clear();
1762 * element_indices_J.clear();
1763 * for (unsigned int k = 0; k < fe.dofs_per_cell; ++k)
1765 * const unsigned int k_group = fe.system_to_base_index(k).first.first;
1766 * if (k_group == u_block)
1767 * element_indices_u.push_back(k);
1768 * else if (k_group == p_block)
1769 * element_indices_p.push_back(k);
1770 * else if (k_group == J_block)
1771 * element_indices_J.push_back(k);
1774 * Assert(k_group <= J_block, ExcInternalError());
1778 * template <int dim>
1779 * void Solid<dim>::setup_qph()
1781 * pcout << "Setting up quadrature point data..." << std::endl;
1782 * quadrature_point_history.initialize(triangulation.begin_active(),
1783 * triangulation.end(),
1785 * FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
1786 * cell (IteratorFilters::SubdomainEqualTo(this_mpi_process),
1787 * dof_handler.begin_active()),
1788 * endc (IteratorFilters::SubdomainEqualTo(this_mpi_process),
1789 * dof_handler.end());
1790 * for (; cell!=endc; ++cell)
1792 * Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());
1793 * const std::vector<std::shared_ptr<PointHistory<dim> > > lqph =
1794 * quadrature_point_history.get_data(cell);
1795 * Assert(lqph.size() == n_q_points, ExcInternalError());
1796 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1797 * lqph[q_point]->setup_lqp(parameters, time);
1800 * template <int dim>
1802 * Solid<dim>::solve_nonlinear_timestep(LA::MPI::BlockVector &solution_delta)
1804 * pcout << std::endl
1805 * << "Timestep " << time.get_timestep() << " @ "
1806 * << time.current() << "s of "
1807 * << time.end() << "s" << std::endl;
1808 * LA::MPI::BlockVector newton_update(locally_owned_partitioning,
1809 * mpi_communicator);
1810 * error_residual.reset();
1811 * error_residual_0.reset();
1812 * error_residual_norm.reset();
1813 * error_update.reset();
1814 * error_update_0.reset();
1815 * error_update_norm.reset();
1816 * print_conv_header();
1817 * unsigned int newton_iteration = 0;
1818 * for (; newton_iteration < parameters.max_iterations_NR;
1819 * ++newton_iteration)
1821 * pcout << " " << std::setw(2) << newton_iteration << " " << std::flush;
1822 * make_constraints(newton_iteration);
1823 * assemble_system(solution_delta);
1824 * get_error_residual(error_residual);
1825 * if (newton_iteration == 0)
1826 * error_residual_0 = error_residual;
1827 * error_residual_norm = error_residual;
1828 * error_residual_norm.normalise(error_residual_0);
1829 * if (newton_iteration > 0 &&
1830 * (error_update_norm.u <= parameters.tol_u &&
1831 * error_residual_norm.u <= parameters.tol_f) )
1833 * pcout << " CONVERGED! " << std::endl;
1834 * print_conv_footer(solution_delta);
1837 * const std::pair<unsigned int, double>
1838 * lin_solver_output = solve_linear_system(newton_update);
1839 * get_error_update(newton_update, error_update);
1840 * if (newton_iteration == 0)
1841 * error_update_0 = error_update;
1842 * error_update_norm = error_update;
1843 * error_update_norm.normalise(error_update_0);
1844 * solution_delta += newton_update;
1845 * newton_update = 0.0;
1846 * pcout << " | " << std::fixed << std::setprecision(3) << std::setw(7)
1847 * << std::scientific << lin_solver_output.first << " "
1848 * << lin_solver_output.second << " " << error_residual_norm.norm
1849 * << " " << error_residual_norm.u << " "
1850 * << error_residual_norm.p << " " << error_residual_norm.J
1851 * << " " << error_update_norm.norm << " " << error_update_norm.u
1852 * << " " << error_update_norm.p << " " << error_update_norm.J
1853 * << " " << std::endl;
1855 * AssertThrow (newton_iteration <= parameters.max_iterations_NR,
1856 * ExcMessage("No convergence in nonlinear solver!"));
1858 * template <int dim>
1859 * void Solid<dim>::print_conv_header()
1861 * pcout << std::string(132,'_
') << std::endl;
1862 * pcout << " SOLVER STEP "
1863 * << " | LIN_IT LIN_RES RES_NORM "
1864 * << " RES_U RES_P RES_J NU_NORM "
1865 * << " NU_U NU_P NU_J " << std::endl;
1866 * pcout << std::string(132,'_
') << std::endl;
1868 * template <int dim>
1869 * void Solid<dim>::print_conv_footer(const LA::MPI::BlockVector &solution_delta)
1871 * pcout << std::string(132,'_
') << std::endl;
1872 * const std::pair<double,std::pair<double,double> > error_dil = get_error_dilation(get_solution_total(solution_delta));
1873 * pcout << "Relative errors:" << std::endl
1874 * << "Displacement:\t" << error_update.u / error_update_0.u << std::endl
1875 * << "Force: \t\t" << error_residual.u / error_residual_0.u << std::endl
1876 * << "Dilatation:\t" << error_dil.first << std::endl
1877 * << "v / V_0:\t" << error_dil.second.second << " / " << error_dil.second.first
1878 * << " = " << (error_dil.second.second/error_dil.second.first) << std::endl;
1880 * template <int dim>
1881 * std::pair<double,std::pair<double,double> >
1882 * Solid<dim>::get_error_dilation(const LA::MPI::BlockVector &solution_total) const
1884 * double vol_reference = 0.0;
1885 * double vol_current = 0.0;
1886 * double dil_L2_error = 0.0;
1887 * FEValues<dim> fe_values_ref(fe, qf_cell,
1888 * update_values | update_gradients | update_JxW_values);
1889 * std::vector<Tensor<2, dim> > solution_grads_u_total (qf_cell.size());
1890 * std::vector<double> solution_values_J_total (qf_cell.size());
1891 * FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
1892 * cell (IteratorFilters::SubdomainEqualTo(this_mpi_process),
1893 * dof_handler.begin_active()),
1894 * endc (IteratorFilters::SubdomainEqualTo(this_mpi_process),
1895 * dof_handler.end());
1896 * for (; cell != endc; ++cell)
1898 * Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());
1899 * fe_values_ref.reinit(cell);
1900 * fe_values_ref[u_fe].get_function_gradients(solution_total,
1901 * solution_grads_u_total);
1902 * fe_values_ref[J_fe].get_function_values(solution_total,
1903 * solution_values_J_total);
1904 * const std::vector<std::shared_ptr<const PointHistory<dim> > > lqph =
1905 * quadrature_point_history.get_data(cell);
1906 * Assert(lqph.size() == n_q_points, ExcInternalError());
1907 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1909 * const double det_F_qp = determinant(Physics::Elasticity::Kinematics::F(solution_grads_u_total[q_point]));
1910 * const double J_tilde_qp = solution_values_J_total[q_point];
1911 * const double the_error_qp_squared = std::pow((det_F_qp - J_tilde_qp),
1913 * const double JxW = fe_values_ref.JxW(q_point);
1914 * dil_L2_error += the_error_qp_squared * JxW;
1915 * vol_reference += JxW;
1916 * vol_current += det_F_qp * JxW;
1919 * Assert(vol_current > 0.0, ExcInternalError());
1922 * Sum across all porcessors
1925 * dil_L2_error = Utilities::MPI::sum(dil_L2_error,mpi_communicator);
1926 * vol_reference = Utilities::MPI::sum(vol_reference,mpi_communicator);
1927 * vol_current = Utilities::MPI::sum(vol_current,mpi_communicator);
1929 * return std::make_pair(std::sqrt(dil_L2_error),
1930 * std::make_pair(vol_reference,vol_current));
1932 * template <int dim>
1933 * void Solid<dim>::get_error_residual(Errors &error_residual)
1937 * Construct a residual vector that has the values for all of its
1938 * constrained DoFs set to zero.
1941 * LA::MPI::BlockVector error_res (system_rhs);
1942 * constraints.set_zero(error_res);
1943 * error_residual.norm = error_res.l2_norm();
1944 * error_residual.u = error_res.block(u_block).l2_norm();
1945 * error_residual.p = error_res.block(p_block).l2_norm();
1946 * error_residual.J = error_res.block(J_block).l2_norm();
1948 * template <int dim>
1949 * void Solid<dim>::get_error_update(const LA::MPI::BlockVector &newton_update,
1950 * Errors &error_update)
1954 * Construct a update vector that has the values for all of its
1955 * constrained DoFs set to zero.
1958 * LA::MPI::BlockVector error_ud (newton_update);
1959 * constraints.set_zero(error_ud);
1960 * error_update.norm = error_ud.l2_norm();
1961 * error_update.u = error_ud.block(u_block).l2_norm();
1962 * error_update.p = error_ud.block(p_block).l2_norm();
1963 * error_update.J = error_ud.block(J_block).l2_norm();
1965 * template <int dim>
1966 * LA::MPI::BlockVector
1967 * Solid<dim>::get_solution_total(const LA::MPI::BlockVector &solution_delta) const
1971 * Cell interpolation -> Ghosted vector
1974 * LA::MPI::BlockVector solution_total (locally_owned_partitioning,
1975 * locally_relevant_partitioning,
1977 * /*vector_writable = */ false);
1978 * LA::MPI::BlockVector tmp (solution_total);
1979 * solution_total = solution_n;
1980 * tmp = solution_delta;
1981 * solution_total += tmp;
1982 * return solution_total;
1984 * template <int dim>
1985 * void Solid<dim>::assemble_system(const LA::MPI::BlockVector &solution_delta)
1987 * timer.enter_subsection("Assemble system");
1988 * pcout << " ASM_SYS " << std::flush;
1989 * tangent_matrix = 0.0;
1991 * const LA::MPI::BlockVector solution_total(get_solution_total(solution_delta));
1992 * const UpdateFlags uf_cell(update_values |
1993 * update_gradients |
1994 * update_JxW_values);
1995 * const UpdateFlags uf_face(update_values |
1996 * update_normal_vectors |
1997 * update_JxW_values);
1998 * PerTaskData_ASM per_task_data(dofs_per_cell);
1999 * ScratchData_ASM scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face, solution_total);
2001 * FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
2002 * cell (IteratorFilters::SubdomainEqualTo(this_mpi_process),
2003 * dof_handler.begin_active()),
2004 * endc (IteratorFilters::SubdomainEqualTo(this_mpi_process),
2005 * dof_handler.end());
2006 * for (; cell != endc; ++cell)
2008 * Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());
2009 * assemble_system_one_cell(cell, scratch_data, per_task_data);
2010 * copy_local_to_global_system(per_task_data);
2012 * tangent_matrix.compress(VectorOperation::add);
2013 * system_rhs.compress(VectorOperation::add);
2014 * timer.leave_subsection();
2016 * template <int dim>
2017 * void Solid<dim>::copy_local_to_global_system(const PerTaskData_ASM &data)
2019 * constraints.distribute_local_to_global(data.cell_matrix, data.cell_rhs,
2020 * data.local_dof_indices,
2021 * tangent_matrix, system_rhs);
2023 * template <int dim>
2025 * Solid<dim>::assemble_system_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2026 * ScratchData_ASM &scratch,
2027 * PerTaskData_ASM &data) const
2029 * Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());
2033 * scratch.fe_values_ref.reinit(cell);
2034 * cell->get_dof_indices(data.local_dof_indices);
2035 * const std::vector<std::shared_ptr<const PointHistory<dim> > > lqph =
2036 * quadrature_point_history.get_data(cell);
2037 * Assert(lqph.size() == n_q_points, ExcInternalError());
2041 * Update quadrature point solution
2044 * scratch.fe_values_ref[u_fe].get_function_gradients(scratch.solution_total,
2045 * scratch.solution_grads_u_total);
2046 * scratch.fe_values_ref[p_fe].get_function_values(scratch.solution_total,
2047 * scratch.solution_values_p_total);
2048 * scratch.fe_values_ref[J_fe].get_function_values(scratch.solution_total,
2049 * scratch.solution_values_J_total);
2053 * Update shape functions and their gradients (push-forward)
2056 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2058 * const Tensor<2, dim> F = Physics::Elasticity::Kinematics::F(scratch.solution_grads_u_total[q_point]);
2059 * const Tensor<2, dim> F_inv = invert(F);
2061 * for (unsigned int k = 0; k < dofs_per_cell; ++k)
2063 * const unsigned int k_group = fe.system_to_base_index(k).first.first;
2064 * if (k_group == u_block)
2066 * scratch.grad_Nx[q_point][k] = scratch.fe_values_ref[u_fe].gradient(k, q_point)
2068 * scratch.symm_grad_Nx[q_point][k] = symmetrize(scratch.grad_Nx[q_point][k]);
2070 * else if (k_group == p_block)
2071 * scratch.Nx[q_point][k] = scratch.fe_values_ref[p_fe].value(k,
2073 * else if (k_group == J_block)
2074 * scratch.Nx[q_point][k] = scratch.fe_values_ref[J_fe].value(k,
2077 * Assert(k_group <= J_block, ExcInternalError());
2080 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2082 * const SymmetricTensor<2, dim> &I = Physics::Elasticity::StandardTensors<dim>::I;
2083 * const Tensor<2, dim> F = Physics::Elasticity::Kinematics::F(scratch.solution_grads_u_total[q_point]);
2084 * const double det_F = determinant(F);
2085 * const double &p_tilde = scratch.solution_values_p_total[q_point];
2086 * const double &J_tilde = scratch.solution_values_J_total[q_point];
2087 * Assert(det_F > 0, ExcInternalError());
2090 * PointHistory<dim> *lqph_q_point_nc = const_cast<PointHistory<dim>*>(lqph[q_point].get());
2091 * lqph_q_point_nc->update_internal_equilibrium(F,p_tilde,J_tilde);
2094 * const SymmetricTensor<2, dim> tau = lqph[q_point]->get_tau(F,p_tilde);
2095 * const Tensor<2, dim> tau_ns (tau);
2096 * const SymmetricTensor<4, dim> Jc = lqph[q_point]->get_Jc(F,p_tilde);
2097 * const double dPsi_vol_dJ = lqph[q_point]->get_dPsi_vol_dJ(J_tilde);
2098 * const double d2Psi_vol_dJ2 = lqph[q_point]->get_d2Psi_vol_dJ2(J_tilde);
2100 * const std::vector<double> &Nx = scratch.Nx[q_point];
2101 * const std::vector<Tensor<2, dim> > &grad_Nx = scratch.grad_Nx[q_point];
2102 * const std::vector<SymmetricTensor<2, dim> > &symm_grad_Nx = scratch.symm_grad_Nx[q_point];
2103 * const double JxW = scratch.fe_values_ref.JxW(q_point);
2105 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2107 * const unsigned int component_i = fe.system_to_component_index(i).first;
2108 * const unsigned int i_group = fe.system_to_base_index(i).first.first;
2109 * if (i_group == u_block)
2110 * data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW;
2111 * else if (i_group == p_block)
2112 * data.cell_rhs(i) -= Nx[i] * (det_F - J_tilde) * JxW;
2113 * else if (i_group == J_block)
2114 * data.cell_rhs(i) -= Nx[i] * (dPsi_vol_dJ - p_tilde) * JxW;
2116 * Assert(i_group <= J_block, ExcInternalError());
2118 * for (unsigned int j = 0; j <= i; ++j)
2120 * const unsigned int component_j = fe.system_to_component_index(j).first;
2121 * const unsigned int j_group = fe.system_to_base_index(j).first.first;
2122 * if ((i_group == u_block) && (j_group == u_block))
2124 * data.cell_matrix(i, j) += symm_grad_Nx[i] * Jc // The material contribution:
2125 * * symm_grad_Nx[j] * JxW;
2126 * if (component_i == component_j) // geometrical stress contribution
2127 * data.cell_matrix(i, j) += grad_Nx[i][component_i] * tau_ns
2128 * * grad_Nx[j][component_j] * JxW;
2130 * else if ((i_group == u_block) && (j_group == p_block))
2132 * data.cell_matrix(i, j) += (symm_grad_Nx[i] * I)
2136 * else if ((i_group == p_block) && (j_group == u_block))
2138 * data.cell_matrix(i, j) += Nx[i] * det_F
2139 * * (symm_grad_Nx[j] * I)
2142 * else if ((i_group == p_block) && (j_group == J_block))
2143 * data.cell_matrix(i, j) -= Nx[i] * Nx[j] * JxW;
2144 * else if ((i_group == J_block) && (j_group == p_block))
2145 * data.cell_matrix(i, j) -= Nx[i] * Nx[j] * JxW;
2146 * else if ((i_group == J_block) && (j_group == J_block))
2147 * data.cell_matrix(i, j) += Nx[i] * d2Psi_vol_dJ2 * Nx[j] * JxW;
2149 * Assert((i_group <= J_block) && (j_group <= J_block),
2150 * ExcInternalError());
2155 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2156 * for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
2157 * data.cell_matrix(i, j) = data.cell_matrix(j, i);
2159 * if (parameters.driver == "Neumann")
2160 * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
2162 * if (cell->face(face)->at_boundary() == true
2163 * && cell->face(face)->boundary_id() == parameters.boundary_id_plus_Y)
2165 * scratch.fe_face_values_ref.reinit(cell, face);
2166 * for (unsigned int f_q_point = 0; f_q_point < n_q_points_f;
2169 * const Tensor<1, dim> &N =
2170 * scratch.fe_face_values_ref.normal_vector(f_q_point);
2171 * static const double pressure_nom = parameters.pressure
2172 * / (parameters.scale * parameters.scale);
2173 * const double time_ramp = (time.current() < parameters.load_time ?
2174 * time.current() / parameters.load_time : 1.0);
2175 * const double pressure = -pressure_nom * time_ramp;
2176 * const Tensor<1, dim> traction = pressure * N;
2177 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2179 * const unsigned int i_group =
2180 * fe.system_to_base_index(i).first.first;
2181 * if (i_group == u_block)
2183 * const unsigned int component_i =
2184 * fe.system_to_component_index(i).first;
2186 * scratch.fe_face_values_ref.shape_value(i,
2188 * const double JxW = scratch.fe_face_values_ref.JxW(
2190 * data.cell_rhs(i) += (Ni * traction[component_i])
2197 * template <int dim>
2198 * void Solid<dim>::make_constraints(const int &it_nr)
2200 * pcout << " CST " << std::flush;
2203 * constraints.clear();
2204 * const bool apply_dirichlet_bc = (it_nr == 0);
2205 * const FEValuesExtractors::Scalar x_displacement(0);
2206 * const FEValuesExtractors::Scalar y_displacement(1);
2208 * const int boundary_id = parameters.boundary_id_minus_X;
2209 * if (apply_dirichlet_bc == true)
2210 * VectorTools::interpolate_boundary_values(dof_handler,
2212 * ZeroFunction<dim>(n_components),
2214 * fe.component_mask(x_displacement));
2216 * VectorTools::interpolate_boundary_values(dof_handler,
2218 * ZeroFunction<dim>(n_components),
2220 * fe.component_mask(x_displacement));
2223 * const int boundary_id = parameters.boundary_id_minus_Y;
2224 * if (apply_dirichlet_bc == true)
2225 * VectorTools::interpolate_boundary_values(dof_handler,
2227 * ZeroFunction<dim>(n_components),
2229 * fe.component_mask(y_displacement));
2231 * VectorTools::interpolate_boundary_values(dof_handler,
2233 * ZeroFunction<dim>(n_components),
2235 * fe.component_mask(y_displacement));
2239 * const FEValuesExtractors::Scalar z_displacement(2);
2241 * const int boundary_id = parameters.boundary_id_minus_Z;
2242 * if (apply_dirichlet_bc == true)
2243 * VectorTools::interpolate_boundary_values(dof_handler,
2245 * ZeroFunction<dim>(n_components),
2247 * fe.component_mask(z_displacement));
2249 * VectorTools::interpolate_boundary_values(dof_handler,
2251 * ZeroFunction<dim>(n_components),
2253 * fe.component_mask(z_displacement));
2256 * const int boundary_id = parameters.boundary_id_plus_Z;
2257 * if (apply_dirichlet_bc == true)
2258 * VectorTools::interpolate_boundary_values(dof_handler,
2260 * ZeroFunction<dim>(n_components),
2262 * fe.component_mask(z_displacement));
2264 * VectorTools::interpolate_boundary_values(dof_handler,
2266 * ZeroFunction<dim>(n_components),
2268 * fe.component_mask(z_displacement));
2271 * if (parameters.driver == "Dirichlet")
2273 * const int boundary_id = parameters.boundary_id_plus_Y;
2274 * if (apply_dirichlet_bc == true)
2277 * if (time.current() < parameters.load_time+0.01*time.get_delta_t())
2279 * const double delta_length = parameters.length*(parameters.stretch - 1.0)*parameters.scale;
2280 * const unsigned int n_stretch_steps = parameters.load_time/time.get_delta_t();
2281 * const double delta_u_y = delta_length/2.0/n_stretch_steps;
2282 * VectorTools::interpolate_boundary_values(dof_handler,
2284 * ConstantFunction<dim>(delta_u_y,n_components),
2286 * fe.component_mask(y_displacement));
2289 * VectorTools::interpolate_boundary_values(dof_handler,
2291 * ZeroFunction<dim>(n_components),
2293 * fe.component_mask(y_displacement));
2296 * VectorTools::interpolate_boundary_values(dof_handler,
2298 * ZeroFunction<dim>(n_components),
2300 * fe.component_mask(y_displacement));
2302 * constraints.close();
2304 * template <int dim>
2305 * std::pair<unsigned int, double>
2306 * Solid<dim>::solve_linear_system(LA::MPI::BlockVector &newton_update)
2308 * unsigned int lin_it = 0;
2309 * double lin_res = 0.0;
2311 * timer.enter_subsection("Linear solver");
2312 * pcout << " SLV " << std::flush;
2314 * const LA::MPI::Vector &f_u = system_rhs.block(u_block);
2315 * const LA::MPI::Vector &f_p = system_rhs.block(p_block);
2316 * const LA::MPI::Vector &f_J = system_rhs.block(J_block);
2317 * LA::MPI::Vector &d_u = newton_update.block(u_block);
2318 * LA::MPI::Vector &d_p = newton_update.block(p_block);
2319 * LA::MPI::Vector &d_J = newton_update.block(J_block);
2320 * const auto K_uu = linear_operator<LA::MPI::Vector>(tangent_matrix.block(u_block, u_block));
2321 * const auto K_up = linear_operator<LA::MPI::Vector>(tangent_matrix.block(u_block, p_block));
2322 * const auto K_pu = linear_operator<LA::MPI::Vector>(tangent_matrix.block(p_block, u_block));
2323 * const auto K_Jp = linear_operator<LA::MPI::Vector>(tangent_matrix.block(J_block, p_block));
2324 * const auto K_JJ = linear_operator<LA::MPI::Vector>(tangent_matrix.block(J_block, J_block));
2326 * LA::PreconditionJacobi preconditioner_K_Jp_inv;
2327 * preconditioner_K_Jp_inv.initialize(
2328 * tangent_matrix.block(J_block, p_block),
2329 * LA::PreconditionJacobi::AdditionalData());
2330 * ReductionControl solver_control_K_Jp_inv (
2331 * tangent_matrix.block(J_block, p_block).m() * parameters.max_iterations_lin,
2333 * ::SolverCG<LA::MPI::Vector> solver_K_Jp_inv (solver_control_K_Jp_inv);
2335 * const auto K_Jp_inv = inverse_operator(K_Jp,
2337 * preconditioner_K_Jp_inv);
2338 * const auto K_pJ_inv = transpose_operator(K_Jp_inv);
2339 * const auto K_pp_bar = K_Jp_inv * K_JJ * K_pJ_inv;
2340 * const auto K_uu_bar_bar = K_up * K_pp_bar * K_pu;
2341 * const auto K_uu_con = K_uu + K_uu_bar_bar;
2343 * LA::PreconditionAMG preconditioner_K_con_inv;
2344 * preconditioner_K_con_inv.initialize(
2345 * tangent_matrix.block(u_block, u_block),
2346 * LA::PreconditionAMG::AdditionalData(
2347 * true /*elliptic*/,
2348 * (parameters.poly_degree > 1 /*higher_order_elements*/)) );
2349 * ReductionControl solver_control_K_con_inv (
2350 * tangent_matrix.block(u_block, u_block).m() * parameters.max_iterations_lin,
2351 * 1.0e-30, parameters.tol_lin);
2352 * ::SolverSelector<LA::MPI::Vector> solver_K_con_inv;
2353 * solver_K_con_inv.select(parameters.type_lin);
2354 * solver_K_con_inv.set_control(solver_control_K_con_inv);
2355 * const auto K_uu_con_inv = inverse_operator(K_uu_con,
2357 * preconditioner_K_con_inv);
2359 * d_u = K_uu_con_inv*(f_u - K_up*(K_Jp_inv*f_J - K_pp_bar*f_p));
2360 * lin_it = solver_control_K_con_inv.last_step();
2361 * lin_res = solver_control_K_con_inv.last_value();
2362 * timer.leave_subsection();
2364 * timer.enter_subsection("Linear solver postprocessing");
2365 * d_J = K_pJ_inv*(f_p - K_pu*d_u);
2366 * d_p = K_Jp_inv*(f_J - K_JJ*d_J);
2367 * timer.leave_subsection();
2369 * constraints.distribute(newton_update);
2370 * return std::make_pair(lin_it, lin_res);
2372 * template <int dim>
2374 * Solid<dim>::update_end_timestep ()
2376 * FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
2377 * cell (IteratorFilters::SubdomainEqualTo(this_mpi_process),
2378 * dof_handler.begin_active()),
2379 * endc (IteratorFilters::SubdomainEqualTo(this_mpi_process),
2380 * dof_handler.end());
2381 * for (; cell != endc; ++cell)
2383 * Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());
2384 * const std::vector<std::shared_ptr<PointHistory<dim> > > lqph =
2385 * quadrature_point_history.get_data(cell);
2386 * Assert(lqph.size() == n_q_points, ExcInternalError());
2387 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2388 * lqph[q_point]->update_end_timestep();
2392 * template<int dim, class DH=DoFHandler<dim> >
2393 * class FilteredDataOut : public DataOut<dim, DH>
2396 * FilteredDataOut (const unsigned int subdomain_id)
2398 * subdomain_id (subdomain_id)
2401 * virtual ~FilteredDataOut() {}
2403 * virtual typename DataOut<dim, DH>::cell_iterator
2406 * typename DataOut<dim, DH>::active_cell_iterator
2407 * cell = this->dofs->begin_active();
2408 * while ((cell != this->dofs->end()) &&
2409 * (cell->subdomain_id() != subdomain_id))
2414 * virtual typename DataOut<dim, DH>::cell_iterator
2415 * next_cell (const typename DataOut<dim, DH>::cell_iterator &old_cell)
2417 * if (old_cell != this->dofs->end())
2419 * const IteratorFilters::SubdomainEqualTo predicate(subdomain_id);
2421 * ++(FilteredIterator
2422 * <typename DataOut<dim, DH>::active_cell_iterator>
2423 * (predicate,old_cell));
2430 * const unsigned int subdomain_id;
2433 * template <int dim>
2434 * void Solid<dim>::output_results(const unsigned int timestep,
2435 * const double current_time) const
2439 * Output -> Ghosted vector
2442 * LA::MPI::BlockVector solution_total (locally_owned_partitioning,
2443 * locally_relevant_partitioning,
2445 * /*vector_writable = */ false);
2446 * LA::MPI::BlockVector residual (locally_owned_partitioning,
2447 * locally_relevant_partitioning,
2449 * /*vector_writable = */ false);
2450 * solution_total = solution_n;
2451 * residual = system_rhs;
2456 * --- Additional data ---
2459 * Vector<double> material_id;
2460 * Vector<double> polynomial_order;
2461 * material_id.reinit(triangulation.n_active_cells());
2462 * polynomial_order.reinit(triangulation.n_active_cells());
2463 * std::vector<types::subdomain_id> partition_int (triangulation.n_active_cells());
2465 * FilteredDataOut<dim> data_out(this_mpi_process);
2466 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
2467 * data_component_interpretation(dim,
2468 * DataComponentInterpretation::component_is_part_of_vector);
2469 * data_component_interpretation.push_back(DataComponentInterpretation::component_is_scalar);
2470 * data_component_interpretation.push_back(DataComponentInterpretation::component_is_scalar);
2472 * GridTools::get_subdomain_association (triangulation, partition_int);
2476 * Can't use filtered iterators here because the cell
2477 * count
"c" is incorrect
for the
parallel case
2480 *
unsigned int c = 0;
2483 * endc = dof_handler.
end();
2484 *
for (; cell!=endc; ++cell, ++c)
2488 *
material_id(c) =
static_cast<int>(cell->material_id());
2491 * std::vector<std::string> solution_name(
n_components,
"solution_");
2492 * std::vector<std::string> residual_name(
n_components,
"residual_");
2495 *
if (block_component[c] == u_block)
2497 * solution_name[c] +=
"u";
2498 * residual_name[c] +=
"u";
2500 *
else if (block_component[c] == p_block)
2502 * solution_name[c] +=
"p";
2503 * residual_name[c] +=
"p";
2505 *
else if (block_component[c] == J_block)
2507 * solution_name[c] +=
"J";
2508 * residual_name[c] +=
"J";
2516 * data_out.attach_dof_handler(dof_handler);
2517 * data_out.add_data_vector(solution_total,
2520 * data_component_interpretation);
2521 * data_out.add_data_vector(residual,
2524 * data_component_interpretation);
2526 * partition_int.end());
2527 * data_out.add_data_vector (
material_id,
"material_id");
2528 * data_out.add_data_vector (partitioning,
"partitioning");
2529 * data_out.build_patches(degree);
2533 *
static std::string get_filename_vtu (
unsigned int process,
2534 *
unsigned int timestep,
2535 *
const unsigned int n_digits = 4)
2537 * std::ostringstream filename_vtu;
2546 *
return filename_vtu.str();
2549 *
static std::string get_filename_pvtu (
unsigned int timestep,
2550 *
const unsigned int n_digits = 4)
2552 * std::ostringstream filename_vtu;
2559 *
return filename_vtu.str();
2562 *
static std::string get_filename_pvd (
void)
2564 * std::ostringstream filename_vtu;
2569 *
return filename_vtu.str();
2575 * Write out main data file
2578 *
const std::string filename_vtu = Filename::get_filename_vtu(
this_mpi_process, timestep);
2579 * std::ofstream output(filename_vtu.c_str());
2580 * data_out.write_vtu(output);
2584 * Collection of files written in
parallel
2585 * This next
set of steps should only be performed
2593 * List of all files written out at
this timestep by all processors
2596 * std::vector<std::string> parallel_filenames_vtu;
2599 * parallel_filenames_vtu.push_back(Filename::get_filename_vtu(p, timestep));
2602 *
const std::string filename_pvtu (Filename::get_filename_pvtu(timestep));
2603 * std::ofstream pvtu_master(filename_pvtu.c_str());
2604 * data_out.write_pvtu_record(pvtu_master,
2605 * parallel_filenames_vtu);
2609 * Time dependent data master file
2612 *
static std::vector<std::pair<double,std::string> > time_and_name_history;
2613 * time_and_name_history.push_back (std::make_pair (current_time,
2615 *
const std::string filename_pvd (Filename::get_filename_pvd());
2616 * std::ofstream pvd_output (filename_pvd.c_str());
2620 *
template <
int dim>
2621 *
void Solid<dim>::compute_vertex_positions(std::vector<double> &real_time,
2622 * std::vector<std::vector<
Point<dim> > > &tracked_vertices,
2625 * real_time.push_back(time.current());
2627 * std::vector<bool> vertex_found (tracked_vertices.size(),
false);
2628 * std::vector<Tensor<1,dim> > vertex_update (tracked_vertices.size());
2634 * dof_handler.
end());
2635 *
for (; cell != endc; ++cell)
2638 *
for (
unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
2640 *
for (
unsigned int p=0; p<tracked_vertices.size(); ++p)
2642 *
if (vertex_found[p] ==
true)
continue;
2644 *
const Point<dim> pt_ref = tracked_vertices[p][0];
2645 *
if (cell->vertex(v).distance(pt_ref) < 1
e-6*parameters.scale)
2647 *
for (
unsigned int d=0;
d<dim; ++
d)
2648 * vertex_update[p][
d] = solution_total(cell->vertex_dof_index(v,u_block+
d));
2650 * vertex_found[p] =
true;
2656 *
for (
unsigned int p=0; p<tracked_vertices.size(); ++p)
2658 *
const int found_on_n_processes =
Utilities::MPI::sum(
int(vertex_found[p]), mpi_communicator);
2659 *
Assert(found_on_n_processes>0,
ExcMessage(
"Vertex not found on any processor"));
2661 *
for (
unsigned int d=0;
d<dim; ++
d)
2663 * update /= found_on_n_processes;
2664 * tracked_vertices[p].push_back(tracked_vertices[p][0] + update);
2669 *
int main (
int argc,
char *argv[])
2671 *
using namespace dealii;
2672 *
using namespace ViscoElasStripHole;
2678 *
const unsigned int dim = 2;
2679 * Solid<dim> solid(
"parameters.prm");
2682 *
catch (std::exception &exc)
2686 * std::cerr << std::endl << std::endl
2687 * <<
"----------------------------------------------------"
2689 * std::cerr <<
"Exception on processing: " << std::endl << exc.what()
2690 * << std::endl <<
"Aborting!" << std::endl
2691 * <<
"----------------------------------------------------"
2700 * std::cerr << std::endl << std::endl
2701 * <<
"----------------------------------------------------"
2703 * std::cerr <<
"Unknown exception!" << std::endl <<
"Aborting!"
2705 * <<
"----------------------------------------------------"