318 *
prm.enter_subsection(
"Boundary conditions");
320 *
prm.declare_entry(
"Driver",
"Dirichlet",
322 *
"Driver boundary condition for the problem");
323 *
prm.declare_entry(
"Final stretch",
"2.0",
325 *
"Positive stretch applied length-ways to the strip");
326 *
prm.declare_entry(
"Applied pressure",
"0.0",
328 *
"Hydrostatic pressure applied (in the referential configuration) to the interior surface of the hole");
329 *
prm.declare_entry(
"Load time",
"2.5",
331 *
"Total time over which the stretch/pressure is ramped up");
333 *
prm.leave_subsection();
337 *
prm.enter_subsection(
"Boundary conditions");
339 *
driver = prm.get(
"Driver");
340 *
stretch = prm.get_double(
"Final stretch");
341 *
pressure = prm.get_double(
"Applied pressure");
342 *
load_time = prm.get_double(
"Load time");
344 *
prm.leave_subsection();
357 *
prm.enter_subsection(
"Finite element system");
359 *
prm.declare_entry(
"Polynomial degree",
"2",
361 *
"Displacement system polynomial order");
362 *
prm.declare_entry(
"Quadrature order",
"3",
364 *
"Gauss quadrature order");
366 *
prm.leave_subsection();
370 *
prm.enter_subsection(
"Finite element system");
372 *
poly_degree = prm.get_integer(
"Polynomial degree");
373 *
quad_order = prm.get_integer(
"Quadrature order");
375 *
prm.leave_subsection();
395 *
prm.enter_subsection(
"Geometry");
397 *
prm.declare_entry(
"Length",
"100.0",
399 *
"Total sample length");
400 *
prm.declare_entry(
"Width",
"50.0",
402 *
"Total sample width");
403 *
prm.declare_entry(
"Thickness",
"5.0",
405 *
"Total sample thickness");
406 *
prm.declare_entry(
"Hole diameter",
"20.0",
409 *
prm.declare_entry(
"Hole division fraction",
"0.5",
411 *
"A geometric factor affecting the discretisation near the hole");
412 *
prm.declare_entry(
"Number of subdivisions in cross-section",
"2",
414 *
"A factor defining the number of initial grid subdivisions in the cross-section");
415 *
prm.declare_entry(
"Number of subdivisions thickness",
"6",
417 *
"A factor defining the number of initial grid subdivisions through the thickness");
418 *
prm.declare_entry(
"Global refinement",
"2",
420 *
"Global refinement level");
421 *
prm.declare_entry(
"Grid scale",
"1e-3",
423 *
"Global grid scaling factor");
425 *
prm.leave_subsection();
429 *
prm.enter_subsection(
"Geometry");
431 *
length = prm.get_double(
"Length");
432 *
width = prm.get_double(
"Width");
433 *
thickness = prm.get_double(
"Thickness");
436 *
n_repetitions_xy = prm.get_integer(
"Number of subdivisions in cross-section");
437 *
n_repetitions_z = prm.get_integer(
"Number of subdivisions thickness");
439 *
scale = prm.get_double(
"Grid scale");
441 *
prm.leave_subsection();
456 *
prm.enter_subsection(
"Material properties");
458 *
prm.declare_entry(
"Poisson's ratio",
"0.4999",
460 *
"Poisson's ratio");
461 *
prm.declare_entry(
"Elastic shear modulus",
"80.194e6",
463 *
"Elastic shear modulus");
464 *
prm.declare_entry(
"Viscous shear modulus",
"80.194e6",
466 *
"Viscous shear modulus");
467 *
prm.declare_entry(
"Viscous relaxation time",
"2.0",
469 *
"Viscous relaxation time");
471 *
prm.leave_subsection();
475 *
prm.enter_subsection(
"Material properties");
477 *
nu_e = prm.get_double(
"Poisson's ratio");
478 *
mu_e = prm.get_double(
"Elastic shear modulus");
479 *
mu_v = prm.get_double(
"Viscous shear modulus");
480 *
tau_v = prm.get_double(
"Viscous relaxation time");
482 *
prm.leave_subsection();
496 *
prm.enter_subsection(
"Linear solver");
498 *
prm.declare_entry(
"Solver type",
"cg",
500 *
"Type of solver used to solve the linear system");
501 *
prm.declare_entry(
"Residual",
"1e-6",
503 *
"Linear solver residual (scaled by residual norm)");
504 *
prm.declare_entry(
"Max iteration multiplier",
"1",
506 *
"Linear solver iterations (multiples of the system matrix size)");
508 *
prm.leave_subsection();
512 *
prm.enter_subsection(
"Linear solver");
514 *
type_lin = prm.get(
"Solver type");
515 *
tol_lin = prm.get_double(
"Residual");
518 *
prm.leave_subsection();
520 *
struct NonlinearSolver
532 *
prm.enter_subsection(
"Nonlinear solver");
534 *
prm.declare_entry(
"Max iterations Newton-Raphson",
"10",
536 *
"Number of Newton-Raphson iterations allowed");
537 *
prm.declare_entry(
"Tolerance displacement",
"1.0e-6",
539 *
"Displacement error tolerance");
540 *
prm.declare_entry(
"Tolerance force",
"1.0e-9",
542 *
"Force residual tolerance");
544 *
prm.leave_subsection();
548 *
prm.enter_subsection(
"Nonlinear solver");
551 *
tol_f = prm.get_double(
"Tolerance force");
552 *
tol_u = prm.get_double(
"Tolerance displacement");
554 *
prm.leave_subsection();
567 *
prm.enter_subsection(
"Time");
569 *
prm.declare_entry(
"End time",
"1",
572 *
prm.declare_entry(
"Time step size",
"0.1",
576 *
prm.leave_subsection();
580 *
prm.enter_subsection(
"Time");
582 *
end_time = prm.get_double(
"End time");
583 *
delta_t = prm.get_double(
"Time step size");
585 *
prm.leave_subsection();
593 *
public NonlinearSolver,
605 *
declare_parameters(prm);
607 *
parse_parameters(prm);
611 *
BoundaryConditions::declare_parameters(prm);
612 *
FESystem::declare_parameters(prm);
613 *
Geometry::declare_parameters(prm);
614 *
Materials::declare_parameters(prm);
615 *
LinearSolver::declare_parameters(prm);
616 *
NonlinearSolver::declare_parameters(prm);
617 *
Time::declare_parameters(prm);
621 *
BoundaryConditions::parse_parameters(prm);
622 *
FESystem::parse_parameters(prm);
623 *
Geometry::parse_parameters(prm);
624 *
Materials::parse_parameters(prm);
625 *
LinearSolver::parse_parameters(prm);
626 *
NonlinearSolver::parse_parameters(prm);
627 *
Time::parse_parameters(prm);
634 *
const double delta_t)
668 *
const double delta_t;
677 *
const double tau_v,
736 *
const double kappa;
739 *
const double tau_v;
819 *
parameters.mu_e, parameters.nu_e,
820 *
parameters.mu_v, parameters.tau_v,
859 *
std::shared_ptr< Material_Compressible_Three_Field_Linear_Viscoelastic<dim> >
material;
900 *
std::pair<unsigned int, double>
902 *
LA::MPI::BlockVector
908 *
const double current_time)
const;
930 *
const unsigned int degree;
933 *
const unsigned int dofs_per_cell;
937 *
static const unsigned int n_blocks = 3;
938 *
static const unsigned int n_components = dim + 2;
971 *
const unsigned int n_q_points;
981 *
norm(1.0),
u(1.0), p(1.0), J(1.0)
1001 *
double norm,
u, p, J;
1010 *
std::pair<double, std::pair<double,double> >
1017 *
template <
int dim>
1023 *
pcout(std::cout, this_mpi_process == 0),
1026 *
time(parameters.end_time, parameters.delta_t),
1027 *
timer(mpi_communicator,
1031 *
degree(parameters.poly_degree),
1032 *
fe(
FE_Q<dim>(parameters.poly_degree), dim,
1036 *
dofs_per_cell (fe.dofs_per_cell),
1041 *
qf_cell(parameters.quad_order),
1042 *
qf_face(parameters.quad_order),
1043 *
n_q_points (
qf_cell.size()),
1046 *
Assert(dim==2 || dim==3,
ExcMessage(
"This problem only works in 2 or 3 space dimensions."));
1049 *
template <
int dim>
1052 *
dof_handler.
clear();
1054 *
template <
int dim>
1063 *
constraints.close();
1085 *
p[1] = parameters.length/2.0;
1090 *
p[1] = parameters.hole_diameter/2.0;
1095 *
p[0] = parameters.hole_diameter/2.0;
1100 *
p[0] = parameters.width/2.0;
1104 *
while (time.current() < time.end()+0.01*time.get_delta_t())
1117 *
pcout <<
"\n\n*** Spatial position history for tracked vertices ***" << std::endl;
1118 *
for (
unsigned int t=0; t<
real_time.size(); ++t)
1125 *
for (
unsigned int d=0;
d<dim; ++
d)
1127 *
pcout <<
"Point " << p <<
" [" <<
d <<
"]";
1132 *
pcout << std::endl;
1135 *
pcout << std::setprecision(6);
1140 *
ExcMessage(
"Vertex not tracked at each timestep"));
1141 *
for (
unsigned int d=0;
d<dim; ++
d)
1148 *
pcout << std::endl;
1151 *
template <
int dim>
1156 *
std::vector<types::global_dof_index> local_dof_indices;
1161 *
local_dof_indices(dofs_per_cell)
1169 *
template <
int dim>
1196 *
std::vector<std::vector<double> >
Nx;
1197 *
std::vector<std::vector<Tensor<2, dim> > >
grad_Nx;
1198 *
std::vector<std::vector<SymmetricTensor<2, dim> > >
symm_grad_Nx;
1212 *
std::vector<double>(
fe_cell.dofs_per_cell)),
1223 *
rhs.fe_values_ref.get_quadrature(),
1224 *
rhs.fe_values_ref.get_update_flags()),
1226 *
rhs.fe_face_values_ref.get_quadrature(),
1227 *
rhs.fe_face_values_ref.get_update_flags()),
1238 *
const unsigned int n_dofs_per_cell =
Nx[0].size();
1264 *
for (
unsigned int k = 0;
k < n_dofs_per_cell; ++
k)
1276 *
const int dim = 2;
1277 *
const double tol = 1
e-12;
1279 *
parameters.length/2.0,
1280 *
parameters.width/2.0,
1281 *
parameters.hole_diameter/2.0,
1282 *
parameters.n_repetitions_xy,
1283 *
parameters.hole_division_fraction);
1290 * for (typename Triangulation<dim>::active_cell_iterator
1291 * cell = triangulation.begin_active();
1292 * cell != triangulation.end(); ++cell)
1294 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1295 * if (cell->face(face)->at_boundary())
1297 * cell->face(face)->set_all_boundary_ids(0);
1303 * Set boundary IDs and and manifolds
1306 * const Point<dim> centre (0,0);
1307 * for (typename Triangulation<dim>::active_cell_iterator
1308 * cell = triangulation.begin_active();
1309 * cell != triangulation.end(); ++cell)
1311 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1312 * if (cell->face(face)->at_boundary())
1319 * if (std::abs(cell->face(face)->center()[0] - 0.0) < tol)
1321 * cell->face(face)->set_boundary_id(parameters.boundary_id_minus_X);
1323 * else if (std::abs(cell->face(face)->center()[0] - parameters.width/2.0) < tol)
1325 * cell->face(face)->set_boundary_id(parameters.boundary_id_plus_X);
1327 * else if (std::abs(cell->face(face)->center()[1] - 0.0) < tol)
1329 * cell->face(face)->set_boundary_id(parameters.boundary_id_minus_Y);
1331 * else if (std::abs(cell->face(face)->center()[1] - parameters.length/2.0) < tol)
1333 * cell->face(face)->set_boundary_id(parameters.boundary_id_plus_Y);
1337 * for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
1338 * if (std::abs(cell->vertex(vertex).distance(centre) - parameters.hole_diameter/2.0) < tol)
1340 * cell->face(face)->set_boundary_id(parameters.boundary_id_hole);
1350 * for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
1351 * if (std::abs(cell->vertex(vertex).distance(centre) - parameters.hole_diameter/2.0) < tol)
1353 * cell->face(face)->set_manifold_id(parameters.manifold_id_hole);
1358 * static SphericalManifold<dim> spherical_manifold (centre);
1359 * triangulation.set_manifold(parameters.manifold_id_hole,spherical_manifold);
1360 * triangulation.refine_global(parameters.global_refinement);
1361 * GridTools::scale(parameters.scale,triangulation);
1364 * void Solid<3>::make_grid()
1366 * const int dim = 3;
1367 * const double tol = 1e-12;
1368 * Triangulation<2> tria_2d;
1369 * make_2d_quarter_plate_with_hole(tria_2d,
1370 * parameters.length/2.0,
1371 * parameters.width/2.0,
1372 * parameters.hole_diameter/2.0,
1373 * parameters.n_repetitions_xy,
1374 * parameters.hole_division_fraction);
1375 * GridGenerator::extrude_triangulation(tria_2d,
1376 * parameters.n_repetitions_z+1,
1377 * parameters.thickness/2.0,
1382 * Clear boundary ID's
1390 *
if (cell->face(face)->at_boundary())
1392 *
cell->face(face)->set_all_boundary_ids(0);
1408 *
if (cell->face(face)->at_boundary())
1415 *
if (
std::abs(cell->face(face)->center()[0] - 0.0) < tol)
1417 *
cell->face(face)->set_boundary_id(parameters.boundary_id_minus_X);
1419 *
else if (
std::abs(cell->face(face)->center()[0] - parameters.width/2.0) < tol)
1421 *
cell->face(face)->set_boundary_id(parameters.boundary_id_plus_X);
1423 *
else if (
std::abs(cell->face(face)->center()[1] - 0.0) < tol)
1425 *
cell->face(face)->set_boundary_id(parameters.boundary_id_minus_Y);
1427 *
else if (
std::abs(cell->face(face)->center()[1] - parameters.length/2.0) < tol)
1429 *
cell->face(face)->set_boundary_id(parameters.boundary_id_plus_Y);
1431 *
else if (
std::abs(cell->face(face)->center()[2] - 0.0) < tol)
1433 *
cell->face(face)->set_boundary_id(parameters.boundary_id_minus_Z);
1435 *
else if (
std::abs(cell->face(face)->center()[2] - parameters.thickness/2.0) < tol)
1437 *
cell->face(face)->set_boundary_id(parameters.boundary_id_plus_Z);
1453 *
cell->face(face)->set_boundary_id(parameters.boundary_id_hole);
1481 *
cell->face(face)->set_all_manifold_ids(parameters.manifold_id_hole);
1489 *
triangulation.refine_global(parameters.global_refinement);
1492 *
template <
int dim>
1512 *
std::set<typename Triangulation<2>::active_cell_iterator >
cells_to_remove;
1522 *
if (cell->center()[0] < 0.0 || cell->center()[1] < 0.0)
1539 * with cells with a decent aspect ratio
1542 * std::vector<std::vector<double> > step_sizes;
1544 * std::vector<double> subdivision_width;
1545 * subdivision_width.push_back(internal_width/2.0);
1546 * const double width_remaining = (width - internal_width)/2.0;
1547 * const unsigned int n_subs = static_cast<unsigned int>(std::max(1.0,std::ceil(width_remaining/(internal_width/2.0))));
1548 * Assert(n_subs>0, ExcInternalError());
1549 * for (unsigned int s=0; s<n_subs; ++s)
1550 * subdivision_width.push_back(width_remaining/n_subs);
1551 * step_sizes.push_back(subdivision_width);
1553 * const double sum_half_width = std::accumulate(subdivision_width.begin(), subdivision_width.end(), 0.0);
1554 * (void)sum_half_width;
1555 * Assert(std::abs(sum_half_width-width/2.0) < 1e-12, ExcInternalError());
1558 * std::vector<double> subdivision_length;
1559 * subdivision_length.push_back(internal_width/2.0);
1560 * const double length_remaining = (length - internal_width)/2.0;
1561 * const unsigned int n_subs = static_cast<unsigned int>(std::max(1.0,std::ceil(length_remaining/(internal_width/2.0))));
1562 * Assert(n_subs>0, ExcInternalError());
1563 * for (unsigned int s=0; s<n_subs; ++s)
1564 * subdivision_length.push_back(length_remaining/n_subs);
1565 * step_sizes.push_back(subdivision_length);
1567 * const double sum_half_length = std::accumulate(subdivision_length.begin(), subdivision_length.end(), 0.0);
1568 * (void)sum_half_length;
1569 * Assert(std::abs(sum_half_length-length/2.0) < 1e-12, ExcInternalError());
1572 * GridGenerator::subdivided_hyper_rectangle(tria_plate,
1574 * Point<2>(0.0, 0.0),
1575 * Point<2>(width/2.0, length/2.0));
1577 * std::set<typename Triangulation<2>::active_cell_iterator > cells_to_remove;
1578 * for (typename Triangulation<2>::active_cell_iterator
1579 * cell = tria_plate.begin_active();
1580 * cell != tria_plate.end(); ++cell)
1584 * Remove all cells that are in the first quadrant
1587 * if (cell->center()[0] < internal_width/2.0 && cell->center()[1] < internal_width/2.0)
1588 * cells_to_remove.insert(cell);
1590 * Assert(cells_to_remove.size() > 0, ExcInternalError());
1591 * Assert(cells_to_remove.size() != tria_plate.n_active_cells(), ExcInternalError());
1592 * GridGenerator::create_triangulation_with_removed_cells(tria_plate,cells_to_remove,tria_cut_plate);
1595 * Triangulation<2> tria_2d_not_flat;
1596 * GridGenerator::merge_triangulations(tria_quarter_plate_hole,
1598 * tria_2d_not_flat);
1602 * Attach a manifold to the curved boundary and refine
1603 * Note: We can only guarantee that the vertices sit on
1604 * the curve, so we must test with their position instead
1605 * of the cell centre.
1608 * const Point<2> centre_2d (0,0);
1609 * for (typename Triangulation<2>::active_cell_iterator
1610 * cell = tria_2d_not_flat.begin_active();
1611 * cell != tria_2d_not_flat.end(); ++cell)
1613 * for (unsigned int face=0; face<GeometryInfo<2>::faces_per_cell; ++face)
1614 * if (cell->face(face)->at_boundary())
1615 * for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_face; ++vertex)
1616 * if (std::abs(cell->vertex(vertex).distance(centre_2d) - hole_diameter/2.0) < 1e-12)
1618 * cell->face(face)->set_manifold_id(10);
1622 * SphericalManifold<2> spherical_manifold_2d (centre_2d);
1623 * tria_2d_not_flat.set_manifold(10,spherical_manifold_2d);
1624 * tria_2d_not_flat.refine_global(std::max (1U, n_repetitions_xy));
1625 * tria_2d_not_flat.reset_manifold(10); // Clear manifold
1627 * GridGenerator::flatten_triangulation(tria_2d_not_flat,tria_2d);
1629 * template <int dim>
1630 * void Solid<dim>::setup_system(LA::MPI::BlockVector &solution_delta)
1632 * timer.enter_subsection("Setup system");
1633 * pcout << "Setting up linear system..." << std::endl;
1637 * Partition triangulation
1640 * GridTools::partition_triangulation (n_mpi_processes,
1643 * block_component = std::vector<unsigned int> (n_components, u_block); // Displacement
1644 * block_component[p_component] = p_block; // Pressure
1645 * block_component[J_component] = J_block; // Dilatation
1646 * dof_handler.distribute_dofs(fe);
1647 * DoFRenumbering::Cuthill_McKee(dof_handler);
1648 * DoFRenumbering::component_wise(dof_handler, block_component);
1652 * Count DoFs in each block
1655 * dofs_per_block = DoFTools::count_dofs_per_fe_block(dof_handler, block_component);
1657 * all_locally_owned_dofs = DoFTools::locally_owned_dofs_per_subdomain (dof_handler);
1658 * std::vector<IndexSet> all_locally_relevant_dofs
1659 * = DoFTools::locally_relevant_dofs_per_subdomain (dof_handler);
1661 * locally_owned_dofs.clear();
1662 * locally_owned_partitioning.clear();
1663 * Assert(all_locally_owned_dofs.size() > this_mpi_process, ExcInternalError());
1664 * locally_owned_dofs = all_locally_owned_dofs[this_mpi_process];
1666 * locally_relevant_dofs.clear();
1667 * locally_relevant_partitioning.clear();
1668 * Assert(all_locally_relevant_dofs.size() > this_mpi_process, ExcInternalError());
1669 * locally_relevant_dofs = all_locally_relevant_dofs[this_mpi_process];
1671 * locally_owned_partitioning.reserve(n_blocks);
1672 * locally_relevant_partitioning.reserve(n_blocks);
1673 * for (unsigned int b=0; b<n_blocks; ++b)
1675 * const types::global_dof_index idx_begin
1676 * = std::accumulate(dofs_per_block.begin(),
1677 * std::next(dofs_per_block.begin(),b), 0);
1678 * const types::global_dof_index idx_end
1679 * = std::accumulate(dofs_per_block.begin(),
1680 * std::next(dofs_per_block.begin(),b+1), 0);
1681 * locally_owned_partitioning.push_back(locally_owned_dofs.get_view(idx_begin, idx_end));
1682 * locally_relevant_partitioning.push_back(locally_relevant_dofs.get_view(idx_begin, idx_end));
1686 * << " Number of active cells: " << triangulation.n_active_cells()
1687 * << " (by partition:";
1688 * for (unsigned int p=0; p<n_mpi_processes; ++p)
1690 * << (p==0 ? ' ' : '+
')
1691 * << (GridTools::count_cells_with_subdomain_association (triangulation,p));
1692 * pcout << ")" << std::endl;
1695 * << " Number of degrees of freedom: " << dof_handler.n_dofs()
1696 * << " (by partition:";
1697 * for (unsigned int p=0; p<n_mpi_processes; ++p)
1699 * << (p==0 ? ' ' : '+
')
1700 * << (DoFTools::count_dofs_with_subdomain_association (dof_handler,p));
1701 * pcout << ")" << std::endl;
1703 * << " Number of degrees of freedom per block: "
1704 * << "[n_u, n_p, n_J] = ["
1705 * << dofs_per_block[u_block] << ", "
1706 * << dofs_per_block[p_block] << ", "
1707 * << dofs_per_block[J_block] << "]"
1711 * Table<2, DoFTools::Coupling> coupling(n_components, n_components);
1712 * for (unsigned int ii = 0; ii < n_components; ++ii)
1713 * for (unsigned int jj = 0; jj < n_components; ++jj)
1714 * if (((ii < p_component) && (jj == J_component))
1715 * || ((ii == J_component) && (jj < p_component))
1716 * || ((ii == p_component) && (jj == p_component)))
1717 * coupling[ii][jj] = DoFTools::none;
1719 * coupling[ii][jj] = DoFTools::always;
1721 * TrilinosWrappers::BlockSparsityPattern bsp (locally_owned_partitioning,
1722 * locally_owned_partitioning,
1723 * locally_relevant_partitioning,
1724 * mpi_communicator);
1725 * DoFTools::make_sparsity_pattern (dof_handler, bsp,
1726 * constraints, false,
1727 * this_mpi_process);
1729 * tangent_matrix.reinit (bsp);
1733 * We then set up storage vectors
1736 * system_rhs.reinit(locally_owned_partitioning,
1737 * mpi_communicator);
1738 * solution_n.reinit(locally_owned_partitioning,
1739 * mpi_communicator);
1740 * solution_delta.reinit(locally_owned_partitioning,
1741 * mpi_communicator);
1743 * timer.leave_subsection();
1745 * template <int dim>
1747 * Solid<dim>::determine_component_extractors()
1749 * element_indices_u.clear();
1750 * element_indices_p.clear();
1751 * element_indices_J.clear();
1752 * for (unsigned int k = 0; k < fe.dofs_per_cell; ++k)
1754 * const unsigned int k_group = fe.system_to_base_index(k).first.first;
1755 * if (k_group == u_block)
1756 * element_indices_u.push_back(k);
1757 * else if (k_group == p_block)
1758 * element_indices_p.push_back(k);
1759 * else if (k_group == J_block)
1760 * element_indices_J.push_back(k);
1763 * Assert(k_group <= J_block, ExcInternalError());
1767 * template <int dim>
1768 * void Solid<dim>::setup_qph()
1770 * pcout << "Setting up quadrature point data..." << std::endl;
1771 * quadrature_point_history.initialize(triangulation.begin_active(),
1772 * triangulation.end(),
1774 * FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
1775 * cell (IteratorFilters::SubdomainEqualTo(this_mpi_process),
1776 * dof_handler.begin_active()),
1777 * endc (IteratorFilters::SubdomainEqualTo(this_mpi_process),
1778 * dof_handler.end());
1779 * for (; cell!=endc; ++cell)
1781 * Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());
1782 * const std::vector<std::shared_ptr<PointHistory<dim> > > lqph =
1783 * quadrature_point_history.get_data(cell);
1784 * Assert(lqph.size() == n_q_points, ExcInternalError());
1785 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1786 * lqph[q_point]->setup_lqp(parameters, time);
1789 * template <int dim>
1791 * Solid<dim>::solve_nonlinear_timestep(LA::MPI::BlockVector &solution_delta)
1793 * pcout << std::endl
1794 * << "Timestep " << time.get_timestep() << " @ "
1795 * << time.current() << "s of "
1796 * << time.end() << "s" << std::endl;
1797 * LA::MPI::BlockVector newton_update(locally_owned_partitioning,
1798 * mpi_communicator);
1799 * error_residual.reset();
1800 * error_residual_0.reset();
1801 * error_residual_norm.reset();
1802 * error_update.reset();
1803 * error_update_0.reset();
1804 * error_update_norm.reset();
1805 * print_conv_header();
1806 * unsigned int newton_iteration = 0;
1807 * for (; newton_iteration < parameters.max_iterations_NR;
1808 * ++newton_iteration)
1810 * pcout << " " << std::setw(2) << newton_iteration << " " << std::flush;
1811 * make_constraints(newton_iteration);
1812 * assemble_system(solution_delta);
1813 * get_error_residual(error_residual);
1814 * if (newton_iteration == 0)
1815 * error_residual_0 = error_residual;
1816 * error_residual_norm = error_residual;
1817 * error_residual_norm.normalise(error_residual_0);
1818 * if (newton_iteration > 0 &&
1819 * (error_update_norm.u <= parameters.tol_u &&
1820 * error_residual_norm.u <= parameters.tol_f) )
1822 * pcout << " CONVERGED! " << std::endl;
1823 * print_conv_footer(solution_delta);
1826 * const std::pair<unsigned int, double>
1827 * lin_solver_output = solve_linear_system(newton_update);
1828 * get_error_update(newton_update, error_update);
1829 * if (newton_iteration == 0)
1830 * error_update_0 = error_update;
1831 * error_update_norm = error_update;
1832 * error_update_norm.normalise(error_update_0);
1833 * solution_delta += newton_update;
1834 * newton_update = 0.0;
1835 * pcout << " | " << std::fixed << std::setprecision(3) << std::setw(7)
1836 * << std::scientific << lin_solver_output.first << " "
1837 * << lin_solver_output.second << " " << error_residual_norm.norm
1838 * << " " << error_residual_norm.u << " "
1839 * << error_residual_norm.p << " " << error_residual_norm.J
1840 * << " " << error_update_norm.norm << " " << error_update_norm.u
1841 * << " " << error_update_norm.p << " " << error_update_norm.J
1842 * << " " << std::endl;
1844 * AssertThrow (newton_iteration <= parameters.max_iterations_NR,
1845 * ExcMessage("No convergence in nonlinear solver!"));
1847 * template <int dim>
1848 * void Solid<dim>::print_conv_header()
1850 * pcout << std::string(132,'_') << std::endl;
1851 * pcout << " SOLVER STEP "
1852 * << " | LIN_IT LIN_RES RES_NORM "
1853 * << " RES_U RES_P RES_J NU_NORM "
1854 * << " NU_U NU_P NU_J " << std::endl;
1855 * pcout << std::string(132,'_') << std::endl;
1857 * template <int dim>
1858 * void Solid<dim>::print_conv_footer(const LA::MPI::BlockVector &solution_delta)
1860 * pcout << std::string(132,'_') << std::endl;
1861 * const std::pair<double,std::pair<double,double> > error_dil = get_error_dilation(get_solution_total(solution_delta));
1862 * pcout << "Relative errors:" << std::endl
1863 * << "Displacement:\t" << error_update.u / error_update_0.u << std::endl
1864 * << "Force: \t\t" << error_residual.u / error_residual_0.u << std::endl
1865 * << "Dilatation:\t" << error_dil.first << std::endl
1866 * << "v / V_0:\t" << error_dil.second.second << " / " << error_dil.second.first
1867 * << " = " << (error_dil.second.second/error_dil.second.first) << std::endl;
1869 * template <int dim>
1870 * std::pair<double,std::pair<double,double> >
1871 * Solid<dim>::get_error_dilation(const LA::MPI::BlockVector &solution_total) const
1873 * double vol_reference = 0.0;
1874 * double vol_current = 0.0;
1875 * double dil_L2_error = 0.0;
1876 * FEValues<dim> fe_values_ref(fe, qf_cell,
1877 * update_values | update_gradients | update_JxW_values);
1878 * std::vector<Tensor<2, dim> > solution_grads_u_total (qf_cell.size());
1879 * std::vector<double> solution_values_J_total (qf_cell.size());
1880 * FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
1881 * cell (IteratorFilters::SubdomainEqualTo(this_mpi_process),
1882 * dof_handler.begin_active()),
1883 * endc (IteratorFilters::SubdomainEqualTo(this_mpi_process),
1884 * dof_handler.end());
1885 * for (; cell != endc; ++cell)
1887 * Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());
1888 * fe_values_ref.reinit(cell);
1889 * fe_values_ref[u_fe].get_function_gradients(solution_total,
1890 * solution_grads_u_total);
1891 * fe_values_ref[J_fe].get_function_values(solution_total,
1892 * solution_values_J_total);
1893 * const std::vector<std::shared_ptr<const PointHistory<dim> > > lqph =
1894 * quadrature_point_history.get_data(cell);
1895 * Assert(lqph.size() == n_q_points, ExcInternalError());
1896 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1898 * const double det_F_qp = determinant(Physics::Elasticity::Kinematics::F(solution_grads_u_total[q_point]));
1899 * const double J_tilde_qp = solution_values_J_total[q_point];
1900 * const double the_error_qp_squared = std::pow((det_F_qp - J_tilde_qp),
1902 * const double JxW = fe_values_ref.JxW(q_point);
1903 * dil_L2_error += the_error_qp_squared * JxW;
1904 * vol_reference += JxW;
1905 * vol_current += det_F_qp * JxW;
1908 * Assert(vol_current > 0.0, ExcInternalError());
1911 * Sum across all processors
1914 * dil_L2_error = Utilities::MPI::sum(dil_L2_error,mpi_communicator);
1915 * vol_reference = Utilities::MPI::sum(vol_reference,mpi_communicator);
1916 * vol_current = Utilities::MPI::sum(vol_current,mpi_communicator);
1918 * return std::make_pair(std::sqrt(dil_L2_error),
1919 * std::make_pair(vol_reference,vol_current));
1921 * template <int dim>
1922 * void Solid<dim>::get_error_residual(Errors &error_residual)
1926 * Construct a residual vector that has the values for all of its
1927 * constrained DoFs set to zero.
1930 * LA::MPI::BlockVector error_res (system_rhs);
1931 * constraints.set_zero(error_res);
1932 * error_residual.norm = error_res.l2_norm();
1933 * error_residual.u = error_res.block(u_block).l2_norm();
1934 * error_residual.p = error_res.block(p_block).l2_norm();
1935 * error_residual.J = error_res.block(J_block).l2_norm();
1937 * template <int dim>
1938 * void Solid<dim>::get_error_update(const LA::MPI::BlockVector &newton_update,
1939 * Errors &error_update)
1943 * Construct a update vector that has the values for all of its
1944 * constrained DoFs set to zero.
1947 * LA::MPI::BlockVector error_ud (newton_update);
1948 * constraints.set_zero(error_ud);
1949 * error_update.norm = error_ud.l2_norm();
1950 * error_update.u = error_ud.block(u_block).l2_norm();
1951 * error_update.p = error_ud.block(p_block).l2_norm();
1952 * error_update.J = error_ud.block(J_block).l2_norm();
1954 * template <int dim>
1955 * LA::MPI::BlockVector
1956 * Solid<dim>::get_solution_total(const LA::MPI::BlockVector &solution_delta) const
1960 * Cell interpolation -> Ghosted vector
1963 * LA::MPI::BlockVector solution_total (locally_owned_partitioning,
1964 * locally_relevant_partitioning,
1966 * /*vector_writable = */ false);
1967 * LA::MPI::BlockVector tmp (solution_total);
1968 * solution_total = solution_n;
1969 * tmp = solution_delta;
1970 * solution_total += tmp;
1971 * return solution_total;
1973 * template <int dim>
1974 * void Solid<dim>::assemble_system(const LA::MPI::BlockVector &solution_delta)
1976 * timer.enter_subsection("Assemble system");
1977 * pcout << " ASM_SYS " << std::flush;
1978 * tangent_matrix = 0.0;
1980 * const LA::MPI::BlockVector solution_total(get_solution_total(solution_delta));
1981 * const UpdateFlags uf_cell(update_values |
1982 * update_gradients |
1983 * update_JxW_values);
1984 * const UpdateFlags uf_face(update_values |
1985 * update_normal_vectors |
1986 * update_JxW_values);
1987 * PerTaskData_ASM per_task_data(dofs_per_cell);
1988 * ScratchData_ASM scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face, solution_total);
1990 * FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
1991 * cell (IteratorFilters::SubdomainEqualTo(this_mpi_process),
1992 * dof_handler.begin_active()),
1993 * endc (IteratorFilters::SubdomainEqualTo(this_mpi_process),
1994 * dof_handler.end());
1995 * for (; cell != endc; ++cell)
1997 * Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());
1998 * assemble_system_one_cell(cell, scratch_data, per_task_data);
1999 * copy_local_to_global_system(per_task_data);
2001 * tangent_matrix.compress(VectorOperation::add);
2002 * system_rhs.compress(VectorOperation::add);
2003 * timer.leave_subsection();
2005 * template <int dim>
2006 * void Solid<dim>::copy_local_to_global_system(const PerTaskData_ASM &data)
2008 * constraints.distribute_local_to_global(data.cell_matrix, data.cell_rhs,
2009 * data.local_dof_indices,
2010 * tangent_matrix, system_rhs);
2012 * template <int dim>
2014 * Solid<dim>::assemble_system_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2015 * ScratchData_ASM &scratch,
2016 * PerTaskData_ASM &data) const
2018 * Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());
2022 * scratch.fe_values_ref.reinit(cell);
2023 * cell->get_dof_indices(data.local_dof_indices);
2024 * const std::vector<std::shared_ptr<const PointHistory<dim> > > lqph =
2025 * quadrature_point_history.get_data(cell);
2026 * Assert(lqph.size() == n_q_points, ExcInternalError());
2030 * Update quadrature point solution
2033 * scratch.fe_values_ref[u_fe].get_function_gradients(scratch.solution_total,
2034 * scratch.solution_grads_u_total);
2035 * scratch.fe_values_ref[p_fe].get_function_values(scratch.solution_total,
2036 * scratch.solution_values_p_total);
2037 * scratch.fe_values_ref[J_fe].get_function_values(scratch.solution_total,
2038 * scratch.solution_values_J_total);
2042 * Update shape functions and their gradients (push-forward)
2045 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2047 * const Tensor<2, dim> F = Physics::Elasticity::Kinematics::F(scratch.solution_grads_u_total[q_point]);
2048 * const Tensor<2, dim> F_inv = invert(F);
2050 * for (unsigned int k = 0; k < dofs_per_cell; ++k)
2052 * const unsigned int k_group = fe.system_to_base_index(k).first.first;
2053 * if (k_group == u_block)
2055 * scratch.grad_Nx[q_point][k] = scratch.fe_values_ref[u_fe].gradient(k, q_point)
2057 * scratch.symm_grad_Nx[q_point][k] = symmetrize(scratch.grad_Nx[q_point][k]);
2059 * else if (k_group == p_block)
2060 * scratch.Nx[q_point][k] = scratch.fe_values_ref[p_fe].value(k,
2062 * else if (k_group == J_block)
2063 * scratch.Nx[q_point][k] = scratch.fe_values_ref[J_fe].value(k,
2066 * Assert(k_group <= J_block, ExcInternalError());
2069 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2071 * const SymmetricTensor<2, dim> &I = Physics::Elasticity::StandardTensors<dim>::I;
2072 * const Tensor<2, dim> F = Physics::Elasticity::Kinematics::F(scratch.solution_grads_u_total[q_point]);
2073 * const double det_F = determinant(F);
2074 * const double &p_tilde = scratch.solution_values_p_total[q_point];
2075 * const double &J_tilde = scratch.solution_values_J_total[q_point];
2076 * Assert(det_F > 0, ExcInternalError());
2079 * PointHistory<dim> *lqph_q_point_nc = const_cast<PointHistory<dim>*>(lqph[q_point].get());
2080 * lqph_q_point_nc->update_internal_equilibrium(F,p_tilde,J_tilde);
2083 * const SymmetricTensor<2, dim> tau = lqph[q_point]->get_tau(F,p_tilde);
2084 * const Tensor<2, dim> tau_ns (tau);
2085 * const SymmetricTensor<4, dim> Jc = lqph[q_point]->get_Jc(F,p_tilde);
2086 * const double dPsi_vol_dJ = lqph[q_point]->get_dPsi_vol_dJ(J_tilde);
2087 * const double d2Psi_vol_dJ2 = lqph[q_point]->get_d2Psi_vol_dJ2(J_tilde);
2089 * const std::vector<double> &Nx = scratch.Nx[q_point];
2090 * const std::vector<Tensor<2, dim> > &grad_Nx = scratch.grad_Nx[q_point];
2091 * const std::vector<SymmetricTensor<2, dim> > &symm_grad_Nx = scratch.symm_grad_Nx[q_point];
2092 * const double JxW = scratch.fe_values_ref.JxW(q_point);
2094 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2096 * const unsigned int component_i = fe.system_to_component_index(i).first;
2097 * const unsigned int i_group = fe.system_to_base_index(i).first.first;
2098 * if (i_group == u_block)
2099 * data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW;
2100 * else if (i_group == p_block)
2101 * data.cell_rhs(i) -= Nx[i] * (det_F - J_tilde) * JxW;
2102 * else if (i_group == J_block)
2103 * data.cell_rhs(i) -= Nx[i] * (dPsi_vol_dJ - p_tilde) * JxW;
2105 * Assert(i_group <= J_block, ExcInternalError());
2107 * for (unsigned int j = 0; j <= i; ++j)
2109 * const unsigned int component_j = fe.system_to_component_index(j).first;
2110 * const unsigned int j_group = fe.system_to_base_index(j).first.first;
2111 * if ((i_group == u_block) && (j_group == u_block))
2113 * data.cell_matrix(i, j) += symm_grad_Nx[i] * Jc // The material contribution:
2114 * * symm_grad_Nx[j] * JxW;
2115 * if (component_i == component_j) // geometrical stress contribution
2116 * data.cell_matrix(i, j) += grad_Nx[i][component_i] * tau_ns
2117 * * grad_Nx[j][component_j] * JxW;
2119 * else if ((i_group == u_block) && (j_group == p_block))
2121 * data.cell_matrix(i, j) += (symm_grad_Nx[i] * I)
2125 * else if ((i_group == p_block) && (j_group == u_block))
2127 * data.cell_matrix(i, j) += Nx[i] * det_F
2128 * * (symm_grad_Nx[j] * I)
2131 * else if ((i_group == p_block) && (j_group == J_block))
2132 * data.cell_matrix(i, j) -= Nx[i] * Nx[j] * JxW;
2133 * else if ((i_group == J_block) && (j_group == p_block))
2134 * data.cell_matrix(i, j) -= Nx[i] * Nx[j] * JxW;
2135 * else if ((i_group == J_block) && (j_group == J_block))
2136 * data.cell_matrix(i, j) += Nx[i] * d2Psi_vol_dJ2 * Nx[j] * JxW;
2138 * Assert((i_group <= J_block) && (j_group <= J_block),
2139 * ExcInternalError());
2144 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2145 * for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
2146 * data.cell_matrix(i, j) = data.cell_matrix(j, i);
2148 * if (parameters.driver == "Neumann")
2149 * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
2151 * if (cell->face(face)->at_boundary() == true
2152 * && cell->face(face)->boundary_id() == parameters.boundary_id_plus_Y)
2154 * scratch.fe_face_values_ref.reinit(cell, face);
2155 * for (unsigned int f_q_point = 0; f_q_point < n_q_points_f;
2158 * const Tensor<1, dim> &N =
2159 * scratch.fe_face_values_ref.normal_vector(f_q_point);
2160 * static const double pressure_nom = parameters.pressure
2161 * / (parameters.scale * parameters.scale);
2162 * const double time_ramp = (time.current() < parameters.load_time ?
2163 * time.current() / parameters.load_time : 1.0);
2164 * const double pressure = -pressure_nom * time_ramp;
2165 * const Tensor<1, dim> traction = pressure * N;
2166 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2168 * const unsigned int i_group =
2169 * fe.system_to_base_index(i).first.first;
2170 * if (i_group == u_block)
2172 * const unsigned int component_i =
2173 * fe.system_to_component_index(i).first;
2175 * scratch.fe_face_values_ref.shape_value(i,
2177 * const double JxW = scratch.fe_face_values_ref.JxW(
2179 * data.cell_rhs(i) += (Ni * traction[component_i])
2186 * template <int dim>
2187 * void Solid<dim>::make_constraints(const int &it_nr)
2189 * pcout << " CST " << std::flush;
2192 * constraints.clear();
2193 * const bool apply_dirichlet_bc = (it_nr == 0);
2194 * const FEValuesExtractors::Scalar x_displacement(0);
2195 * const FEValuesExtractors::Scalar y_displacement(1);
2197 * const int boundary_id = parameters.boundary_id_minus_X;
2198 * if (apply_dirichlet_bc == true)
2199 * VectorTools::interpolate_boundary_values(dof_handler,
2201 * Functions::ZeroFunction<dim>(n_components),
2203 * fe.component_mask(x_displacement));
2205 * VectorTools::interpolate_boundary_values(dof_handler,
2207 * Functions::ZeroFunction<dim>(n_components),
2209 * fe.component_mask(x_displacement));
2212 * const int boundary_id = parameters.boundary_id_minus_Y;
2213 * if (apply_dirichlet_bc == true)
2214 * VectorTools::interpolate_boundary_values(dof_handler,
2216 * Functions::ZeroFunction<dim>(n_components),
2218 * fe.component_mask(y_displacement));
2220 * VectorTools::interpolate_boundary_values(dof_handler,
2222 * Functions::ZeroFunction<dim>(n_components),
2224 * fe.component_mask(y_displacement));
2228 * const FEValuesExtractors::Scalar z_displacement(2);
2230 * const int boundary_id = parameters.boundary_id_minus_Z;
2231 * if (apply_dirichlet_bc == true)
2232 * VectorTools::interpolate_boundary_values(dof_handler,
2234 * Functions::ZeroFunction<dim>(n_components),
2236 * fe.component_mask(z_displacement));
2238 * VectorTools::interpolate_boundary_values(dof_handler,
2240 * Functions::ZeroFunction<dim>(n_components),
2242 * fe.component_mask(z_displacement));
2245 * const int boundary_id = parameters.boundary_id_plus_Z;
2246 * if (apply_dirichlet_bc == true)
2247 * VectorTools::interpolate_boundary_values(dof_handler,
2249 * Functions::ZeroFunction<dim>(n_components),
2251 * fe.component_mask(z_displacement));
2253 * VectorTools::interpolate_boundary_values(dof_handler,
2255 * Functions::ZeroFunction<dim>(n_components),
2257 * fe.component_mask(z_displacement));
2260 * if (parameters.driver == "Dirichlet")
2262 * const int boundary_id = parameters.boundary_id_plus_Y;
2263 * if (apply_dirichlet_bc == true)
2266 * if (time.current() < parameters.load_time+0.01*time.get_delta_t())
2268 * const double delta_length = parameters.length*(parameters.stretch - 1.0)*parameters.scale;
2269 * const unsigned int n_stretch_steps = static_cast<unsigned int>(parameters.load_time/time.get_delta_t());
2270 * const double delta_u_y = delta_length/2.0/n_stretch_steps;
2271 * VectorTools::interpolate_boundary_values(dof_handler,
2273 * Functions::ConstantFunction<dim>(delta_u_y,n_components),
2275 * fe.component_mask(y_displacement));
2278 * VectorTools::interpolate_boundary_values(dof_handler,
2280 * Functions::ZeroFunction<dim>(n_components),
2282 * fe.component_mask(y_displacement));
2285 * VectorTools::interpolate_boundary_values(dof_handler,
2287 * Functions::ZeroFunction<dim>(n_components),
2289 * fe.component_mask(y_displacement));
2291 * constraints.close();
2293 * template <int dim>
2294 * std::pair<unsigned int, double>
2295 * Solid<dim>::solve_linear_system(LA::MPI::BlockVector &newton_update)
2297 * unsigned int lin_it = 0;
2298 * double lin_res = 0.0;
2300 * timer.enter_subsection("Linear solver");
2301 * pcout << " SLV " << std::flush;
2303 * const LA::MPI::Vector &f_u = system_rhs.block(u_block);
2304 * const LA::MPI::Vector &f_p = system_rhs.block(p_block);
2305 * const LA::MPI::Vector &f_J = system_rhs.block(J_block);
2306 * LA::MPI::Vector &d_u = newton_update.block(u_block);
2307 * LA::MPI::Vector &d_p = newton_update.block(p_block);
2308 * LA::MPI::Vector &d_J = newton_update.block(J_block);
2309 * const auto K_uu = linear_operator<LA::MPI::Vector>(tangent_matrix.block(u_block, u_block));
2310 * const auto K_up = linear_operator<LA::MPI::Vector>(tangent_matrix.block(u_block, p_block));
2311 * const auto K_pu = linear_operator<LA::MPI::Vector>(tangent_matrix.block(p_block, u_block));
2312 * const auto K_Jp = linear_operator<LA::MPI::Vector>(tangent_matrix.block(J_block, p_block));
2313 * const auto K_JJ = linear_operator<LA::MPI::Vector>(tangent_matrix.block(J_block, J_block));
2315 * LA::PreconditionJacobi preconditioner_K_Jp_inv;
2316 * preconditioner_K_Jp_inv.initialize(
2317 * tangent_matrix.block(J_block, p_block),
2318 * LA::PreconditionJacobi::AdditionalData());
2319 * ReductionControl solver_control_K_Jp_inv (
2320 * static_cast<unsigned int>(tangent_matrix.block(J_block, p_block).m()
2321 * * parameters.max_iterations_lin),
2323 * ::SolverCG<LA::MPI::Vector> solver_K_Jp_inv (solver_control_K_Jp_inv);
2325 * const auto K_Jp_inv = inverse_operator(K_Jp,
2327 * preconditioner_K_Jp_inv);
2328 * const auto K_pJ_inv = transpose_operator(K_Jp_inv);
2329 * const auto K_pp_bar = K_Jp_inv * K_JJ * K_pJ_inv;
2330 * const auto K_uu_bar_bar = K_up * K_pp_bar * K_pu;
2331 * const auto K_uu_con = K_uu + K_uu_bar_bar;
2333 * LA::PreconditionAMG preconditioner_K_con_inv;
2334 * preconditioner_K_con_inv.initialize(
2335 * tangent_matrix.block(u_block, u_block),
2336 * LA::PreconditionAMG::AdditionalData(
2337 * true /*elliptic*/,
2338 * (parameters.poly_degree > 1 /*higher_order_elements*/)) );
2339 * ReductionControl solver_control_K_con_inv (
2340 * static_cast<unsigned int>(tangent_matrix.block(u_block, u_block).m()
2341 * * parameters.max_iterations_lin),
2342 * 1.0e-30, parameters.tol_lin);
2343 * ::SolverSelector<LA::MPI::Vector> solver_K_con_inv;
2344 * solver_K_con_inv.select(parameters.type_lin);
2345 * solver_K_con_inv.set_control(solver_control_K_con_inv);
2346 * const auto K_uu_con_inv = inverse_operator(K_uu_con,
2348 * preconditioner_K_con_inv);
2350 * d_u = K_uu_con_inv*(f_u - K_up*(K_Jp_inv*f_J - K_pp_bar*f_p));
2351 * lin_it = solver_control_K_con_inv.last_step();
2352 * lin_res = solver_control_K_con_inv.last_value();
2353 * timer.leave_subsection();
2355 * timer.enter_subsection("Linear solver postprocessing");
2356 * d_J = K_pJ_inv*(f_p - K_pu*d_u);
2357 * d_p = K_Jp_inv*(f_J - K_JJ*d_J);
2358 * timer.leave_subsection();
2360 * constraints.distribute(newton_update);
2361 * return std::make_pair(lin_it, lin_res);
2363 * template <int dim>
2365 * Solid<dim>::update_end_timestep ()
2367 * FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
2368 * cell (IteratorFilters::SubdomainEqualTo(this_mpi_process),
2369 * dof_handler.begin_active()),
2370 * endc (IteratorFilters::SubdomainEqualTo(this_mpi_process),
2371 * dof_handler.end());
2372 * for (; cell != endc; ++cell)
2374 * Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());
2375 * const std::vector<std::shared_ptr<PointHistory<dim> > > lqph =
2376 * quadrature_point_history.get_data(cell);
2377 * Assert(lqph.size() == n_q_points, ExcInternalError());
2378 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2379 * lqph[q_point]->update_end_timestep();
2384 * class FilteredDataOut : public DataOut<dim>
2387 * FilteredDataOut (const unsigned int subdomain_id)
2389 * subdomain_id (subdomain_id)
2392 * virtual ~FilteredDataOut() {}
2394 * virtual typename DataOut<dim>::cell_iterator
2397 * auto cell = this->dofs->begin_active();
2398 * while ((cell != this->dofs->end()) &&
2399 * (cell->subdomain_id() != subdomain_id))
2404 * virtual typename DataOut<dim>::cell_iterator
2405 * next_cell (const typename DataOut<dim>::cell_iterator &old_cell)
2407 * if (old_cell != this->dofs->end())
2409 * const IteratorFilters::SubdomainEqualTo predicate(subdomain_id);
2411 * ++(FilteredIterator
2412 * <typename DataOut<dim>::cell_iterator>
2413 * (predicate,old_cell));
2420 * const unsigned int subdomain_id;
2423 * template <int dim>
2424 * void Solid<dim>::output_results(const unsigned int timestep,
2425 * const double current_time) const
2429 * Output -> Ghosted vector
2432 * LA::MPI::BlockVector solution_total (locally_owned_partitioning,
2433 * locally_relevant_partitioning,
2435 * /*vector_writable = */ false);
2436 * LA::MPI::BlockVector residual (locally_owned_partitioning,
2437 * locally_relevant_partitioning,
2439 * /*vector_writable = */ false);
2440 * solution_total = solution_n;
2441 * residual = system_rhs;
2446 * --- Additional data ---
2449 * Vector<double> material_id;
2450 * Vector<double> polynomial_order;
2451 * material_id.reinit(triangulation.n_active_cells());
2452 * polynomial_order.reinit(triangulation.n_active_cells());
2453 * std::vector<types::subdomain_id> partition_int (triangulation.n_active_cells());
2455 * FilteredDataOut<dim> data_out(this_mpi_process);
2456 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
2457 * data_component_interpretation(dim,
2458 * DataComponentInterpretation::component_is_part_of_vector);
2459 * data_component_interpretation.push_back(DataComponentInterpretation::component_is_scalar);
2460 * data_component_interpretation.push_back(DataComponentInterpretation::component_is_scalar);
2462 * GridTools::get_subdomain_association (triangulation, partition_int);
2470 *
unsigned int c = 0;
2472 *
cell = dof_handler.begin_active(),
2473 *
endc = dof_handler.end();
2474 *
for (; cell!=
endc; ++cell, ++c)
2478 *
material_id(c) =
static_cast<int>(cell->material_id());
2481 *
std::vector<std::string>
solution_name(n_components,
"solution_");
2482 *
std::vector<std::string>
residual_name(n_components,
"residual_");
2483 *
for (
unsigned int c=0; c<n_components; ++c)
2506 *
data_out.attach_dof_handler(dof_handler);
2510 *
data_component_interpretation);
2511 *
data_out.add_data_vector(residual,
2514 *
data_component_interpretation);
2517 *
data_out.add_data_vector (material_id,
"material_id");
2518 *
data_out.add_data_vector (
partitioning,
"partitioning");
2519 *
data_out.build_patches(degree);
2530 *
<< (std::to_string(dim) +
"d")
2545 *
<< (std::to_string(dim) +
"d")
2557 *
<< (std::to_string(dim) +
"d")
2570 *
data_out.write_vtu(output);
2579 *
if (this_mpi_process == 0)
2610 *
template <
int dim>
2622 *
dof_handler.begin_active()),
2624 *
dof_handler.end());
2625 *
for (; cell !=
endc; ++cell)
2635 *
if (cell->vertex(v).distance(
pt_ref) < 1e-6*parameters.scale)
2637 *
for (
unsigned int d=0;
d<dim; ++
d)
2651 *
for (
unsigned int d=0;
d<dim; ++
d)
2661 *
using namespace dealii;
2668 *
const unsigned int dim = 2;
2672 *
catch (std::exception &exc)
2676 *
std::cerr << std::endl << std::endl
2677 *
<<
"----------------------------------------------------"
2679 *
std::cerr <<
"Exception on processing: " << std::endl << exc.what()
2680 *
<< std::endl <<
"Aborting!" << std::endl
2681 *
<<
"----------------------------------------------------"
2690 *
std::cerr << std::endl << std::endl
2691 *
<<
"----------------------------------------------------"
2693 *
std::cerr <<
"Unknown exception!" << std::endl <<
"Aborting!"
2695 *
<<
"----------------------------------------------------"
numbers::NumberTraits< Number >::real_type norm() const
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::active_cell_iterator active_cell_iterator
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &, const Matrix &)
std::vector< index_type > data
void write_pvd_record(std::ostream &out, const std::vector< std::pair< double, std::string > > ×_and_names)
void hyper_cube_with_cylindrical_hole(Triangulation< dim, spacedim > &triangulation, const double inner_radius=.25, const double outer_radius=.5, const double L=.5, const unsigned int repetitions=1, const bool colorize=false)
void cylinder(Triangulation< dim > &tria, const double radius=1., const double half_length=1.)
void create_triangulation_with_removed_cells(const Triangulation< dim, spacedim > &input_triangulation, const std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > &cells_to_remove, Triangulation< dim, spacedim > &result)
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
VectorType::value_type * end(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
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)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)