simulation.
-<h3>Possible directions for extensions</h3>
-
-The program as is does not really solve an equation that has many applications
-in practice: quasi-static material deformation based on a purely elastic law
-is almost boring. However, the program may serve as the starting point for
-more interesting experiments, and that indeed was the initial motivation for
-writing it. Here are some suggestions of what the program is missing and in
-what direction it may be extended:
-
-<h5>Plasticity models</h5>
-
- The most obvious extension is to use a more
-realistic material model for large-scale quasistatic deformation. The natural
-choice for this would be plasticity, in which a nonlinear relationship between
-stress and strain replaces equation <a href="#step-18.stress-strain">[stress-strain]</a>. Plasticity
-models are usually rather complicated to program since the stress-strain
-dependence is generally non-smooth. The material can be thought of being able
-to withstand only a maximal stress (the yield stress) after which it starts to
-“flow”. A mathematical description to this can be given in the form of a
-variational inequality, which alternatively can be treated as minimizing the
-elastic energy
-@f[
- E(\mathbf{u}) =
- (\varepsilon(\mathbf{u}), C\varepsilon(\mathbf{u}))_{\Omega}
- - (\mathbf{f}, \mathbf{u})_{\Omega} - (\mathbf{b}, \mathbf{u})_{\Gamma_N},
-@f]
-subject to the constraint
-@f[
- f(\sigma(\mathbf{u})) \le 0
-@f]
-on the stress. This extension makes the problem to be solved in each time step
-nonlinear, so we need another loop within each time step.
-
-Without going into further details of this model, we refer to the excellent
-book by Simo and Hughes on “Computational Inelasticity” for a
-comprehensive overview of computational strategies for solving plastic
-models. Alternatively, a brief but concise description of an algorithm for
-plasticity is given in an article by S. Commend, A. Truty, and Th. Zimmermann,
-titled “Stabilized finite elements applied to
-elastoplasticity: I. Mixed displacement-pressure formulation”
-(Computer Methods in Applied Mechanics and Engineering, vol. 193,
-pp. 3559-3586, 2004).
-
-
-<h5>Stabilization issues</h5>
-
-The formulation we have chosen, i.e. using
-piecewise (bi-, tri-)linear elements for all components of the displacement
-vector, and treating the stress as a variable dependent on the displacement is
-appropriate for most materials. However, this so-called displacement-based
-formulation becomes unstable and exhibits spurious modes for incompressible or
-nearly-incompressible materials. While fluids are usually not elastic (in most
-cases, the stress depends on velocity gradients, not displacement gradients,
-although there are exceptions such as electro-rheologic fluids), there are a
-few solids that are nearly incompressible, for example rubber. Another case is
-that many plasticity models ultimately let the material become incompressible,
-although this is outside the scope of the present program.
-
-Incompressibility is characterized by Poisson's ratio
-@f[
- \nu = \frac{\lambda}{2(\lambda+\mu)},
-@f]
-where $\lambda,\mu$ are the Lam\'e constants of the material.
-Physical constraints indicate that $-1\le \nu\le \frac 12$ (the condition
-also follows from mathematical stability considerations). If $\nu$
-approaches $\frac 12$, then the material becomes incompressible. In that
-case, pure displacement-based formulations are no longer appropriate for the
-solution of such problems, and stabilization techniques have to be employed
-for a stable and accurate solution. The book and paper cited above give
-indications as to how to do this, but there is also a large volume of
-literature on this subject; a good start to get an overview of the topic can
-be found in the references of the paper by
-H.-Y. Duan and Q. Lin on “Mixed finite elements of least-squares type for
-elasticity” (Computer Methods in Applied Mechanics and Engineering, vol. 194,
-pp. 1093-1112, 2005).
-
-
-<h5>Refinement during timesteps</h5>
-
-In the present form, the program
-only refines the initial mesh a number of times, but then never again. For any
-kind of realistic simulation, one would want to extend this so that the mesh
-is refined and coarsened every few time steps instead. This is not hard to do,
-in fact, but has been left for future tutorial programs or as an exercise, if
-you wish. The main complication one has to overcome is that one has to
-transfer the data that is stored in the quadrature points of the cells of the
-old mesh to the new mesh, preferably by some sort of projection scheme. This
-is only slightly messy in the sequential case; in fact, the functions
-FETools::get_projection_from_quadrature_points_matrix will do
-the projection, and the FiniteElement::get_restriction_matrix and
-FiniteElement::get_prolongation_matrix functions will do the
-transfer between mother and child cells. Alternatively, one can use
-FETools::get_projection_from_quadrature_points_matrix to convert data only
-defined in quadrature points to a finite element field, use the
-SolutionTransfer class to move the field to a refined/coarsened version of the
-mesh, and then use FETools::compute_interpolation_to_quadrature_points_matrix
-to get the information back into the quadrature points of the new mesh.
-
-However, whichever way we do it, it becomes complicated
-once we run the program in parallel, since then each process only stores this
-data for the cells it owned on the old mesh, and it may need to know the
-values of the quadrature point data on other cells if the corresponding cells
-on the new mesh are assigned to this process after subdividing the new mesh. A
-global communication of these data elements is therefore necessary, making the
-entire process a little more unpleasant. The step-28 tutorial
-program shows how to work with different meshes at the same time, albeit for a
-different kind of problem, giving indications on how to approach the problem
-for time-dependent adaptivity as well.
-
-
-<h5>Ensuring mesh regularity</h5>
-
-At present, the program makes no attempt
-to make sure that a cell, after moving its vertices at the end of the time
-step, still has a valid geometry (i.e. that its Jacobian determinant is
-positive and bounded away from zero everywhere). It is, in fact, not very hard
-to set boundary values and forcing terms in such a way that one gets distorted
-and inverted cells rather quickly. Certainly, in some cases of large
-deformation, this is unavoidable with a mesh of finite mesh size, but in some
-other cases this should be preventable by appropriate mesh refinement and/or a
-reduction of the time step size. The program does not do that, but a more
-sophisticated version definitely should employ some sort of heuristic defining
-what amount of deformation of cells is acceptable, and what isn't.
-
-
@note PETSc appears not to co-operate well when using
threads and the program crashes when deal.II is compiled in its usual mode
supporting multithreading. The program can therefore only be compiled and run
itself during deformation. Without such a formulation we cannot expect anything
to make physical sense, even if it produces nice pictures!
+
+<h3>Possible directions for extensions</h3>
+
+The program as is does not really solve an equation that has many applications
+in practice: quasi-static material deformation based on a purely elastic law
+is almost boring. However, the program may serve as the starting point for
+more interesting experiments, and that indeed was the initial motivation for
+writing it. Here are some suggestions of what the program is missing and in
+what direction it may be extended:
+
+<h5>Plasticity models</h5>
+
+ The most obvious extension is to use a more
+realistic material model for large-scale quasistatic deformation. The natural
+choice for this would be plasticity, in which a nonlinear relationship between
+stress and strain replaces equation <a href="#step-18.stress-strain">[stress-strain]</a>. Plasticity
+models are usually rather complicated to program since the stress-strain
+dependence is generally non-smooth. The material can be thought of being able
+to withstand only a maximal stress (the yield stress) after which it starts to
+“flow”. A mathematical description to this can be given in the form of a
+variational inequality, which alternatively can be treated as minimizing the
+elastic energy
+@f[
+ E(\mathbf{u}) =
+ (\varepsilon(\mathbf{u}), C\varepsilon(\mathbf{u}))_{\Omega}
+ - (\mathbf{f}, \mathbf{u})_{\Omega} - (\mathbf{b}, \mathbf{u})_{\Gamma_N},
+@f]
+subject to the constraint
+@f[
+ f(\sigma(\mathbf{u})) \le 0
+@f]
+on the stress. This extension makes the problem to be solved in each time step
+nonlinear, so we need another loop within each time step.
+
+Without going into further details of this model, we refer to the excellent
+book by Simo and Hughes on “Computational Inelasticity” for a
+comprehensive overview of computational strategies for solving plastic
+models. Alternatively, a brief but concise description of an algorithm for
+plasticity is given in an article by S. Commend, A. Truty, and Th. Zimmermann,
+titled “Stabilized finite elements applied to
+elastoplasticity: I. Mixed displacement-pressure formulation”
+(Computer Methods in Applied Mechanics and Engineering, vol. 193,
+pp. 3559-3586, 2004).
+
+
+<h5>Stabilization issues</h5>
+
+The formulation we have chosen, i.e. using
+piecewise (bi-, tri-)linear elements for all components of the displacement
+vector, and treating the stress as a variable dependent on the displacement is
+appropriate for most materials. However, this so-called displacement-based
+formulation becomes unstable and exhibits spurious modes for incompressible or
+nearly-incompressible materials. While fluids are usually not elastic (in most
+cases, the stress depends on velocity gradients, not displacement gradients,
+although there are exceptions such as electro-rheologic fluids), there are a
+few solids that are nearly incompressible, for example rubber. Another case is
+that many plasticity models ultimately let the material become incompressible,
+although this is outside the scope of the present program.
+
+Incompressibility is characterized by Poisson's ratio
+@f[
+ \nu = \frac{\lambda}{2(\lambda+\mu)},
+@f]
+where $\lambda,\mu$ are the Lam\'e constants of the material.
+Physical constraints indicate that $-1\le \nu\le \frac 12$ (the condition
+also follows from mathematical stability considerations). If $\nu$
+approaches $\frac 12$, then the material becomes incompressible. In that
+case, pure displacement-based formulations are no longer appropriate for the
+solution of such problems, and stabilization techniques have to be employed
+for a stable and accurate solution. The book and paper cited above give
+indications as to how to do this, but there is also a large volume of
+literature on this subject; a good start to get an overview of the topic can
+be found in the references of the paper by
+H.-Y. Duan and Q. Lin on “Mixed finite elements of least-squares type for
+elasticity” (Computer Methods in Applied Mechanics and Engineering, vol. 194,
+pp. 1093-1112, 2005).
+
+
+<h5>Refinement during timesteps</h5>
+
+In the present form, the program
+only refines the initial mesh a number of times, but then never again. For any
+kind of realistic simulation, one would want to extend this so that the mesh
+is refined and coarsened every few time steps instead. This is not hard to do,
+in fact, but has been left for future tutorial programs or as an exercise, if
+you wish.
+
+The main complication one has to overcome is that one has to
+transfer the data that is stored in the quadrature points of the cells of the
+old mesh to the new mesh, preferably by some sort of projection scheme. The
+general approach to this would go like this:
+
+- At the beginning, the data is only available in the quadrature points of
+ individual cells, not as a finite element field that is defined everywhere.
+
+- So let us find a finite element field that <i>is</i> defined everywhere so
+ that we can later interpolate it to the quadrature points of the new
+ mesh. In general, it will be difficult to find a continuous finite element
+ field that matches the values in the quadrature points exactly because the
+ number of degrees of freedom of these fields does not match the number of
+ quadrature points there are, and the nodal values of this global field will
+ either be over- or underdetermined. But it is usually not very difficult to
+ find a discontinuous field that matches the values in the quadrature points;
+ for example, if you have a QGauss(2) quadrature formula (i.e. 4 points per
+ cell in 2d, 8 points in 3d), then one would use a finite element of kind
+ FE_DGQ(1), i.e. bi-/tri-linear functions as these have 4 degrees of freedom
+ per cell in 2d and 8 in 3d. There are functions that can make this
+ conversion from individual points to a global field simpler; the following
+ piece of pseudo-code should help if you use a QGauss(2) quadrature formula
+ (the prefix <code>history_</code> indicates that we work with quantities
+ related to the history variables defined in the quadrature points):
+ @code
+ FE_DGQ<dim> history_fe (1);
+ DoFHandler<dim> history_dof_handler (triangulation);
+ history_dof_handler.distribute_dofs (history_fe);
+
+ Vector<double> history_field (history_dof_handler.n_dofs());
+
+ FullMatrix<double> qpoint_to_dof_matrix (history_fe.dofs_per_cell,
+ quadrature.size());
+ FETools::compute_projection_from_quadrature_points_matrix
+ (history_fe,
+ quadrature, quadrature,
+ qpoint_to_dof_matrix);
+
+ Vector<double> local_history_values_at_qpoints (quadrature.size());
+ Vector<double> local_history_fe_values (fe.dofs_per_cell);
+ for (cell=...)
+ {
+ ...collect values from quadrature points into
+ local_history_values_at_qpoints...
+ qpoint_to_dof_matrix.vmult (local_history_fe_values,
+ local_history_values_at_qpoints);
+ cell->set_dof_values (local_history_fe_values,
+ history_field);
+ }
+ @endcode
+
+- Now that we have a global field, we can refine the mesh and transfer the
+ history_field vector as usual using the SolutionTransfer class. This will
+ interpolate everything from the old to the new mesh.
+
+- In a final step, we have to get the data back from the now interpolated
+ global field to the quadrature points on the new mesh. The following code
+ will do that:
+ @code
+ FullMatrix<double> dof_to_qpoint_matrix (quadrature.size(),
+ history_fe.dofs_per_cell);
+ FETools::compute_interpolation_to_quadrature_points_matrix
+ (history_fe,
+ quadrature,
+ dof_to_qpoint_matrix);
+
+ Vector<double> local_history_values_at_qpoints (quadrature.size());
+ Vector<double> local_history_fe_values (fe.dofs_per_cell);
+ for (cell=...)
+ {
+ cell->set_get_values (history_field,
+ local_history_fe_values);
+ dof_to_qpoint_matrix.vmult (local_history_values_at_qpoints,
+ local_history_fe_values);
+ ...put values back from local_history_values_at_qpoints
+ quadrature points into...
+ }
+ @endcode
+
+It becomes a bit more complicated once we run the program in parallel, since
+then each process only stores this data for the cells it owned on the old
+mesh. That said, using a parallel vector for <code>history_field</code> will
+do the trick if you put a call to <code>compress</code> after the transfer
+from quadrature points into the global vector.
+
+
+<h5>Ensuring mesh regularity</h5>
+
+At present, the program makes no attempt
+to make sure that a cell, after moving its vertices at the end of the time
+step, still has a valid geometry (i.e. that its Jacobian determinant is
+positive and bounded away from zero everywhere). It is, in fact, not very hard
+to set boundary values and forcing terms in such a way that one gets distorted
+and inverted cells rather quickly. Certainly, in some cases of large
+deformation, this is unavoidable with a mesh of finite mesh size, but in some
+other cases this should be preventable by appropriate mesh refinement and/or a
+reduction of the time step size. The program does not do that, but a more
+sophisticated version definitely should employ some sort of heuristic defining
+what amount of deformation of cells is acceptable, and what isn't.
+