1438 *
const double g = 9.81;
1439 *
const double rho = 7700;
1442 *
values(dim - 1) = -rho * g;
1447 *
template <
int dim>
1448 *
void BodyForce<dim>::vector_value_list(
1452 *
const unsigned int n_points = points.size();
1454 *
Assert(value_list.size() == n_points,
1457 *
for (
unsigned int p = 0; p < n_points; ++p)
1458 * BodyForce<dim>::vector_value(points[p], value_list[p]);
1466 * <a name=
"ThecodeIncrementalBoundaryValuecodeclass"></a>
1467 * <h3>The <code>IncrementalBoundaryValue</code>
class</h3>
1471 * In addition to body forces, movement can be induced by boundary forces
1472 * and forced boundary displacement. The latter
case is equivalent to forces
1473 * being chosen in such a way that they induce certain displacement.
1477 * For quasistatic displacement, typical boundary forces would be pressure
1478 * on a body, or tangential friction against another body. We chose a
1479 * somewhat simpler
case here: we prescribe a certain movement of (parts of)
1480 * the boundary, or at least of certain components of the displacement
1481 * vector. We describe
this by another vector-valued function that,
for a
1482 * given
point on the boundary, returns the prescribed displacement.
1486 * Since we have a time-dependent problem, the displacement increment of the
1487 * boundary equals the displacement accumulated during the length of the
1488 * timestep. The
class therefore has to know both the present time and the
1489 * length of the present time step, and can then
approximate the incremental
1490 * displacement as the present velocity times the present timestep.
1494 * For the purposes of
this program, we choose a simple form of boundary
1495 * displacement: we displace the top boundary with constant velocity
1496 * downwards. The rest of the boundary is either going to be fixed (and is
1497 * then described
using an
object of type
1499 *
nothing special has to be done). The implementation of the
class
1500 * describing the constant downward motion should then be obvious
using the
1501 * knowledge we gained through all the previous example programs:
1504 *
template <
int dim>
1505 *
class IncrementalBoundaryValues :
public Function<dim>
1508 * IncrementalBoundaryValues(
const double present_time,
1509 *
const double present_timestep);
1519 *
const double velocity;
1520 *
const double present_time;
1521 *
const double present_timestep;
1525 *
template <
int dim>
1526 * IncrementalBoundaryValues<dim>::IncrementalBoundaryValues(
1527 *
const double present_time,
1528 *
const double present_timestep)
1531 * , present_time(present_time)
1532 * , present_timestep(present_timestep)
1536 *
template <
int dim>
1538 * IncrementalBoundaryValues<dim>::vector_value(
const Point<dim> & ,
1544 *
values(2) = -present_timestep * velocity;
1549 *
template <
int dim>
1550 *
void IncrementalBoundaryValues<dim>::vector_value_list(
1554 *
const unsigned int n_points = points.size();
1556 *
Assert(value_list.size() == n_points,
1559 *
for (
unsigned int p = 0; p < n_points; ++p)
1560 * IncrementalBoundaryValues<dim>::vector_value(points[p], value_list[p]);
1568 * <a name=
"ImplementationofthecodeTopLevelcodeclass"></a>
1569 * <h3>Implementation of the <code>TopLevel</code>
class</h3>
1573 * Now
for the implementation of the main
class. First, we initialize the
1574 * stress-strain tensor, which we have declared as a
static const
1575 * variable. We chose Lamé constants that are appropriate
for steel:
1578 *
template <
int dim>
1580 * get_stress_strain_tensor<dim>( 9.695e10,
1588 * <a name=
"Thepublicinterface"></a>
1589 * <h4>The
public interface</h4>
1593 * The next step is the definition of constructors and destructors. There
1594 * are no surprises here: we choose linear and continuous finite elements
1595 *
for each of the <code>dim</code> vector components of the solution, and a
1596 * Gaussian quadrature formula with 2 points in each coordinate
1597 * direction. The destructor should be obvious:
1600 *
template <
int dim>
1601 * TopLevel<dim>::TopLevel()
1605 * , quadrature_formula(fe.degree + 1)
1606 * , present_time(0.0)
1607 * , present_timestep(1.0)
1610 * , mpi_communicator(MPI_COMM_WORLD)
1618 *
template <
int dim>
1619 * TopLevel<dim>::~TopLevel()
1621 * dof_handler.clear();
1628 * The last of the
public functions is the
one that directs all the work,
1629 * <code>
run()</code>. It initializes the variables that describe where in
1630 * time we presently are, then runs the
first time step, then loops over all
1631 * the other time steps. Note that
for simplicity we use a fixed time step,
1632 * whereas a more sophisticated program would of course have to choose it in
1633 * some more reasonable way adaptively:
1636 *
template <
int dim>
1639 * do_initial_timestep();
1641 *
while (present_time < end_time)
1649 * <a name=
"TopLevelcreate_coarse_grid"></a>
1650 * <h4>TopLevel::create_coarse_grid</h4>
1654 * The next function in the order in which they were declared above is the
1655 *
one that creates the coarse grid from which we start. For
this example
1656 * program, we want to compute the deformation of a
cylinder under axial
1657 * compression. The
first step therefore is to generate a mesh
for a
1658 *
cylinder of length 3 and with inner and outer radii of 0.8 and 1,
1659 * respectively. Fortunately, there is a library function
for such a mesh.
1663 * In a
second step, we have to associated boundary conditions with the
1664 * upper and lower faces of the
cylinder. We choose a boundary indicator of
1665 * 0
for the boundary faces that are characterized by their midpoints having
1666 * z-coordinates of either 0 (bottom face), an indicator of 1
for z=3 (top
1667 * face);
finally, we use boundary indicator 2
for all faces on the inside
1668 * of the
cylinder shell, and 3
for the outside.
1671 *
template <
int dim>
1672 *
void TopLevel<dim>::create_coarse_grid()
1674 *
const double inner_radius = 0.8, outer_radius = 1;
1676 *
for (
const auto &cell :
triangulation.active_cell_iterators())
1677 *
for (
const auto &face : cell->face_iterators())
1678 *
if (face->at_boundary())
1680 *
const Point<dim> face_center = face->center();
1682 *
if (face_center[2] == 0)
1683 * face->set_boundary_id(0);
1684 *
else if (face_center[2] == 3)
1685 * face->set_boundary_id(1);
1686 *
else if (
std::sqrt(face_center[0] * face_center[0] +
1687 * face_center[1] * face_center[1]) <
1688 * (inner_radius + outer_radius) / 2)
1689 * face->set_boundary_id(2);
1691 * face->set_boundary_id(3);
1696 * Once all
this is done, we can
refine the mesh once globally:
1703 * As the
final step, we need to
set up a clean state of the data that we
1704 * store in the quadrature points on all cells that are treated on the
1705 * present processor.
1708 * setup_quadrature_point_history();
1716 * <a name=
"TopLevelsetup_system"></a>
1717 * <h4>TopLevel::setup_system</h4>
1721 * The next function is the
one that sets up the data structures
for a given
1722 * mesh. This is done in most the same way as in @ref step_17
"step-17": distribute the
1723 * degrees of freedom, then sort these degrees of freedom in such a way that
1724 * each processor gets a contiguous chunk of them. Note that subdivisions into
1725 * chunks
for each processor is handled in the
functions that create or
1726 *
refine grids, unlike in the previous example program (the
point where
1727 *
this happens is mostly a matter of taste; here, we chose to
do it when
1728 * grids are created since in the <code>do_initial_timestep</code> and
1729 * <code>do_timestep</code>
functions we want to output the number of cells
1730 * on each processor at a
point where we haven
't called the present function
1734 * template <int dim>
1735 * void TopLevel<dim>::setup_system()
1737 * dof_handler.distribute_dofs(fe);
1738 * locally_owned_dofs = dof_handler.locally_owned_dofs();
1739 * DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
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
1755 * which we simply used PETSc's ability to just know about the size of the
1756 *
matrix and later allocate those
nonzero elements that are being written
1757 * to. While
this works just fine from a correctness viewpoint, it is not
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
1799 * and we realize that we don't need all that much memory, but at time it
1800 * is already too late:
for large problems, the temporary allocation of
1801 * too much memory can lead to out-of-memory situations.
1805 * In order to avoid
this, we resort to the
1807 * not require any up-front estimate on the number of
nonzero entries per
1808 * row. It therefore only ever allocates as much memory as it needs at any
1809 * given time, and we can build it even
for large 3
d problems.
1813 * It is also worth noting that due to the specifics of
1815 * global, i.e. comprises all degrees of freedom whether they will be
1816 * owned by the processor we are on or another
one (in
case this program
1817 * is
run in %
parallel via MPI). This of course is not optimal -- it
1818 * limits the size of the problems we can solve, since storing the entire
1819 * sparsity pattern (even
if only
for a
short time) on each processor does
1820 * not
scale well. However, there are several more places in the program
1821 * in which we
do this,
for example we
always keep the global
1822 *
triangulation and DoF handler objects around, even
if we only work on
1823 * part of them. At present, deal.II does not have the necessary
1824 * facilities to completely distribute these objects (a task that, indeed,
1825 * is very hard to achieve with adaptive meshes, since well-balanced
1826 * subdivisions of a domain tend to become unbalanced as the mesh is
1827 * adaptively refined).
1831 * With
this data structure, we can then go to the PETSc sparse
matrix and
1832 * tell it to preallocate all the entries we will later want to write to:
1835 * system_matrix.reinit(locally_owned_dofs,
1836 * locally_owned_dofs,
1838 * mpi_communicator);
1841 * After
this point, no further
explicit knowledge of the sparsity pattern
1842 * is required any more and we can let the <code>sparsity_pattern</code>
1843 * variable go out of scope without any problem.
1847 * The last task in
this function is then only to reset the right hand
1848 * side vector as well as the solution vector to its correct size;
1849 * remember that the solution vector is a local
one, unlike the right hand
1850 * side that is a distributed %
parallel one and therefore needs to know
1851 * the MPI communicator over which it is supposed to transmit messages:
1854 * system_rhs.reinit(locally_owned_dofs, mpi_communicator);
1855 * incremental_displacement.reinit(dof_handler.n_dofs());
1863 * <a name=
"TopLevelassemble_system"></a>
1864 * <h4>TopLevel::assemble_system</h4>
1868 * Again, assembling the system
matrix and right hand side follows the same
1869 * structure as in many example programs before. In particular, it is mostly
1870 * equivalent to @ref step_17
"step-17", except
for the different right hand side that now
1871 * only has to take into account
internal stresses. In addition, assembling
1872 * the
matrix is made significantly more transparent by
using the
1873 * <code>
SymmetricTensor</code>
class: note the elegance of forming the
1874 * scalar products of
symmetric tensors of rank 2 and 4. The implementation
1875 * is also more
general since it is independent of the fact that we may or
1876 * may not be
using an isotropic elasticity tensor.
1880 * The
first part of the assembly routine is as
always:
1883 *
template <
int dim>
1884 *
void TopLevel<dim>::assemble_system()
1887 * system_matrix = 0;
1890 * quadrature_formula,
1894 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1895 *
const unsigned int n_q_points = quadrature_formula.size();
1900 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1902 * BodyForce<dim> body_force;
1903 * std::vector<Vector<double>> body_force_values(n_q_points,
1908 * As in @ref step_17
"step-17", we only need to
loop over all cells that belong to the
1909 * present processor:
1912 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1913 *
if (cell->is_locally_owned())
1918 * fe_values.reinit(cell);
1922 * Then
loop over all indices i,j and quadrature points and
assemble
1923 * the system
matrix contributions from
this cell. Note how we
1925 * at a given quadrature
point from the <code>
FEValues</code>
1926 * object, and the elegance with which we form the triple
1927 * contraction <code>eps_phi_i :
C : eps_phi_j</code>; the latter
1928 * needs to be compared to the clumsy computations needed in
1929 * @ref step_17
"step-17", both in the introduction as well as in the respective
1930 * place in the program:
1933 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1934 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1935 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1938 * eps_phi_i = get_strain(fe_values, i, q_point),
1939 * eps_phi_j = get_strain(fe_values, j, q_point);
1942 * stress_strain_tensor *
1945 * fe_values.JxW(q_point);
1951 * Then also
assemble the local right hand side contributions. For
1952 *
this, we need to access the prior stress
value in
this quadrature
1953 *
point. To get it, we use the user pointer of
this cell that
1954 * points into the global array to the quadrature
point data
1955 * corresponding to the
first quadrature
point of the present cell,
1956 * and then add an offset corresponding to the index of the
1957 * quadrature
point we presently consider:
1960 *
const PointHistory<dim> *local_quadrature_points_data =
1961 *
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
1964 * In addition, we need the
values of the external body forces at
1965 * the quadrature points on
this cell:
1968 * body_force.vector_value_list(fe_values.get_quadrature_points(),
1969 * body_force_values);
1972 * Then we can
loop over all degrees of freedom on
this cell and
1973 * compute local contributions to the right hand side:
1976 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1978 *
const unsigned int component_i =
1979 * fe.system_to_component_index(i).first;
1981 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1984 * local_quadrature_points_data[q_point].old_stress;
1987 * (body_force_values[q_point](component_i) *
1988 * fe_values.shape_value(i, q_point) -
1989 * old_stress * get_strain(fe_values, i, q_point)) *
1990 * fe_values.JxW(q_point);
1996 * Now that we have the local contributions to the linear system, we
1997 * need to transfer it into the global objects. This is done exactly
1998 * as in @ref step_17
"step-17":
2001 * cell->get_dof_indices(local_dof_indices);
2003 * hanging_node_constraints.distribute_local_to_global(
cell_matrix,
2005 * local_dof_indices,
2021 * The last step is to again fix up boundary
values, just as we already
2022 * did in previous programs.
A slight complication is that the
2024 * vector compatible with the
matrix and right hand side (i.e. here a
2025 * distributed %
parallel vector, rather than the sequential vector we use
2026 * in
this program) in order to preset the entries of the solution vector
2027 * with the correct boundary
values. We provide such a compatible vector
2028 * in the form of a temporary vector which we then
copy into the
2033 * We make up
for this complication by showing how boundary
values can be
2034 * used flexibly: following the way we create the
triangulation, there are
2035 * three distinct boundary indicators used to describe the domain,
2036 * corresponding to the bottom and top faces, as well as the inner/outer
2037 * surfaces. We would like to impose boundary conditions of the following
2038 * type: The inner and outer
cylinder surfaces are
free of external
2039 * forces, a fact that corresponds to natural (Neumann-type) boundary
2040 * conditions
for which we don
't have to do anything. At the bottom, we
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 * 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="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="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
2127 * PETSc's solver routines. To
this end, we declare a local temporary
2128 * variable
for the distributed vector and initialize it with the contents
2129 * of the local variable (remember that the
2130 * <code>apply_boundary_values</code> function called in
2131 * <code>assemble_system</code> preset the
values of boundary nodes in
this
2132 * vector), solve with it, and at the
end of the function
copy it again into
2133 * the complete local vector that we declared as a member variable. Hanging
2134 * node constraints are then distributed only on the local
copy,
2135 * i.e. independently of each other on each of the processors:
2138 *
template <
int dim>
2139 *
unsigned int TopLevel<dim>::solve_linear_problem()
2142 * locally_owned_dofs, mpi_communicator);
2143 * distributed_incremental_displacement = incremental_displacement;
2146 * 1
e-16 * system_rhs.l2_norm());
2152 * cg.solve(system_matrix,
2153 * distributed_incremental_displacement,
2157 * incremental_displacement = distributed_incremental_displacement;
2159 * hanging_node_constraints.distribute(incremental_displacement);
2161 *
return solver_control.last_step();
2169 * <a name=
"TopLeveloutput_results"></a>
2170 * <h4>TopLevel::output_results</h4>
2174 * This function generates the graphical output in .vtu format as explained
2175 * in the introduction. Each process will only work on the cells it owns,
2176 * and then write the result into a file of its own. Additionally, processor
2177 * 0 will write the record files the reference all the .vtu files.
2181 * The crucial part of
this function is to give the <code>
DataOut</code>
2182 *
class a way to only work on the cells that the present process owns.
2188 *
template <
int dim>
2189 *
void TopLevel<dim>::output_results() const
2196 * Then, just as in @ref step_17
"step-17", define the names of solution variables (which
2197 * here are the displacement increments) and queue the solution vector
for
2198 * output. Note in the following
switch how we make sure that
if the space
2199 * dimension should be unhandled that we
throw an exception saying that we
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 * Assert(false, ExcNotImplemented());
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
2233 * stresses
for all cells. However, remember that our
class derived from
2234 * <code>
DataOut</code> only iterates over those cells that actually
do
2235 * belong to the present processor, i.e. we don
't have to compute anything
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.push_back(
2354 * std::pair<double, std::string>(present_time, pvtu_filename));
2355 * std::ofstream pvd_output("solution.pvd");
2356 * DataOutBase::write_pvd_record(pvd_output, times_and_names);
2365 * <a name="TopLeveldo_initial_timestep"></a>
2366 * <h4>TopLevel::do_initial_timestep</h4>
2370 * This and the next function handle the overall structure of the first and
2371 * following timesteps, respectively. The first timestep is slightly more
2372 * involved because we want to compute it multiple times on successively
2373 * refined meshes, each time starting from a clean state. At the end of
2374 * these computations, in which we compute the incremental displacements
2375 * each time, we use the last results obtained for the incremental
2376 * displacements to compute the resulting stress updates and move the mesh
2377 * accordingly. On this new mesh, we then output the solution and any
2378 * additional data we consider important.
2382 * All this is interspersed by generating output to the console to update
2383 * the person watching the screen on what is going on. As in @ref step_17 "step-17", the
2384 * use of <code>pcout</code> instead of <code>std::cout</code> makes sure
2385 * that only one of the parallel processes is actually writing to the
2386 * console, without having to explicitly code an if-statement in each place
2387 * where we generate output:
2390 * template <int dim>
2391 * void TopLevel<dim>::do_initial_timestep()
2393 * present_time += present_timestep;
2395 * pcout << "Timestep " << timestep_no << " at time " << present_time
2398 * for (unsigned int cycle = 0; cycle < 2; ++cycle)
2400 * pcout << " Cycle " << cycle << ':
' << std::endl;
2403 * create_coarse_grid();
2405 * refine_initial_grid();
2407 * pcout << " Number of active cells: "
2408 * << triangulation.n_active_cells() << " (by partition:";
2409 * for (unsigned int p = 0; p < n_mpi_processes; ++p)
2410 * pcout << (p == 0 ? ' ' : '+
')
2411 * << (GridTools::count_cells_with_subdomain_association(
2412 * triangulation, p));
2413 * pcout << ")" << std::endl;
2417 * pcout << " Number of degrees of freedom: " << dof_handler.n_dofs()
2418 * << " (by partition:";
2419 * for (unsigned int p = 0; p < n_mpi_processes; ++p)
2420 * pcout << (p == 0 ? ' ' : '+
')
2421 * << (DoFTools::count_dofs_with_subdomain_association(dof_handler,
2423 * pcout << ")" << std::endl;
2431 * pcout << std::endl;
2439 * <a name="TopLeveldo_timestep"></a>
2440 * <h4>TopLevel::do_timestep</h4>
2444 * Subsequent timesteps are simpler, and probably do not require any more
2445 * documentation given the explanations for the previous function above:
2448 * template <int dim>
2449 * void TopLevel<dim>::do_timestep()
2451 * present_time += present_timestep;
2453 * pcout << "Timestep " << timestep_no << " at time " << present_time
2455 * if (present_time > end_time)
2457 * present_timestep -= (present_time - end_time);
2458 * present_time = end_time;
2467 * pcout << std::endl;
2474 * <a name="TopLevelrefine_initial_grid"></a>
2475 * <h4>TopLevel::refine_initial_grid</h4>
2479 * The following function is called when solving the first time step on
2480 * successively refined meshes. After each iteration, it computes a
2481 * refinement criterion, refines the mesh, and sets up the history variables
2482 * in each quadrature point again to a clean state.
2485 * template <int dim>
2486 * void TopLevel<dim>::refine_initial_grid()
2490 * First, let each process compute error indicators for the cells it owns:
2493 * Vector<float> error_per_cell(triangulation.n_active_cells());
2494 * KellyErrorEstimator<dim>::estimate(
2496 * QGauss<dim - 1>(fe.degree + 1),
2497 * std::map<types::boundary_id, const Function<dim> *>(),
2498 * incremental_displacement,
2502 * MultithreadInfo::n_threads(),
2503 * this_mpi_process);
2507 * Then set up a global vector into which we merge the local indicators
2508 * from each of the %parallel processes:
2511 * const unsigned int n_local_cells =
2512 * triangulation.n_locally_owned_active_cells();
2514 * PETScWrappers::MPI::Vector distributed_error_per_cell(
2515 * mpi_communicator, triangulation.n_active_cells(), n_local_cells);
2517 * for (unsigned int i = 0; i < error_per_cell.size(); ++i)
2518 * if (error_per_cell(i) != 0)
2519 * distributed_error_per_cell(i) = error_per_cell(i);
2520 * distributed_error_per_cell.compress(VectorOperation::insert);
2524 * Once we have that, copy it back into local copies on all processors and
2525 * refine the mesh accordingly:
2528 * error_per_cell = distributed_error_per_cell;
2529 * GridRefinement::refine_and_coarsen_fixed_number(triangulation,
2533 * triangulation.execute_coarsening_and_refinement();
2537 * Finally, set up quadrature point data again on the new mesh, and only
2538 * on those cells that we have determined to be ours:
2541 * setup_quadrature_point_history();
2549 * <a name="TopLevelmove_mesh"></a>
2550 * <h4>TopLevel::move_mesh</h4>
2554 * At the end of each time step, we move the nodes of the mesh according to
2555 * the incremental displacements computed in this time step. To do this, we
2556 * keep a vector of flags that indicate for each vertex whether we have
2557 * already moved it around, and then loop over all cells and move those
2558 * vertices of the cell that have not been moved yet. It is worth noting
2559 * that it does not matter from which of the cells adjacent to a vertex we
2560 * move this vertex: since we compute the displacement using a continuous
2561 * finite element, the displacement field is continuous as well and we can
2562 * compute the displacement of a given vertex from each of the adjacent
2563 * cells. We only have to make sure that we move each node exactly once,
2564 * which is why we keep the vector of flags.
2568 * There are two noteworthy things in this function. First, how we get the
2569 * displacement field at a given vertex using the
2570 * <code>cell-@>vertex_dof_index(v,d)</code> function that returns the index
2571 * of the <code>d</code>th degree of freedom at vertex <code>v</code> of the
2572 * given cell. In the present case, displacement in the k-th coordinate
2573 * direction corresponds to the k-th component of the finite element. Using a
2574 * function like this bears a certain risk, because it uses knowledge of the
2575 * order of elements that we have taken together for this program in the
2576 * <code>FESystem</code> element. If we decided to add an additional
2577 * variable, for example a pressure variable for stabilization, and happened
2578 * to insert it as the first variable of the element, then the computation
2579 * below will start to produce nonsensical results. In addition, this
2580 * computation rests on other assumptions: first, that the element we use
2581 * has, indeed, degrees of freedom that are associated with vertices. This
2582 * is indeed the case for the present Q1 element, as would be for all Qp
2583 * elements of polynomial order <code>p</code>. However, it would not hold
2584 * for discontinuous elements, or elements for mixed formulations. Secondly,
2585 * it also rests on the assumption that the displacement at a vertex is
2586 * determined solely by the value of the degree of freedom associated with
2587 * this vertex; in other words, all shape functions corresponding to other
2588 * degrees of freedom are zero at this particular vertex. Again, this is the
2589 * case for the present element, but is not so for all elements that are
2590 * presently available in deal.II. Despite its risks, we choose to use this
2591 * way in order to present a way to query individual degrees of freedom
2592 * associated with vertices.
2596 * In this context, it is instructive to point out what a more general way
2597 * would be. For general finite elements, the way to go would be to take a
2598 * quadrature formula with the quadrature points in the vertices of a
2599 * cell. The <code>QTrapezoid</code> formula for the trapezoidal rule does
2600 * exactly this. With this quadrature formula, we would then initialize an
2601 * <code>FEValues</code> object in each cell, and use the
2602 * <code>FEValues::get_function_values</code> function to obtain the values
2603 * of the solution function in the quadrature points, i.e. the vertices of
2604 * the cell. These are the only values that we really need, i.e. we are not
2605 * at all interested in the weights (or the <code>JxW</code> values)
2606 * associated with this particular quadrature formula, and this can be
2607 * specified as the last argument in the constructor to
2608 * <code>FEValues</code>. The only point of minor inconvenience in this
2609 * scheme is that we have to figure out which quadrature point corresponds
2610 * to the vertex we consider at present, as they may or may not be ordered
2611 * in the same order.
2615 * This inconvenience could be avoided if finite elements have support
2616 * points on vertices (which the one here has; for the concept of support
2617 * points, see @ref GlossSupport "support points"). For such a case, one
2618 * could construct a custom quadrature rule using
2619 * FiniteElement::get_unit_support_points(). The first
2620 * <code>cell->n_vertices()*fe.dofs_per_vertex</code>
2621 * quadrature points will then correspond to the vertices of the cell and
2622 * are ordered consistent with <code>cell-@>vertex(i)</code>, taking into
2623 * account that support points for vector elements will be duplicated
2624 * <code>fe.dofs_per_vertex</code> times.
2628 * Another point worth explaining about this short function is the way in
2629 * which the triangulation class exports information about its vertices:
2630 * through the <code>Triangulation::n_vertices</code> function, it
2631 * advertises how many vertices there are in the triangulation. Not all of
2632 * them are actually in use all the time -- some are left-overs from cells
2633 * that have been coarsened previously and remain in existence since deal.II
2634 * never changes the number of a vertex once it has come into existence,
2635 * even if vertices with lower number go away. Secondly, the location
2636 * returned by <code>cell-@>vertex(v)</code> is not only a read-only object
2637 * of type <code>Point@<dim@></code>, but in fact a reference that can be
2638 * written to. This allows to move around the nodes of a mesh with relative
2639 * ease, but it is worth pointing out that it is the responsibility of an
2640 * application program using this feature to make sure that the resulting
2641 * cells are still useful, i.e. are not distorted so much that the cell is
2642 * degenerated (indicated, for example, by negative Jacobians). Note that we
2643 * do not have any provisions in this function to actually ensure this, we
2648 * After this lengthy introduction, here are the full 20 or so lines of
2652 * template <int dim>
2653 * void TopLevel<dim>::move_mesh()
2655 * pcout << " Moving mesh..." << std::endl;
2657 * std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
2658 * for (auto &cell : dof_handler.active_cell_iterators())
2659 * for (const auto v : cell->vertex_indices())
2660 * if (vertex_touched[cell->vertex_index(v)] == false)
2662 * vertex_touched[cell->vertex_index(v)] = true;
2664 * Point<dim> vertex_displacement;
2665 * for (unsigned int d = 0; d < dim; ++d)
2666 * vertex_displacement[d] =
2667 * incremental_displacement(cell->vertex_dof_index(v, d));
2669 * cell->vertex(v) += vertex_displacement;
2677 * <a name="TopLevelsetup_quadrature_point_history"></a>
2678 * <h4>TopLevel::setup_quadrature_point_history</h4>
2682 * At the beginning of our computations, we needed to set up initial values
2683 * of the history variables, such as the existing stresses in the material,
2684 * that we store in each quadrature point. As mentioned above, we use the
2685 * <code>user_pointer</code> for this that is available in each cell.
2689 * To put this into larger perspective, we note that if we had previously
2690 * available stresses in our model (which we assume do not exist for the
2691 * purpose of this program), then we would need to interpolate the field of
2692 * preexisting stresses to the quadrature points. Likewise, if we were to
2693 * simulate elasto-plastic materials with hardening/softening, then we would
2694 * have to store additional history variables like the present yield stress
2695 * of the accumulated plastic strains in each quadrature
2696 * points. Pre-existing hardening or weakening would then be implemented by
2697 * interpolating these variables in the present function as well.
2700 * template <int dim>
2701 * void TopLevel<dim>::setup_quadrature_point_history()
2705 * For good measure, we set all user pointers of all cells, whether
2706 * ours of not, to the null pointer. This way, if we ever access the user
2707 * pointer of a cell which we should not have accessed, a segmentation
2708 * fault will let us know that this should not have happened:
2714 * triangulation.clear_user_data();
2718 * Next, allocate the quadrature objects that are within the responsibility
2719 * of this processor. This, of course, equals the number of cells that
2720 * belong to this processor times the number of quadrature points our
2721 * quadrature formula has on each cell. Since the `resize()` function does
2722 * not actually shrink the amount of allocated memory if the requested new
2723 * size is smaller than the old size, we resort to a trick to first free all
2724 * memory, and then reallocate it: we declare an empty vector as a temporary
2725 * variable and then swap the contents of the old vector and this temporary
2726 * variable. This makes sure that the `quadrature_point_history` is now
2727 * really empty, and we can let the temporary variable that now holds the
2728 * previous contents of the vector go out of scope and be destroyed. In the
2729 * next step we can then re-allocate as many elements as we need, with the
2730 * vector default-initializing the `PointHistory` objects, which includes
2731 * setting the stress variables to zero.
2735 * std::vector<PointHistory<dim>> tmp;
2736 * quadrature_point_history.swap(tmp);
2738 * quadrature_point_history.resize(
2739 * triangulation.n_locally_owned_active_cells() * quadrature_formula.size());
2743 * Finally loop over all cells again and set the user pointers from the
2744 * cells that belong to the present processor to point to the first
2745 * quadrature point objects corresponding to this cell in the vector of
2749 * unsigned int history_index = 0;
2750 * for (auto &cell : triangulation.active_cell_iterators())
2751 * if (cell->is_locally_owned())
2753 * cell->set_user_pointer(&quadrature_point_history[history_index]);
2754 * history_index += quadrature_formula.size();
2759 * At the end, for good measure make sure that our count of elements was
2760 * correct and that we have both used up all objects we allocated
2761 * previously, and not point to any objects beyond the end of the
2762 * vector. Such defensive programming strategies are always good checks to
2763 * avoid accidental errors and to guard against future changes to this
2764 * function that forget to update all uses of a variable at the same
2765 * time. Recall that constructs using the <code>Assert</code> macro are
2766 * optimized away in optimized mode, so do not affect the run time of
2770 * Assert(history_index == quadrature_point_history.size(),
2771 * ExcInternalError());
2779 * <a name="TopLevelupdate_quadrature_point_history"></a>
2780 * <h4>TopLevel::update_quadrature_point_history</h4>
2784 * At the end of each time step, we should have computed an incremental
2785 * displacement update so that the material in its new configuration
2786 * accommodates for the difference between the external body and boundary
2787 * forces applied during this time step minus the forces exerted through
2788 * preexisting internal stresses. In order to have the preexisting
2789 * stresses available at the next time step, we therefore have to update the
2790 * preexisting stresses with the stresses due to the incremental
2791 * displacement computed during the present time step. Ideally, the
2792 * resulting sum of internal stresses would exactly counter all external
2793 * forces. Indeed, a simple experiment can make sure that this is so: if we
2794 * choose boundary conditions and body forces to be time independent, then
2795 * the forcing terms (the sum of external forces and internal stresses)
2796 * should be exactly zero. If you make this experiment, you will realize
2797 * from the output of the norm of the right hand side in each time step that
2798 * this is almost the case: it is not exactly zero, since in the first time
2799 * step the incremental displacement and stress updates were computed
2800 * relative to the undeformed mesh, which was then deformed. In the second
2801 * time step, we again compute displacement and stress updates, but this
2802 * time in the deformed mesh -- there, the resulting updates are very small
2803 * but not quite zero. This can be iterated, and in each such iteration the
2804 * residual, i.e. the norm of the right hand side vector, is reduced; if one
2805 * makes this little experiment, one realizes that the norm of this residual
2806 * decays exponentially with the number of iterations, and after an initial
2807 * very rapid decline is reduced by roughly a factor of about 3.5 in each
2808 * iteration (for one testcase I looked at, other testcases, and other
2809 * numbers of unknowns change the factor, but not the exponential decay).
2813 * In a sense, this can then be considered as a quasi-timestepping scheme to
2814 * resolve the nonlinear problem of solving large-deformation elasticity on
2815 * a mesh that is moved along in a Lagrangian manner.
2819 * Another complication is that the existing (old) stresses are defined on
2820 * the old mesh, which we will move around after updating the stresses. If
2821 * this mesh update involves rotations of the cell, then we need to also
2822 * rotate the updated stress, since it was computed relative to the
2823 * coordinate system of the old cell.
2827 * Thus, what we need is the following: on each cell which the present
2828 * processor owns, we need to extract the old stress from the data stored
2829 * with each quadrature point, compute the stress update, add the two
2830 * together, and then rotate the result together with the incremental
2831 * rotation computed from the incremental displacement at the present
2832 * quadrature point. We will detail these steps below:
2835 * template <int dim>
2836 * void TopLevel<dim>::update_quadrature_point_history()
2840 * First, set up an <code>FEValues</code> object by which we will evaluate
2841 * the incremental displacements and the gradients thereof at the
2842 * quadrature points, together with a vector that will hold this
2846 * FEValues<dim> fe_values(fe,
2847 * quadrature_formula,
2848 * update_values | update_gradients);
2850 * std::vector<std::vector<Tensor<1, dim>>> displacement_increment_grads(
2851 * quadrature_formula.size(), std::vector<Tensor<1, dim>>(dim));
2855 * Then loop over all cells and do the job in the cells that belong to our
2859 * for (auto &cell : dof_handler.active_cell_iterators())
2860 * if (cell->is_locally_owned())
2864 * Next, get a pointer to the quadrature point history data local to
2865 * the present cell, and, as a defensive measure, make sure that
2866 * this pointer is within the bounds of the global array:
2869 * PointHistory<dim> *local_quadrature_points_history =
2870 * reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
2871 * Assert(local_quadrature_points_history >=
2872 * &quadrature_point_history.front(),
2873 * ExcInternalError());
2874 * Assert(local_quadrature_points_history <=
2875 * &quadrature_point_history.back(),
2876 * ExcInternalError());
2880 * Then initialize the <code>FEValues</code> object on the present
2881 * cell, and extract the gradients of the displacement at the
2882 * quadrature points for later computation of the strains
2885 * fe_values.reinit(cell);
2886 * fe_values.get_function_gradients(incremental_displacement,
2887 * displacement_increment_grads);
2891 * Then loop over the quadrature points of this cell:
2894 * for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
2898 * On each quadrature point, compute the strain increment from
2899 * the gradients, and multiply it by the stress-strain tensor to
2900 * get the stress update. Then add this update to the already
2901 * existing strain at this point:
2904 * const SymmetricTensor<2, dim> new_stress =
2905 * (local_quadrature_points_history[q].old_stress +
2906 * (stress_strain_tensor *
2907 * get_strain(displacement_increment_grads[q])));
2911 * Finally, we have to rotate the result. For this, we first
2912 * have to compute a rotation matrix at the present quadrature
2913 * point from the incremental displacements. In fact, it can be
2914 * computed from the gradients, and we already have a function
2918 * const Tensor<2, dim> rotation =
2919 * get_rotation_matrix(displacement_increment_grads[q]);
2922 * Note that the result, a rotation matrix, is in general an
2923 * antisymmetric tensor of rank 2, so we must store it as a full
2928 * With this rotation matrix, we can compute the rotated tensor
2929 * by contraction from the left and right, after we expand the
2930 * symmetric tensor <code>new_stress</code> into a full tensor:
2933 * const SymmetricTensor<2, dim> rotated_new_stress =
2934 * symmetrize(transpose(rotation) *
2935 * static_cast<Tensor<2, dim>>(new_stress) * rotation);
2938 * Note that while the result of the multiplication of these
2939 * three matrices should be symmetric, it is not due to floating
2940 * point round off: we get an asymmetry on the order of 1e-16 of
2941 * the off-diagonal elements of the result. When assigning the
2942 * result to a <code>SymmetricTensor</code>, the constructor of
2943 * that class checks the symmetry and realizes that it isn't
2944 * exactly
symmetric; it will then raise an exception. To avoid
2945 * that, we explicitly
symmetrize the result to make it exactly
2950 * The result of all these operations is then written back into
2951 * the original place:
2954 * local_quadrature_points_history[q].old_stress =
2955 * rotated_new_stress;
2962 * This ends the
project specific namespace <code>Step18</code>. The rest is
2963 * as usual and as already shown in @ref step_17
"step-17": A <code>main()</code> function
2964 * that initializes and terminates PETSc, calls the classes that
do the
2965 * actual work, and makes sure that we
catch all exceptions that propagate
2972 *
int main(
int argc,
char **argv)
2976 *
using namespace dealii;
2977 *
using namespace Step18;
2981 * TopLevel<3> elastic_problem;
2982 * elastic_problem.run();
2984 *
catch (std::exception &exc)
2986 * std::cerr << std::endl
2988 * <<
"----------------------------------------------------"
2990 * std::cerr <<
"Exception on processing: " << std::endl
2991 * << exc.what() << std::endl
2992 * <<
"Aborting!" << std::endl
2993 * <<
"----------------------------------------------------"
3000 * std::cerr << std::endl
3002 * <<
"----------------------------------------------------"
3004 * std::cerr <<
"Unknown exception!" << std::endl
3005 * <<
"Aborting!" << std::endl
3006 * <<
"----------------------------------------------------"
3014<a name=
"Results"></a><h1>Results</h1>
3018Running the program takes a good
while if one uses debug mode; it takes about
3019eleven minutes on my i7 desktop. Fortunately, the version compiled with
3020optimizations is much faster; the program only takes about a minute and a half
3021after recompiling with the command <tt>make release</tt> on the same machine, a
3022much more reasonable time.
3025If
run, the program prints the following output, explaining what it is
3026doing during all that time:
3029[ 66%] Built target step-18
3030[100%] Run step-18 with Release configuration
3033 Number of active cells: 3712 (by
partition: 3712)
3034 Number of degrees of freedom: 17226 (by
partition: 17226)
3035 Assembling system...
norm of rhs is 1.88062e+10
3036 Solver converged in 103 iterations.
3037 Updating quadrature
point data...
3039 Number of active cells: 12812 (by
partition: 12812)
3040 Number of degrees of freedom: 51738 (by
partition: 51738)
3041 Assembling system...
norm of rhs is 1.86145e+10
3042 Solver converged in 121 iterations.
3043 Updating quadrature
point data...
3047 Assembling system...
norm of rhs is 1.84169e+10
3048 Solver converged in 122 iterations.
3049 Updating quadrature
point data...
3053 Assembling system...
norm of rhs is 1.82355e+10
3054 Solver converged in 122 iterations.
3055 Updating quadrature
point data...
3059 Assembling system...
norm of rhs is 1.80728e+10
3060 Solver converged in 117 iterations.
3061 Updating quadrature
point data...
3065 Assembling system...
norm of rhs is 1.79318e+10
3066 Solver converged in 116 iterations.
3067 Updating quadrature
point data...
3071 Assembling system...
norm of rhs is 1.78171e+10
3072 Solver converged in 115 iterations.
3073 Updating quadrature
point data...
3077 Assembling system...
norm of rhs is 1.7737e+10
3078 Solver converged in 112 iterations.
3079 Updating quadrature
point data...
3083 Assembling system...
norm of rhs is 1.77127e+10
3084 Solver converged in 111 iterations.
3085 Updating quadrature
point data...
3089 Assembling system...
norm of rhs is 1.78207e+10
3090 Solver converged in 113 iterations.
3091 Updating quadrature
point data...
3094Timestep 10 at time 10
3095 Assembling system...
norm of rhs is 1.83544e+10
3096 Solver converged in 115 iterations.
3097 Updating quadrature
point data...
3100[100%] Built target
run
3101make
run 176.82s user 0.15s system 198% cpu 1:28.94 total
3103In other words, it is computing on 12,000 cells and with some 52,000
3104unknowns. Not a whole lot, but enough
for a coupled three-dimensional
3105problem to keep a computer busy
for a
while. At the
end of the day,
3106this is what we have
for output:
3109-rw-r--r-- 1 drwells users 1706059 Feb 13 19:36 solution-0010.000.vtu
3110-rw-r--r-- 1 drwells users 761 Feb 13 19:36 solution-0010.pvtu
3111-rw-r--r-- 1 drwells users 33 Feb 13 19:36 solution-0010.visit
3112-rw-r--r-- 1 drwells users 1707907 Feb 13 19:36 solution-0009.000.vtu
3113-rw-r--r-- 1 drwells users 761 Feb 13 19:36 solution-0009.pvtu
3114-rw-r--r-- 1 drwells users 33 Feb 13 19:36 solution-0009.visit
3115-rw-r--r-- 1 drwells users 1703771 Feb 13 19:35 solution-0008.000.vtu
3116-rw-r--r-- 1 drwells users 761 Feb 13 19:35 solution-0008.pvtu
3117-rw-r--r-- 1 drwells users 33 Feb 13 19:35 solution-0008.visit
3118-rw-r--r-- 1 drwells users 1693671 Feb 13 19:35 solution-0007.000.vtu
3119-rw-r--r-- 1 drwells users 761 Feb 13 19:35 solution-0007.pvtu
3120-rw-r--r-- 1 drwells users 33 Feb 13 19:35 solution-0007.visit
3121-rw-r--r-- 1 drwells users 1681847 Feb 13 19:35 solution-0006.000.vtu
3122-rw-r--r-- 1 drwells users 761 Feb 13 19:35 solution-0006.pvtu
3123-rw-r--r-- 1 drwells users 33 Feb 13 19:35 solution-0006.visit
3124-rw-r--r-- 1 drwells users 1670115 Feb 13 19:35 solution-0005.000.vtu
3125-rw-r--r-- 1 drwells users 761 Feb 13 19:35 solution-0005.pvtu
3126-rw-r--r-- 1 drwells users 33 Feb 13 19:35 solution-0005.visit
3127-rw-r--r-- 1 drwells users 1658559 Feb 13 19:35 solution-0004.000.vtu
3128-rw-r--r-- 1 drwells users 761 Feb 13 19:35 solution-0004.pvtu
3129-rw-r--r-- 1 drwells users 33 Feb 13 19:35 solution-0004.visit
3130-rw-r--r-- 1 drwells users 1639983 Feb 13 19:35 solution-0003.000.vtu
3131-rw-r--r-- 1 drwells users 761 Feb 13 19:35 solution-0003.pvtu
3132-rw-r--r-- 1 drwells users 33 Feb 13 19:35 solution-0003.visit
3133-rw-r--r-- 1 drwells users 1625851 Feb 13 19:35 solution-0002.000.vtu
3134-rw-r--r-- 1 drwells users 761 Feb 13 19:35 solution-0002.pvtu
3135-rw-r--r-- 1 drwells users 33 Feb 13 19:35 solution-0002.visit
3136-rw-r--r-- 1 drwells users 1616035 Feb 13 19:34 solution-0001.000.vtu
3137-rw-r--r-- 1 drwells users 761 Feb 13 19:34 solution-0001.pvtu
3138-rw-r--r-- 1 drwells users 33 Feb 13 19:34 solution-0001.visit
3142If we visualize these files with VisIt or Paraview, we get to see the full picture
3143of the disaster our forced compression wreaks on the
cylinder (colors in the
3144images encode the
norm of the stress in the material):
3147<div class=
"threecolumn" style=
"width: 80%">
3148 <div class=
"parent">
3149 <div class=
"img" align=
"center">
3150 <img src=
"https://www.dealii.org/images/steps/developer/step-18.sequential-0002.0000.png"
3154 <div class=
"text" align=
"center">
3158 <div class=
"parent">
3159 <div class=
"img" align=
"center">
3160 <img src=
"https://www.dealii.org/images/steps/developer/step-18.sequential-0005.0000.png"
3164 <div class=
"text" align=
"center">
3168 <div class=
"parent">
3169 <div class=
"img" align=
"center">
3170 <img src=
"https://www.dealii.org/images/steps/developer/step-18.sequential-0007.0000.png"
3174 <div class=
"text" align=
"center">
3181<div class=
"threecolumn" style=
"width: 80%">
3182 <div class=
"parent">
3183 <div class=
"img" align=
"center">
3184 <img src=
"https://www.dealii.org/images/steps/developer/step-18.sequential-0008.0000.png"
3188 <div class=
"text" align=
"center">
3192 <div class=
"parent">
3193 <div class=
"img" align=
"center">
3194 <img src=
"https://www.dealii.org/images/steps/developer/step-18.sequential-0009.0000.png"
3198 <div class=
"text" align=
"center">
3202 <div class=
"parent">
3203 <div class=
"img" align=
"center">
3204 <img src=
"https://www.dealii.org/images/steps/developer/step-18.sequential-0010.0000.png"
3208 <div class=
"text" align=
"center">
3215As is clearly visible, as we keep compressing the
cylinder, it starts
3216to bow out near the fully constrained bottom surface and, after about eight
3217time units, buckle in an azimuthally
symmetric manner.
3220Although the result appears plausible for the
symmetric geometry and loading,
3221it is yet to be established whether or not the computation is fully converged.
3222In order to see whether it is, we ran the program again with
one more global
3223refinement at the beginning and with the time step halved. This would have
3224taken a very long time on a single machine, so we used a proper workstation and
3225ran it on 16 processors in
parallel. The beginning of the output now looks like
3228Timestep 1 at time 0.5
3230 Number of active cells: 29696 (by
partition: 1808+1802+1894+1881+1870+1840+1884+1810+1876+1818+1870+1884+1854+1903+1816+1886)
3231 Number of degrees of freedom: 113100 (by
partition: 6936+6930+7305+7116+7326+6869+7331+6786+7193+6829+7093+7162+6920+7280+6843+7181)
3232 Assembling system...
norm of rhs is 1.10765e+10
3233 Solver converged in 209 iterations.
3234 Updating quadrature
point data...
3236 Number of active cells: 102034 (by
partition: 6387+6202+6421+6341+6408+6201+6428+6428+6385+6294+6506+6244+6417+6527+6299+6546)
3237 Number of degrees of freedom: 359337 (by
partition: 23255+21308+24774+24019+22304+21415+22430+22184+22298+21796+22396+21592+22325+22553+21977+22711)
3238 Assembling system...
norm of rhs is 1.35759e+10
3239 Solver converged in 268 iterations.
3240 Updating quadrature
point data...
3244 Assembling system...
norm of rhs is 1.34674e+10
3245 Solver converged in 267 iterations.
3246 Updating quadrature
point data...
3249Timestep 3 at time 1.5
3250 Assembling system...
norm of rhs is 1.33607e+10
3251 Solver converged in 265 iterations.
3252 Updating quadrature
point data...
3256 Assembling system...
norm of rhs is 1.32558e+10
3257 Solver converged in 263 iterations.
3258 Updating quadrature
point data...
3263Timestep 20 at time 10
3264 Assembling system...
norm of rhs is 1.47755e+10
3265 Solver converged in 425 iterations.
3266 Updating quadrature
point data...
3269That
's quite a good number of unknowns, given that we are in 3d. The output of
3270this program are 16 files for each time step:
3272$ ls -l solution-0001*
3273-rw-r--r-- 1 wellsd2 user 761065 Feb 13 21:09 solution-0001.000.vtu
3274-rw-r--r-- 1 wellsd2 user 759277 Feb 13 21:09 solution-0001.001.vtu
3275-rw-r--r-- 1 wellsd2 user 761217 Feb 13 21:09 solution-0001.002.vtu
3276-rw-r--r-- 1 wellsd2 user 761605 Feb 13 21:09 solution-0001.003.vtu
3277-rw-r--r-- 1 wellsd2 user 756917 Feb 13 21:09 solution-0001.004.vtu
3278-rw-r--r-- 1 wellsd2 user 752669 Feb 13 21:09 solution-0001.005.vtu
3279-rw-r--r-- 1 wellsd2 user 735217 Feb 13 21:09 solution-0001.006.vtu
3280-rw-r--r-- 1 wellsd2 user 750065 Feb 13 21:09 solution-0001.007.vtu
3281-rw-r--r-- 1 wellsd2 user 760273 Feb 13 21:09 solution-0001.008.vtu
3282-rw-r--r-- 1 wellsd2 user 777265 Feb 13 21:09 solution-0001.009.vtu
3283-rw-r--r-- 1 wellsd2 user 772469 Feb 13 21:09 solution-0001.010.vtu
3284-rw-r--r-- 1 wellsd2 user 760833 Feb 13 21:09 solution-0001.011.vtu
3285-rw-r--r-- 1 wellsd2 user 782241 Feb 13 21:09 solution-0001.012.vtu
3286-rw-r--r-- 1 wellsd2 user 748905 Feb 13 21:09 solution-0001.013.vtu
3287-rw-r--r-- 1 wellsd2 user 738413 Feb 13 21:09 solution-0001.014.vtu
3288-rw-r--r-- 1 wellsd2 user 762133 Feb 13 21:09 solution-0001.015.vtu
3289-rw-r--r-- 1 wellsd2 user 1421 Feb 13 21:09 solution-0001.pvtu
3290-rw-r--r-- 1 wellsd2 user 364 Feb 13 21:09 solution-0001.visit
3294Here are first the mesh on which we compute as well as the partitioning
3295for the 16 processors:
3298<div class="twocolumn" style="width: 80%">
3299 <div class="parent">
3300 <div class="img" align="center">
3301 <img src="https://www.dealii.org/images/steps/developer/step-18.parallel-000mesh.png"
3302 alt="Discretization"
3305 <div class="text" align="center">
3309 <div class="parent">
3310 <div class="img" align="center">
3311 <img src="https://www.dealii.org/images/steps/developer/step-18.parallel-0002.p.png"
3312 alt="Parallel partitioning"
3315 <div class="text" align="center">
3316 Parallel partitioning
3322Finally, here is the same output as we have shown before for the much smaller
3325<div class="threecolumn" style="width: 80%">
3326 <div class="parent">
3327 <div class="img" align="center">
3328 <img src="https://www.dealii.org/images/steps/developer/step-18.parallel-0002.s.png"
3332 <div class="text" align="center">
3336 <div class="parent">
3337 <div class="img" align="center">
3338 <img src="https://www.dealii.org/images/steps/developer/step-18.parallel-0005.s.png"
3342 <div class="text" align="center">
3346 <div class="parent">
3347 <div class="img" align="center">
3348 <img src="https://www.dealii.org/images/steps/developer/step-18.parallel-0007.s.png"
3352 <div class="text" align="center">
3359<div class="threecolumn" style="width: 80%">
3360 <div class="parent">
3361 <div class="img" align="center">
3362 <img src="https://www.dealii.org/images/steps/developer/step-18.parallel-0008.s.png"
3366 <div class="text" align="center">
3370 <div class="parent">
3371 <div class="img" align="center">
3372 <img src="https://www.dealii.org/images/steps/developer/step-18.parallel-0009.s.png"
3376 <div class="text" align="center">
3380 <div class="parent">
3381 <div class="img" align="center">
3382 <img src="https://www.dealii.org/images/steps/developer/step-18.parallel-0010.s.png"
3386 <div class="text" align="center">
3393As before, we observe that at high axial compression the cylinder begins
3394to buckle, but this time ultimately collapses on itself. In contrast to our
3395first run, towards the end of the simulation the deflection pattern becomes
3396nonsymmetric (the central bulge deflects laterally). The model clearly does not
3397provide for this (all our forces and boundary deflections are symmetric) but the
3398effect is probably physically correct anyway: in reality, small inhomogeneities
3399in the body's material properties would lead it to buckle to
one side
3400to evade the forcing; in numerical simulations, small perturbations
3401such as numerical round-off or an inexact solution of a linear system
3402by an iterative solver could have the same effect. Another typical source
for
3403asymmetries in adaptive computations is that only a certain fraction of cells
3404is refined in each step, which may lead to asymmetric meshes even
if the
3408If
one compares
this with the previous
run, the results both qualitatively
3409and quantitatively different. The previous computation was
3410therefore certainly not converged, though we can
't say for sure anything about
3411the present one. One would need an even finer computation to find out. However,
3412the point may be moot: looking at the last picture in detail, it is pretty
3413obvious that not only is the linear small deformation model we chose completely
3414inadequate, but for a realistic simulation we would also need to make sure that
3415the body does not intersect itself during deformation (if we continued
3416compressing the cylinder we would observe some self-intersection).
3417Without such a formulation we cannot expect anything to make physical sense,
3418even if it produces nice pictures!
3421<a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
3424The program as is does not really solve an equation that has many applications
3425in practice: quasi-static material deformation based on a purely elastic law
3426is almost boring. However, the program may serve as the starting point for
3427more interesting experiments, and that indeed was the initial motivation for
3428writing it. Here are some suggestions of what the program is missing and in
3429what direction it may be extended:
3431<a name="Plasticitymodels"></a><h5>Plasticity models</h5>
3434 The most obvious extension is to use a more
3435realistic material model for large-scale quasistatic deformation. The natural
3436choice for this would be plasticity, in which a nonlinear relationship between
3437stress and strain replaces equation <a href="#step_18.stress-strain">[stress-strain]</a>. Plasticity
3438models are usually rather complicated to program since the stress-strain
3439dependence is generally non-smooth. The material can be thought of being able
3440to withstand only a maximal stress (the yield stress) after which it starts to
3441“flow”. A mathematical description to this can be given in the form of a
3442variational inequality, which alternatively can be treated as minimizing the
3446 (\varepsilon(\mathbf{u}), C\varepsilon(\mathbf{u}))_{\Omega}
3447 - (\mathbf{f}, \mathbf{u})_{\Omega} - (\mathbf{b}, \mathbf{u})_{\Gamma_N},
3449subject to the constraint
3451 f(\sigma(\mathbf{u})) \le 0
3453on the stress. This extension makes the problem to be solved in each time step
3454nonlinear, so we need another loop within each time step.
3456Without going into further details of this model, we refer to the excellent
3457book by Simo and Hughes on “Computational Inelasticity” for a
3458comprehensive overview of computational strategies for solving plastic
3459models. Alternatively, a brief but concise description of an algorithm for
3460plasticity is given in an article by S. Commend, A. Truty, and Th. Zimmermann;
3464<a name="Stabilizationissues"></a><h5>Stabilization issues</h5>
3467The formulation we have chosen, i.e. using
3468piecewise (bi-, tri-)linear elements for all components of the displacement
3469vector, and treating the stress as a variable dependent on the displacement is
3470appropriate for most materials. However, this so-called displacement-based
3471formulation becomes unstable and exhibits spurious modes for incompressible or
3472nearly-incompressible materials. While fluids are usually not elastic (in most
3473cases, the stress depends on velocity gradients, not displacement gradients,
3474although there are exceptions such as electro-rheologic fluids), there are a
3475few solids that are nearly incompressible, for example rubber. Another case is
3476that many plasticity models ultimately let the material become incompressible,
3477although this is outside the scope of the present program.
3479Incompressibility is characterized by Poisson's ratio
3481 \nu = \frac{\lambda}{2(\lambda+\mu)},
3483where @f$\lambda,\mu@f$ are the Lamé constants of the material.
3484Physical constraints indicate that @f$-1\le \nu\le \frac 12@f$ (the condition
3485also follows from mathematical stability considerations). If @f$\nu@f$
3486approaches @f$\frac 12@f$, then the material becomes incompressible. In that
3487case, pure displacement-based formulations are no longer appropriate
for the
3488solution of such problems, and stabilization techniques have to be employed
3489for a stable and accurate solution. The book and paper cited above give
3490indications as to how to
do this, but there is also a large
volume of
3491literature on
this subject; a good start to get an overview of the topic can
3492be found in the references of the paper by H.-Y. Duan and Q. Lin; @cite DL05.
3495<a name=
"Refinementduringtimesteps"></a><h5>Refinement during timesteps</h5>
3498In the present form, the program
3499only refines the
initial mesh a number of times, but then never again. For any
3500kind of realistic simulation,
one would want to extend
this so that the mesh
3501is refined and coarsened every few time steps instead. This is not hard to
do,
3502in fact, but has been left
for future tutorial programs or as an exercise,
if
3505The main complication
one has to overcome is that
one has to
3506transfer the data that is stored in the quadrature points of the cells of the
3507old mesh to the
new mesh, preferably by some sort of projection scheme. The
3508general approach to
this would go like
this:
3510- At the beginning, the data is only available in the quadrature points of
3511 individual cells, not as a finite element field that is defined everywhere.
3513- So let us find a finite element field that <i>is</i> defined everywhere so
3514 that we can later
interpolate it to the quadrature points of the
new
3515 mesh. In
general, it will be difficult to find a continuous finite element
3516 field that matches the
values in the quadrature points exactly because the
3517 number of degrees of freedom of these fields does not match the number of
3518 quadrature points there are, and the nodal
values of
this global field will
3519 either be over- or underdetermined. But it is usually not very difficult to
3520 find a discontinuous field that matches the
values in the quadrature points;
3521 for example,
if you have a
QGauss(2) quadrature formula (i.
e. 4 points per
3522 cell in 2
d, 8 points in 3
d), then
one would use a finite element of kind
3523 FE_DGQ(1), i.
e. bi-/tri-linear
functions as these have 4 degrees of freedom
3524 per cell in 2
d and 8 in 3
d.
3526- There are
functions that can make this conversion from individual points to
3527 a global field simpler. The following piece of pseudo-code should help if
3528 you use a
QGauss(2) quadrature formula. Note that the multiplication by the
3529 projection
matrix below takes a vector of scalar components, i.
e., we can only
3530 convert
one set of scalars at a time from the quadrature points to the degrees
3531 of freedom and vice versa. So we need to store each component of stress separately,
3532 which requires <code>dim*dim</code> vectors. We'll store this
set of vectors in a 2D array to
3533 make it easier to read off components in the same way you would the stress tensor.
3534 Thus, we'll
loop over the components of stress on each cell and store
3535 these
values in the global history field. (The prefix <code>history_</code>
3536 indicates that we work with quantities related to the history variables defined
3537 in the quadrature points.)
3539 FE_DGQ<dim> history_fe (1);
3541 history_dof_handler.distribute_dofs (history_fe);
3544 history_field (dim,
std::vector<
Vector<
double> >(dim)),
3545 local_history_values_at_qpoints (dim,
std::vector<
Vector<
double> >(dim)),
3546 local_history_fe_values (dim,
std::vector<
Vector<
double> >(dim));
3548 for (
unsigned int i=0; i<dim; i++)
3549 for (
unsigned int j=0; j<dim; j++)
3551 history_field[i][j].
reinit(history_dof_handler.n_dofs());
3552 local_history_values_at_qpoints[i][j].reinit(quadrature.size());
3553 local_history_fe_values[i][j].reinit(history_fe.n_dofs_per_cell());
3560 quadrature, quadrature,
3561 qpoint_to_dof_matrix);
3564 endc = dof_handler.end(),
3565 dg_cell = history_dof_handler.begin_active();
3567 for (; cell!=endc; ++cell, ++dg_cell)
3570 PointHistory<dim> *local_quadrature_points_history
3571 =
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
3573 Assert (local_quadrature_points_history >= &quadrature_point_history.front(),
3575 Assert (local_quadrature_points_history < &quadrature_point_history.back(),
3578 for (
unsigned int i=0; i<dim; i++)
3579 for (
unsigned int j=0; j<dim; j++)
3581 for (
unsigned int q=0; q<quadrature.size(); ++q)
3582 local_history_values_at_qpoints[i][j](q)
3583 = local_quadrature_points_history[q].old_stress[i][j];
3585 qpoint_to_dof_matrix.vmult (local_history_fe_values[i][j],
3586 local_history_values_at_qpoints[i][j]);
3588 dg_cell->set_dof_values (local_history_fe_values[i][j],
3589 history_field[i][j]);
3594- Now that we have a global field, we can
refine the mesh and transfer the
3596 interpolate everything from the old to the
new mesh.
3598- In a
final step, we have to get the data back from the now interpolated
3599 global field to the quadrature points on the
new mesh. The following code
3603 history_fe.n_dofs_per_cell());
3607 dof_to_qpoint_matrix);
3610 endc = dof_handler.end(),
3611 dg_cell = history_dof_handler.begin_active();
3613 for (; cell != endc; ++cell, ++dg_cell)
3615 PointHistory<dim> *local_quadrature_points_history
3616 =
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
3618 Assert (local_quadrature_points_history >= &quadrature_point_history.front(),
3620 Assert (local_quadrature_points_history < &quadrature_point_history.back(),
3623 for (
unsigned int i=0; i<dim; i++)
3624 for (
unsigned int j=0; j<dim; j++)
3626 dg_cell->get_dof_values (history_field[i][j],
3627 local_history_fe_values[i][j]);
3629 dof_to_qpoint_matrix.vmult (local_history_values_at_qpoints[i][j],
3630 local_history_fe_values[i][j]);
3632 for (
unsigned int q=0; q<quadrature.size(); ++q)
3633 local_quadrature_points_history[q].old_stress[i][j]
3634 = local_history_values_at_qpoints[i][j](q);
3638It becomes a bit more complicated once we
run the program in
parallel, since
3639then each process only stores
this data
for the cells it owned on the old
3640mesh. That said,
using a
parallel vector for <code>history_field</code> will
3641do the trick
if you put a
call to <code>
compress</code> after the transfer
3642from quadrature points into the global vector.
3645<a name=
"Ensuringmeshregularity"></a><h5>Ensuring mesh regularity</h5>
3648At present, the program makes no attempt
3649to make sure that a cell, after moving its
vertices at the
end of the time
3651positive and bounded away from
zero everywhere). It is, in fact, not very hard
3652to
set boundary
values and forcing terms in such a way that
one gets distorted
3653and inverted cells rather quickly. Certainly, in some cases of large
3654deformation,
this is unavoidable with a mesh of finite mesh size, but in some
3655other cases
this should be preventable by appropriate mesh refinement and/or a
3656reduction of the time step size. The program does not
do that, but a more
3657sophisticated version definitely should employ some sort of heuristic defining
3658what amount of deformation of cells is acceptable, and what isn
't.
3661<a name="PlainProg"></a>
3662<h1> The plain program</h1>
3663@include "step-18.cc"
void attach_dof_handler(const DoFHandlerType &)
active_cell_iterator begin_active(const unsigned int level=0) const
virtual void vector_value_list(const std::vector< Point< dim > > &points, std::vector< Vector< RangeNumberType > > &values) const
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &values) const
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
__global__ void reduction(Number *result, const Number *v, const size_type N)
__global__ void set(Number *val, const Number s, const size_type N)
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
void approximate(SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > const &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)
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.
static const types::blas_int zero
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
@ general
No special properties.
static const 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.)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
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)
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)
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
VectorType::value_type * end(VectorType &V)
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
std::string compress(const std::string &input)
void run(const Iterator &begin, const typename identity< Iterator >::type &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)
TriangulationBase< dim, spacedim > Triangulation
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation