1439 *
const double g = 9.81;
1440 *
const double rho = 7700;
1448 *
template <
int dim>
1453 *
const unsigned int n_points = points.size();
1457 *
for (
unsigned int p = 0; p < n_points; ++p)
1466 * <a name=
"step_18-ThecodeIncrementalBoundaryValuecodeclass"></a>
1504 *
template <
int dim>
1511 *
virtual void vector_value(
const Point<dim> &p,
1515 *
vector_value_list(
const std::vector<
Point<dim>> &points,
1525 *
template <
int dim>
1536 *
template <
int dim>
1549 *
template <
int dim>
1554 *
const unsigned int n_points = points.size();
1558 *
for (
unsigned int p = 0; p < n_points; ++p)
1567 * <a name=
"step_18-ImplementationofthecodeTopLevelcodeclass"></a>
1577 *
template <
int dim>
1587 * <a name=
"step_18-Thepublicinterface"></a>
1599 *
template <
int dim>
1612 *
,
pcout(std::cout, this_mpi_process == 0)
1617 *
template <
int dim>
1620 *
dof_handler.
clear();
1635 *
template <
int dim>
1648 * <a name=
"step_18-TopLevelcreate_coarse_grid"></a>
1670 *
template <
int dim>
1675 *
for (
const auto &cell :
triangulation.active_cell_iterators())
1677 *
if (face->at_boundary())
1682 *
face->set_boundary_id(0);
1684 *
face->set_boundary_id(1);
1688 *
face->set_boundary_id(2);
1690 *
face->set_boundary_id(3);
1715 * <a name=
"step_18-TopLevelsetup_system"></a>
1733 * template <int dim>
1734 * void TopLevel<dim>::setup_system()
1736 * dof_handler.distribute_dofs(fe);
1737 * locally_owned_dofs = dof_handler.locally_owned_dofs();
1738 * locally_relevant_dofs =
1739 * DoFTools::extract_locally_relevant_dofs(dof_handler);
1743 * The next step is to set up constraints due to hanging nodes. This has
1744 * been handled many times before:
1747 * hanging_node_constraints.clear();
1748 * DoFTools::make_hanging_node_constraints(dof_handler,
1749 * hanging_node_constraints);
1750 * hanging_node_constraints.close();
1754 * And then we have to set up the matrix. Here we deviate from @ref step_17 "step-17", in
1758 * at all
efficient:
if we don't give PETSc a clue as to which elements
1759 * are written to, it is (at least at the time of this writing) unbearably
1760 * slow when we set the elements in the matrix for the first time (i.e. in
1761 * the first timestep). Later on, when the elements have been allocated,
1762 * everything is much faster. In experiments we made, the first timestep
1763 * can be accelerated by almost two orders of magnitude if we instruct
1764 * PETSc which elements will be used and which are not.
1768 * To do so, we first generate the sparsity pattern of the matrix we are
1769 * going to work with, and make sure that the condensation of hanging node
1770 * constraints add the necessary additional entries in the sparsity
1774 * DynamicSparsityPattern sparsity_pattern(locally_relevant_dofs);
1775 * DoFTools::make_sparsity_pattern(dof_handler,
1777 * hanging_node_constraints,
1778 * /*keep constrained dofs*/ false);
1779 * SparsityTools::distribute_sparsity_pattern(sparsity_pattern,
1780 * locally_owned_dofs,
1782 * locally_relevant_dofs);
1785 * Note that we have used the <code>DynamicSparsityPattern</code> class
1786 * here that was already introduced in @ref step_11 "step-11", rather than the
1787 * <code>SparsityPattern</code> class that we have used in all other
1788 * cases. The reason for this is that for the latter class to work we have
1789 * to give an initial upper bound for the number of entries in each row, a
1790 * task that is traditionally done by
1791 * <code>DoFHandler::max_couplings_between_dofs()</code>. However, this
1792 * function suffers from a serious problem: it has to compute an upper
1793 * bound to the number of nonzero entries in each row, and this is a
1794 * rather complicated task, in particular in 3d. In effect, while it is
1795 * quite accurate in 2d, it often comes up with much too large a number in
1796 * 3d, and in that case the <code>SparsityPattern</code> allocates much
1797 * too much memory at first, often several 100 MBs. This is later
1798 * corrected when <code>DoFTools::make_sparsity_pattern</code> is called
1835 *
system_matrix.reinit(locally_owned_dofs,
1836 *
locally_owned_dofs,
1838 *
mpi_communicator);
1854 *
system_rhs.reinit(locally_owned_dofs, mpi_communicator);
1863 * <a name=
"step_18-TopLevelassemble_system"></a>
1883 *
template <
int dim>
1887 *
system_matrix = 0;
1894 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1900 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1912 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1913 *
if (cell->is_locally_owned())
1918 *
fe_values.reinit(cell);
1933 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1934 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
1965 *
the quadrature points
on this cell:
1968 *
body_force.vector_value_list(fe_values.get_quadrature_points(),
1976 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1979 *
fe.system_to_component_index(i).first;
1988 *
fe_values.shape_value(i,
q_point) -
1998 *
as in @
ref step_17
"step-17":
2001 *
cell->get_dof_indices(local_dof_indices);
2005 *
local_dof_indices,
2041 * want no movement at all, corresponding to the cylinder being clamped or
2042 * cemented in at this part of the boundary. At the top, however, we want
2043 * a prescribed vertical downward motion compressing the cylinder; in
2044 * addition, we only want to restrict the vertical movement, but not the
2045 * horizontal ones -- one can think of this situation as a well-greased
2046 * plate sitting on top of the cylinder pushing it downwards: the atoms of
2047 * the cylinder are forced to move downward, but they are free to slide
2048 * horizontally along the plate.
2052 * The way to describe this is as follows: for boundary indicator zero
2053 * (bottom face) we use a dim-dimensional zero function representing no
2054 * motion in any coordinate direction. For the boundary with indicator 1
2055 * (top surface), we use the <code>IncrementalBoundaryValues</code> class,
2056 * but we specify an additional argument to the
2057 * <code>VectorTools::interpolate_boundary_values</code> function denoting
2058 * which vector components it should apply to; this is a vector of bools
2059 * for each vector component and because we only want to restrict vertical
2060 * motion, it has only its last component set:
2063 * const FEValuesExtractors::Scalar z_component(dim - 1);
2064 * std::map<types::global_dof_index, double> boundary_values;
2065 * VectorTools::interpolate_boundary_values(dof_handler,
2067 * Functions::ZeroFunction<dim>(dim),
2069 * VectorTools::interpolate_boundary_values(
2072 * IncrementalBoundaryValues<dim>(present_time, present_timestep),
2074 * fe.component_mask(z_component));
2076 * PETScWrappers::MPI::Vector tmp(locally_owned_dofs, mpi_communicator);
2077 * MatrixTools::apply_boundary_values(
2078 * boundary_values, system_matrix, tmp, system_rhs, false);
2079 * incremental_displacement = tmp;
2087 * <a name="step_18-TopLevelsolve_timestep"></a>
2088 * <h4>TopLevel::solve_timestep</h4>
2092 * The next function is the one that controls what all has to happen within
2093 * a timestep. The order of things should be relatively self-explanatory
2094 * from the function names:
2097 * template <int dim>
2098 * void TopLevel<dim>::solve_timestep()
2100 * pcout << " Assembling system..." << std::flush;
2101 * assemble_system();
2102 * pcout << " norm of rhs is " << system_rhs.l2_norm() << std::endl;
2104 * const unsigned int n_iterations = solve_linear_problem();
2106 * pcout << " Solver converged in " << n_iterations << " iterations."
2109 * pcout << " Updating quadrature point data..." << std::flush;
2110 * update_quadrature_point_history();
2111 * pcout << std::endl;
2119 * <a name="step_18-TopLevelsolve_linear_problem"></a>
2120 * <h4>TopLevel::solve_linear_problem</h4>
2124 * Solving the linear system again works mostly as before. The only
2125 * difference is that we want to only keep a complete local copy of the
2126 * solution vector instead of the distributed one that we get as output from
2138 *
template <
int dim>
2142 *
locally_owned_dofs, mpi_communicator);
2152 *
cg.solve(system_matrix,
2161 *
return solver_control.last_step();
2169 * <a name=
"step_18-TopLeveloutput_results"></a>
2188 *
template <
int dim>
2192 *
data_out.attach_dof_handler(dof_handler);
2200 *
haven't implemented this case yet (another case of defensive
2204 * std::vector<std::string> solution_names;
2208 * solution_names.emplace_back("delta_x");
2211 * solution_names.emplace_back("delta_x");
2212 * solution_names.emplace_back("delta_y");
2215 * solution_names.emplace_back("delta_x");
2216 * solution_names.emplace_back("delta_y");
2217 * solution_names.emplace_back("delta_z");
2220 * DEAL_II_NOT_IMPLEMENTED();
2223 * data_out.add_data_vector(incremental_displacement, solution_names);
2228 * The next thing is that we wanted to output something like the average
2229 * norm of the stresses that we have stored in each cell. This may seem
2230 * complicated, since on the present processor we only store the stresses
2231 * in quadrature points on those cells that actually belong to the present
2232 * process. In other words, it seems as if we can't compute
the average
2236 * for all the other cells as this information would not be touched. The
2237 * following little loop does this. We enclose the entire block into a
2238 * pair of braces to make sure that the iterator variables do not remain
2239 * accidentally visible beyond the end of the block in which they are
2243 * Vector<double> norm_of_stress(triangulation.n_active_cells());
2247 * Loop over all the cells...
2250 * for (auto &cell : triangulation.active_cell_iterators())
2251 * if (cell->is_locally_owned())
2255 * On these cells, add up the stresses over all quadrature
2259 * SymmetricTensor<2, dim> accumulated_stress;
2260 * for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
2261 * accumulated_stress +=
2262 * reinterpret_cast<PointHistory<dim> *>(cell->user_pointer())[q]
2267 * ...then write the norm of the average to their destination:
2270 * norm_of_stress(cell->active_cell_index()) =
2271 * (accumulated_stress / quadrature_formula.size()).norm();
2275 * And on the cells that we are not interested in, set the respective
2276 * value in the vector to a bogus value (norms must be positive, and a
2277 * large negative value should catch your eye) in order to make sure
2278 * that if we were somehow wrong about our assumption that these
2279 * elements would not appear in the output file, that we would find out
2280 * by looking at the graphical output:
2284 * norm_of_stress(cell->active_cell_index()) = -1e+20;
2288 * Finally attach this vector as well to be treated for output:
2291 * data_out.add_data_vector(norm_of_stress, "norm_of_stress");
2295 * As a last piece of data, let us also add the partitioning of the domain
2296 * into subdomains associated with the processors if this is a parallel
2297 * job. This works in the exact same way as in the @ref step_17 "step-17" program:
2300 * std::vector<types::subdomain_id> partition_int(
2301 * triangulation.n_active_cells());
2302 * GridTools::get_subdomain_association(triangulation, partition_int);
2303 * const Vector<double> partitioning(partition_int.begin(),
2304 * partition_int.end());
2305 * data_out.add_data_vector(partitioning, "partitioning");
2309 * Finally, with all this data, we can instruct deal.II to munge the
2310 * information and produce some intermediate data structures that contain
2311 * all these solution and other data vectors:
2314 * data_out.build_patches();
2318 * Let us call a function that opens the necessary output files and writes
2319 * the data we have generated into them. The function automatically
2320 * constructs the file names from the given directory name (the first
2321 * argument) and file name base (second argument). It augments the resulting
2322 * string by pieces that result from the time step number and a "piece
2323 * number" that corresponds to a part of the overall domain that can consist
2324 * of one or more subdomains.
2328 * The function also writes a record files (with suffix `.pvd`) for Paraview
2329 * that describes how all of these output files combine into the data for
2330 * this single time step:
2333 * const std::string pvtu_filename = data_out.write_vtu_with_pvtu_record(
2334 * "./", "solution", timestep_no, mpi_communicator, 4);
2338 * The record files must be written only once and not by each processor,
2339 * so we do this on processor 0:
2342 * if (this_mpi_process == 0)
2346 * Finally, we write the paraview record, that references all .pvtu
2347 * files and their respective time. Note that the variable
2348 * times_and_names is declared static, so it will retain the entries
2349 * from the previous timesteps.
2352 * static std::vector<std::pair<double, std::string>> times_and_names;
2353 * times_and_names.emplace_back(present_time, pvtu_filename);
2354 * std::ofstream pvd_output("solution.pvd");
2355 * DataOutBase::write_pvd_record(pvd_output, times_and_names);
2364 * <a name="step_18-TopLeveldo_initial_timestep"></a>
2365 * <h4>TopLevel::do_initial_timestep</h4>
2369 * This and the next function handle the overall structure of the first and
2370 * following timesteps, respectively. The first timestep is slightly more
2371 * involved because we want to compute it multiple times on successively
2372 * refined meshes, each time starting from a clean state. At the end of
2373 * these computations, in which we compute the incremental displacements
2374 * each time, we use the last results obtained for the incremental
2375 * displacements to compute the resulting stress updates and move the mesh
2376 * accordingly. On this new mesh, we then output the solution and any
2377 * additional data we consider important.
2381 * All this is interspersed by generating output to the console to update
2382 * the person watching the screen on what is going on. As in @ref step_17 "step-17", the
2383 * use of <code>pcout</code> instead of <code>std::cout</code> makes sure
2384 * that only one of the parallel processes is actually writing to the
2385 * console, without having to explicitly code an if-statement in each place
2386 * where we generate output:
2389 * template <int dim>
2390 * void TopLevel<dim>::do_initial_timestep()
2392 * present_time += present_timestep;
2394 * pcout << "Timestep " << timestep_no << " at time " << present_time
2397 * for (unsigned int cycle = 0; cycle < 2; ++cycle)
2399 * pcout << " Cycle " << cycle << ':
' << std::endl;
2402 * create_coarse_grid();
2404 * refine_initial_grid();
2406 * pcout << " Number of active cells: "
2407 * << triangulation.n_active_cells() << " (by partition:";
2408 * for (unsigned int p = 0; p < n_mpi_processes; ++p)
2409 * pcout << (p == 0 ? ' ' : '+
')
2410 * << (GridTools::count_cells_with_subdomain_association(
2411 * triangulation, p));
2412 * pcout << ')
' << std::endl;
2416 * pcout << " Number of degrees of freedom: " << dof_handler.n_dofs()
2417 * << " (by partition:";
2418 * for (unsigned int p = 0; p < n_mpi_processes; ++p)
2419 * pcout << (p == 0 ? ' ' : '+
')
2420 * << (DoFTools::count_dofs_with_subdomain_association(dof_handler,
2422 * pcout << ')
' << std::endl;
2430 * pcout << std::endl;
2438 * <a name="step_18-TopLeveldo_timestep"></a>
2439 * <h4>TopLevel::do_timestep</h4>
2443 * Subsequent timesteps are simpler, and probably do not require any more
2444 * documentation given the explanations for the previous function above:
2447 * template <int dim>
2448 * void TopLevel<dim>::do_timestep()
2450 * present_time += present_timestep;
2452 * pcout << "Timestep " << timestep_no << " at time " << present_time
2454 * if (present_time > end_time)
2456 * present_timestep -= (present_time - end_time);
2457 * present_time = end_time;
2466 * pcout << std::endl;
2473 * <a name="step_18-TopLevelrefine_initial_grid"></a>
2474 * <h4>TopLevel::refine_initial_grid</h4>
2478 * The following function is called when solving the first time step on
2479 * successively refined meshes. After each iteration, it computes a
2480 * refinement criterion, refines the mesh, and sets up the history variables
2481 * in each quadrature point again to a clean state.
2484 * template <int dim>
2485 * void TopLevel<dim>::refine_initial_grid()
2489 * First, let each process compute error indicators for the cells it owns:
2492 * Vector<float> error_per_cell(triangulation.n_active_cells());
2493 * KellyErrorEstimator<dim>::estimate(
2495 * QGauss<dim - 1>(fe.degree + 1),
2496 * std::map<types::boundary_id, const Function<dim> *>(),
2497 * incremental_displacement,
2501 * MultithreadInfo::n_threads(),
2502 * this_mpi_process);
2506 * Then set up a global vector into which we merge the local indicators
2507 * from each of the %parallel processes:
2510 * const unsigned int n_local_cells =
2511 * triangulation.n_locally_owned_active_cells();
2513 * PETScWrappers::MPI::Vector distributed_error_per_cell(
2514 * mpi_communicator, triangulation.n_active_cells(), n_local_cells);
2516 * for (unsigned int i = 0; i < error_per_cell.size(); ++i)
2517 * if (error_per_cell(i) != 0)
2518 * distributed_error_per_cell(i) = error_per_cell(i);
2519 * distributed_error_per_cell.compress(VectorOperation::insert);
2523 * Once we have that, copy it back into local copies on all processors and
2524 * refine the mesh accordingly:
2527 * error_per_cell = distributed_error_per_cell;
2528 * GridRefinement::refine_and_coarsen_fixed_number(triangulation,
2532 * triangulation.execute_coarsening_and_refinement();
2536 * Finally, set up quadrature point data again on the new mesh, and only
2537 * on those cells that we have determined to be ours:
2540 * setup_quadrature_point_history();
2548 * <a name="step_18-TopLevelmove_mesh"></a>
2549 * <h4>TopLevel::move_mesh</h4>
2553 * At the end of each time step, we move the nodes of the mesh according to
2554 * the incremental displacements computed in this time step. To do this, we
2555 * keep a vector of flags that indicate for each vertex whether we have
2556 * already moved it around, and then loop over all cells and move those
2557 * vertices of the cell that have not been moved yet. It is worth noting
2558 * that it does not matter from which of the cells adjacent to a vertex we
2559 * move this vertex: since we compute the displacement using a continuous
2560 * finite element, the displacement field is continuous as well and we can
2561 * compute the displacement of a given vertex from each of the adjacent
2562 * cells. We only have to make sure that we move each node exactly once,
2563 * which is why we keep the vector of flags.
2567 * There are two noteworthy things in this function. First, how we get the
2568 * displacement field at a given vertex using the
2569 * <code>cell-@>vertex_dof_index(v,d)</code> function that returns the index
2570 * of the <code>d</code>th degree of freedom at vertex <code>v</code> of the
2571 * given cell. In the present case, displacement in the k-th coordinate
2572 * direction corresponds to the k-th component of the finite element. Using a
2573 * function like this bears a certain risk, because it uses knowledge of the
2574 * order of elements that we have taken together for this program in the
2575 * <code>FESystem</code> element. If we decided to add an additional
2576 * variable, for example a pressure variable for stabilization, and happened
2577 * to insert it as the first variable of the element, then the computation
2578 * below will start to produce nonsensical results. In addition, this
2579 * computation rests on other assumptions: first, that the element we use
2580 * has, indeed, degrees of freedom that are associated with vertices. This
2581 * is indeed the case for the present Q1 element, as would be for all Qp
2582 * elements of polynomial order <code>p</code>. However, it would not hold
2583 * for discontinuous elements, or elements for mixed formulations. Secondly,
2584 * it also rests on the assumption that the displacement at a vertex is
2585 * determined solely by the value of the degree of freedom associated with
2586 * this vertex; in other words, all shape functions corresponding to other
2587 * degrees of freedom are zero at this particular vertex. Again, this is the
2588 * case for the present element, but is not so for all elements that are
2589 * presently available in deal.II. Despite its risks, we choose to use this
2590 * way in order to present a way to query individual degrees of freedom
2591 * associated with vertices.
2595 * In this context, it is instructive to point out what a more general way
2596 * would be. For general finite elements, the way to go would be to take a
2597 * quadrature formula with the quadrature points in the vertices of a
2598 * cell. The <code>QTrapezoid</code> formula for the trapezoidal rule does
2599 * exactly this. With this quadrature formula, we would then initialize an
2600 * <code>FEValues</code> object in each cell, and use the
2601 * <code>FEValues::get_function_values</code> function to obtain the values
2602 * of the solution function in the quadrature points, i.e. the vertices of
2603 * the cell. These are the only values that we really need, i.e. we are not
2604 * at all interested in the weights (or the <code>JxW</code> values)
2605 * associated with this particular quadrature formula, and this can be
2606 * specified as the last argument in the constructor to
2607 * <code>FEValues</code>. The only point of minor inconvenience in this
2608 * scheme is that we have to figure out which quadrature point corresponds
2609 * to the vertex we consider at present, as they may or may not be ordered
2610 * in the same order.
2614 * This inconvenience could be avoided if finite elements have support
2615 * points on vertices (which the one here has; for the concept of support
2616 * points, see @ref GlossSupport "support points"). For such a case, one
2617 * could construct a custom quadrature rule using
2618 * FiniteElement::get_unit_support_points(). The first
2619 * <code>cell->n_vertices()*fe.dofs_per_vertex</code>
2620 * quadrature points will then correspond to the vertices of the cell and
2621 * are ordered consistent with <code>cell-@>vertex(i)</code>, taking into
2622 * account that support points for vector elements will be duplicated
2623 * <code>fe.dofs_per_vertex</code> times.
2627 * Another point worth explaining about this short function is the way in
2628 * which the triangulation class exports information about its vertices:
2629 * through the <code>Triangulation::n_vertices</code> function, it
2630 * advertises how many vertices there are in the triangulation. Not all of
2631 * them are actually in use all the time -- some are left-overs from cells
2632 * that have been coarsened previously and remain in existence since deal.II
2633 * never changes the number of a vertex once it has come into existence,
2634 * even if vertices with lower number go away. Secondly, the location
2635 * returned by <code>cell-@>vertex(v)</code> is not only a read-only object
2636 * of type <code>Point@<dim@></code>, but in fact a reference that can be
2637 * written to. This allows to move around the nodes of a mesh with relative
2638 * ease, but it is worth pointing out that it is the responsibility of an
2639 * application program using this feature to make sure that the resulting
2640 * cells are still useful, i.e. are not distorted so much that the cell is
2641 * degenerated (indicated, for example, by negative Jacobians). Note that we
2642 * do not have any provisions in this function to actually ensure this, we
2647 * After this lengthy introduction, here are the full 20 or so lines of
2651 * template <int dim>
2652 * void TopLevel<dim>::move_mesh()
2654 * pcout << " Moving mesh..." << std::endl;
2656 * std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
2657 * for (auto &cell : dof_handler.active_cell_iterators())
2658 * for (const auto v : cell->vertex_indices())
2659 * if (vertex_touched[cell->vertex_index(v)] == false)
2661 * vertex_touched[cell->vertex_index(v)] = true;
2663 * Point<dim> vertex_displacement;
2664 * for (unsigned int d = 0; d < dim; ++d)
2665 * vertex_displacement[d] =
2666 * incremental_displacement(cell->vertex_dof_index(v, d));
2668 * cell->vertex(v) += vertex_displacement;
2676 * <a name="step_18-TopLevelsetup_quadrature_point_history"></a>
2677 * <h4>TopLevel::setup_quadrature_point_history</h4>
2681 * At the beginning of our computations, we needed to set up initial values
2682 * of the history variables, such as the existing stresses in the material,
2683 * that we store in each quadrature point. As mentioned above, we use the
2684 * <code>user_pointer</code> for this that is available in each cell.
2688 * To put this into larger perspective, we note that if we had previously
2689 * available stresses in our model (which we assume do not exist for the
2690 * purpose of this program), then we would need to interpolate the field of
2691 * preexisting stresses to the quadrature points. Likewise, if we were to
2692 * simulate elasto-plastic materials with hardening/softening, then we would
2693 * have to store additional history variables like the present yield stress
2694 * of the accumulated plastic strains in each quadrature
2695 * points. Pre-existing hardening or weakening would then be implemented by
2696 * interpolating these variables in the present function as well.
2699 * template <int dim>
2700 * void TopLevel<dim>::setup_quadrature_point_history()
2704 * For good measure, we set all user pointers of all cells, whether
2705 * ours of not, to the null pointer. This way, if we ever access the user
2706 * pointer of a cell which we should not have accessed, a segmentation
2707 * fault will let us know that this should not have happened:
2713 * triangulation.clear_user_data();
2717 * Next, allocate the quadrature objects that are within the responsibility
2718 * of this processor. This, of course, equals the number of cells that
2719 * belong to this processor times the number of quadrature points our
2720 * quadrature formula has on each cell. Since the `resize()` function does
2721 * not actually shrink the amount of allocated memory if the requested new
2722 * size is smaller than the old size, we resort to a trick to first free all
2723 * memory, and then reallocate it: we declare an empty vector as a temporary
2724 * variable and then swap the contents of the old vector and this temporary
2725 * variable. This makes sure that the `quadrature_point_history` is now
2726 * really empty, and we can let the temporary variable that now holds the
2727 * previous contents of the vector go out of scope and be destroyed. In the
2728 * next step we can then re-allocate as many elements as we need, with the
2729 * vector default-initializing the `PointHistory` objects, which includes
2730 * setting the stress variables to zero.
2734 * std::vector<PointHistory<dim>> tmp;
2735 * quadrature_point_history.swap(tmp);
2737 * quadrature_point_history.resize(
2738 * triangulation.n_locally_owned_active_cells() * quadrature_formula.size());
2742 * Finally loop over all cells again and set the user pointers from the
2743 * cells that belong to the present processor to point to the first
2744 * quadrature point objects corresponding to this cell in the vector of
2748 * unsigned int history_index = 0;
2749 * for (auto &cell : triangulation.active_cell_iterators())
2750 * if (cell->is_locally_owned())
2752 * cell->set_user_pointer(&quadrature_point_history[history_index]);
2753 * history_index += quadrature_formula.size();
2758 * At the end, for good measure make sure that our count of elements was
2759 * correct and that we have both used up all objects we allocated
2760 * previously, and not point to any objects beyond the end of the
2761 * vector. Such defensive programming strategies are always good checks to
2762 * avoid accidental errors and to guard against future changes to this
2763 * function that forget to update all uses of a variable at the same
2764 * time. Recall that constructs using the <code>Assert</code> macro are
2765 * optimized away in optimized mode, so do not affect the run time of
2769 * Assert(history_index == quadrature_point_history.size(),
2770 * ExcInternalError());
2778 * <a name="step_18-TopLevelupdate_quadrature_point_history"></a>
2779 * <h4>TopLevel::update_quadrature_point_history</h4>
2783 * At the end of each time step, we should have computed an incremental
2784 * displacement update so that the material in its new configuration
2785 * accommodates for the difference between the external body and boundary
2786 * forces applied during this time step minus the forces exerted through
2787 * preexisting internal stresses. In order to have the preexisting
2788 * stresses available at the next time step, we therefore have to update the
2789 * preexisting stresses with the stresses due to the incremental
2790 * displacement computed during the present time step. Ideally, the
2791 * resulting sum of internal stresses would exactly counter all external
2792 * forces. Indeed, a simple experiment can make sure that this is so: if we
2793 * choose boundary conditions and body forces to be time independent, then
2794 * the forcing terms (the sum of external forces and internal stresses)
2795 * should be exactly zero. If you make this experiment, you will realize
2796 * from the output of the norm of the right hand side in each time step that
2797 * this is almost the case: it is not exactly zero, since in the first time
2798 * step the incremental displacement and stress updates were computed
2799 * relative to the undeformed mesh, which was then deformed. In the second
2800 * time step, we again compute displacement and stress updates, but this
2801 * time in the deformed mesh -- there, the resulting updates are very small
2802 * but not quite zero. This can be iterated, and in each such iteration the
2803 * residual, i.e. the norm of the right hand side vector, is reduced; if one
2804 * makes this little experiment, one realizes that the norm of this residual
2805 * decays exponentially with the number of iterations, and after an initial
2806 * very rapid decline is reduced by roughly a factor of about 3.5 in each
2807 * iteration (for one testcase I looked at, other testcases, and other
2808 * numbers of unknowns change the factor, but not the exponential decay).
2812 * In a sense, this can then be considered as a quasi-timestepping scheme to
2813 * resolve the nonlinear problem of solving large-deformation elasticity on
2814 * a mesh that is moved along in a Lagrangian manner.
2818 * Another complication is that the existing (old) stresses are defined on
2819 * the old mesh, which we will move around after updating the stresses. If
2820 * this mesh update involves rotations of the cell, then we need to also
2821 * rotate the updated stress, since it was computed relative to the
2822 * coordinate system of the old cell.
2826 * Thus, what we need is the following: on each cell which the present
2827 * processor owns, we need to extract the old stress from the data stored
2828 * with each quadrature point, compute the stress update, add the two
2829 * together, and then rotate the result together with the incremental
2830 * rotation computed from the incremental displacement at the present
2831 * quadrature point. We will detail these steps below:
2834 * template <int dim>
2835 * void TopLevel<dim>::update_quadrature_point_history()
2839 * First, set up an <code>FEValues</code> object by which we will evaluate
2840 * the incremental displacements and the gradients thereof at the
2841 * quadrature points, together with a vector that will hold this
2845 * FEValues<dim> fe_values(fe,
2846 * quadrature_formula,
2847 * update_values | update_gradients);
2849 * std::vector<std::vector<Tensor<1, dim>>> displacement_increment_grads(
2850 * quadrature_formula.size(), std::vector<Tensor<1, dim>>(dim));
2854 * Then loop over all cells and do the job in the cells that belong to our
2858 * for (auto &cell : dof_handler.active_cell_iterators())
2859 * if (cell->is_locally_owned())
2863 * Next, get a pointer to the quadrature point history data local to
2864 * the present cell, and, as a defensive measure, make sure that
2865 * this pointer is within the bounds of the global array:
2868 * PointHistory<dim> *local_quadrature_points_history =
2869 * reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
2870 * Assert(local_quadrature_points_history >=
2871 * &quadrature_point_history.front(),
2872 * ExcInternalError());
2873 * Assert(local_quadrature_points_history <=
2874 * &quadrature_point_history.back(),
2875 * ExcInternalError());
2879 * Then initialize the <code>FEValues</code> object on the present
2880 * cell, and extract the gradients of the displacement at the
2881 * quadrature points for later computation of the strains
2884 * fe_values.reinit(cell);
2885 * fe_values.get_function_gradients(incremental_displacement,
2886 * displacement_increment_grads);
2890 * Then loop over the quadrature points of this cell:
2893 * for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
2897 * On each quadrature point, compute the strain increment from
2898 * the gradients, and multiply it by the stress-strain tensor to
2899 * get the stress update. Then add this update to the already
2900 * existing strain at this point:
2903 * const SymmetricTensor<2, dim> new_stress =
2904 * (local_quadrature_points_history[q].old_stress +
2905 * (stress_strain_tensor *
2906 * get_strain(displacement_increment_grads[q])));
2910 * Finally, we have to rotate the result. For this, we first
2911 * have to compute a rotation matrix at the present quadrature
2912 * point from the incremental displacements. In fact, it can be
2913 * computed from the gradients, and we already have a function
2917 * const Tensor<2, dim> rotation =
2918 * get_rotation_matrix(displacement_increment_grads[q]);
2921 * Note that the result, a rotation matrix, is in general an
2922 * antisymmetric tensor of rank 2, so we must store it as a full
2927 * With this rotation matrix, we can compute the rotated tensor
2928 * by contraction from the left and right, after we expand the
2929 * symmetric tensor <code>new_stress</code> into a full tensor:
2932 * const SymmetricTensor<2, dim> rotated_new_stress =
2933 * symmetrize(transpose(rotation) *
2934 * static_cast<Tensor<2, dim>>(new_stress) * rotation);
2937 * Note that while the result of the multiplication of these
2938 * three matrices should be symmetric, it is not due to floating
2939 * point round off: we get an asymmetry on the order of 1e-16 of
2940 * the off-diagonal elements of the result. When assigning the
2941 * result to a <code>SymmetricTensor</code>, the constructor of
2942 * that class checks the symmetry and realizes that it isn't
2975 *
using namespace dealii;
2976 *
using namespace Step18;
2983 *
catch (std::exception &exc)
2985 *
std::cerr << std::endl
2987 *
<<
"----------------------------------------------------"
2989 *
std::cerr <<
"Exception on processing: " << std::endl
2990 *
<< exc.what() << std::endl
2991 *
<<
"Aborting!" << std::endl
2992 *
<<
"----------------------------------------------------"
2999 *
std::cerr << std::endl
3001 *
<<
"----------------------------------------------------"
3003 *
std::cerr <<
"Unknown exception!" << std::endl
3004 *
<<
"Aborting!" << std::endl
3005 *
<<
"----------------------------------------------------"
3149 <
img src=
"https://www.dealii.org/images/steps/developer/step-18.sequential-0002.0000.png"
3159 <
img src=
"https://www.dealii.org/images/steps/developer/step-18.sequential-0005.0000.png"
3169 <
img src=
"https://www.dealii.org/images/steps/developer/step-18.sequential-0007.0000.png"
3183 <
img src=
"https://www.dealii.org/images/steps/developer/step-18.sequential-0008.0000.png"
3193 <
img src=
"https://www.dealii.org/images/steps/developer/step-18.sequential-0009.0000.png"
3203 <
img src=
"https://www.dealii.org/images/steps/developer/step-18.sequential-0010.0000.png"
3229 Number
of active cells: 29696 (
by partition: 1808+1802+1894+1881+1870+1840+1884+1810+1876+1818+1870+1884+1854+1903+1816+1886)
3230 Number
of degrees
of freedom: 113100 (
by partition: 6936+6930+7305+7116+7326+6869+7331+6786+7193+6829+7093+7162+6920+7280+6843+7181)
3235 Number
of active cells: 102034 (
by partition: 6387+6202+6421+6341+6408+6201+6428+6428+6385+6294+6506+6244+6417+6527+6299+6546)
3236 Number
of degrees
of freedom: 359337 (
by partition: 23255+21308+24774+24019+22304+21415+22430+22184+22298+21796+22396+21592+22325+22553+21977+22711)
3268That's quite a good number of unknowns, given that we are in 3d. The output of
3269this program are 16 files for each time step:
3271\$ ls -l solution-0001*
3272-rw-r--r-- 1 wellsd2 user 761065 Feb 13 21:09 solution-0001.000.vtu
3273-rw-r--r-- 1 wellsd2 user 759277 Feb 13 21:09 solution-0001.001.vtu
3274-rw-r--r-- 1 wellsd2 user 761217 Feb 13 21:09 solution-0001.002.vtu
3275-rw-r--r-- 1 wellsd2 user 761605 Feb 13 21:09 solution-0001.003.vtu
3276-rw-r--r-- 1 wellsd2 user 756917 Feb 13 21:09 solution-0001.004.vtu
3277-rw-r--r-- 1 wellsd2 user 752669 Feb 13 21:09 solution-0001.005.vtu
3278-rw-r--r-- 1 wellsd2 user 735217 Feb 13 21:09 solution-0001.006.vtu
3279-rw-r--r-- 1 wellsd2 user 750065 Feb 13 21:09 solution-0001.007.vtu
3280-rw-r--r-- 1 wellsd2 user 760273 Feb 13 21:09 solution-0001.008.vtu
3281-rw-r--r-- 1 wellsd2 user 777265 Feb 13 21:09 solution-0001.009.vtu
3282-rw-r--r-- 1 wellsd2 user 772469 Feb 13 21:09 solution-0001.010.vtu
3283-rw-r--r-- 1 wellsd2 user 760833 Feb 13 21:09 solution-0001.011.vtu
3284-rw-r--r-- 1 wellsd2 user 782241 Feb 13 21:09 solution-0001.012.vtu
3285-rw-r--r-- 1 wellsd2 user 748905 Feb 13 21:09 solution-0001.013.vtu
3286-rw-r--r-- 1 wellsd2 user 738413 Feb 13 21:09 solution-0001.014.vtu
3287-rw-r--r-- 1 wellsd2 user 762133 Feb 13 21:09 solution-0001.015.vtu
3288-rw-r--r-- 1 wellsd2 user 1421 Feb 13 21:09 solution-0001.pvtu
3289-rw-r--r-- 1 wellsd2 user 364 Feb 13 21:09 solution-0001.visit
3293Here are first the mesh on which we compute as well as the partitioning
3294for the 16 processors:
3297<div class="twocolumn" style="width: 80%">
3298 <div class="parent">
3299 <div class="img" align="center">
3300 <img src="https://www.dealii.org/images/steps/developer/step-18.parallel-000mesh.png"
3301 alt="Discretization"
3304 <div class="text" align="center">
3308 <div class="parent">
3309 <div class="img" align="center">
3310 <img src="https://www.dealii.org/images/steps/developer/step-18.parallel-0002.p.png"
3311 alt="Parallel partitioning"
3314 <div class="text" align="center">
3315 Parallel partitioning
3321Finally, here is the same output as we have shown before for the much smaller
3324<div class="threecolumn" style="width: 80%">
3325 <div class="parent">
3326 <div class="img" align="center">
3327 <img src="https://www.dealii.org/images/steps/developer/step-18.parallel-0002.s.png"
3331 <div class="text" align="center">
3335 <div class="parent">
3336 <div class="img" align="center">
3337 <img src="https://www.dealii.org/images/steps/developer/step-18.parallel-0005.s.png"
3341 <div class="text" align="center">
3345 <div class="parent">
3346 <div class="img" align="center">
3347 <img src="https://www.dealii.org/images/steps/developer/step-18.parallel-0007.s.png"
3351 <div class="text" align="center">
3358<div class="threecolumn" style="width: 80%">
3359 <div class="parent">
3360 <div class="img" align="center">
3361 <img src="https://www.dealii.org/images/steps/developer/step-18.parallel-0008.s.png"
3365 <div class="text" align="center">
3369 <div class="parent">
3370 <div class="img" align="center">
3371 <img src="https://www.dealii.org/images/steps/developer/step-18.parallel-0009.s.png"
3375 <div class="text" align="center">
3379 <div class="parent">
3380 <div class="img" align="center">
3381 <img src="https://www.dealii.org/images/steps/developer/step-18.parallel-0010.s.png"
3385 <div class="text" align="center">
3392As before, we observe that at high axial compression the cylinder begins
3393to buckle, but this time ultimately collapses on itself. In contrast to our
3394first run, towards the end of the simulation the deflection pattern becomes
3395nonsymmetric (the central bulge deflects laterally). The model clearly does not
3396provide for this (all our forces and boundary deflections are symmetric) but the
3397effect is probably physically correct anyway: in reality, small inhomogeneities
3410the present one. One would need an even finer computation to find out. However,
3411the point may be moot: looking at the last picture in detail, it is pretty
3412obvious that not only is the linear small deformation model we chose completely
3413inadequate, but for a realistic simulation we would also need to make sure that
3414the body does not intersect itself during deformation (if we continued
3415compressing the cylinder we would observe some self-intersection).
3416Without such a formulation we cannot expect anything to make physical sense,
3417even if it produces nice pictures!
3420<a name="step_18-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
3423The program as is does not really solve an equation that has many applications
3424in practice: quasi-static material deformation based on a purely elastic law
3425is almost boring. However, the program may serve as the starting point for
3426more interesting experiments, and that indeed was the initial motivation for
3427writing it. Here are some suggestions of what the program is missing and in
3428what direction it may be extended:
3430<a name="step_18-Plasticitymodels"></a><h5>Plasticity models</h5>
3433 The most obvious extension is to use a more
3434realistic material model for large-scale quasistatic deformation. The natural
3435choice for this would be plasticity, in which a nonlinear relationship between
3436stress and strain replaces equation <a href="#step_18.stress-strain">[stress-strain]</a>. Plasticity
3437models are usually rather complicated to program since the stress-strain
3438dependence is generally non-smooth. The material can be thought of being able
3439to withstand only a maximal stress (the yield stress) after which it starts to
3440“flow”. A mathematical description to this can be given in the form of a
3441variational inequality, which alternatively can be treated as minimizing the
3445 (\varepsilon(\mathbf{u}), C\varepsilon(\mathbf{u}))_{\Omega}
3446 - (\mathbf{f}, \mathbf{u})_{\Omega} - (\mathbf{b}, \mathbf{u})_{\Gamma_N},
3448subject to the constraint
3450 f(\sigma(\mathbf{u})) \le 0
3452on the stress. This extension makes the problem to be solved in each time step
3453nonlinear, so we need another loop within each time step.
3455Without going into further details of this model, we refer to the excellent
3456book by Simo and Hughes on “Computational Inelasticity” for a
3457comprehensive overview of computational strategies for solving plastic
3458models. Alternatively, a brief but concise description of an algorithm for
3459plasticity is given in an article by S. Commend, A. Truty, and Th. Zimmermann;
3463<a name="step_18-Stabilizationissues"></a><h5>Stabilization issues</h5>
3466The formulation we have chosen, i.e. using
3467piecewise (bi-, tri-)linear elements for all components of the displacement
3468vector, and treating the stress as a variable dependent on the displacement is
3469appropriate for most materials. However, this so-called displacement-based
3470formulation becomes unstable and exhibits spurious modes for incompressible or
3471nearly-incompressible materials. While fluids are usually not elastic (in most
3472cases, the stress depends on velocity gradients, not displacement gradients,
3473although there are exceptions such as electro-rheologic fluids), there are a
3474few solids that are nearly incompressible, for example rubber. Another case is
3475that many plasticity models ultimately let the material become incompressible,
3476although this is outside the scope of the present program.
3478Incompressibility is characterized by Poisson's
ratio
3494<a name=
"step_18-Refinementduringtimesteps"></a><
h5>Refinement
during timesteps</
h5>
3523 per cell in 2d
and 8 in 3d.
3536 in
the quadrature points.)
3547 for (
unsigned int i=0; i<dim; ++i)
3548 for (
unsigned int j=0;
j<dim; ++
j)
3559 quadrature, quadrature,
3563 endc = dof_handler.end(),
3577 for (
unsigned int i=0; i<dim; ++i)
3578 for (
unsigned int j=0;
j<dim; ++
j)
3580 for (
unsigned int q=0;
q<quadrature.size(); ++
q)
3609 endc = dof_handler.end(),
3622 for (
unsigned int i=0; i<dim; ++i)
3623 for (
unsigned int j=0;
j<dim; ++
j)
3631 for (
unsigned int q=0;
q<quadrature.size(); ++
q)
3644<a name=
"step_18-Ensuringmeshregularity"></a><
h5>
Ensuring mesh regularity</
h5>
3660<a name="step_18-PlainProg"></a>
3661<h1> The plain program</h1>
3662@include "step-18.cc"
active_cell_iterator begin_active(const unsigned int level=0) const
static constexpr unsigned int rank
numbers::NumberTraits< Number >::real_type norm() const
static constexpr unsigned int dimension
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
void cylinder(Triangulation< dim > &tria, const double radius=1., const double half_length=1.)
void cylinder_shell(Triangulation< dim > &tria, const double length, const double inner_radius, const double outer_radius, const unsigned int n_radial_cells=0, const unsigned int n_axial_cells=0, const bool colorize=false)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
@ valid
Iterator points to a valid object.
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
@ general
No special properties.
constexpr types::blas_int one
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.)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
VectorType::value_type * end(VectorType &V)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
std::string compress(const std::string &input)
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)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
::VectorizedArray< Number, width > sqrt(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)