]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Expand discussion and move possibilities for extensions to the end, where they usuall...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 3 Jun 2011 03:00:38 +0000 (03:00 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 3 Jun 2011 03:00:38 +0000 (03:00 +0000)
git-svn-id: https://svn.dealii.org/trunk@23784 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-18/doc/intro.dox
deal.II/examples/step-18/doc/results.dox

index 10940ee4cdba615a60aab1afa836489863655230..d912f3b072d5752127061b3e8c9cf5d16a665913 100644 (file)
@@ -555,131 +555,6 @@ the results section will show what displacement pattern comes out of this
 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
-&ldquo;flow&rdquo;. 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 &ldquo;Computational Inelasticity&rdquo; 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 &ldquo;Stabilized finite elements applied to
-elastoplasticity: I. Mixed displacement-pressure formulation&rdquo;
-(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 &ldquo;Mixed finite elements of least-squares type for
-elasticity&rdquo; (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
index 4a80e0ad4aae0cbd7bad00fb832abae572ff29ef..8bb54ed7c80686777a3942f1a7a386aa87bf70de 100644 (file)
@@ -334,3 +334,190 @@ simulation we would also need to make sure that the body does not intersect
 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
+&ldquo;flow&rdquo;. 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 &ldquo;Computational Inelasticity&rdquo; 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 &ldquo;Stabilized finite elements applied to
+elastoplasticity: I. Mixed displacement-pressure formulation&rdquo;
+(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 &ldquo;Mixed finite elements of least-squares type for
+elasticity&rdquo; (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.
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.